Seismic tomography of the asthenosphere beneath mid-ocean ridges has produced images of wave speed and anisotropy that are asymmetric across the ridge axis. These features have been interpreted as resulting from an asymmetric distribution of upwelling and melting. Using computational models of coupled magma/mantle dynamics beneath mid-ocean ridges, I show that such asymmetry should be expected if buoyancy forces contribute to mantle upwelling beneath ridges. The sole source of buoyancy considered here is the dynamic retention of less dense magma within the pores of the mantle matrix. Through a scaling analysis and comparison with a suite of simulations, I derive a quantitative prediction of the contribution of such buoyancy to upwelling; this prediction of convective vigor is based on parameters that, for the Earth, can be constrained through natural observations and experiments. I show how the width of the melting region and the crustal thickness, as well as the susceptibility to asymmetric upwelling, are related to convective vigor. I consider three causes of symmetry breaking: gradients in mantle potential temperature and composition and ridge migration. I also report that in numerical experiments performed for this study, the fluid dynamical instability associated with porosity/shear band formation is not observed to occur.

1. Introduction

The global mid-ocean ridge (MOR) system is greater than 50,000 km in length and spans all of the ocean basins on the planet [Ito and Dunn, 2009]. Physical and chemical properties of the ridge axis, and their geographical variation, provide fundamental constraints on the dynamics and petrology of the otherwise inaccessible mantle below. A primary variable in the geographical diversity of the ridge system is the rate of divergence of adjacent plates; the full-spreading rate varies from ultraslow (∼1 cm/yr) to fast (∼15 cm/yr). It has long been a goal of MOR studies to correlate the spreading rate with other features that are observed to vary geographically. These efforts have met with partial success [e.g., Bown and White, 1994], but it is clear that much of the diversity of the ridge system cannot be explained by spreading rate alone [e.g., Small, 1994; Rubin and Sinton, 2007]. Mantle heterogeneity must play a role too, and Langmuir et al. [1992] have shown that geographic variation in mantle potential temperature can account for significant differences between segments of the ridge. Despite such progress, and despite the simplicity of MOR volcanism relative to other volcanic settings, many observations of MORs remain unexplained. These observations may help us to refine models of the fluid dynamics and petrology of the asthenosphere beneath ridges.

Past work on subridge mantle dynamics has focused on two basic models: passive and active flow. Passive mantle flow is driven by external forcing, imposed as boundary conditions on models. Tectonic plate divergence is the primary example of such forcing, and when applied as a boundary condition to viscous flow models of mantle dynamics, it gives rise to upwelling flow that turns a corner and moves laterally away from the ridge axis [Batchelor, 1967; McKenzie, 1969; Lachenbruch, 1976]. In contrast, active flow is calculated by accounting for buoyancy forces as well as boundary conditions. In a model that accounts for buoyancy, density gradients can give rise to time-dependent, convective flows that are not present in passive models. Previous authors have argued that the effects of active flow are absent from MOR data, and therefore the passive model (i.e., neglecting buoyancy forces) is appropriate [e.g., Spiegelman and Reynolds, 1999]. The purpose of the present work is to reconsider the effects of buoyancy on MOR dynamics, and show that it may help to explain the asymmetric distribution of melt observed in the subridge mantle by seismic tomography.

Because most magma produced in the MOR environment is focused to the ridge axis, the spatial distribution of mantle melting cannot typically be inferred from petrologic observations. Hence petrologic studies of mid-ocean ridge basalt (MORB) formation [e.g., McKenzie and Bickle, 1988; Langmuir et al., 1992] have tended to favor the simpler theory of passive flow of the mantle beneath ridges. Petrology tells us that melting occurs in the region where the mantle is upwelling at temperatures above the mantle solidus. For passive flow, this region is found to be approximately triangular in cross section; it is centered beneath the ridge axis and symmetric across it.

Observations of the mantle asthenosphere made with seismic tomography are a challenge to this simple picture of symmetric, plate-driven flow. Tomographic studies exploit the sensitivity of seismic wave speed to the presence of magma at mantle grain boundaries [e.g., Blackman and Kendall, 1997; Yang et al., 2007] to resolve the outlines of the melting region. Wang et al. [2009] applied such techniques to produce a survey of seismic wave speeds in the mantle beneath the Gulf of California. The Gulf hosts a geometrically complex ridge system that creates oceanic crust between Baja California and mainland Mexico; the overall ridge trend is oblique to the spreading direction. Nonetheless, perpendicularly spreading segments are discernible along the ridge, and a model invoking passive mantle upwelling would predict melting and anomalously slow wave speeds beneath these segments. In contrast, Wang et al. [2009] find that the slow mantle seismic anomalies are not geographically associated with spreading segments, but are instead located at tens to hundreds of kilometers off axis. The authors interpret their findings as evidence of buoyancy-driven mantle upwelling. The geometric complexity of the ridge, and the close juxtaposition of continental crust make the Gulf a setting that requires three-dimensional flow models to elucidate the convective dynamics.

The MELT experiment [Forsyth et al., 1998a] was conducted in a simpler tectonic setting of the southern East Pacific Rise (EPR). It revealed a broad, seismically slow region of mantle and interpreted it as a zone of partial melting, with porosity ∼1% [Toomey et al., 1998]. Various indicators of asymmetry of melt production across the ridge axis were also discovered by the MELT experiment. Scheirer et al. [1998] analyzed shipboard measurements of gravity, seamount distribution, and subsidence rate in the MELT region and found consistent asymmetry among these indicators. Forsyth et al. [1998b] used Rayleigh wave tomography to image the shear wave velocity structure of the mantle between 20 and 70 km depth in the same geographic area. Their results exhibited a broad region of anomalously slow velocities, centered to the west of the ridge axis, that they attribute to the presence of partial melt. A similar pattern of asymmetry was observed with teleseismic body wave tomography [Toomey et al., 1998; Hammond and Toomey, 2003], shear wave splitting [Wolfe and Solomon, 1998], and a magnetotelluric survey [Evans et al., 1999]. Dunn and Forsyth [2003] used short-period Love waves to probe S wave velocity to a depth of ∼50 km for up to ∼100 km on each side of the axis and found strong cross-ridge asymmetry in the slow anomaly. Together these studies provide a clear indication that in the MELT region at present, there is more mantle melting and residual mantle porosity on the west side of the southern EPR.

Two studies, Conder et al. [2002] and Toomey et al. [2002], have sought to model the observed asymmetry at the MELT region by imposing asymmetric boundary conditions on a passive flow model. Both invoke a combination of anomalously hot mantle beneath the Pacific plate, and a large-scale pressure gradient driving mantle from west to east. The authors justify inclusion of these special circumstances by the proximity of the EPR to hot spots of the Pacific Superswell. To reproduce the observed tomographic asymmetry, however, the models require either very rapid, pressure-driven mantle flow across the ridge (∼30 cm/yr), or a very large temperature anomaly (>100 K) that is confined to the mantle beneath the Pacific plate. In the latter case, the modeled asymmetry is critically dependent on the precise depth extent of the imposed thermal anomaly on the inflow boundary. Though the MELT region could be a special case, one expects a more general explanation for the observed asymmetry that requires only moderate external forcing. Furthermore, lateral temperature and melting rate gradients of the magnitude invoked by Conder et al. [2002] and Toomey et al. [2002] would yield sharp gradients in mantle density, and hence would likely cause buoyancy-driven flow, a scenario not considered by those studies.

Chemical and thermal buoyancy are known to be significant in the balance of stresses in the deeper mantle; they are the main driving forces for mantle convection and plate tectonics [e.g., Hager and O'Connell, 1981; Faccenna and Becker, 2010]. The role of buoyancy forces beneath mid-ocean ridges remains controversial, however. Rabinowicz et al. [1984] presented an early model of convection beneath ridges that was inspired by observations of flow structures in ophiolites [Nicolas and Violette, 1982]. They reasoned that because magma has a lower density than the solid mantle, partially molten mantle is buoyant relative to unmelted mantle. Rabinowicz et al. [1984] found that this buoyancy, associated with the residual porosity of the mantle, can drive enhanced upwelling and melting, which can result in a self-reinforcing convective feedback. Subsequent investigations included explicit calculation of coupled mantle flow, melting, and magmatic segregation and explored the dynamics of active upwelling beneath MORs [e.g., Buck and Su, 1989; Scott and Stevenson, 1989; Spiegelman, 1996] and extending lithosphere [Hernlund et al., 2008a, 2008b]. Tackley and Stevenson [1993] employed these concepts to produce a model of self-perpetuating volcanism on terrestrial planets. A selection of these contributions are reviewed in more detail below.

All previous studies of active flow beneath MORs have explicitly assumed symmetry across the ridge axis. However, buoyancy-driven flows are known to produce symmetry-breaking behavior in fluid systems [e.g., Turner, 1973]. As an illustration of this, consider the end-member case where buoyancy forces in a partially molten aggregate are much greater than the stresses imposed by plate motion: upwelling will be centered on the regions of maximum buoyancy, rather than aligned with the zone of plate divergence. If buoyancy created by melting is important, we should expect that small, across-axis gradients in mantle temperature, fertility, or dynamic pressure could lead to the development of strong asymmetry in mantle flow. This expectation is borne out by new calculations presented here, and may provide a physical framework in which to better explain tomographic observations.

2. Previous Geodynamic Models of Active Flow at MORs

Previous models of active flow beneath mid-ocean ridges can be divided into two categories: single- and two-phase models. In single-phase models, equations derived from conservation principles are used to compute the creeping, incompressible flow of the solid mantle only. In some cases, single-phase models use a parametrization of magmatic transport to compute the rate of crust formation and the residual mantle porosity [e.g., McKenzie, 1985]. Two-phase models, in contrast, solve conservation equations for the interpenetrating flows of the creeping mantle and liquid magma. The magma is assumed to flow through the pores spaces between grains of the mantle. Both phases are independently incompressible, but within a volume of mantle containing liquid and solid, the liquid-filled pore network can exhale or inhale (by compacting or dilating) such that there is a negative or positive divergence of solid flow at the continuum scale. Both single- and two-phase models can compute a spatially variable bulk density based on the local porosity, temperature, and composition, however only two-phase models are consistent representations of the dynamic porosity in a natural magma/mantle system.

The first computational study of active convection at MORs, by Rabinowicz et al. [1984], used a 2-D, single-phase formulation. They imposed an aggregate density decrease in the melting region assuming 5% retained porosity; the melting region was defined with one of two recipes. They concluded that for mantle viscosity of 10¹⁸ Pa s, active upwelling makes a significant but not dominant contribution, whereas at 10¹⁷ Pa s, active convection dominates and upwelling is focused to a narrow zone. This conclusion is roughly consistent with all subsequent work on the subject. Sotin and Parmentier [1989] also used parameterized melt transport, but allowed for spatially variable porosity, and incorporated chemical depletion buoyancy into calculations. Parmentier and Phipps Morgan [1990] and Sparks and Parmentier [1993] computed a single-phase model in three dimensions, prescribed instantaneous and complete melt extraction, and considered compositional and thermal buoyancy, but not buoyancy from residual porosity. These studies used constant (or piecewise constant) mantle viscosity. Recent single-phase models have used a temperature-, pressure-, composition-, and porosity-dependent, non-Newtonian viscosity to model features such as off-axis volcanic chains [Ballmer et al., 2009].

Two-phase models are more challenging to compute, but are more realistic than single-phase models. They capture magmatic flow and predict the distribution of porosity based on a balance between melting and magmatic segregation. This allows the aggregate density to be computed consistently throughout the domain. Buck and Su [1989] introduced temperature- and porosity-dependent mantle viscosity into a 2-D model with Darcian melt segregation. They obtained a sharply localized upwelling beneath the ridge axis, and mantle recirculation at its flanks. Their choice of parameters yielded a porosity of ∼25% in the upwelling region. Scott and Stevenson [1989] took a similar approach, but applied a piecewise constant mantle viscosity. They explored a larger parameter space than did Buck and Su [1989], and presented a reduced analytical model of active flow. Spiegelman [1993] assumed constant mantle viscosity in order to decouple the incompressible and irrotational components of the mantle flow field. He obtained solutions for the irrotational component that, in contrast to previous work, incorporate viscous compaction stresses. He derived an inhomogeneous biharmonic equation for the incompressible part, and rescaled it to obtain a dimensionless number analogous to the Rayleigh number, but for buoyancy arising from residual porosity. Spiegelman [1993, 1996] presented end-member solutions for large and small values of this number. Choblet and Parmentier [2001] constructed a three-dimensional model of active flow with viscosity and bulk density both dependent on porosity and depletion. That study incorporated Darcian magmatic flow, but neglect compaction and associated stresses. Other authors have considered the two-phase dynamics of diapirs and plumes away from ridges [Tackley and Stevenson, 1993; Ghods and Arkani-Hamed, 2002; Schmeling, 2010].

All previous models of active flow at MORs have assumed symmetry of flow across the ridge axis. The authors performed computations on a domain with one boundary aligned with a vertical line (or plane) directly beneath the ridge axis, where a reflection boundary condition was applied. While their results show that significant active upwelling is expected for mantle viscosities below ∼10¹⁹ Pa s, they do not address the seismic observations of an asymmetrical distribution of mantle flow and porosity that are reviewed above.

Work presented here explores the development of asymmetry in active flow models of coupled magma/mantle dynamics with a domain that straddles the ridge axis. These models account for buoyancy arising from residual porosity only, but represent a fully consistent solution of the magma-dynamics equations of McKenzie [1984]. An earlier version of this model system was presented by Katz [2008]; that work incorporated the Enthalpy Method to transparently account for melting/freezing in a two-phase, two-component, energy-conserving system in thermodynamic equilibrium. The models presented here incorporate temperature- and porosity-dependent viscosity. This and other differences from Katz [2008] are described in 3. 9 extends the scaling argument of Spiegelman [1993] to obtain a dimensionless ratio analogous to the Rayleigh number. 10 presents new results on convective vigor and crustal production for comparison with past studies. It then explains the model response to asymmetrical forcing by lateral temperature and compositional gradients, as well as ridge migration. 17 discusses present results in the context of observations and previous work, and speculates on how additional model complexity would modify results.

3. Model Description

The models presented below are an extension of those described by Katz [2008]. They are the result of numerical solution of a set of coupled partial differential equations, representing conservation of mass, momentum, energy, and composition in a two-phase, two-component medium. The equations, presented in Appendix A, are a reformulated version of the set derived by McKenzie [1984], and subsequently by other authors with minor differences [e.g., Fowler, 1985; Ribe, 1985b; Scott and Stevenson, 1986; Bercovici et al., 2001]. Use of the Enthalpy Method [e.g., Alexiades and Solomon, 1993] as a means for computing the evolution of thermochemical variables in magma-dynamic simulations was introduced by Katz [2008] and is continued here. Numerical solutions are obtained using the Portable Extensible Toolkit for Scientific Computation; the discretization, solution strategy, and algorithms are outlined in 27.

There are several important differences between the current models and past work. In this manuscript I examine the effects of buoyancy forces on the flow of the matrix phase, whereas Katz [2008] considered only passively driven mantle flow. Mathematical details of the buoyancy terms in the governing equations can be found in Appendix A. To model active mantle convection in a physically reasonable way (i.e., respecting the rigidity of the lithosphere), I incorporate a temperature-dependent shear viscosity, representing deformation by diffusion creep [e.g., Karato and Wu, 1993]. The shear viscosity in the present models also depends on porosity [e.g., Hirth and Kohlstedt, 2003]. Models described by Katz [2008] used constant shear viscosity. Constitutive relations employed in the present models are detailed below. This is preceded by an explanation of the highly simplified petrological relations parameterized as a two-component phase diagram, and followed by a discussion of the means of melt extraction to the model crust.

Two-phase mid-ocean ridge models by Katz [2008] as well those by previous workers [Buck and Su, 1989; Scott and Stevenson, 1989; Spiegelman, 1996; Ghods and Arkani-Hamed, 2000] have assumed mirror symmetry across the vertical boundary beneath the ridge axis. Motivated by an interest in observed asymmetries in tomographic images [e.g., Forsyth et al., 1998b; Dunn and Forsyth, 2003], I have extended the simulation code described by Katz [2008] to optionally model both sides of the ridge axis, without any assumption regarding symmetry. Indeed, this allows for symmetry-breaking boundary conditions such as ridge migration and large-scale temperature gradients; these are explored in detail below.

3.1. Petrology

The petrological model is based on that employed by Katz [2008]: temperature-pressure-composition relations in an equilibrium, two-component system are parameterized with a phase diagram. The solidus and liquidus temperatures increase linearly with pressure according the Clapeyron slope. In the work by Katz [2008], the phase compositions were described with cubic polynomials of temperature; in the current models these relations are linearized about the ambient mantle concentration C₀, as shown in Figure 1b. The surfaces are defined by

C_S and C_L define the compositions of the solid and liquid, respectively, where they coexist in thermodynamic equilibrium; T₀ is the reference solidus temperature at zero pressure and ambient mantle composition C₀; ΔC is the reference composition difference between solidus and liquidus surfaces at temperature T₀; P_f is the pressure of the fluid; and γ = dP/dT∣_C is the Clapeyron slope (assumed equal for both phases). The change in temperature with composition at fixed pressure is given by M_S and M_L for the solidus and the liquidus, respectively.

Details are in the caption following the image — **Figure 1**
Open in figure viewer PowerPoint

Phase diagrams for the two-phase, two-component system. The black mesh represents the solidus temperature as a function of pressure and composition; the blue mesh represents the liquidus temperature. Solid red lines highlight the liquidus andsolidus at one particular value of pressure. Along lines of constant composition, the change in the liquidus and solidus surfaces is given by the Clausius-Clapeyron parameter, dP/dT∣_C = γ. (a) Solubility loop phase diagram as implemented by *Katz* [2008]. (b) Liquidus and solidus surfaces linearized about C₀. The slopes at constant pressure are given by M_L and M_S for the liquidus and solidus, respectively.

Katz [2008] chose C₀, the concentration of the less fusible phase in the ambient mantle, to be 12 percent by mass. In the current model I choose C₀ to be 85 percent; I then set the solidus temperature T₀ for composition C₀ and mantle potential temperature equation image _m so as to impose the bottom of the melting regime at a depth of ∼60 km. The different values of C₀ make no difference to the dynamics of the system, though the current choice is consistent with an asthenosphere that is composed mainly of olivine, which is less fusible than other mantle minerals.

It would be a mistake, however, to quantitatively map the compositional space of the current models to actual phase relations under mantle melting. The purpose of this simple petrological system is to provide a basis for identifying “enriched” (C smaller) and “depleted” (C larger) compositions, and to hence enable the model to capture a minimal aspect of chemical variability. Furthermore, the assumption of thermodynamic equilibrium made here is problematic in terms of fidelity to the natural mantle system. This assumption is justified by the aim of the current work to study the coupled fluid dynamics of the magma/mantle system, which requires a petrologically reasonable, energy conservative, internally consistent, and simple approach to mantle melting/freezing.

The solidus equation (1a) can be combined with the linearized expression for the adiabatic temperature gradient, T(z) ≈ equation image

_m (1 + αρz/c_P), to give the initial depth of melting. For a homogeneous mantle of composition C₀ and potential temperature equation image

_m, the result is

where ρ is the mean density of a the mantle column, α is the coefficient of thermal expansivity, and c_P is the specific heat capacity. Using γ = 60⁻¹ GPa/K, equation image

_m = 1648 K, T₀ = 1565 K, ρ = 3000 kg/m³, g = 10 m/s², α = 3 × 10⁻⁵ K⁻¹, and c_P = 1200 J/kg/K, I obtain z_m = 60 km.

3.2. Rheology

Experimental studies of mantle rheology indicate that it varies with parameters including temperature, pressure, strain rate, melt fraction, and water content [e.g., Hirth and Kohlstedt, 2003]. These dependencies themselves vary as a function of the mode of deformation. New evidence from analogue experiments [Takei, 2010] and homogenization theory [Takei and Holtzman, 2009a, 2009b, 2009c] suggest that the mantle is viscously anisotropic. While these complexities could have significant consequences for coupled magma/mantle dynamics, to keep models relatively simple here, I capture the only the leading-order variation in mantle viscosity, which is isotropic.

3.2.1. Shear Viscosity

The shear viscosity is taken to vary with temperature and porosity according to

where E* is the activation energy, R is the universal gas constant, T is the temperature in Kelvin, λ is a positive constant, and ϕ is the porosity. equation image

is a constant, set equal to the temperature at which the mantle first crosses the solidus, equation image

_m exp(αgz_m/c_P), with z_m from equation (2). Hence η₀ is the viscosity just below the melting region. For numerical efficiency, variation of viscosity in the present models is limited by an upper bound of 10²⁵ Pa s; this limit is sufficiently high to have a negligible effect on the mantle flow field. No lower bound is imposed.

Equation (3) can be thought of as a simplified parametrization of diffusion creep viscosity. The temperature dependence allows for the important effect of rigid lithospheric plates. By incorporating a porosity dependence, the emergence of localized porosity structures, as observed in experiments and theory [Kohlstedt and Holtzman, 2009], are not excluded a priori. As we shall see in 19, however, under the current rheological formulation, a mechanical porosity-banding instability is not predicted by models.

3.2.2. Bulk Viscosity

The bulk viscosity is taken to vary with temperature and porosity according to

where ζ₀ is a reference value of bulk viscosity. The Arrhenius factor controls the shear viscosity of the solid grains, and hence appears in the relation for bulk viscosity as well. The form of (4) is consistent with theoretical predictions [e.g., Batchelor, 1967; Schmeling, 2000; Bercovici et al., 2001; Hewitt and Fowler, 2008; Simpson et al., 2010a] and experimental constraints [Cooper, 1990] on the variation of bulk viscosity. These studies indicate that ζ₀ should be equal to η₀ times an order 1 constant.

Katz [2008] neglected the temperature dependence of viscosity and showed how the magnitude of the bulk viscosity is a control on sublithospheric melt focusing to the ridge axis [Sparks and Parmentier, 1991]; this is confirmed by current models with the rheology as given above. In the calculations presented below, unless otherwise indicated, I have taken ζ₀ = 10¹⁹ Pa s. This choice may be inconsistent with theoretical predictions [e.g., Simpson et al., 2010b], however it yields compaction lengths within the partially molten region that can be resolved numerically, and it ensures a high focusing efficiency by minimizing thermomagmatic erosion into the base of the lithosphere [see Katz, 2008, Figure 4]. Efficient focusing promotes the attainment of a steady state in numerical solutions; hence it facilitates investigation of the convective dynamics of the subridge mantle.

3.3. Melt Extraction at the Ridge

The present models also differ from those of Katz [2008] in the way that melt is extracted at the ridge axis. In past work, magma escaped through the upper boundary of the domain via a narrow, imposed “window” at the ridge axis. This approach, though simple, has various problems; in particular, it precludes the use of slower spreading rates and variable viscosity. In the present models, magma is extracted through a “dike.” In the computational domain, this is implemented as a boundary condition between grid cells, directly beneath the ridge axis. This internal boundary begins between the shallowest cells with nonzero porosity, and extends a fixed number of cells downward. The choice to impose the extent of the dike in terms of a number of grid cells rather than a physical length is based on technical rather than physical grounds. Three or four grid cells, independent of the grid spacing, are the minimum number required such that entry into the dike is not a rate-limiting factor for melt extraction.

The grid cells that are adjacent to the dike experience a fluid pressure gradient that is proportional to the difference between the lithostatic pressure outside the dike and the magma-static pressure within it

where

is a fraction between zero and unity, Δx/2 is the x distance between the center of a grid cell and its edge. The case of equation image

= 1 applies if the dike is continuously open all the way to the surface; assuming dynamic forces are negligible, this case provides an upper bound on the pressure gradient. Reducing equation image

leads to larger porosity and smaller magmatic flow speed in the immediate neighborhood of the dike. I use equation image

= 0.005 to obtain solutions reported below; this value is large enough to maintain low porosity beneath the ridge axis, but small enough to avoid a suction effect being felt beyond the immediate neighborhood of the dike.

Magma that enters the dike (by crossing the internal boundary between cells) is immediately extracted from the domain and pooled, along with the composition and enthalpy that it carries. The rate of extraction and the composition of the extracted melt is recorded at each time step.

4. Scaling Analysis

If gradients of buoyancy in the mantle beneath a mid-ocean ridges are of the same order as gradients in the viscous stresses arising from plate-induced mantle flow, then there will be a component of active mantle convection beneath the ridge axis [Rabinowicz et al., 1987]. Buoyancy-driven convection is typically studied in terms of the dimensionless Rayleigh number, which is the ratio of factors driving convection to those resisting it [e.g., Turner, 1973]. A parameter similar to the Rayleigh number was derived by Spiegelman [1993] for the mantle beneath a mid-ocean ridge, where buoyancy forces arise from the presence of melt. Here I extend that analysis in order to obtain the scaling of active upwelling velocity with problem parameters.

In the absence of melt, shearing flow of the crystalline mantle on tectonic length and time scales is governed by the Stokes equation. Appendix A reviews the derivation of the corresponding equation for the case of a two-phase system, with small amounts of melt present in the pore spaces between mantle grains. In terms of dimensional variables, this equation becomes

where P is a dynamic pressure, η is the shear viscosity of the two-phase aggregate, v_m is the mantle velocity, Δρ is the difference in density between the mantle and the magma, and g is the gravity vector. This equation states that pressure gradients are balanced by viscous stresses and a body force arising from buoyancy. Assuming the phase densities are constant in space and time, the buoyancy force on the two-phase aggregate is due solely to the porosity.

To simplify, I further assume that beneath the lithosphere, where mantle temperatures are approximately adiabatic, the mantle viscosity is constant. This approximation, applied in this section only, neglects the effects temperature and porosity, which are included in simulations, as well as lithostatic pressure and water content [e.g., Hirth and Kohlstedt, 2003]. While the former may be a good assumption in the shallow asthenosphere, the latter two are questionable; even small amounts of water in crystals has a large effect on their viscosity. Choblet and Parmentier [2001] found that the rheological effects of water on convection beneath ridges are significant. Neglecting viscosity variation due to porosity is less troublesome; simulations show that it does not lead to large variations within the melting region. There are other limitations associated with these approximations, and I do not expect this scaling analysis to be strictly accurate for the mantle. It will provide a basis, however, for understanding current simulations (that employ temperature- and porosity-dependent viscosity throughout the domain) and may give some insight into the physical regime relevant for the Earth.

For constant η = η₀, I take the curl of equation (6) to obtain, within the adiabatic asthenosphere,

where ω =

× v_m and ∂_x indicates a partial derivative in the x direction. This is the vorticity equation and states that lateral gradients in buoyancy are an internal source of vorticity. Vorticity is also generated by the boundary conditions, the spreading tectonic plates. To estimate the importance of buoyancy versus spreading-induced vorticity, following Spiegelman [1993] I rescale variables with characteristic scales:

The characteristic scale for distance is h₀, the height of the melting column. The characteristic scale for porosity is ϕ₀, which is a representative value for the column of rock beneath the ridge axis. The characteristic scale for vorticity is U₀/h₀, with U₀ the half-spreading rate. Applying these scales and dropping primes from dimensionless variables gives

where

is a dimensionless number analogous to the Rayleigh number in thermal convection. This result was obtained by Spiegelman [1993]; I extend it by constraining the characteristic porosity ϕ₀ using a result from Ribe [1985a]. That work considered an upwelling, melting column of mantle rock and neglected compaction stresses to obtain

where μ is the viscosity of the magma, F(z) is the degree of melting, W₀ represents the rate of mantle upwelling at the base of the column, and K₀ and n are associated with the permeability-porosity relation, K = K₀ϕⁿ. Since the contribution of buoyancy to upwelling is a priori unknown, I will take the upwelling rate W₀ equal to U₀ for convenience, and remain cognizant that this will be an underestimate (though the exponent 1/n reduces the impact of this error). Furthermore, since I seek an estimate of the characteristic porosity, I will take a representative value for degree of melting, F₀, equal to the degree of melting at the top of the column.

Incorporating these modifications, taking n = 3, and substituting equation (10) into equation (9) gives

The dimensionless number equation image

predicts how the relative importance of active upwelling scales with problem parameters. It is more informative than the version presented by Spiegelman [1993] because it expands the dependence of equation image

₀ to parameters that can be estimated from observations or experiments. When equation image

≪ 1, equation (11) predicts a passive mantle flow regime, whereas when &#55349;&#56497; ≫ 1, buoyancy-driven flow will dominate. As we shall see in the results 10, the vigor parameter can provide a quantitative prediction of upwelling rate relative to the spreading rate.

Inspection of equation (11) shows that the melting column height and mantle viscosity are the primary controls on the vigor parameter. The spreading rate appears explicitly as U₀^−2/3, indicating a decrease in the relative importance of active convection for faster spreading. However, for slow and ultraslow spreading centers, the column height h₀ decreases sharply with decreasing spreading rate, building in an implicit dependence on U₀.

5. Results

Figure 2 shows output at t = 1 Ma from a representative simulation with buoyancy-driven mantle flow. The calculation was initialized following the steps made by Katz [2008] and Ghods and Arkani-Hamed [2000]: (1) the steady state mantle flow, temperature, and porosity are calculated assuming passive upwelling and zero permeability; (2) the enthalpy field is reduced by the latent heat at each point [ equation image ← ( − ρLϕ)]; (3) the bulk composition is set equal to the mantle composition [C ← C_m]; and (4) the porosity is set to zero everywhere [ϕ ← 0]. At later times, the history of this initial condition is reflected in the sharp gradient in the degree of melting F, contoured in Figure 2a; this feature marks what was the edge of the melting region at t = 0.

A solution to the governing equations comprises mantle and magmatic flow fields, dynamic and compaction pressure, temperature, porosity, and the bulk, magma, and mantle compositions; some of these outputs are shown in Figure 2 (see caption for details). Melting rate can be obtained from solutions at two consecutive time steps as Γ = ∂ϕ/∂t − ∇ · (1 − ϕ)v_m and is shown in Figure 2b.

Streamlines in Figure 2b show that for a background shear viscosity of η₀ = 10¹⁸ Pa s and residual mantle porosity of ∼1%, mantle flow is not drastically different from the typical pattern of passive upwelling. The melting region appears roughly triangular, however its width is diminished relative to the case of passive upwelling only. The choice of ζ₀ = 10¹⁹ Pa s ensures that the focusing efficiency is high [Katz, 2008]; it is also responsible for the large values of compaction length shown in Figure 2d. The compaction length is given by

when ϕ ≪ 1 and ζ₀ ≳ 4η₀/3. For the parameters used in the calculation shown in Figure 2, ϕ ∼ 0.01 and T ∼ T₀, the compaction length is on the order of 10 km.

In 10–12, I describe results from suites of models, where each model is a variation of the one in Figure 2. 10 quantifies the relationship between the vigor parameter equation image and simulation results on the rate of mantle upwelling. 11 considers variation in crustal production with problem parameters. 12 considers the symmetry breaking by active convection beneath a mid-ocean ridge. In 17, these results are interpreted in terms of geophysical predictions and consistency with observations.

5.1. Convective Vigor

The contribution of active convection to upwelling beneath the ridge should be related to the importance of the buoyancy term in equation (6), hence it should be possible to predict the rate of upwelling using the vigor parameter equation image . This parameter was introduced in a simpler form by Spiegelman [1993], but he presented mantle flow for only two end-member cases, ≪ 1 and ≫ 1. A natural system may exist anywhere along this spectrum; I test for a quantitative relationship between &#55349;&#56497; and upwelling rate using a suite of numerical models.

Calculations of passive, variable viscosity mantle flow produce a relationship between half-spreading rate and maximum mantle upwelling rate given by W_max ≈ 1.4U₀. I therefore seek to quantify the contribution of active upwelling by considering the ratio W_max/(1.4U₀) for different values of equation image . This data from a suite of simulations is shown in Figure 3. Figure 3a demonstrates that there is a clear relationship between the vigor parameter and the contribution of buoyancy to upwelling. The colored points represent model results for half-spreading rates of 1.5, 3, 4, and 6 cm/yr, with viscosity parameter η₀ ranging from 2 × 10¹⁷ to 10²¹ Pa s, permeability parameter K₀ ranging from 10⁻⁸ to 10⁻⁶ m², and other parameter variations as given in Table B1. The y values are spread over more than an order of magnitude; by plotting against equation image , they collapse onto a line. Figure 3b shows that the data collapse is equally tight if plotted against the dimensionless number _W = &#55349;&#56497;(W₀/U₀)^1/n. This is a modified vigor parameter that includes the mantle upwelling velocity extracted from each simulation result (i.e., it does not assume that W₀ ≈ U₀).

The solid lines in Figure 3 are best fit regressions of the form W_max/(1.4U₀) = 1 + a equation image ^b, where a and b are fit parameters. Values for these parameters were obtained by nonlinear least squares error minimization (see caption of Figure 3). Assuming that the components of are independently constrained by observations, experiments, and thermodynamic calculations, then the lines of best fit in Figure 3 can be used to predict the contribution to mantle upwelling of the buoyancy from residual mantle porosity (explicitly for the fit from Figure 3a, implicitly for that of Figure 3b). More generally, however, these “empirical” relationships provide a means for assessing the dynamical consequences of hypotheses on mantle properties and processes. It is interesting to note that the relationship between W_max/(1.4U₀) and equation image _W is nearly linear in the regime where active upwelling is important (Figure 3b); this suggests that the assumptions made in deriving the scaling parameter are appropriate in the context of current simulations.

The data collapse obtained by plotting against 1 + a equation image ^b suggests that the parameter space of possible active upwelling scenarios is smaller than anticipated: instead of having dimensions of the number of ill-constrained model parameters, Figure 3 suggests that it is a one-dimensional space described by (or _W). Within the context of the physics that is captured by these models, this dimensionless number quantifies how parameter values can trade off. For example, the effect of doubling the mantle viscosity parameter η₀ is equivalent to the effect of a factor 2ⁿ change in permeability parameter K₀, where n is the exponent in the permeability-porosity relation. I use this and other trade-offs quantified by equation image to explore the behavioral space of the simulations without varying all of the parameters; in 11 and 12, only the parameter η₀ is varied between simulations, but the results can be approximately mapped to variations in the vigor of active convection.

Figure 4 shows the half-width of the melting regime from the same set of simulations that was used to construct Figure 3. The data are plotted against 1 − a equation image ^b, which approximates the ratio of the maximum upwelling rate to that expected for passive flow. For x → 1 (passive limit), the data trend toward different y values, as expected based on spreading rate control of lithospheric thermal structure. With increasing convective vigor this spread decreases, indicating that the width of the melting regime is here controlled by the horizontal length scale of active convection.

5.2. Crustal Generation

How does vigor of active convection affect the thickness of crust produced at the ridge axis? If we assume thermodynamic reversibility and neglect the details of melt segregation, the melt productivity dF/dz∣_S of adiabatically upwelling lherzolite is independent of upwelling rate [e.g., Asimow et al., 1997]; one therefore expects that the maximum degree of melting at the top of a melting column, equation image dF/dz∣_S dz, does not depend on upwelling rate. A faster upwelling rate should produce a proportionally larger melt production rate Γ beneath the ridge axis. Passive mantle flow, driven by divergence of tectonic plates at the ridge axis, produces a rate of upwelling that scales with the spreading rate. As described in 10, buoyancy forces arising from residual porosity can lead to upwelling at greater rates; if these are substantial relative to the rate of passive flow, buoyancy-driven flow will lead to enhanced melting, and thicker crust.

Figure 5 shows the results of a suite of models with different relative contributions of buoyancy-driven upwelling. Each line corresponds to a single value of η₀, and each point on the line corresponds to a different spreading rate. Figure 5a shows crustal thickness as a function of half-spreading rate; the data points are reproduced from White et al. [2001], and represent seismic measurements of crustal thickness from spreading ridges around the globe. The black curve in Figure 5a represents models with background viscosity η₀ = 8 × 10¹⁸ Pa s; in these calculations, the pattern of upwelling is almost indistinguishable from passive flow. With decreasing η₀ (increasing equation image ), subridge mantle upwelling rates increase, and hence crustal thickness increases at any value of half-spreading rate. This effect is pronounced for η₀ less than 2 × 10¹⁸ Pa s.

Figure 5b shows the mean degree of melting that produces the crust, based on the composition of the crust. This is calculated using F_crust = (C_crust + ΔC − C₀)/ΔC, where C_crust(t) is the composition of aggregated melts entering the dike at time t and ΔC is the reference compositional difference between the liquidus and solidus. This plot shows that the mean degree of melting is nearly independent of η₀, and hence nearly independent of the component of active upwelling. This is unsurprising, given that F is controlled by the column height of upwelling mantle at temperatures above the solidus [Langmuir et al., 1992]. At slow spreading rates, the current simulations indicate that active upwelling can modify the vertical balance of advection and diffusion of heat below the ridge axis. This causes a thinning of the thermal boundary layer, an increase in height of the melting column, and an increase in the maximum and the mean degree of melting. At U₀ = 0.5 cm/yr, the maximum difference in melting column height between two simulations from Figure 5 is 35 − 29 = 6 km. The simulations produce a mantle-upwelling melt productivity dF/dz that is nearly constant at 0.35 %/km; multiplying this by 6 km gives a difference in degree of melting of about 2 percent. Note, however, that Ito and Mahoney [2005] used a petrological model with variable adiabatic productivity and found that the mean degree of melting varies with the spatial distribution of melting rate, even for a fixed melting column.

A comparison between the data points and model curves in Figure 5a would seem to indicate that a nonnegligible contribution from active upwelling is required to fit the data. It is important to note, however, that the amplitude of these curves is sensitive to a range of ill-constrained model parameters. For example, an increase in mantle potential temperature would lead to thicker crust for the same spreading rate and mantle viscosity structure. A change of parameters in the two-component petrological model used here could have the same effect. Hence in the context of the present model, crustal thickness alone does not permit us to distinguish between a range of active upwelling scenarios. This topic is examined in more detail in 17.

The variation of crustal thickness with spreading rate at intermediate to fast rates may be useful for making this distinction. The array of data in Figure 5a seems to trend toward smaller values of crustal thickness with increasing spreading rate, above 2 cm/yr. A similar trend is seen in the red and pink curves, representing lower values of viscosity. Figure 5c indicates that in the models, this trend arises from the decreasing contribution of active upwelling relative to passive upwelling. This can be understood by considering that F_max is relatively constant over this range of increasing spreading rates, and hence equation image is decreasing. If the decreasing trend in the data is robust, it may provide evidence for a contribution of active upwelling.

Figure 6 shows that at a fixed spreading rate, crustal thickness increases linearly with the contribution of active upwelling. The slope of the lines in Figure 6 are smaller than might be expected, however. At W_max/(1.4U₀) = 2, the maximum upwelling rate is double its expected value for passive upwelling, and hence one might expect the crustal thickness to have increased by a factor of 2. This doubling does not occur because the distribution of upwelling and melting changes with increasing contribution of buoyancy: the melting region narrows (Figure 4) and upwelling is localized around the zone of maximum buoyancy. As shown in 12, this zone is not necessarily located directly beneath the ridge axis.

5.3. Symmetry Breaking by Convection

Plate spreading is clearly an essential aspect of the localization of magma genesis to the region beneath and around mid-ocean ridges. It provides a baseline mantle-upwelling flow that leads to melt production and residual porosity. As 10 and 11 show, in a model where residual porosity is a source of matrix buoyancy, plate spreading can initiate and maintain a component of buoyancy-driven mantle convection. In the extreme model scenario where the rate of buoyancy-driven upwelling is much larger than that caused by plate spreading, mantle convective flow is self-sustaining, and would continue even if plate spreading ceased. Furthermore, in this case, the upwelling/melting regime is a dynamic feature that is centered on the zone of greatest buoyancy, even if this is not directly beneath the ridge axis. Spatial variation of mantle temperature or composition might lead to movement of the melting regime away from its position beneath the axis.

Such an extreme scenario seems unlikely for typical MORs on the Earth, but as an end-member, it provides insight into the results presented in this section. It demonstrates that convection can lead to symmetry breaking with respect to the distribution of the melting regime relative to the ridge axis. In 10 and 11 of the current work, and in previous work [e.g., Scott and Stevenson, 1989], symmetry was forced on the problem by the choice of boundary conditions. In this section I generalize the model by considering a full-ridge domain in 2-D, allowing for a linear gradient in temperature or composition at the bottom (inflow) boundary, and for ridge migration.

5.3.1. Temperature Forcing

A representative simulation result is shown in Figure 7 for U₀ = 3 cm/yr and η₀ = 10¹⁸ Pa s. A horizontal gradient of potential temperature of 0.09 K/km is imposed at 120 km depth, on the bottom boundary. There is a temperature difference of only 3.0 K in Figure 7a from the left side (−75 km) to the right side (+40 km) of the melting region at 60 km depth. Evidently, under the conditions of this model, a small lateral temperature gradient is sufficient to cause a large asymmetry in the distribution of upwelling, melting and porosity. I quantify the extent of porosity asymmetry by defining an asymmetry parameter

In this equation, Ω represents the entire domain, and Ω₊ represents the right half of the domain, where x > 0. Equation (13) states that Ψ = 1 when all of the porosity is found to the right of the ridge axis, Ψ = 0 when porosity is distributed evenly across the axis, and Ψ = −1 when all the porosity is on the left side.

Figure 8 shows how Ψ varies with spreading rate, mantle viscosity, and the magnitude of temperature gradient forcing. In Figure 8a, Ψ is plotted against the viscosity parameter η₀; curves for different spreading rates are nearly superimposed on each other. Figure 8b shows Ψ plotted against the predicted ratio of maximum upwelling rate to that expected for passive spreading. In contrast to Figure 8a, lines representing different spreading rates can be clearly distinguished. Although the data is sparse, there is a roughly linear relation between the vigor of active convection and the degree of asymmetry. The slope of this line increases with spreading rate. Further consideration of this observation is given in 17.

5.3.2. Compositional Forcing

Mantle heterogeneity in the Earth is composed of spatial variation in both temperature and composition. Figure 9a shows how Ψ varies with mantle viscosity and the magnitude of composition gradient forcing dC_m/dx. This gradient is imposed on the domain's bottom boundary, so that inflowing mantle has composition C₀ at the center of the domain, and deviates from this with distance from that point according to xdC_m/dx. For reasons that are not well understood, these model calculations are more computationally challenging than the calculations for gradients in mantle temperature. The number of model runs, and the magnitude of the applied forcing, were therefore limited relative to Figure 8. It was also necessary to use a larger bulk viscosity (ζ₀ = 4 × 10¹⁹ Pa s) to obtain convergence in this case. Despite these differences, results shown in Figure 9 tell the same story as those of Figure 8: active convection beneath a mid-ocean ridge makes the upwelling and melting regime sensitive to modest asymmetry in forcing; for the same forcing, the melting regime produced under passive flow is less asymmetric.

Figure 9b shows that the scaling of asymmetry Ψ is roughly linear with 1 + a equation image ^b for a fixed compositional gradient forcing, similar to the case of temperature-driven asymmetry in Figure 8b.

5.3.3. Ridge Migration

Symmetry of sub-MOR mantle flow can also be broken by ridge migration. For a low-viscosity asthenosphere, it is reasonable to assume that mantle beneath the asthenosphere has zero horizontal velocity in the hot spot reference frame [e.g., Lenardic et al., 2006]. Taking the frame of reference of the migrating ridge axis, the mantle flow induced by ridge migration can then be calculated by imposing the ridge migration velocity as a boundary condition on the bottom of the domain, assumed to coincide with the base of the asthenosphere. Previous work has considered the effect of ridge migration on passive mantle flow and melting in two [Davis and Karsten, 1986; Schouten et al., 1987; Conder et al., 2002; Toomey et al., 2002; Katz et al., 2004] and three dimensions [Weatherley and Katz, 2010].

Figure 10 shows the development of asymmetry of the melting region with time for a ridge that is spreading with a half rate of 3 cm/yr and migrating to the left. Results from numerical experiments for a migration rate of U_m = 1.5, 3, 6 cm/yr and η₀ = 0.5, 1, 2, 4, 8 × 10¹⁸ Pa s are plotted and show two behavioral regimes. For a more viscous mantle, ridge migration leads to enhanced melting on the leading side of the ridge (“upwind”) and a corresponding negative value of Ψ. For η₀ = 0.5 × 10¹⁸ Pa s (red lines), ridge migration leads to positive (“downwind”) asymmetry in Figure 10. Model results shown by red lines in Figure 10 indicate that the low-viscosity runs did not reach a steady state after two million years of model time.

Figure 11 gives a physical picture of the perturbation to mantle upwelling for both asymmetrical regimes. Enhanced upwelling on the upwind side of the ridge, corresponding to Ψ < 0 and high η₀ (Figure 11a), is the result of the mantle shear that arises from ridge migration; this perturbation flow must conform to the curved bottom boundary of the lithosphere [e.g., Katz et al., 2004; Conrad et al., 2010]. Enhanced upwelling on the downwind side of the ridge, corresponding to Ψ > 0 and low η₀ (Figure 11b), is the result of the convection cell that is rooted in the mantle asthenosphere and is lagging behind the leftward migrating ridge. The transition between these two regimes occurs with decreasing η₀: active upwelling makes an increasing contribution, the melting region narrows and hence the kinematic effect of the lithospheric boundary on upwelling is diminished. At the same time, the dynamic effect of buoyancy is augmented and the ridge can move away laterally, upwind from the center of buoyancy-driven upwelling.

In both classes of asymmetry development, Figure 10 shows that more rapid ridge migration relative to the half-spreading rate leads to larger ∣Ψ∣. An assessment of global MOR data by Small and Danyushevsky [2003] indicates that the magnitude of ridge-perpendicular migration velocity clusters tightly around the half spreading rate, U₀. The magnitude of asymmetry also scales with the domain depth; work on passive flow, single-phase models demonstrates that ∣Ψ∣ ∝ h_A⁻¹, where h_A is the depth of the base of the asthenosphere [Weatherley and Katz, 2010]. The domain height (and hence asthenospheric depth) for the models presented in this section and 10 is 120 km, which is probably an underestimate by at least a factor of 2 of the actual asthenospheric thickness [e.g., Craig and McKenzie, 1986]. The magnitude of Ψ shown in Figure 10 is small relative to that arising from thermal or compositional gradients, and it would even smaller for larger h_A.

6. Discussion

In 3 (and in Appendix A) I presented a modified version of the two-phase, two-component model by Katz [2008] of coupled magma/mantle dynamics at mid-ocean ridges. The modifications incorporate temperature- and porosity-dependent mantle viscosity, a dike-type internal boundary for draining magma to the crust, and the capability to model buoyancy-driven mantle flow within a domain that straddles the ridge axis. A suite of half-ridge simulations performed with different model parameters was presented in 10; results from these simulations were shown to be consistent with the scaling relation presented in 9. 11 presented the dependence of crustal thickness and degree of melting on spreading rate and mantle viscosity. 12 considered asymmetrical mantle flow, melting, and porosity beneath the ridge axis. Aspects of these results are discussed in more detail below.

The vigor parameter has significant predictive power in the context of the assumptions used to develop the models presented here, in particular, the assumptions of a mantle buoyancy that derives solely from the residual porosity, and a fixed density difference between phases. Figure 3 shows that a simple function of equation image can be used to reasonably estimate the contribution of buoyancy forces to upwelling, over a range of half-spreading rates. It is not surprising, then, that the same function correlates with crustal thickness and susceptibility to melting zone asymmetry. Figures 6, 8, and 9, however, show that the correlation is different for different spreading rates. This may be due to the varying shape of the lithosphere-asthenosphere boundary as a function of spreading rate. At slower spreading rates, the thermal and rheological boundary descends steeply on either side of the ridge, and laterally confines the low-viscosity region where active convection occurs. This confinement inhibits off-axis migration of the upwelling center under asymmetric buoyancy forcing, and helps explain the fanning of lines of equal forcing magnitude in Figure 8b.

Nevertheless, 12 shows that the amount of asymmetry developed for a given forcing is sensitive to the contribution of active convection. The maximum magnitudes of temperature and composition gradient forcing chosen for the present study are modest. In the mantle, one might expect larger gradients sustained over shorter distances; such heterogeneity is not considered here. Sharp gradients in mantle temperature were considered by Conder et al. [2002] and Toomey et al. [2002], but those studies calculated passive mantle flow only. Accounting for buoyancy forces would certainly have modified the mantle flows they computed and would have produced increased asymmetry in the melting region for equal forcing. The key question then regards the relative contribution of buoyancy-driven upwelling to MOR melting regimes.

6.1. Observational Constraints on Active Flow

There are a variety of constraints that bear on the contribution of active flow, though none of them are conclusive. Global variation in oceanic crustal thickness is an example, and White et al. [2001] have produced a compilation of such data. Figure 5 shows that the data do not rule out a significant contribution from active convection. In fact, the downward trend of the data at U₀ > 2 may signify a diminishing (but still important) contribution of active upwelling. The curves computed for Figure 5a derive their amplitude from a variety of ill-constrained parameters in addition to buoyancy, and thus it may be impossible to distinguish between curves on the basis of a fit to the amplitude of data. For example, past studies [e.g., Sotin and Parmentier, 1989] have computed crustal thickness assuming complete melt extraction from the mantle to the ridge axis. If melt focusing is inefficient, for example, or if its efficiency depends on the style of mantle flow, then the computed curves may contain a systematic offset. The present models assume values for parameters that control the focusing efficiency (especially ζ₀) but do not prescribe melt focusing [Katz, 2008].

Unlike passive flow, however, active convection may yield 3-D variations in upwelling away from transform faults [Sparks and Parmentier, 1993; Jha et al., 1994; Barnouin-Jha et al., 1997; Choblet and Parmentier, 2001]. Observed along-axis variations in crustal thickness lead Dick et al. [2003] to suggest that buoyancy-driven convection is important beneath the ultraslow spreading Gakkel ridge.

The predicted width of the melting regime beneath a MOR varies with the contribution of buoyancy (Figure 4). Observations of this feature are indirect, however, and highly uncertain. The MELT experiment [Forsyth et al., 1998a] imaged a ∼ 500 km wide region of low seismic velocity that the authors attributed to the presence of mantle melt. They imaged this anomaly to extend to depths of ∼150 km, a large depth that may be indicative of the role of water in the melting regime [e.g., Asimow and Langmuir, 2003]. If a small amount of hydrous flux melting was occurring, it would explain the large width of the melting region without excluding a contribution of active upwelling near the axis. The off-axis distance of active seamounts on the Pacific plate, for example, might be another indicator of the width of the melting region. Scheirer et al. [1998] report that most of the recent, off-axis lava flows occur within 60 km of the axis, however on the west side of the ridge some fresh flows can be found more than 100 km away. In contrast to the MELT results, tomography of the Gulf of California [Wang et al., 2009] suggests melting regions that are approximately 100 km in width.

Geochemical observations may provide a means to constrain the contribution of buoyancy-driven flow to the melting regime. Figure 5b, however, shows that in the context of models reported here, active upwelling causes only a small change to the degree of melting represented by the crust, and does so only for slow to ultraslow spreading. Spiegelman [1996] showed that the distribution of trace elements transported by magma may provide a more powerful constraint. He computed two end-member, isoviscous models, one with equation image 1 (passive), and one with &#55349;&#56497; ≫ 1 (active). In the passive flow model, melt trajectories from a broad region converged to the ridge axis, delivering enriched melts there. This focusing of magmatic flow was produced by a dynamic pressure gradient [Spiegelman and McKenzie, 1987], not a sublithospheric decompaction channel [Sparks and Parmentier, 1991], although the resulting melt flow trajectories are similar. The active flow calculation by Spiegelman [1996] exhibited strong mantle recirculation, which injected previously depleted mantle back into the melting region (however this depleted rock was not allowed to melt). This recirculating flow regime also yielded laterally convergent streamlines of undepleted mantle just above the base of the melting region. Early melts of this undepleted source were highly enriched, while later melts, produced directly below the axis, were depleted. The magmatic streamlines of the early, enriched melts diverged away from the ridge axis, leaving only the depleted melts to arrive on axis. This combination of convergent solid flow and divergent melt flow resulted in the geochemical distinction between active and passive flow in the work of Spiegelman [1996] and Spiegelman and Reynolds [1999].

Active flow models in the present work exhibit neither convergent solid flow nor divergent melt flow in the melting region. Figure 12 shows mantle streamlines and the logarithm of porosity for a simulation with η₀ = 5 × 10¹⁷ Pa s ( equation image = 190, W_max/(1.4U₀) = 2.7). A domain size of 300 × 300 km was used for this model to minimize the influence of boundary conditions on flow in the melting region. It is evident that there is no mantle recirculation, and within the melting region there is no lateral convergence of solid streamlines. This difference from the active flow end-member of Spiegelman [1996] may be attributable to the larger domain size and the T-dependent viscosity that yields a rigid lithosphere. Magmatic streamlines (not shown in Figure 12) are focused to the ridge axis. One might therefore expect, in contrast to the prediction of Spiegelman [1996], that the trace element enrichment of on-axis lavas would not vary strongly with &#55349;&#56497;. Clearly this hypothesis requires further testing with calculations of trace element transport.

A broad range of different observations can be applied to constrain the components of the vigor parameter equation image . The mantle viscosity, for example, could be as large as 10²¹ Pa s, according to studies of postglacial rebound [Peltier, 1998], or as low as 5 × 10¹⁷, according to studies of postseismic relaxation [e.g., Panet et al., 2010], or it may be both, depending on the location and time scale of the relevant flow. Constraints on mantle permeability come from experiments and models [e.g., Wark and Watson, 1998; Faul, 2001; Zhu and Hirth, 2003], as well as from inverse calculations based on geochemical transport [e.g., McKenzie, 2000; Stracke et al., 2006] and magmatic flow [e.g., Maclennan et al., 2002], however these different approaches yield vastly different estimates for mantle permeability. The depth of the melting region, which exerts an important influence on the vigor parameter, depends on the water content of the mantle [e.g., Katz et al., 2003], which is thought to be about 125 ± 75 ppm [e.g., Hirth and Kohlstedt, 1996]. Clearly, putting tight bounds on the value of equation image is not yet possible.

Given the considerations above, it is clear that existing constraints on mantle flow do not preclude a significant contribution of buoyancy-driven upwelling at mid-ocean ridges. A trend of slowly decreasing crustal thickness for U₀ > 2 cm/yr, if it exists in the observational data, would indicate a significant role for active convection. Such a role would help to explain the asymmetry observed by seismic tomography at the MELT region of the EPR [e.g., Forsyth et al., 1998a] and beneath the Gulf of California [e.g., Wang et al., 2009], and would have other implications for the dynamics of a thermally and chemically heterogeneous asthenosphere.

6.2. Thermal and Compositional Density Variations

Simulations presented here consider density variations arising only from residual porosity; previous studies have also investigated the effects of density variations arising from depletion and temperature. Sotin and Parmentier [1989] considered the effect of melt depletion buoyancy beneath a ridge where the mantle surrounding the melting region was assumed to be undepleted. They concluded that depletion is a source of buoyancy driving active flow. Beneath a ridge, chemically buoyant, depleted mantle created in the melting region must spread laterally with time, however, and create a horizontal layer of depletion beneath the crust; this represents a stable stratification [Raddick et al., 2002]. In a homogeneous mantle, thermal buoyancy would be unlikely to substantially modify any convective upwelling near the ridge axis. Ito et al. [1996] showed that for the melting region within a mantle plume, the loss of thermal buoyancy due to transfer of energy to latent heat roughly balances the gain of depletion buoyancy due to melt extraction. Farther from the ridge axis, where the thermal boundary layer has thickened substantially, thermal buoyancy may be important for triggering convective rolls aligned with the spreading direction [e.g., Barnouin-Jha et al., 1997; Dumoulin et al., 2008; Ballmer et al., 2009].

In models of a heterogeneous mantle, thermal and depletion effects on buoyancy may have greater significance. For example, the calculated asymmetry in upwelling due to imposed gradients in mantle potential temperature and composition would certainly be modified. Thermal buoyancy would reinforce the asymmetric porosity distribution, while compositional buoyancy would suppress it. This suppression would occur because fertile, denser mantle melts deeper, and to a greater extent, than depleted, less dense mantle. Such issues will be reconsidered in future work.

6.3. On the Absence of Porosity/Shear Bands

Although it is not the main purpose of the present work, simulations described here offer an opportunity to investigate the relevance of a porosity-localizing mechanical instability beneath ridges. Experimental [e.g., Holtzman et al., 2003; King et al., 2010] and theoretical [e.g., Stevenson, 1989; Spiegelman, 2003; Katz et al., 2006] studies have shown that partially molten mantle materials subject to shear can develop bands of locally increased porosity and deformation. These rheologically weak bands tend to emerge at a particular angle to the shear plane under simple shear deformation. Katz et al. [2006] proposed that they might provide a mechanism for melt focusing to the ridge axis. Butler [2009] showed that simulations that account for buoyancy-driven magmatic segregation in addition to large-scale mantle shear can produce porosity bands. His work demonstrated that the orientation of bands can be modified by the action of buoyancy forces; he quantified the importance of buoyancy by introducing a dimensionless ratio B_u = (Δρgδ_c)/[(ζ + 4η/3) equation image ].

Localized bands of high porosity did not appear in the present simulations, despite inclusion of the porosity-weakening viscosity that is required to produce them. Most simulations did, of course, predict a high-porosity decompaction channel at the base of the lithosphere, but this feature is independent of the mechanical instability associated with porosity bands [Sparks and Parmentier, 1991]. Emergence of porosity bands in simulations of experimental conditions by Katz et al. [2006] were sensitive to the amplitude of the initial noise that was introduced; simulations described above contained no initial noise. To test the importance of initial noise for porosity band emergence in MOR models, I ran the two simulations shown in Figure 13. These two runs (for η₀ = 10¹⁹ Pa s and 10¹⁸ Pa s) started with an initial noise in the bulk composition of ±0.6% of C₀. They were run with enhanced porosity weakening: λ = 90 in equation (3), rather than the usual value of 27.

The compaction length in the center of the melting region of Figure 13 is ∼15 km, and so it is well resolved by the 1 km grid spacing. Various studies have shown that the spacing of porosity bands arising from mechanical instability is expected to be less than or equal to the compaction length [e.g., Spiegelman, 2003]. Numerical models performed by Katz et al. [2006] suggest that within this range, it is the wavenumber spectrum of initial noise that controls the spacing of bands. The simulations presented in Figure 13 use an initial noise field with spectral power at scales smaller than the compaction length, and hence should enable growth of porosity bands at that scale.

Variations in the porosity field in Figures 13a and 13d indicate the presence of compositional heterogeneity, but no shear bands are evident after two million years of simulated time. The second invariant of the strain rate tensor, shown in Figures 13b and 13e, is enhanced along the decompaction channel at the base of the lithosphere, but is otherwise smooth. Figures 13c and 13f show the point-wise value of log₁₀B_u. The size of this parameter within the melting region is small relative to the prediction of Butler [2009, Figure 1]. According to his analysis, at small values of B_u, buoyancy does not affect the orientation of shear bands.

Porosity band emergence and orientation was considered by Takei and Holtzman [2009c], who introduced an anisotropic viscosity tensor based on the homogenized contiguity between mantle grains. Their work predicted larger growth rates of instabilities than did Katz et al. [2006], and obtained band orientations consistent with experiments [Holtzman et al., 2003] without invoking non-Newtonian viscosity. It is possible that their proposed rheology and/or a non-Newtonian viscosity, and/or a much lower value of ζ₀ would give a different result in terms of mechanical instability and the emergence of porosity bands. A network of high-permeability melt channels, regardless of their origin, might change the melt retention properties of the melting region and hence affect the vigor of active convection, relative to diffuse porous flow.

6.4. Toward Models of 3-D Magma/Mantle Flow and Ridge Relocation

Much past work, and in particular that of Barnouin-Jha et al. [1997] and Choblet and Parmentier [2001], showed that active flow can lead to three-dimensional upwelling and melting structure beneath an idealized MOR without ridge offsets. The natural ridge system, of course, imposes a 3-D structure on asthenospheric flow due to the presence of transform faults [Dumoulin et al., 2008], overlapping spreading centers, and other forms of offset. These offsets are associated with asymmetries in ridge properties [Carbotte et al., 2004] that are indicative of complexity in the mantle flow and melting regime beneath. Present two-phase models have not been extended to three dimensions, though this is a goal for future work. Significant 3-D structure is expected based on past work with numerical models by other authors, and is consistent with seismic tomography [e.g., Wang et al., 2009]. Most models, however, have not explicitly incorporated magmatic flow. While this omission reduces computational cost and complexity, it precludes the models from addressing questions of 3-D magmatic transport near ridge offsets, and along ridge segments. A 3-D extension of present models will be a powerful tool for investigating these dynamics.

Mittelstaedt et al. [2008] investigated the interaction of a mantle plume with a MOR and showed that thermomagmatic erosion of the lithosphere by plume-derived magma can produce dynamic relocation of the ridge axis. Key to this process is the transport and deposition of latent heat by magma. Two-phase models can quantify this heat transport and model thermomagmatic erosion of the lithosphere. Combining dynamic ridge relocation [e.g., Gerya, 2010] with three dimensional models of coupled magma/mantle dynamics will enable a study of how magmatism affects the morphological evolution of the MOR system. In particular, it will be possible to investigate the influence of an asymmetric melting region on the location and kinematics of the ridge axis.

7. Summary

In the present article I introduced a modified version of the two-dimensional, two-phase, two-component model of magma/mantle interaction beneath mid-ocean ridges developed by Katz [2008]. The thermodynamics of the system are captured with the Enthalpy Method, and with a simple, binary phase diagram. The components are not actual chemical species, but instead represent “enriched” and “depleted” end-members in a one dimensional composition space. The fluid mechanics of the system are governed by a set of equations derived by McKenzie [1984] and reconsidered by many other authors.

The inclusion of constitutive laws controlling the temperature and porosity dependence of mantle rheology allowed simulations to capture the rigid oceanic lithosphere. The equations account for buoyancy forces that arise from the presence of a less dense magma in the pores of the mantle matrix.

I obtained numerical solutions to the governing equations in domains that are oriented vertically, perpendicular to the ridge axis. Some solutions were computed on half-ridge domains, while others were computed in domains that straddle the ridge axis. Solutions were analyzed to understand the characteristics and consequences of buoyancy-driven mantle flow on the mid-ocean ridge melting regime. In particular, motivated by observations of asymmetric seismic wave speed distributions in tomographic surveys of the subridge asthenosphere, I investigated the development of asymmetric upwelling and melting in the presence of buoyancy-driven mantle flow.

The key findings of this study are as follows:

1. The effect of buoyancy due to residual porosity on the flow and melting regime beneath MORs (i.e., the vigor of active convection) can be characterized by a simple function of a nondimensional parameter that is analogous to the Rayleigh number.

2. The width of the melting region is a decreasing function of convective vigor. Crustal thickness increases linearly with vigor, while the degree of melting inferred from the composition of the crust is independent of vigor.

3. The response of the system to an asymmetric forcing (e.g., ridge migration or a gradient in mantle potential temperature or composition) is amplified when buoyancy forces are important. Even under vigorous convection, however, ridge migration has a smaller effect than asymmetric thermal or compositional gradients.

4. Although the viscosity law used in simulations provides the necessary conditions for mechanical instability and emergence of porosity/shear bands, this instability was not observed to occur. Even tests with compositional noise and enhanced porosity weakening did not exhibit the emergence of porosity/shear bands.

The current models are based on simplifications that limit a quantitative mapping between results obtained here and observations of natural systems. First, convection is an inherently three dimensional process; 2-D models of convection cannot properly capture the scaling properties of the real system. Second, current models neglect thermal and compositional density variation. And finally, the symmetry-breaking thermal and compositional gradients are a highly idealized representation of mantle heterogeneity. Despite these limitations, the results presented here indicate that buoyancy-driven mantle flow increases the sensitivity of the MOR melting regime to the development of asymmetry, and may help to explain the unexpectedly asymmetric wave speed distributions in the subridge asthenosphere, as observed by seismic tomography.

Acknowledgments

I am grateful to G. Ito and B. Kaus for insightful reviews, M. Spiegelman and C. Langmuir for scientific discussions, R. S. White for providing crustal thickness data, Y. Sun at the Numerical Algorithms Group for technical support, and the Oxford Supercomputing Centre and the UK National Supercomputing Service for access to well-maintained clusters. This work was supported by an academic fellowship from the Research Councils UK and by NERC grant NE/H00081X/1.

Appendix A:: Governing Equations

The equations used in the present work to model the fluid mechanics of coupled magma/mantle flow are well established. Their first derivation was by McKenzie [1984] but other authors have derived equivalent or more general versions [Fowler, 1985; Ribe, 1985b; Scott and Stevenson, 1986; Spiegelman, 1993; Bercovici et al., 2001]. Katz et al. [2007] introduced a pressure decomposition and reformulated the equations in a manner that is well suited to numerical solution; this decomposition is used in the present work. McKenzie [1984] derived a statement of conservation of energy for the two-phase system and recent work has clarified, extended, and developed solutions for this equation [e.g., Šrámek et al., 2007; Hewitt and Fowler, 2008].

The equations and numerical approach used here are related to the work of all of these authors, but are most closely based on Katz [2008]. In fact, the differences between that work and the current modeling approach are small. They are (1) inclusion of buoyancy terms in the governing equations, (2) consideration of the temperature dependence of bulk and shear viscosity and the porosity dependence of shear viscosity, and (3) linearization of the petrological phase boundaries used as closure conditions for the Enthalpy Method. Like the model proposed by Katz [2008], the current simulations use a two-component thermodynamic system, assume thermodynamic equilibrium, and apply the Enthalpy Method to derive closure conditions.

A1. Conservation of Mass and Momentum

To expose these differences and for clarity of presentation, I review the full set of conservation equations and constitutive laws below. I begin with the conservation of mass and momentum equations derived by McKenzie [1984]:

where subscripts f and m represent the fluid and matrix phases, respectively. ρ is density; ϕ is the porosity or volume fraction of fluid phase; v is velocity; Γ is the mass transfer rate between phases (e.g., melting); P_f is the fluid pressure; K is the permeability; and μ, η, and ζ are the fluid, matrix shear, and matrix bulk viscosities, respectively; g is the gravity vector. For any variable or product of variables, Q, defined for each phase, equation image

= ϕQ_f + (1 − ϕ)Q_m. The first two equations represent conservation of mass for each phase, the second two equations represent force balance for each phase. For further details, see any of the references above.

A2. Phase Densities, Pressure Decomposition, and Reformulation

To account for variations in density due to temperature and composition, I introduce linearized expressions for the densities

where α and β are the coefficients of thermal and compositional expansion, respectively; T is the temperature; C is the concentration of the less dense and less fusible component; and the subscript i is replaced with either the fluid f or matrix m phase. The constants ρ_f0, ρ_m0, T₀, and C₀ are reference values. B can be thought of as the relative buoyancy of each phase. Below I make the Boussinesq approximation, taking ρ_f = ρ_m = ρ as a constant in all nonbuoyancy terms.

Following Katz et al. [2007], I decompose pressure into three parts

where P_l = ρ_m0gz is the reference lithostatic pressure, equation image

= (ζ − 2η/3)∇ · v_m is the compaction pressure, and P is the dynamic pressure.

Using equations (A2) and (A3) and applying the Boussinesq approximation, I reformulate the conservation equations (A1c) and (A1d) as

where ξ = (ζ − 2η/3) and Δρ = ρ_m0 − ρ_f0. The first of this set is the Compaction equation, governing volumetric deformation of the matrix; the second is the Stokes equation.

In addition to the model for density variations in (A2), the system of equations requires closure relations in the form of constitutive laws for permeability, shear viscosity, and bulk viscosity. These were presented in the main text, above, and are not repeated here.

A3. Conservation of Energy and Species Mass

Following [Katz, 2008], conservation of energy is cast as an equation for enthalpy equation image

where c_P is the specific heat capacity, equation image

is the potential temperature in Kelvin, L is the latent heat per kilogram, and k is the thermal conductivity (assumed equal between phases). This equation states that changes in enthalpy at a point are due to advection of sensible heat, advection of latent heat, and thermal diffusion. The conservation of bulk species mass used by Katz [2008] is

where

is the chemical diffusivity. This equation states that changes in bulk composition are due to advection by each phase and diffusion through the fluid phase.

Equations (A5) and (A6) contain variables equation image

, C_m, and C_f. Closure requires that I express the four variables as functions of the enthalpy ℋ and bulk composition C. Katz [2008] showed that the Enthalpy Method [Alexiades and Solomon, 1993] is suited for this task. Assuming thermal equilibrium in a two-phase, two-component system and neglecting variations in partial specific entropy, I write the following four equations:

The first of these is the definition of bulk enthalpy, and says that the enthalpy equals the sum of the latent and sensible energy; the second is the definition of bulk composition; the third and fourth are the solidus and liquidus equations, respectively, that specify the composition of the coexisting solid and liquid phase in equilibrium at a given pressure and temperature. These can be combined to give an equation for the porosity,

The solution of (A8) can be substituted back into (A7a) to obtain T, which can then be used to compute the phase compositions from (A7c) and (A7d). Katz [2008] generalized this method to two-phase systems with more than two components; I do not pursue that extension here.

A4. Nondimensionalization

To nondimensionalize the above equations, I employ the following variables

where ΔT = M_SΔC (the definitions of M_S and ΔC are given in 3). I also define the following characteristic scales

where H is a characteristic domain size, w₀ is a characteristic Darcian fluid velocity, K₀ is a proportionality constant for permeability (that should vary with grain size, dihedral angle, etc. but is held constant here for simplicity), and η₀ is a reference value for viscosity.

Using these scales, dropping primes and rearranging gives

where

= g/∣g∣ is a unit vector pointing in the direction of gravity and

are dimensionless buoyancy parameters. This set of equations contains nondimensional ratios δ, α*, β*, equation image

, Pe_T, and Pe_C; these are defined and explained in Table A1.

Table A1. A Summary of Nondimensional Parameters^a

Parameter	Typical Size	Comment
δ = ^1/2	3	Compaction length
α* =	7 × 10⁻³	Thermal expansivity
β* =	0	Compositional expansivity
	0.03	Adiabatic parameter
	8	Stefan number
Pe_T =	6 × 10⁷	Thermal Peclet number
Pe_C =	6 × 10⁹	Compositional Peclet number
G =	9 × 10⁸	Inverse Clapeyron slope
R_M =	1	Ratio of solidus to liquidus slope

a Typical sizes are calculated using a domain height H = 120 km. These parameters are for nondimensionalization of the problem and are not necessarily indicative of physical balances of terms in the governing equations.

The linearized equations for the solidus (1a) and liquidus (1b) are nondimensionalized according the same scheme and become

where P_f is the nondimensional pressure, G is a nondimensional inverse Clapeyron slope, and R_M is the ratio of solidus to liquidus slope (see Table A1).

Using (A12a) and (A12b) with nondimensionalized definitions of bulk enthalpy and composition, equation (A8) can be rewritten in nondimensional form as

which is a quadratic equation that can be solved for the porosity as a function of enthalpy and bulk composition, ϕ = ϕ( equation image

, Θ).

A5. Numerical Solution

The governing equations (A10) and constitutive laws are discretized using a finite volume scheme on a staggered mesh with 1 × 1 km cell size. At each time step, the full system of implicit, nonlinear algebraic equations is separated into two parts that are solved iteratively. The first part comprises equations governing the mechanics (A10a)–(A10c), while the second part is built from the equations governing thermochemistry, (A10d) and (A10e). Each block of equations is solved in the PETSc framework using an Incomplete LU preconditioned Newton-Krylov (Generalized Minimal Residual) method with an additive Schwartz scheme of domain decomposition. Details and references are given by Katz et al. [2007].

The only algorithmic difference from the approach presented by Katz et al. [2007] is that the current simulations are accelerated with staggered time stepping. Thermochemical variables equation image and Θ (along with ϕ, θ, Θ_m and Θ_f) are updated at each time step, Δt. Also, the compaction equation (A10b) is solved at each time step. The computationally expensive Stokes equation (A10a) and associated continuity equation (A10c) are solved at a larger time interval, however, because the tectonic-scale pattern of mantle flow evolves very slowly relative to the porosity or composition. The interval for Stokes updates is set at 5Δt for the models described in this article. Comparison between models runs with and without staggered time stepping confirms that the loss of accuracy is insignificant at a lag of 5Δt.

Appendix B:

Table B1 gives a full list of the simulations and their parameter values that are plotted in Figure 3. Parameter values not listed in Table B1 or its caption are given in Table B2.

Table B1. Model Runs, Parameters, and Results Used in Figure 3^a

Run ID	U₀ (cm/yr)	η₀ (10¹⁸ Pa s)	K₀ (10⁻⁷ m²)	Δρ (kg/m³)	f		_W	W_max/(1.4U₀)
mr5674	1.5	0.125	1	500	0.1	1062.02	3418.70	23.83
mr6020	1.5	0.2	1	500	0.1	635.98	1764.70	15.26
mr5659	1.5	0.25	1	500	0.1	505.88	1306.09	12.29
mr5675	1.5	0.35	1	500	0.1	358.70	805.40	8.09
mr5658	1.5	0.5	1	500	0.1	239.30	466.25	5.28
mr5655	1.5	1	1	500	0.1	114.06	173.10	2.50
mr5660	1.5	1	0.1	500	0.1	281.41	553.30	5.43
mr5656	1.5	2	1	500	0.1	56.88	73.15	1.52
mr5661	1.5	2	0.1	500	0.1	134.17	199.44	2.35
mr5657	1.5	4	1	500	0.1	28.39	34.19	1.25
mr5662	1.5	4	0.1	500	0.1	64.04	81.93	1.50
mr5663	1.5	8	0.1	500	0.1	31.95	38.37	1.24
mr5665	1.5	100	1	500	0.1	1.13	1.29	1.07
mr5681	1.5	1000	1	500	0.1	0.11	0.13	1.04
mr5643	3	0.25	10	500	0.1	186.07	312.58	3.39
mr6021	3	0.3	1	500	0.1	334.79	708.67	6.77
mr5679	3	0.35	1	500	0.1	286.22	572.46	5.71
mr4741	3	0.5	1	500	0.1	200.13	351.73	3.88
mr5644	3	0.5	10	500	0.1	92.97	129.02	1.91
mr4740	3	1	1	500	0.1	99.84	139.08	1.93
mr5641	3	1	0.1	500	0.1	233.12	413.17	3.98
mr5645	3	1	10	500	0.1	44.68	54.08	1.27
mr5683	3	1	1	500	0.01	103.87	145.07	1.95
mr5684	3	1	1	500	0.1	103.60	138.86	1.72
mr5685	3	1	1	250	0.1	62.87	78.87	1.41
mr5686	3	1	1	1000	0.1	158.61	255.20	2.97
mr4739	3	2	1	500	0.1	47.93	57.97	1.26
mr5640	3	2	0.1	500	0.1	112.11	153.63	1.84
mr5636	3	4	1	500	0.1	23.91	27.49	1.09
mr5639	3	4	0.1	500	0.1	53.86	64.99	1.25
mr5637	3	8	1	500	0.1	11.91	13.50	1.04
mr5638	3	8	0.1	500	0.1	26.86	30.74	1.07
mr5666	3	100	1	500	0.1	0.95	1.07	1.02
mr5668	3	1000	1	500	0.1	0.10	0.11	1.00
mr3370	4	1	0.0625	500	0.1	236.17	422.20	4.08
mr3375	4	1	4	500	0.1	56.50	68.73	1.29
mr3491	4	1	1	500	0.1	89.73	118.07	1.63
mr3492	4	1	0.25	500	0.1	147.92	219.74	2.34
mr3494	4	2	1	500	0.1	44.86	52.73	1.16
mr3495	4	2	0.25	500	0.1	73.87	91.66	1.36
mr3496	4	2	0.0625	500	0.1	117.86	160.95	1.82
mr3497	4	2	0.015625	500	0.1	197.79	363.13	4.42
mr5647	6	0.5	1	500	0.1	143.86	225.67	2.76
mr5646	6	1	1	500	0.1	71.82	91.84	1.49
mr5650	6	1	0.1	500	0.1	160.89	262.95	3.12
mr5649	6	2	0.1	500	0.1	80.36	103.38	1.52
mr5651	6	4	0.1	500	0.1	40.17	45.98	1.07
mr5653	6	8	0.1	500	0.1	20.08	22.77	1.04
mr5654	6	16	0.1	500	0.1	10.04	11.40	1.04
mr5667	6	100	1	500	0.1	0.74	0.83	0.98

a The bulk viscosity parameter ζ₀ = 4 × 10¹⁹ Pa s for all of these simulations. This value was chosen to avoid magmatic pooling along the solidus boundary, which leads to significant slowdown of the simulation code.

Table B2. Dimensional Parameters Used in the Model and Their Preferred Values

Quantity	Symbol	Definition	Range Considered	Preferred Value	Units
Half-spreading rate	U₀		0.5–7	3	cm/yr
Ridge migration rate	U_m		0.5–7	U₀	cm/yr
Permeability const.	K₀	K = K₀ϕⁿ	10⁻⁸–10⁻⁶	10⁻⁷	m²
Permeability exponent	n	ibid		3
Mantle shear visc. const.	η₀	equation (3)	5 × 10¹⁷–8 × 10¹⁹	10¹⁸	Pa s
Reference viscous temp.		equation (3)		see note^a	°C
Activation energy	E*	equation (3)		3 × 10⁵	J mol⁻¹
Viscosity constant	λ	equation (3)		27
Mantle bulk visc. const.	ζ₀	equation (4)	10¹⁸–4 × 10¹⁹	2 × 10¹⁹	Pa s
Magma viscosity	μ			1	Pa s
Boussinesq density	ρ			3000	kg m⁻³
Boussinesq density diff.	Δρ			500	kg m⁻³
Spec. heat capacity	c_P			1200	J kg⁻¹ K⁻¹
Thermal diffusivity	κ			10⁻⁶	m² s⁻¹
Coef. thermal exp.^b	α			3 × 10⁻⁵	K⁻¹
Latent heat	L			4 × 10⁵	J kg⁻¹
Chemical diffusivity				10⁻⁸	m² s⁻¹
Potential temperature	_m			1375	°C
Reference melting temp.	T₀	equation (1)		1292	°C
Clapeyron slope	γ	equation (1)		60⁻¹	GPa K⁻¹
Composition diff.	ΔC	equation (1)		0.1	wt frac.
Solidus, liquidus slope	M_S, M_L	equation (1)		400	K (wt frac.)⁻¹
Dike underpressure coef.		equation (5)	0.001–0.1	0.005

a The value of this parameter is chosen such that η = η₀ for potential temperature _m, composition C₀, and zero porosity, at the depth given by equation (2).
b The coefficient of thermal expansion is used to calculate the adiabatic temperature gradient but is not used in calculating variations in density for buoyancy terms.

References

Alexiades, V., and A. D. Solomon (1993), Mathematical Modeling of Melting and Freezing Processes, Hemisphere, Washington, D. C.
10.1115/1.2930032
Web of Science®Google Scholar
Asimow, P. D., and C. H. Langmuir (2003), The importance of water to oceanic mantle melting regimes, Nature, 421(6925), 815–820, doi:10.1038/nature01429.
10.1038/nature01429
CASADSPubMedWeb of Science®Google Scholar
Asimow, P. D., M. M. Hirschmann, and E. M. Stolper (1997), An analysis of variations in isentropic melt productivity, Philos. Trans. R. Soc. London A, 355(1723), 255–281.
10.1098/rsta.1997.0009
ADSWeb of Science®Google Scholar
Ballmer, M. D., J. van Hunen, G. Ito, T. A. Bianco, and P. J. Tackley (2009), Intraplate volcanism with complex age-distance patterns: A case for small-scale sublithospheric convection, Geochem. Geophys. Geosyst., 10, Q06015, doi:10.1029/2009GC002386.
10.1029/2009GC002386
ADSWeb of Science®Google Scholar
Barnouin-Jha, K., E. M. Parmentier, and D. W. Sparks (1997), Buoyant mantle upwelling and crustal production at oceanic spreading centers: On-axis segmentation and off-axis melting, J. Geophys. Res., 102(B6), 11,979–11,989.
10.1029/96JB03807
ADSWeb of Science®Google Scholar
Batchelor, G. K. (1967), An Introduction to Fluid Mechanics, Cambridge Univ. Press, Cambridge, U. K.

Google Scholar
Bercovici, D., Y. Ricard, and G. Schubert (2001), A two-phase model for compaction and damage: 1. General theory, J. Geophys. Res., 106, 8887–8906.
10.1029/2000JB900430
ADSWeb of Science®Google Scholar
Blackman, D. K., and J. M. Kendall (1997), Sensitivity of teleseismic body waves to mineral texture and melt in the mantle beneath a mid-ocean ridge, Philos. Trans. R. Soc. London A, 355(1723), 217–231.
10.1098/rsta.1997.0007
ADSWeb of Science®Google Scholar
Bown, J. W., and R. S. White (1994), Variation with spreading rate of oceanic crustal thickness and geochemistry, Earth Planet. Sci. Lett., 121, 435–449.
10.1016/0012-821X(94)90082-5
CASADSWeb of Science®Google Scholar
Buck, W. R., and W. S. Su (1989), Focused mantle upwelling below mid-ocean ridges due to feedback between viscosity and melting, Geophys. Res. Lett., 16(7), 641–644.
10.1029/GL016i007p00641
ADSWeb of Science®Google Scholar
Butler, S. L. (2009), The effects of buoyancy on shear-induced melt bands in a compacting porous medium, Phys. Earth Planet. Inter., 173(1–2), 51–59, doi:10.1016/j.pepi.2008.10.022.
10.1016/j.pepi.2008.10.022
ADSWeb of Science®Google Scholar
Carbotte, S. M., C. Small, and K. Donnelly (2004), The influence of ridge migration on the magmatic segmentation of mid-ocean ridges, Nature, 429, 743–746.
10.1038/nature02652
CASADSPubMedWeb of Science®Google Scholar
Choblet, G., and E. M. Parmentier (2001), Mantle upwelling and melting beneath slow spreading centers: Effects of variable rheology and melt productivity, Earth Planet. Sci. Lett., 184, 589–604.
10.1016/S0012-821X(00)00330-7
CASADSWeb of Science®Google Scholar
Conder, J. A., D. W. Forsyth, and E. M. Parmentier (2002), Asthenospheric flow and asymmetry of the East Pacific Rise, MELT area, J. Geophys. Res., 107(B12), 2344, doi:10.1029/2001JB000807.
10.1029/2001JB000807
ADSWeb of Science®Google Scholar
Conrad, C. P., B. Wu, E. I. Smith, T. A. Bianco, and A. Tibbetts (2010), Shear-driven upwrelling induced by lateral viscosity variations and asthenospheric shear: A mechanism for intraplate volcanism, Phys. Earth Planet. Inter., 178(3–4), 162–175, doi:10.1016/j.pepi.2009.10.001.
10.1016/j.pepi.2009.10.001
ADSWeb of Science®Google Scholar
Cooper, R. F. (1990), Differential stress-induced melt migration: An experimental approach, J. Geophys. Res., 95(B5), 6979–6992.
10.1029/JB095iB05p06979
ADSWeb of Science®Google Scholar
Craig, C. H., and D. McKenzie (1986), The existence of a thin low-viscosity layer beneath the lithosphere, Earth Planet. Sci. Lett., 78(4), 420–426.
10.1016/0012-821X(86)90008-7
ADSWeb of Science®Google Scholar
Davis, E. E., and J. L. Karsten (1986), On the cause of the asymmetric distribution of seamounts about the Juan De Fuca ridge—Ridge-crest migration over a heterogeneous asthenosphere, Earth Planet. Sci. Lett., 79(3–4), 385–396.
10.1016/0012-821X(86)90194-9
ADSWeb of Science®Google Scholar
Dick, H. J. B., J. Lin, and H. Schouten (2003), An ultraslow-spreading class of ocean ridge, Nature, 426(6965), 405–412, doi:10.1038/nature02128.
10.1038/nature02128
CASADSPubMedWeb of Science®Google Scholar
Dumoulin, C., G. Choblet, and M. P. Doin (2008), Convective interactions between oceanic lithosphere and asthenosphere: Influence of a transform fault, Earth Planet. Sci. Lett., 274(3–4), 301–309, doi:10.1016/j.epsl.2008.07.017.
10.1016/j.epsl.2008.07.017
CASADSWeb of Science®Google Scholar
Dunn, R. A., and D. W. Forsyth (2003), Imaging the transition between the region of mantle melt generation and the crustal magma chamber beneath the southern East Pacific Rise with short-period Love waves, J. Geophys. Res., 108(B7), 2352, doi:10.1029/2002JB002217.
10.1029/2002JB002217
ADSPubMedWeb of Science®Google Scholar
Evans, R. L., P. Tarits, A. D. Chave, A. White, G. Heinson, J. H. Filloux, H. Toh, N. Seama, H. Utada, J. R. Booker, and M. J. Unsworth (1999), Asymmetric electrical structure in the mantle beneath the East Pacific Rise at 17°S, Science, 286(5440), 752–756.
10.1126/science.286.5440.752
CASPubMedWeb of Science®Google Scholar
Faccenna, C., and T. W. Becker (2010), Shaping mobile belts by small-scale convection, Nature, 465(7298), 602–605, doi:10.1038/nature09064.
10.1038/nature09064
CASADSPubMedWeb of Science®Google Scholar
Faul, U. H. (2001), Melt retention and segregation beneath mid-ocean ridges, Nature, 410, 920–923.
10.1038/35073556
CASADSPubMedWeb of Science®Google Scholar
Forsyth, D. W., et al. (1998a), Imaging the deep seismic structure beneath a mid-ocean ridge: The melt experiment, Science, 280, 1215–1218.
10.1126/science.280.5367.1215
ADSPubMedWeb of Science®Google Scholar
Forsyth, D. W., S. C. Webb, L. M. Dorman, and Y. Shen (1998b), Phase velocities of Rayleigh waves in the melt experiment on the East Pacific Rise, Science, 280(5367), 1235–1238.
10.1126/science.280.5367.1235
CASADSPubMedWeb of Science®Google Scholar
Fowler, A. (1985), A mathematical model of magma transport in the asthenosphere, Geophys. Astrophys. Fluid Dyn., 33, 63–96.
10.1080/03091928508245423
ADSWeb of Science®Google Scholar
Gerya, T. (2010), Dynamical instability produces transform faults at mid-ocean ridges, Science, 329(5995), 1047–1050, doi:10.1126/science.1191349.
10.1126/science.1191349
CASADSPubMedWeb of Science®Google Scholar
Ghods, A., and J. Arkani-Hamed (2000), Melt migration beneath mid-ocean ridges, Geophys. J. Int., 140, 687–697.
10.1046/j.1365-246X.2000.00032.x
ADSWeb of Science®Google Scholar
Ghods, A., and J. Arkani-Hamed (2002), Effects of melt migration on the dynamics and melt generation of diapirs ascending through asthenosphere, J. Geophys. Res., 107(B1), 2026, doi:10.1029/2000JB000070.
10.1029/2000JB000070
CASWeb of Science®Google Scholar
Hager, B. H., and R. J. O'Connell (1981), A simple global-model of plate dynamics and mantle convection, J. Geophys. Res., 86(B6), 4843–4867.
10.1029/JB086iB06p04843
ADSWeb of Science®Google Scholar
Hammond, W. C., and D. R. Toomey (2003), Seismic velocity anisotropy and heterogeneity beneath the Mantle Electromagnetic and Tomography Experiment (MELT) region of the East Pacific Rise from analysis of P and S body waves, J. Geophys. Res., 108(B4), 2176, doi:10.1029/2002JB001789.
10.1029/2002JB001789
ADSPubMedWeb of Science®Google Scholar
Hernlund, J. W., D. J. Stevenson, and P. J. Tackley (2008a), Buoyant melting instabilities beneath extending lithosphere: 2. Linear analysis, J. Geophys. Res., 113, B04406, doi:10.1029/2006JB004863.
10.1029/2006JB004863
Google Scholar
Hernlund, J. W., P. J. Tackley, and D. J. Stevenson (2008b), Buoyant melting instabilities beneath extending lithosphere: 1. Numerical models, J. Geophys. Res., 113, B04405, doi:10.1029/2006JB004862.
10.1029/2006JB004862
Google Scholar
Hewitt, I. J., and A. C. Fowler (2008), Partial melting in an upwelling mantle column, Philos. Trans. R. Soc. A, 464, 2467–2491, doi:10.1098/rspa.2008.0045.
10.1098/rspa.2008.0045
CASWeb of Science®Google Scholar
Hirth, G., and D. L. Kohlstedt (1996), Water in the oceanic upper mantle: Implications for rheology, melt extraction and the evolution of the lithosphere, Earth Planet. Sci. Lett., 144(1–2), 93–108.
10.1016/0012-821X(96)00154-9
CASADSWeb of Science®Google Scholar
Hirth, G., and D. Kohlstedt (2003), Rheology of the upper mantle and the mantle wedge: A view from the experimentalists, in Inside the Subduction Factory, Geophys. Monogr. Ser., vol. 138, edited by J. Eiler, pp. 83–105, AGU, Washington, D. C.
10.1029/138GM06
Google Scholar
Holtzman, B. K., N. J. Groebner, M. E. Zimmerman, S. B. Ginsberg, and D. L. Kohlstedt (2003), Stress-driven melt segregation in partially molten rocks, Geochem. Geophys. Geosyst., 4(5), 8607, doi:10.1029/2001GC000258.
10.1029/2001GC000258
ADSWeb of Science®Google Scholar
Ito, G., and R. A. Dunn (2009), Mid-ocean ridges: Mantle convection and formation of the lithosphere, in Encyclopedia of Ocean Sciences, edited by J. H. Steele, K. K. Turekian, and S. A. Thorpe, pp. 867–880, doi:10.1016/B978-012374473-9.00654-8, Academic, Oxford, U. K.
10.1016/B978-012374473-9.00654-8
Web of Science®Google Scholar
Ito, G., and J. J. Mahoney (2005), Flow and melting of a heterogeneous mantle: 1. Method and importance to the geochemistry of ocean island and mid-ocean ridge basalts, Earth Planet. Sci. Lett., 230(1–2), 29–46, doi:10.1016/j.epsl.2004.10.035.
10.1016/j.epsl.2004.10.035
CASADSWeb of Science®Google Scholar
Ito, G., J. Lin, and C. W. Gable (1996), Dynamics of mantle flow and melting at a ridge-centered hotspot: Iceland and the Mid-Atlantic Ridge, Earth Planet. Sci. Lett., 144(1–2), 53–74.
10.1016/0012-821X(96)00151-3
CASADSWeb of Science®Google Scholar
Jha, K., E. M. Parmentier, and J. P. Morgan (1994), The role of mantle-depletion and melt-retention buoyancy in spreading-center segmentation, Earth Planet. Sci. Lett., 125(1–4), 221–234.
10.1016/0012-821X(94)90217-8
ADSWeb of Science®Google Scholar
Karato, S., and P. Wu (1993), Rheology of the upper mantle—A synthesis, Science, 260, 771–778.
10.1126/science.260.5109.771
CASADSPubMedWeb of Science®Google Scholar
Katz, R. F. (2008), Magma dynamics with the Enthalpy Method: Benchmark solutions and magmatic focusing at mid-ocean ridges, J. Petrol., 49, 2099–2121, doi:10.1093/petrology/egn058.
10.1093/petrology/egn058
CASADSWeb of Science®Google Scholar
Katz, R. F., M. Spiegelman, and C. H. Langmuir (2003), A new parameterization of hydrous mantle melting, Geochem. Geophys. Geosyst., 4(9), 1073, doi:10.1029/2002GC000433.
10.1029/2002GC000433
CASADSWeb of Science®Google Scholar
Katz, R. F., M. Spiegelman, and S. M. Carbotte (2004), Ridge migration, asthenospheric flow and the origin of magmatic segmentation in the global mid-ocean ridge system, Geophys. Res. Lett., 31, L15605, doi:10.1029/2004GL020388.
10.1029/2004GL020388
ADSWeb of Science®Google Scholar
Katz, R. F., M. Spiegelman, and B. Holtzman (2006), The dynamics of melt and shear localization in partially molten aggregates, Nature, 442, 676–679, doi:10.1038/nature05039.
10.1038/nature05039
CASADSPubMedWeb of Science®Google Scholar
Katz, R. F., M. G. Knepley, B. Smith, M. Spiegelman, and E. T. Coon (2007), Numerical simulation of geodynamic processes with the Portable Extensible Toolkit for Scientific Computation, Phys. Earth Planet. Inter., 163, 52–68, doi:10.1016/j.pepi.2007.04.016.
10.1016/j.pepi.2007.04.016
ADSWeb of Science®Google Scholar
King, D. S. H., M. E. Zimmerman, and D. L. Kohlstedt (2010), Stress-driven melt segregation in partially molten olivine-rich rocks deformed in torsion, J. Petrol., 51, 21–42, doi:10.1093/petrology/egp062.
10.1093/petrology/egp062
CASADSWeb of Science®Google Scholar
Kohlstedt, D. L., and B. K. Holtzman (2009), Shearing melt out of the Earth: An experimentalist's perspective on the influence of deformation on melt extraction, Annu. Rev. Earth Planet. Sci., 37, 561–593, doi:10.1146/annurev.earth.031208.100104.
10.1146/annurev.earth.031208.100104
CASADSWeb of Science®Google Scholar
Lachenbruch, A. H. (1976), Dynamics of a passive spreading center, J. Geophys. Res., 81(11), 1883–1902.
10.1029/JB081i011p01883
ADSWeb of Science®Google Scholar
Langmuir, C. H., E. Klein, and T. Plank (1992), Petrological systematics of mid-oceanic ridge basalts: Constraints on melt generation beneath ocean ridges, in Mantle Flow and Melt Generation at Mid-ocean Ridges, Geophys. Monogr. Ser., vol. 71, edited by J. Phipps Morgan, D. K. Blackman, and J. M. Sinton, pp. 183–280, AGU, Washington, D. C.
10.1029/GM071p0183
Web of Science®Google Scholar
Lenardic, A., M. A. Richards, and F. H. Busse (2006), Depth-dependent rheology and the horizontal length scale of mantle convection, J. Geophys. Res., 111, B07404, doi:10.1029/2005JB003639.
10.1029/2005JB003639
ADSPubMedWeb of Science®Google Scholar
Maclennan, J., M. Jull, D. McKenzie, L. Slater, and K. Grönvold (2002), The link between volcanism and deglaciation in Iceland, Geochem. Geophys. Geosyst., 3(11), 1062, doi:10.1029/2001GC000282.
10.1029/2001GC000282
CASADSWeb of Science®Google Scholar
McKenzie, D. (1969), Speculations on the consequences and causes of plate motion, Geophys. J. R. Astron. Soc., 18(1), 1–32, doi:10.1111/j.1365-246X.1969.tb00259.x.
10.1111/j.1365-246X.1969.tb00259.x
ADSWeb of Science®Google Scholar
McKenzie, D. (1984), The generation and compaction of partially molten rock, J. Petrol., 25(3), 713–765.
10.1093/petrology/25.3.713
CASADSWeb of Science®Google Scholar
McKenzie, D. (1985), ²³⁰Th–²³⁸U disequilibrium and the melting processes beneath ridge axes, Earth Planet. Sci. Lett., 72(2–3), 149–157, doi:10.1016/0012-821X(85)90001-9.
10.1016/0012‐821X(85)90001‐9
CASADSWeb of Science®Google Scholar
McKenzie, D. (2000), Constraints on melt generation and transport from U-series activity ratios, Chem. Geol., 162(2), 81–94.
10.1016/S0009-2541(99)00126-6
CASADSWeb of Science®Google Scholar
McKenzie, D., and M. J. Bickle (1988), The volume and composition of melt generated by extension of the lithosphere, J. Petrol., 29(3), 625–679.
10.1093/petrology/29.3.625
CASADSWeb of Science®Google Scholar
Mittelstaedt, E., G. Ito, and M. D. Behn (2008), Mid-ocean ridge jumps associated with hotspot magmatism, Earth Planet. Sci. Lett., 266(3–4), 256–270, doi:10.1016/j.epsl.2007.10.055.
10.1016/j.epsl.2007.10.055
CASADSWeb of Science®Google Scholar
Nicolas, A., and J. F. Violette (1982), Mantle flow at oceanic spreading centers—Models derived from ophiolites, Tectonophysics, 81(3–4), 319–339.
10.1016/0040-1951(82)90136-6
ADSWeb of Science®Google Scholar
Panet, I., F. Pollitz, V. Mikhailov, M. Diament, P. Banerjee, and K. Grijalva (2010), Upper mantle rheology from GRACE and GPS postseismic deformation after the 2004 Sumatra-Andaman earthquake, Geochem. Geophys. Geosyst., 11, Q06008, doi:10.1029/2009GC002905.
10.1029/2009GC002905
ADSWeb of Science®Google Scholar
Parmentier, E. M., and J. Phipps Morgan (1990), Spreading rate dependence of 3-dimensional structure in oceanic spreading centers, Nature, 348(6299), 325–328.
10.1038/348325a0
ADSWeb of Science®Google Scholar
Peltier, W. R. (1998), Postglacial variations in the level of the sea: Implications for climate dynamics and solid-earth geophysics, Rev. Geophys., 36(4), 603–689.
10.1029/98RG02638
ADSWeb of Science®Google Scholar
Rabinowicz, M., A. Nicolas, and J. L. Vigneresse (1984), A rolling-mill effect in the asthenosphere beneath oceanic spreading centers, Earth Planet. Sci. Lett., 67, 97–108.
10.1016/0012-821X(84)90042-6
ADSWeb of Science®Google Scholar
Rabinowicz, M., G. Ceuleneer, and A. Nicolas (1987), Melt segregation and flow in mantle diapirs below spreading centers: Evidence from the Oman ophiolite, J. Geophys. Res., 92(B5), 3475–3486.
10.1029/JB092iB05p03475
ADSWeb of Science®Google Scholar
Raddick, M. J., E. M. Parmentier, and D. S. Scheirer (2002), Buoyant decompression melting: A possible mechanism for intraplate volcanism, J. Geophys. Res., 107(B10), 2228, doi:10.1029/2001JB000617.
10.1029/2001JB000617
ADSWeb of Science®Google Scholar
Ribe, N. M. (1985a), The generation and composition of partial melts in the Earth's mantle, Earth Planet. Sci. Lett., 73, 361–376.
10.1016/0012-821X(85)90084-6
CASADSWeb of Science®Google Scholar
Ribe, N. M. (1985b), The deformation and compaction of partial molten zones, Geophys. J. R. Astron. Soc., 83, 487–501.
10.1111/j.1365-246X.1985.tb06499.x
ADSWeb of Science®Google Scholar
Rubin, K. H., and J. M. Sinton (2007), Inferences on mid-ocean ridge thermal and magmatic structure from MORB compositions, Earth Planet. Sci. Lett., 260(1–2), 257–276, doi:10.1016/j.epsl.2007.05.035.
10.1016/j.epsl.2007.05.035
CASADSWeb of Science®Google Scholar
Scheirer, D. S., D. W. Forsyth, M. H. Cormier, and K. C. Macdonald (1998), Shipboard geophysical indications of asymmetry and melt production beneath the East Pacific Rise near the MELT experiment, Science, 280(5367), 1221–1224.
10.1126/science.280.5367.1221
CASADSPubMedWeb of Science®Google Scholar
Schmeling, H. (2000), Partial melting and melt segregation in a convecting mantle, in Physics and Chemistry of Partially Molten Rocks, edited by N. Bagdassarov, D. Laporte, and A. B. Thompson, pp. 141–178, Kluwer Acad., Dordrecht, Netherlands.
10.1007/978-94-011-4016-4_5
Google Scholar
Schmeling, H. (2010), Dynamic models of continental rifting with melt generation, Tectonophysics, 480(1–4), 33–47, doi:10.1016/j.tecto.2009.09.005.
10.1016/j.tecto.2009.09.005
ADSWeb of Science®Google Scholar
Schouten, H., H. J. B. Dick, and K. D. Klitgord (1987), Migration of mid-ocean ridge volcanic segments, Nature, 326(6116), 835–839.
10.1038/326835a0
ADSWeb of Science®Google Scholar
Scott, D. R., and D. J. Stevenson (1986), Magma ascent by porous flow, J. Geophys. Res., 91, 9283–9296.
10.1029/JB091iB09p09283
ADSWeb of Science®Google Scholar
Scott, D. R., and D. J. Stevenson (1989), A self-consistent model of melting, magma migration and buoyancy-driven circulation beneath mid-ocean ridges, J. Geophys. Res., 94, 2973–2988.
10.1029/JB094iB03p02973
ADSWeb of Science®Google Scholar
Simpson, G., M. Spiegelman, and M. I. Weinstein (2010a), A multiscale model of partial melts: 1. Effective equations, J. Geophys. Res., 115, B04410, doi:10.1029/2009JB006375.
10.1029/2009JB006375
ADSPubMedWeb of Science®Google Scholar
Simpson, G., M. Spiegelman, and M. I. Weinstein (2010b), A multiscale model of partial melts: 2. Numerical results, J. Geophys. Res., 115, B04411, doi:10.1029/2009JB006376.
10.1029/2009JB006376
PubMedGoogle Scholar
Small, C. (1994), A global analysis of mid-ocean ridge axial topography, Geophys. J. Int., 116(1), 64–84, doi:10.1111/j.1365-246X.1994.tb02128.x.
10.1111/j.1365-246X.1994.tb02128.x
ADSWeb of Science®Google Scholar
Small, C., and L. V. Danyushevsky (2003), Plate-kinematic explanation for mid-oceanic-ridge depth discontinuities, Geology, 31(5), 399–402.
10.1130/0091-7613(2003)031<0399:PEFMDD>2.0.CO;2
ADSWeb of Science®Google Scholar
Sotin, C., and E. M. Parmentier (1989), Dynamical consequences of compositional and thermal density stratification beneath spreading centers, Geophys. Res. Lett., 16(8), 835–838.
10.1029/GL016i008p00835
ADSWeb of Science®Google Scholar
Sparks, D. W., and E. M. Parmentier (1991), Melt extraction from the mantle beneath spreading centers, Earth Planet. Sci. Lett., 105(4), 368–377.
10.1016/0012-821X(91)90178-K
ADSWeb of Science®Google Scholar
Sparks, D. W., and E. M. Parmentier (1993), The structure of 3-dimensional convection beneath oceanic spreading centers, Geophys. J. Int., 112(1), 81–91.
10.1111/j.1365-246X.1993.tb01438.x
ADSWeb of Science®Google Scholar
Spiegelman, M. (1993), Physics of melt extraction: Theory, implications, and applications, Philos. Trans. R. Soc. London A, 342, 23–41.
10.1098/rsta.1993.0002
ADSWeb of Science®Google Scholar
Spiegelman, M. (1996), Geochemical consequences of melt transport in 2-D: The sensitivity of trace elements to mantle dynamics, Earth Planet. Sci. Lett., 139, 115–132.
10.1016/0012-821X(96)00008-8
CASADSWeb of Science®Google Scholar
Spiegelman, M. (2003), Linear analysis of melt band formation by simple shear, Geochem. Geophys. Geosyst., 4(9), 8615, doi:10.1029/2002GC000499.
10.1029/2002GC000499
CASADSWeb of Science®Google Scholar
Spiegelman, M., and D. McKenzie (1987), Simple 2-D models for melt extraction at mid-ocean ridges and island arcs, Earth Planet. Sci. Lett., 83, 137–152.
10.1016/0012-821X(87)90057-4
ADSWeb of Science®Google Scholar
Spiegelman, M., and J. R. Reynolds (1999), Combined dynamic and geochemical evidence for convergent melt flow beneath the East Pacific Rise, Nature, 402, 282–285.
10.1038/46260
CASADSWeb of Science®Google Scholar
Šrámek, O., Y. Ricard, and D. Bercovici (2007), Simultaneous melting and compaction in deformable two-phase media, Geophys. J. Int., 168, 964–982, doi:10.1111/j.1365-246X.2006.03269.x.
10.1111/j.1365-246X.2006.03269.x
ADSWeb of Science®Google Scholar
Stevenson, D. J. (1989), Spontaneous small-scale melt segregation in partial melts undergoing deformation, Geophys. Res. Lett., 16, 1067–1070.
10.1029/GL016i009p01067
ADSWeb of Science®Google Scholar
Stracke, A., B. Bourdon, and D. McKenzie (2006), Melt extraction in the Earth's mantle: Constraints from U-Th-Pa-Ra studies in oceanic basalts, Earth Planet. Sci. Lett., 244, 97–112, doi:10.1016/j.epsl.2006.01.057.
10.1016/j.epsl.2006.01.057
CASADSWeb of Science®Google Scholar
Tackley, P. J., and D. J. Stevenson (1993), A mechanism for spontaneous self-perpetuating volcanism on the terrestrial planets, in Flow and Creep in the Solar System: Observations, Modeling and Theory, NATO ASI Ser., Ser. C, vol. 391, edited by D. B. Stone, and S. K. Runcorn, pp. 307–321, Kluwer, Dordrecht, Netherlands.
10.1007/978-94-015-8206-3_19
Web of Science®Google Scholar
Takei, Y. (2010), Stress-induced anisotropy of partially molten rock analogue deformed under quasi-static loading test, J. Geophys. Res., 115, B03204, doi:10.1029/2009JB006568.
10.1029/2009JB006568
ADSPubMedWeb of Science®Google Scholar
Takei, Y., and B. K. Holtzman (2009a), Viscous constitutive relations of solid-liquid composites in terms of grain boundary contiguity: 1. Grain boundary diffusion control model, J. Geophys. Res., 114, B06205, doi:10.1029/2008JB005850.
10.1029/2008JB005850
PubMedWeb of Science®Google Scholar
Takei, Y., and B. K. Holtzman (2009b), Viscous constitutive relations of solid-liquid composites in terms of grain boundary contiguity: 2. Compositional model for small melt fractions, J. Geophys. Res., 114, B06206, doi:10.1029/2008JB005851.
10.1029/2008JB005851
Web of Science®Google Scholar
Takei, Y., and B. K. Holtzman (2009c), Viscous constitutive relations of solid-liquid composites in terms of grain boundary contiguity: 3. Causes and consequences of viscous anisotropy, J. Geophys. Res., 114, B06207, doi:10.1029/2008JB005852.
10.1029/2008JB005852
Web of Science®Google Scholar
Toomey, D. R., W. S. D. Wilcock, S. C. Solomon, W. C. Hammond, and J. A. Orcutt (1998), Mantle seismic structure beneath the MELT region of the East Pacific Rise from P and S wave tomography, Science, 280(5367), 1224–1227.
10.1126/science.280.5367.1224
CASADSPubMedWeb of Science®Google Scholar
Toomey, D. R., W. S. D. Wilcock, J. A. Conder, D. W. Forsyth, J. D. Blundy, E. M. Parmentier, and W. C. Hammond (2002), Asymmetric mantle dynamics in the MELT region of the East Pacific Rise, Earth Planet. Sci. Lett., 200(3–4), 287–295.
10.1016/S0012-821X(02)00655-6
CASADSWeb of Science®Google Scholar
Turner, J. S. (1973), Buoyancy Effects in Fluids, Cambridge Univ. Press, Cambridge, U. K.
10.1017/CBO9780511608827
Google Scholar
Wang, Y., D. W. Forsyth, and B. Savage (2009), Convective upwelling in the mantle beneath the Gulf of California, Nature, 462(7272), 499–501, doi:10.1038/nature08552.
10.1038/nature08552
CASADSPubMedWeb of Science®Google Scholar
Wark, D. A., and E. B. Watson (1998), Grain-scale permeabilities of texturally equilibrated, monomineralic rocks, Earth Planet. Sci. Lett., 164, 591–605.
10.1016/S0012-821X(98)00252-0
CASADSWeb of Science®Google Scholar
Weatherley, S. M., and R. F. Katz (2010), Plate-driven mantle dynamics and global patterns of mid-ocean ridge bathymetry, Geochem. Geophys. Geosyst., 11, Q10003, doi:10.1029/2010GC003192.
10.1029/2010GC003192
ADSWeb of Science®Google Scholar
White, R. S., T. A. Minshull, M. J. Bickle, and C. J. Robinson (2001), Melt generation at very slow-spreading oceanic ridges: Constraints from geochemical and geophysical data, J. Petrol., 42(6), 1171–1196.
10.1093/petrology/42.6.1171
CASADSWeb of Science®Google Scholar
Wolfe, C. J., and S. C. Solomon (1998), Shear-wave splitting and implications for mantle flow beneath the MELT region of the East Pacific Rise, Science, 280(5367), 1230–1232.
10.1126/science.280.5367.1230
CASADSPubMedWeb of Science®Google Scholar
Yang, Y., D. W. Forsyth, and D. S. Weeraratne (2007), Seismic attenuation near the East Pacific Rise and the origin of the low-velocity zone, Earth Planet. Sci. Lett., 258(1–2), 260–268, doi:10.1016/j.epsl.2007.03.040.
10.1016/j.epsl.2007.03.040
CASADSWeb of Science®Google Scholar
Zhu, W. L., and G. Hirth (2003), A network model for permeability in partially molten rocks, Earth Planet. Sci. Lett., 212, 407–416.
10.1016/S0012-821X(03)00264-4
CASADSWeb of Science®Google Scholar

Citing Literature

Volume11, Issue11

November 2010

This article also appears in:

Geochemical Heterogeneities in Oceanic Island Basalt and Mid-Ocean Ridge Basalt Sources: Implications for Melting Processes and Mantle Dynamics

Porosity-driven convection and asymmetry beneath mid-ocean ridges

Abstract

1. Introduction

2. Previous Geodynamic Models of Active Flow at MORs

3. Model Description

3.1. Petrology

3.2. Rheology

3.2.1. Shear Viscosity

3.2.2. Bulk Viscosity

3.3. Melt Extraction at the Ridge

4. Scaling Analysis

5. Results

5.1. Convective Vigor

5.2. Crustal Generation

5.3. Symmetry Breaking by Convection

5.3.1. Temperature Forcing

5.3.2. Compositional Forcing

5.3.3. Ridge Migration

6. Discussion

6.1. Observational Constraints on Active Flow

6.2. Thermal and Compositional Density Variations

6.3. On the Absence of Porosity/Shear Bands

6.4. Toward Models of 3-D Magma/Mantle Flow and Ridge Relocation

7. Summary

Acknowledgments

Appendix A:: Governing Equations

A1. Conservation of Mass and Momentum

A2. Phase Densities, Pressure Decomposition, and Reformulation

A3. Conservation of Energy and Species Mass

A4. Nondimensionalization

A5. Numerical Solution

Appendix B:

References

Citing Literature

Figures

References

Related

Information

RESOURCES

PUBLICATION INFO