Journal of Petrology

Journal of Petrology Advance Access published online on December 17, 2008
Journal of Petrology, doi:10.1093/petrology/egn058

This Article Abstract FREE Full Text (PDF) All Versions of this Article: 49/12/2099    most recent egn058v1 Alert me when this article is cited Alert me if a correction is posted Services Email this article to a friend Similar articles in this journal Alert me to new issues of the journal Add to My Personal Archive Download to citation manager Request Permissions Google Scholar Articles by Katz, R. F. Search for Related Content GeoRef GeoRef Citation Social Bookmarking          What's this?

© The Author 2008. Published by Oxford University Press. All rights reserved. For Permissions, please e-mail: journals.permissions@oxfordjournals.org

Magma Dynamics with the Enthalpy Method: Benchmark Solutions and Magmatic Focusing at Mid-ocean Ridges

Richard F. Katz^*

Institute for Theoretical Geophysics, Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK

Received May 13, 2008; Revised typescript accepted October 15, 2008

	ABSTRACT

TOP ABSTRACT INTRODUCTION MODEL FORMULATION SOLUTIONS TO THE EQUATIONS DISCUSSION SUMMARY AND CONCLUSIONS Appendix A: Conservation of... Appendix B: Non... Appendix C: The enthalpy... REFERENCES

 
Magma genesis and transport link mantle convection with surface volcanism and hence with the long-term chemical and morphological evolution of the Earth's; crust. Modeling the dynamics of magma–mantle interaction in tectonic settings remains a challenge, however, because of the complexity of multi-component thermodynamics and melt segregation in a permeable, compactible, and actively deforming mantle matrix. Here I describe a flexible approach to formulating the thermochemistry of such models based on the Enthalpy Method, a technique commonly used in simulations of alloy solidification. This approach allows for melting and freezing based on a familiar binary phase diagram, consistent with conservation of energy and two-phase compaction and flow. I present an extension of the Enthalpy Method to more than two thermodynamic components. Simulation of a one-dimensional upwelling and melting column provides a benchmark for the method. Two-dimensional simulations of the melting region that feeds magma to a rapidly spreading mid-ocean ridge demonstrate the utility of the Enthalpy Method. These calculations provide a new estimate of the efficiency of magmatic focusing along the base of the oceanic lithosphere. Modeled focusing efficiency varies with mantle permeability and resistance to compaction. To yield 5–7 km of oceanic crust with 20% melting of a homogeneous, sub-ridge mantle, a focusing efficiency of greater than 70% is required. This, in turn, suggests that matrix permeability and bulk viscosity are at the high end of previously estimated values.

KEY WORDS: mantle; simulation; PETSc; compaction; thermodynamics

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MODEL FORMULATION SOLUTIONS TO THE EQUATIONS DISCUSSION SUMMARY AND CONCLUSIONS Appendix A: Conservation of... Appendix B: Non... Appendix C: The enthalpy... REFERENCES

 
The convective dynamics of the mantle exert a primary yet incompletely understood influence on the surface environment of the Earth. Mantle convection is linked to the physical and chemical characteristics of the planet's; surface by magma genesis and transport. This connection affects the evolution of hotspot chains, explosive volcanoes found above subduction zones, and tectonic seams of the ocean floor, as well as the existence and composition of continents. Understanding the transport of magma through the subsurface is an essential component in our knowledge of the Earth system.

Much previous work has focused on solid-state, creeping convection of the whole mantle (e.g. Schubert et al., 2001) and the origins of plate tectonics (e.g. Tackley, 2000). Relatively less attention has focused on the interaction of fluid magma with the permeable, crystalline mantle in the partially molten zones beneath volcanoes. McKenzie (1984) and subsequent studies (see references below) have derived general equations that aim to describe the conservation of mass, momentum and energy in such settings. Solutions to these equations, plus some treatment of melting, freezing and geochemical transport, can be used to make predictions that are testable using geophysical and geochemical data. The set of predictions that can be made using any formulation of the governing equations depends on the choice of complexities included in the model.

Observations of mid-ocean ridges (MORs) pose a set of fundamental problems in magmatic transport. Of particular interest here is the observed uniformity of crustal thickness (5–7 km) for ridge full-spreading rates above 2 cm/yr. This feature of MORs was demonstrated and modeled by Bown & White (1994). Those researchers showed that crustal thickness can be modeled by decompression melting of a mantle with a potential temperature between 1280 and 1320°C, assuming instantaneous, complete melt extraction. A higher mantle potential temperature would result in a larger melt production rate and would require less efficient extraction, although geochemical constraints on the degree of melting indicate that the potential temperature cannot be significantly higher than 1350°C. This reasoning suggests that melt extraction is uniformly efficient for all oceanic spreading ridges.

Beneath MORs, magma is produced over a volume of mantle that can extend to more than 100 km on either side of the ridge axis (Forsyth et al., 1998); efficient melt extraction requires that this magma be focused laterally toward the ridge axis. The mechanics of magmatic focusing at mid-ocean ridges remains incompletely understood [for a review see Kelemen et al. (1997)]. Models include flow focusing as a result of anisotropic permeability (Phipps Morgan, 1987; Daines & Kohlstedt, 1997; Katz et al., 2006), pressure effects caused by mantle corner flow (Phipps Morgan, 1987; Spiegelman & McKenzie, 1987) and channelized flow along the base of the sloping thermal boundary layer in a high-porosity ‘decompaction channel’ (Sparks & Parmentier, 1991; Spiegelman, 1993c). Of these, the last remains a promising explanation for high-efficiency focusing, although it has not been thoroughly investigated. Sparks & Parmentier (1991) used a semi-analytical analysis of the problem to derive an estimate of focusing efficiency as a function of mantle permeability and magma viscosity. Spiegelman (1993c) developed two-dimensional (2D), isoviscous, numerical simulations of a fixed sloping boundary with a prescribed freezing rate. The work showed that the efficiency of focusing depends on the ratio of the crystallization-region thickness to the local compaction length. In particular, deflection of flow into the decompaction channel occurs only if the crystallizing boundary layer is sharp relative to the compaction length. Ghods & Arkani-Hamed (2000) performed 2D numerical simulations of melting, freezing and magmatic transport at a MOR to better constrain the efficiency of focusing in a ridge setting. Their estimates of efficiency were lower than those of Sparks & Parmentier (1991) and those presented here. More recently, observations of the Oman ophiolite were interpreted by Rabinowicz & Ceuleneer (2005) in terms of the presence of a decompaction layer. That work also described numerical simulations but did not report the efficiency of focusing. The question of the efficiency of magmatic focusing by flow along a sloping decompaction layer thus remains unresolved.

It is the purpose of the present paper to reconsider this problem with the introduction of a new model that incorporates a flexible approach to handling the chemical thermodynamics of magma dynamics simulations. Model complexity is limited to that necessary to address the phenomenon in question. Melting and freezing are required, as are melt segregation and matrix compaction (Spiegelman, 1993c). To establish a consistent ridge thermal structure and melting budget, solution of a conservation equation for energy is also required. This equation should account for heat transport by the solid and fluid, latent heat of melting or freezing, thermal diffusion, and adiabatic temperature changes (Bown & White, 1994). To account for the compositional dependence of mantle melting in the simplest possible way, the system should include at least two thermochemical components such that the melting point of a parcel of solid is at least univariant (at a given pressure). Furthermore, the pressure dependence of the solidus and liquidus is important for calculating the distribution of melting (Asimow et al., 1997). Finally, to properly calculate the compaction length, a physically consistent bulk-viscosity formulation must be considered (Schmelling, 2000; Bercovici et al., 2001). The approach used in the present work combines the magma dynamics formulation of Katz et al. (2007) and the Enthalpy Method (Alexiades & Solomon, 1993). The latter has been frequently used in simulations of solidification of binary alloys (e.g. Oertling & Watts, 2004; Katz & Worster, 2008). This approach accommodates the minimum set of requirements for the problem of magmatic focusing as defined by Sparks & Parmentier (1991); it can be extended in a straightforward manner to more complex constitutive relations and thermodynamic systems. It is described in detail below.

The present work builds on that of previous researchers by incorporating thermodynamic and fluid-mechanical complexity into a 2D, tectonic-scale model. Past work on applications of magma dynamics theory has typically considered a simplified mechanical system with more detailed thermochemistry, or vice versa. Details of thermodynamics and melting have typically been included in 1D ‘melting column’ models. Ribe (1985a) considered a two-component system in thermodynamic equilibrium prescribed by a phase diagram (full solid solution and eutectic). Asimow & Stopler (1999) extended this theory to account for a multi-component system with variations in density and partial specific entropy of multiple phases, using MELTS (Ghioroso, 1994; Ghiorso & Sack, 1995) to resolve the thermodynamic quantities. Other melting column models, by rámek et al. (2007) and Hewitt & Fowler (2008), have considered only a single thermodynamic component and have instead focused on deriving analytical solutions to expose the relevant fluid mechanical processes. Hewitt & Fowler (2008) considered a column capped by a cold boundary that yielded a lithospheric boundary layer and magmatic under-plating beneath it.

The present model adopts a level of thermodynamic complexity equivalent to that of Ribe (1985a) but is most closely associated with previous tectonic-scale 2D and 3D simulations of MORs. Most of these, however, have neglected compaction stresses so as to reduce the momentum relation for the solid to the familiar incompressible Stokes' equation. Spiegelman & McKenzie (1987) considered 2D, constant viscosity, constant porosity models of ridges and arcs to investigate fluid focusing caused by dynamic pressure gradients. Scott & Stevenson (1989) neglected compaction stresses but allowed for variable viscosity in 2D numerical models of melt flow beneath ridges. This was extended to three dimensions and was augmented with wet melting by Choblet & Parmentier (2001). Compaction stresses were accounted for by Sparks & Parmentier (1991) in a study of melt focusing by buoyant flow along the base of the impermeable lithosphere at ridges. Spiegelman (1993c, 1996) included compaction stresses in numerical models and investigated the consequences of active, buoyancy-driven upwelling on melting and melt migration. These studies used a simple melting parameterization in which the melting rate is directly proportional to the upwelling rate and neglected freezing entirely. Ghods & Arkani-Hamed (2000) extended previous models by incorporating temperature-dependent melting or freezing and conservation of energy.

The present model follows most closely from Ghods & Arkani-Hamed (2001). It is a 2D solution of the equations expressing conservation of mass, momentum and energy in the mid-ocean ridge setting. It consistently accounts for melting and freezing using a binary phase diagram—a familiar and transparent thermodynamic parameterization. It incorporates compaction stresses and employs a physically consistent bulk viscosity formulation that has a singularity for zero porosity (Batchelor, 1967; Schmelling, 2000; Bercovici et al., 2001; Hewitt & Fowler, 2008; Simpson, 2008).

Reactive and mechanical localization instabilities are absent from the simulations presented here. Localization caused by matrix deformation and porosity-dependent viscosity was first noted in a linear stability analysis by Stevenson (1989). It has been the subject of both experimental (e.g. Holtzman et al., 2003) and theoretical investigations. Richardson (1998) considered the effects of magmatic buoyancy and Hall & Parmentier (2000) investigated the effects of water on the instability. Spiegelman (2003) performed a stability analysis for simple-shear deformation of a partially molten aggregate with a Newtonian viscosity. This analysis was extended by Katz et al. (2006) to model non-Newtonian viscosity. Magmatic localization as a result of fluid–matrix reactions was considered by Aharonov et al. (1995), who showed that uniform porous flow up a sufficiently strong gradient in solubility would localize into a channelized porosity structure. Spiegelman et al. (2001) presented numerical simulations of reactive, channelized flows with pressure-dependent solubility. Temperature-dependent reactive melting was modeled by Katz (2005) in the context of simplified models of subduction zones.

Questions regarding the formation and behaviour of high-permeability channels beneath mid-ocean ridges have important implications for our understanding of magmatic transport. For example, if excesses of ²²⁶Ra relative to ²³⁰Th observed in young ridge lavas (e.g. Sims et al., 2002; Stracke et al., 2006) originate by fractionation of uranium-series elements in the garnet stability field and are preserved until eruption at the ridge, then magma ascent rates must be high and residual porosity low (McKenzie, 1985). Jull et al. (2002) has shown that rapid, equilibrium melt extraction through high-porosity channels may be sufficient to explain thorium excesses that originate at depths of garnet stability. In this model, radium excesses are generated in the low-porosity inter-channel regions that feed the channels at shallow depths [see also Elliott & Spiegelman (2003)]. If radium excesses are generated by chromatographic effects (without channelization) (Spiegelman & Elliott, 1993) then the constraint on transport rate is relaxed somewhat, although residual porosity in the mantle must still be of the order of the elemental distribution coefficients. Saal & van Orman (2004) proposed an alternative theory in which radium excesses are generated within the magma chamber itself.

In seeking to present a simple exposition of magma dynamics with the Enthalpy Method, I have avoided the conditions that give rise to magmatic localization instabilities. In particular, for the purposes of the present study, the shear viscosity is taken as constant, independent of porosity, to avoid the possibility of mechanical instability. Furthermore, although the pressure and temperature dependence of solubility in the binary phase diagram suggest the possibility of reactive localization, compositional perturbations required to nucleate reactive instabilities have been suppressed. Modification of these attributes may allow for the investigation of mechanical and reactive localization processes leading to channelization of magmatic flux in the ridge melting region using a modified version of the same simulation code that was used for this study.

The present study concentrates on an exposition and benchmark of magma dynamics simulations using the Enthalpy Method, as well as a prediction of the efficiency of melt focusing by buoyant magmatic flow along the base of the lithospheric thermal boundary layer. I explore the sensitivity of focusing efficiency to the magnitude of permeability and bulk viscosity. Given the uncertainty in these parameters, simulations predict a wide range of possible focusing efficiency (see below).

The paper is organized as follows. The next section introduces the relevant theory and equations. The section about solutions to the equations briefly describes the numerical approach taken here and reports on results from 1D and 2D simulations. The following sections provide a discussion of the results and some conclusions. Three appendices detail the derivation of the conservation of enthalpy equation, non-dimensionalization of the governing equations and an extension of the Enthalpy Method to systems with more than two thermochemical components.

	MODEL FORMULATION

TOP ABSTRACT INTRODUCTION MODEL FORMULATION SOLUTIONS TO THE EQUATIONS DISCUSSION SUMMARY AND CONCLUSIONS Appendix A: Conservation of... Appendix B: Non... Appendix C: The enthalpy... REFERENCES

 
The model consists of a set of coupled partial differential equations (PDEs) to describe the essential features of the magma–mantle system, allowing for consistent melting and freezing, segregation of melt from the crystalline mantle matrix, shear and compactive deformation of the matrix, and transport of energy by both fluid and matrix phases. To assemble the appropriate PDEs, I adopt the continuum theory of McKenzie (1984) and develop an approach to the thermodynamics of magma transport based on the Enthalpy Method.

The Enthalpy Method allows the calculation of melting and freezing rates based entirely on an equilibrium phase diagram (Alexiades & Solomon, 1993). It has been commonly used to model solidification problems concerning the formation of a mushy layer (e.g. Oertling & Watts, 2004; Katz & Worster, 2008). Using the Enthalpy Method, the partial differential equation describing the evolution of porosity (the volume fraction of fluid present in a representative volume element of the domain) is replaced with a closure condition between the local bulk enthalpy and bulk composition. This closure condition is derived from the prescribed phase diagram that spans all the relevant thermodynamic components. This approach has several advantages over other melting models. First, it prevents the occurrence of negative porosity values that can appear in numerical solutions of the PDE governing the evolution of porosity. Second, it avoids the need for opaque melting parameterizations. Finally, it reduces the number of coupled equations that must be solved by eliminating the PDE for porosity evolution. The disadvantage of the Enthalpy Method, as it is described here, is that it requires the assumption of thermodynamic equilibrium throughout the domain. Fortunately, its use does not preclude the addition of auxiliary calculations of disequilibrium geochemical transport of trace elements and radiogenic nucleides (e.g. Spiegelman, 1996).

The assumption of thermodynamic equilibrium is valid if reactions toward equilibrium are sufficiently fast. In particular, the mantle–melt system should be in equilibrium if the length-scale over which melt equilibrates with the mantle is of the order of the continuum scale (a few tens of grain diameters). Because the crystals are mostly composed of fusible material (i.e. major elements), reaction is limited by the rate of diffusion of major elements into the melt, away from the crystal–melt interfaces. Aharonov et al. (1995) analyzed this system for a broad range of reasonable parameter values and estimated an equilibration length that is between ångströms and meters. Assuming mantle grains are no larger than a few centimeters in diameter puts an upper bound on the continuum scale for magma dynamics of about 1m, a good match with the maximum equilibration length-scale from Aharonov et al. (1995). Hence it is plausible that the mantle can be described by a model with local thermochemical equilibrium, at least for major elements.

Using that approach, the thermodynamic quantity that must be explicitly conserved is the gravity-compensated enthalpy (Ramberg, 1971). In a unit volume containing both fluid and matrix phases, the magnitude of this quantity is given by – g · x, where is the density, h is the enthalpy per unit mass, g is the gravitational acceleration vector, x is the position vector, and overbars represent volume-averaged bulk quantities (i.e. for some quantity Q with different values for the fluid and matrix phases, = Q_f + (1–)Q_m). The first term, , represents the volumetric enthalpy and the second term, g · x, represents the potential energy; kinetic energy is neglected because the Reynolds number for both phases is much less than unity. For what follows, it will be useful to define the volumetric bulk enthalpy, H = .

Fluid mechanics
Coupled partial differential equations describing the motion of magma in the convecting mantle have been derived by several workers (Ahern & Turcotte, 1979; McKenzie, 1984; Fowler, 1985; Ribe, 1985b; Scott & Stevenson, 1986). A useful review has been provided by Stevenson & Scott (1991). Bercovici et al. (2001) derived a more general version describing two symmetric, immiscible fluids. However it was shown by Bercovici & Ricard (2003) that in the case relevant to magma dynamics of a fluid with viscosity that is much smaller that the matrix viscosity (and neglecting surface tension and damage), the general equations reduce to the set derived by McKenzie (1984). Here, the formalism of McKenzie (1984) is used for consistency with past work.

The equations describing conservation of mass for the fluid and matrix phases, respectively, are

(1)

(2)

where is the porosity, v_f and v_m are the fluid and matrix velocities, and is the rate of mass transfer from the fluid to the solid phase (the melting rate).

The balance of forces in the fluid phase takes the form of a modified Darcy's; law,

(3)

where K is the permeability, µ is the fluid viscosity, and P_f is the fluid pressure. For the matrix phase, the balance of forces is a modified Stokes' law,

(4)

where and are the shear and bulk viscosity of the aggregate, and the superscript T indicates the matrix transpose. The basic properties of these equations have been previously studied by many workers including Barcilon & Lovera (1989) and Spiegelman (1993a, 1993b). Katz et al. (2007) described a reformulation of the equations that is convenient for implicit numerical solution; this form is given below in equations (17)–(21).

It is worth stressing that the formulation of McKenzie (1984) does take into account the pressure difference between fluid and matrix phases as the driving force for compaction. The stress tensor of the matrix phase, , was given by equations. (A14) and (A15) of McKenzie (1984) as

(5)

where _ij is the identity matrix, subscripts represent component indicies into vectors and matricies, and the Einstein summation convention of repeated indicies applies. As usual, the dynamic pressure of the solid is defined as

(6)

where Tr denotes the trace operator. Equation (6) can be rearranged to give the familiar and physically reasonable equation for the pressure difference between phases P_f – P_m = · v_m.

Energy and composition
Energy is conserved for an arbitrary Eulerian volume V containing, in general, both fluid and solid phases. Changes in the total energy within the volume can be related to fluxes of energy across its boundary V. This relation is given by

(7)

where x is the coordinate vector, is the phase-averaged thermal conductivity, T represents temperature, and is an outward pointing unit vector normal to the volume boundary. Terms on the right hand side represent fluxes as a result of advection of enthalpy, advection of potential energy, and diffusion of sensible heat respectively. For clarity of interpretation, viscous dissipation, radiogenic heat production and other irreversible processes have been neglected in writing equation (7). These contributions have been considered by other workers (e.g. McKenzie, 1984; Bercovici & Ricard, 2003; rámek et al., 2007) and may be included in future work.

The differential form of the conservation of energy equation is derived from (7) in Appendix A. This derivation is based on assumptions of thermal equilibrium everywhere within the domain as well as constant, phase-independent material properties (density, specific heat, thermal expansivity and thermal conductivity). The result is

(8)

where c_p is the specific heat, is the thermal expansivity, z is a coordinate representing depth, T = T exp(–gz/c_P) is the potential temperature, k is the thermal diffusivity and g is the magnitude of the gravity vector. Equation (8) states that changes in the volumetric bulk enthalpy are due to advection of sensible heat, advection of latent heat and thermal diffusion. The contribution of potential energy is accounted for implicitly through the use of potential temperature (see Appendix A). As g/c_P << 1, the exponential coefficients in equation (8) could be linearized. However, this would neither clarify the interpretation of the equation nor speed the numerical solution and so the exponentials are left unchanged.

Application of the Enthalpy Method requires that we know both the bulk enthalpy H and the bulk composition C everywhere within the domain. For simplicity, we limit the composition to two thermodynamic components. In that limit, and with _f = _m =, the conservation of bulk composition is governed by the single equation

(9)

where C_f and C_m are the mass concentrations of the less fusible component in the fluid and matrix phases respectively, and is the chemical diffusivity or dispersivity (see Appendix A). This equation states that changes in the bulk composition are due to advection by the fluid and matrix, as well as by diffusion within the fluid phase.

Equations (8) and (9) constrain the bulk enthalpy and composition. However, their solution requires knowledge of four other variables, , T, C_f and C_m. These are provided by the Enthalpy Method, which allows for the derivation of a set of algebraic closure conditions that relate these four unknowns to bulk enthalpy and composition, as detailed in the following subsection.

The Enthalpy Method
The Enthalpy Method is based on the prescription of thermodynamic equilibrium everywhere in the domain. This condition, quantified by a phase diagram, provides closure conditions for porosity, temperature, and the two phase compositions in equations (8) and (9) as a function of pressure or depth. For a two-component system, the simplest phase diagrams are the binary phase loop, which applies in the case of total solid solution of the thermodynamic components, and the eutectic phase diagram. Ribe (1985a) has incorporated both of these diagrams in calculations of 1D melting columns. Here I consider only the phase loop, shown in Fig. 1. This a vast simplification of the full thermodynamic system; however, it has long been known that mantle melting is not eutectic-like; a continuous variation of the solidus and liquidus with extent of melting is a more reasonable model (Asimow et al., 1997). The details of the phase diagram shown in Fig. 1 are ad hoc—they are chosen for mathematical convenience and for consistency with an idealized conception of mantle melting. A more rigorous treatment would use thermodynamic laws and measurements to constrain the shape of the liquidus and solidus. For the present purposes, such an approach is not required because leading-order features of the overall system are insensitive to these details.

View larger version (60K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 1. A binary, pressure-dependent, phase diagram for a two-component where system. The functions used to generate this figure can be written as and TL , where T = T1 – T0, Pl is the lithostatic presure, and is the Clapeyron slope. It should be noted that these functions are the inverse of fs and fl from equations (14) and (15).

 
At each Eulerian grid cell in the domain, the energy available for partition between sensible and latent heat is given by the value of the bulk enthalpy, H = . Neglecting small changes to the pressure within the cell as a result of fluid dynamics (i.e. assuming dP 0), the total differential (A5) can be integrated to give the enthalpy per unit mass for the matrix h_m and fluid h_f phases as

(10)

(11)

where h₀ is a reference enthalpy at the reference temperature T₀, L is the latent heat of the fluid, and T is the temperature. h₀ is assigned to be zero at the minimum melting temperature T₀ over all possible compositions. As above, the specific heat c_p has been taken as constant and equal between phases.

From equations (10) and (11) the bulk volumetric enthalpy can be constructed as

(12)

Unlike Asimow et al. (1997), I neglect variations in the partial specific entropy of the fluid and matrix and assume a constant latent heat of fusion. Experiments on basalt by Bouhifd et al. (2007) over temperatures relevant to mantle melting have shown that L/T 375 J/kg/K. This means that a 100 degree change in temperature changes L by approximately 9% from the value used here (see Table 1).

View this table: [in this window] [in a new window]   Table 1: Parameters and their representative values for the mantle

 
Three additional equations are needed to solve for , T, C_f and C_m; these are given by the definition of bulk composition and the phase diagram as

(13)

(14)

(15)

P_l is the lithostatic pressure, a good approximation of the total pressure at any point within the domain. Equation (13) defines the bulk composition C, equations (14) and (13) are parameterizations of the solidus and liquidus surfaces in Fig. 1 (see figure caption for details). Combining equations (12)–(15) gives an equation for porosity,

(16)

In general, equation (16) must be solved numerically for given values of H and C. The solution can then be used in equations (12)–(15) to obtain the other three required variables, T, C_f and C_m.

Although I have used the binary phase loop here, this formulation is easily adapted to other binary phase diagrams. For example, the binary eutectic system was adopted by Katz & Worster (2008). Appendix C generalizes the Enthalpy Method to N thermochemical components, although it is clearly more difficult to construct functions describing the liquidus and solidus surfaces in this case.

Nondimensional system
Physical considerations of compaction of a two-phase medium suggest that the bulk viscosity must approach infinity as porosity tends toward zero (Batchelor, 1967; Schmelling, 2000; Bercovici et al., 2001; Hewitt & Fowler, 2008; Simpson, 2008). This introduces a singularity in equation (4). Rearranging the equations to put into the denominator facilitates the handling of this singularity. To accomplish this, I adopt the pressure decomposition described by Katz et al. (2007) and in Appendix B. This allows the system of fluid mechanical equations to be reduced to incompressible Stokes' flow when and where the porosity goes to zero.

Nondimensionalization of the governing equations is described in Appendix B. The system of equations can be written in terms of nondimensional variables as

(17)

(18)

(19)

(20)

(21)

where is the compaction pressure (known in rock mechanics as the excess pressure), P is the ‘dynamic’ pressure, =( –2/3) is the viscous resistance to compaction, is a unit vector in the direction of gravity, _t is a partial derivative with respect to time, and is the nondimensional temperature as defined in Appendix B. The adiabatic parameter is the proportional change in temperature as a result of adiabatic decompression over the compaction length. The Stefan number S = L/(c_P T) fixes the importance of latent heat relative to sensible heat in controlling changes in enthalpy. The thermal and compositional Peclet numbers, Pe_T = w₀/ and , characterize the importance of advection relative to diffusion. Representative mantle values for dimensional parameters are given in Table 1. Closure conditions for dimensionless Enthalpy Method variables , , C_f and C_m and the fluid velocity v_f are given in nondimensional form in Appendix B. With these closures, equations (17)–(21) are a closed set of 4+D coupled partial differential equations for 4+D primary variables , where D is the number of spatial dimensions.

It is important to note that the melting rate, , does not appear in the final set of PDEs. The assumption of local thermodynamic equilibrium everywhere in the domain means that is determined implicitly by the local rate of change of conditions that affect melting. Given a solution to the governing equations, can be extracted using the conservation of mass equation (2) in nondimensional form

(22)

Thus the melting rate is the time-rate of change of porosity that is not due to compaction.

The derivation of equations (17)–(21) is based on a number of assumptions. The most important is the assumption of local thermodynamic equilibrium everywhere in the domain. This allows for the local phase fractions, compositions and temperature to be determined from bulk enthalpy and composition using the Enthalpy Method. To simplify the equations to a manageable level of complexity, an extended Boussinesq approximation is applied. Using this approximation, variations in density that are not associated with buoyancy terms are neglected (except in allowing for adiabatic changes in temperature). Moreover, buoyancy is driven by a constant density difference between the phases—thermal and compositional buoyancy are neglected. The system of equations is further simplified by setting material properties such as specific heat and thermal expansivity equal between phases. A more rigorous treatment of the thermodynamics would require consistency between these material properties and the details of the binary phase diagram (e.g. Denbigh, 1981). Finally, irreversible sources of heat including dissipation and radiogenic heating are neglected.

Much complexity remains in these equations, however. As written, the equations allow for variations in viscosity. Such variations arise from, for example, gradients in porosity, temperature, and stress (e.g. Karato & Wu, 1993). Buoyancy caused by the presence of fluid is included in the equations, allowing for modeling of melting-induced solid upwelling (Buck & Su, 1989; Scott & Stevenson, 1989; Cordery & Phipps Morgan, 1992, 1993). The equations allow for segregation of melt and, with it, fluid transport of heat and chemistry. Chemical reactions between fluid and matrix are implicit; they occur as melt segregates and rises along a lithostatic pressure gradient. The inclusion of pressure gradients caused by compaction allows for localization of melt as a result of reactive (e.g. Aharonov et al., 1995; Spiegelman et al., 2001) and mechanical (e.g. Stevenson, 1989; Katz et al., 2006) mechanisms. Because of this complexity, analytical solutions to the full equations do not exist; numerical methods are required and these are discussed below.

Permeability and viscosity
In keeping with the principal goal of this paper, to demonstrate the utility of the Enthalpy Method for problems of magma dynamics, I have chosen to reduce constitutive equations to the simplest reasonable form. For the binary phase diagram (Fig. 1), I have applied the same Clapeyron slope to the entire phase diagram. Dimensional permeability is calculated according to the standard Kozeny–Carmen relationship (Bear, 1972), simplified for small porosity and constant grain size to (Wark & Watson, 1998; Wark et al., 2003)

(23)

where n is a constant exponent typically taken to equal 2 or 3.

The bulk viscosity is taken as a constant for the 1D melting column described below and as being proportional to the inverse porosity in 2D calculations (Batchelor, 1967; Schmelling, 2000; Bercovici et al., 2001; Simpson, 2008). In the latter case, it is given in dimensional form by

(24)

where _R is the ratio of bulk to shear viscosity at the reference porosity ₀. The product _R0 has been estimated experimentally by Cooper (1990) to be about 10 and theoretically by Hewitt & Fowler (2008) and Simpson (2008) to be about unity. The singularity in equation (24) has implications for the behaviour of the compaction length as 0. Ghods & Arkani-Hamed (2000) used a bulk viscosity formulation without a singularity, giving smaller compaction lengths near = 0 than those calculated here. This may have lead to the low focusing efficiency of the decompaction channel in their simulations.

For the purposes of the present research, the shear viscosity is taken to be a constant ₀, independent of temperature, pressure, porosity, stress, etc. This simplification is inconsistent with experiments on mantle deformation (e.g. Karato & Wu, 1993; Hirth & Kohlstedt, 2003) and excludes the possibility of plate-like behaviour of cold thermal boundary layers, as well as localization as a result of mechanical interactions between fluids and solids (e.g. Holtzman et al., 2003). Such behaviour is considered to be important for melt segregation in the Earth (e.g. Kelemen et al., 2002; Katz et al., 2004, 2006). However, the simulations described below are readily extended to more complex rheologies, and this extension is among the goals for future work.

The compaction length _c is the length-scale that emerges from nondimensionalization of the equations of magma dynamics (1)–(4). It is the scale over which perturbations to the compaction pressure decay away (Spiegelman, 1993b). As pointed out by Schmelling (2001), with a bulk viscosity proportional to ^–1 the compaction length can be considerably longer than previous estimates, where the bulk viscosity was taken to be approximately equal to . Applying the relation (24) for the bulk viscosity and assuming that 3_R0 >> 4,

(25)

This equation states that for a porosity of 1% and _R0=1, the compaction length is about a factor of 10 larger than estimated by McKenzie (1984).

	SOLUTIONS TO THE EQUATIONS

TOP ABSTRACT INTRODUCTION MODEL FORMULATION SOLUTIONS TO THE EQUATIONS DISCUSSION SUMMARY AND CONCLUSIONS Appendix A: Conservation of... Appendix B: Non... Appendix C: The enthalpy... REFERENCES

 
The equations (17)–(21) are in a form amenable to numerical solution. Below I describe the discretization approach and the numerical method used to solve the discrete system. Results of 1D simulations are presented and compared with semi-analytical calculations (Ribe, 1985a) as a benchmark of the simulation code and as a means to constrain certain parameters. Two-dimensional simulations were performed for a restricted set of parameters; results of these simulations are presented following the 1D results. Parameter values used in the constitutive and governing equations are given in Table 1.

Discretization and numerical solution
The governing equations are discretized on a staggered, Cartesian grid with matrix and fluid velocities located on cell boundaries and all other variables located at cell centers. The nonlinear system resulting from the discrete mass and momentum equations is separated from that of the enthalpy and bulk composition equations; the two are solved in a Picard iteration loop. Advection terms are handled using an upwind Fromm scheme with second order accuracy (Albers, 2000). Solution of each set of discrete equations is performed using a Newton–Krylov–Schwartz method provided by the Portable, Extensible Toolkit for Scientific Computation (PETSc, Balay et al., 2001, 2004). Details and references for this approach have been provided by Katz et al. (2007) for the equations of magma dynamics and by Katz & Worster (2008) for the Enthalpy Method.

Time-stepping of the enthalpy and composition equations is performed semi-implicitly using a Crank–Nicolson scheme. Although this scheme is unconditionally stable, the time-step size is limited, to preserve accuracy, to be close to the limit prescribed by the Courant–Friedrichs–Levy (CFL) condition. This limit is derived from the velocity of the magma, which can be orders of magnitude larger than that of the mantle matrix, depending on the permeability of the matrix and the buoyancy of the magma. For permeability in the range considered here, a grid spacing of 0.75 km and a domain of about 150 kmx70 km, the simulation of 0.5–2 Ma of model-time required about 40 h of clock-time on eight nodes of a cluster with one Intel® Xeon® (2.4 GHz, 1 GB RAM) processor per node.

Convergence of the simulations to an accurate solution for decreasing grid-spacing and time-step is not rigorously proven here. Katz & Worster (2008) performed simplified benchmark simulations of the Enthalpy Method and of thermal convection in a fixed porous medium and found excellent convergence with analytical or accepted solutions; much of the discretization details and code from that work are reused here. A further benchmark of the Enthalpy Method is performed below with a 1D upwelling column. In two dimensions, comparison of simulations at different spatial resolutions indicates that for grid-spacing smaller than 1 km, integrated results such as focusing distance do not vary systematically with grid-spacing. This result is shown below in Fig. 7a; it gives confidence that the simulations are convergent.

Upwelling column model
One-dimensional simulations of upwelling and melting mantle rock were performed as a benchmark for the simulations code. Ribe (1985a) considered the melting of a two-component mantle with complete solid solution and no thermal or chemical diffusion. He derived a simplified ordinary differential equation for the steady-state profile of temperature. From the temperature profile, Ribe (1985a) calculated the degree of melting, and the approximate porosity and fluid upwelling profiles, all of which can be compared directly with results of simulations. Neglecting diffusion, the steady-state, nondimensional governing equations become

(26)

(27)

(28)

(29)

(30)

where w_f and w_m are the vertical components of fluid and solid velocity, respectively. The matrix-buoyancy term in equation (19) has been dropped because matrix convection is not possible in 1D solutions. Boundary conditions include a fixed potential temperature of 1350°C at the bottom of the domain and a fixed solid upwelling rate, U₀ = 5 cm/yr, at the top of the domain. Incoming mantle is taken to be 12% of the less fusible component and 88% of the more fusible component. Other parameters are as given in Table 1.

Although it is possible to reduce equations (26)–(30) analytically to ordinary differential equations matching those solved by Ribe (1985a), I have solved them numerically to demonstrate the validity of the method and implementation. Figure 2 shows the comparison between results from these simulations and the curves calculated by Ribe (1985a). The correspondence is imperfect for porosity, fluid velocity and bulk composition because Ribe (1985a) neglected pressure gradients arising from compaction in his solution. Temperature and the degree of melting for a steady-state, equilibrium, 1D melting column are independent of the flow parameters, as demonstrated by Ribe (1985a), and as evident in Fig. 2a and b, where curves for all values of K₀ exactly coincide.

View larger version (40K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 2. Results of 1D simulations with the Enthalpy Method (continous curves) compared with semi-analytical calculations after Ribe (1985a) (dashed curves). Darker line colour represents a smaller value of K0; the black curve denotes the line for which K0 = 0. The other six curves represent K0 = 10–11, 10–10, 10–9, 10–8, 10–7 and 10–6 m2. The poor match between simulations and theory for porosity, fluid velocity and bulk composition is the result of the use of the ‘zero-compaction length’ approximation by Ribe (1985a). General agreement validates the numerical approach and implementation.

 
These results demonstrate that the Enthalpy Method and this implementation are valid for simulations of magma dynamics. They also show that the calibration of the phase diagram (i.e. the chosen values of T₀, T₁, and ) gives an amount of melting for 1350°C potential temperature mantle that is consistent with expectations (Langmuir et al., 1992). Figure 2 shows that the porosity and fluid upwelling rate depend on the choice of K₀ [in particular, and (Ribe, 1985a)]. Constraints on the permeability of mantle rock come mainly from experimental measurements of monomineralic samples (e.g. Wark & Watson, 1998), texturally equilibrated rocks (e.g. Faul, 1997), or from grain-scale models (e.g. Zhu & Hirth, 2003; Cheadle et al., 2004). These studies generally used an equation similar to (23) to fit their results. The permeability exponent n is typically estimated to be two or three. Wark et al. (2003) has shown that permeability laws derived for monomineralic systems can be applied to upper mantle assemblages. Given the variability of experimental results and the poor constraints on grain size in the mantle, other workers such as Faul (2001) have taken an inverse approach, seeking a permeability model that is consistent with geophysical and geochemical data. Results from the Mantle ELectromagnetic and Tomography (MELT) experiment (Forsyth et al., 1998) were interpreted by Faul (2001) based on reductions in P-wave speed (Faul et al., 1994) to indicate the presence of 1–2% melt in disk-shaped pores beneath the East Pacific Rise. Observed thorium disequilibria (see the Discussion below) suggest even smaller residual porosities. Simulations shown in Fig. 2 indicate that 10^–9 < K₀ < 10^–6 m² is the range of permeability constants that produces such a range in porosity.

Mid-ocean ridge model
In this section I describe 2D simulations of mid-ocean ridge melting and melt transport. The simulations combine melting, freezing, and porous flow of magma with a consistent determination of the matrix flow field as a result of both compaction and lithospheric motion. The Enthalpy Method is unchanged in higher spatial dimensions. The domain, shown schematically in Fig. 3, contains a vertical section of the ridge perpendicular to the ridge strike. It extends laterally from the ridge axis to a specified width and from the surface to a specified depth. Reflection boundary conditions on the vertical boundary beneath the axis enforce symmetry across the axis and prevent flow through the boundary. The other vertical boundary uses open, outflow conditions to minimize disturbance to the interior of the domain. The bottom boundary has fixed enthalpy corresponding to the desired potential temperature at zero porosity. The top boundary has a specified matrix velocity v_m = (U₀, 0) representing the moving tectonic plate. Away from the ridge, the top boundary has a fixed enthalpy corresponding to a temperature of 0°C. In contrast, within a region less than 6 km from the ridge, the boundary condition on enthalpy is insulating, H/z=0. This allows hot, upwelling mantle to reach the surface and provides a porous conduit for magma to leave the domain. Further details on boundary conditions are given in Fig. 3.

View larger version (17K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 3. Schematic diagram of 2D model domain with boundary conditions. Partial derivatives with respect to x, z are denoted by x, z. The horizontal and vertical components of the matrix velocity vm are denoted U and W, respectively. HC is the enthalpy corresponding to T = 0°C and = 0; HP is the enthalpy corresponding to the prescribed inflow temperature and = 0. C0 is the concentration of the incoming mantle.

 
The segment of the top boundary with an insulating boundary condition in enthalpy allows hot mantle to reach the surface and thus gives magma a route to leave the domain at the ridge. Unfortunately, this approach works only when the spreading rate U₀ is large enough that vertical advective heat transfer dominates lateral diffusive cooling beneath the ridge. In this paper I consider only one half-spreading rate, 5 cm/yr, that is large enough to avoid this problem. In similar calculations, Ghods & Arkani-Hamed (2001) specified a fixed, low temperature over the entire upper boundary and extracted melt from where it pooled below the thermal boundary layer. Such an approach is also possible for the present model but has not been implemented at present.

Below I describe the initial condition, time-evolution and steady state of a typical simulation, as well as the efficiency of melt focusing. Additional simulation results are considered in the Discussion, with predictions for crustal thickness and magmatic transit time as a function of key model parameters.

Initial condition
Because of the complexity of the equations, the lack of a good starting guess, and the possible presence of solitary waves, it has not been possible to solve the steady-state equations directly for the 2D mid-ocean ridge model. Instead, a fully time-dependent solution is computed; this solution may approach steady state. Because such an approach can be computationally expensive, it is important to choose an initial condition that minimizes the transient time in the model. One possible choice is to prescribe sub-solidus enthalpy throughout the domain and allow heat advection by the solid to establish the thermal ridge structure and melting regime. This is impractical, however, because once melting and melt transport begin, the advective time-scale for the magma reduces the time-step size to a small fraction of that determined by the advective time-scale for matrix motion.

A better approach was developed by Ghods & Arkani-Hamed (2000) and has been adapted for use here. To initialize the time-dependent simulation, the matrix velocity v_m is specified by the isoviscous corner-flow solution (Batchelor, 1967) for an incompressible fluid. The prescribed velocity field is then used to solve the energy equation (20) with K₀ = 0 to eliminate melt-segregation effects. This solution provides a map of (non-dimensional) enthalpy, porosity, temperature and phase compositions everywhere in the domain. The enthalpy is then reduced by subtracting off the fraction that is contained in latent heat, . Next, the porosity is set to zero and the bulk composition is set to the composition of the solid phase after the initial melting; temperature is left unchanged. The result is a melting region that is perched precisely on the solidus with a temperature structure very close to steady state.

Steady and quasi-steady solutions
There is an initial transient phase in all simulations during which porosity in the melting region increases from zero and melt begins to percolate upward. Melt produced directly under the ridge ascends vertically and never encounters the base of the cold lithospheric plate. In contrast, off-axis melts rise nearly vertically and accumulate at the depth of the solidus, which is a sloping boundary beneath the lithospheric thermal boundary layer, as shown in Fig. 4a. Pooling of magma beneath this boundary leads to a high-porosity layer that Sparks & Parmentier (1991) termed the decompaction layer. The slope of this layer directs the buoyant flow of magma toward the ridge. Porosity and permeability in the decompaction layer increase with time during the transient phase until a balance exists between (1) the flux of magma up the decompaction layer to the ridge, (2) the flux of magma from the melting region into the decompaction layer, and (3) the flux of magma from the decompaction layer into the cold lithosphere via under-plating (Sparks & Parmentier, 1991).

View larger version (47K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 4. Example simulation result for two simulations with K0 = 10–7 m2 and different values of R. Mesh spacing is 0.75 km. (a) Steady-state conditions for a simulation with bulk to shear viscosity ratio R = 100. The grey scale represents porosity, dashed lines represent instantaneous fluid streamlines, fine continuous lines represent isotherms and the bold continuous line marks the position of the solidus. The box labelled represents the subdomain used for integration in equation (31). (b) Conditions at 900 kyr in a simulation with R = 25. The porosity grey scale is clipped; porosity reaches 9% within the blob at 10 km depth and 25 km from the ridge. Lines as defined in (a). (c) Concentration of less fusible component in the solid for the same simulation as in (a). Lines are streamlines of matrix flow. (d) Concentration in the matrix for the same simulation as in (b). The white blob along the solidus has a concentration of 10 wt % of the less fusible component. White blobs in the top-right corner of (c) and (d) are remnants of the initial condition and should be ignored.

 
If the decompaction layer is morphologically stable, an approximately steady-state solution may emerge after the initial transient. The model time required to establish this steady state decreases with increasing background fluid velocity, w₀, as defined in Appendix B. This means that although simulations with larger magmatic flow velocities are more computationally intensive, such simulations also pass through the initial transient phase in less model time. In general, larger values of the bulk-to-shear viscosity ratio, _R, are associated with a larger compaction length, stability of the decompaction layer, and the existence of a quasi-steady state. Instability of the decompaction layer is considered below.

Figure 4c shows that steady-state melting and melt segregation deplete the solid mantle of the more easily fusible component. The lithosphere is slightly enriched relative to the solid within the compaction layer because of the under-plating that occurs at the top of the compaction layer. This reintroduces a measure of the more easily fusible component into the the solid phase.

Time-dependent solutions
Bulk viscosity plays an important role in determining the stability of the decompaction layer. Figure 4b shows a snapshot of a simulation in which _R=25, a factor of four smaller than that of figure 4a. Magma has ponded on the solidus boundary and formed narrow regions of high porosity up to 9%. Streamlines show that magma rises within the decompaction layer only until it reaches a high-porosity zone where it is trapped. A comparison of the solidus and isotherm contours in Fig. 4a and b demonstrates that the trapped magmatic flux creates a local disturbance to the temperature and bulk composition. The disturbance grows with time as more magma is trapped, modifying the composition of the lithosphere, as shown in Fig. 4c and d.

Magmatic focusing
Tracers are introduced into the flow to determine the path lines of fluid parcels. The particles have zero distribution coefficient and are thus advected at the velocity v_f, which is equal to the matrix velocity where the porosity is zero. The particles do not otherwise affect the simulation. For a simulation with K₀=10^–7 m² and _R=100, a swarm of tracers was introduced into the bottom of the domain after the simulation had reached steady state. Figure 5a is a histogram of the lateral distance from the ridge axis where these tracers cross into the bottom of the melting region. Black bars represent particles that arrive at the ridge axis and white bars represent particles that are frozen into the oceanic lithosphere. It is evident that for this simulation the focusing distance is 60 km. For tracers that are advected to the ridge, Fig. 5b shows the elapsed time for the tracer to travel from a given depth near the bottom of the melting region. The elapsed time is longer for particles that travel through the decompaction channel to reach the ridge axis.

View larger version (25K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 5. Results of tracer transport calculations in a simulation with K0 = 10–7 m2 and R = 100, after that simulation had reached approximate steady state (Fig. 4a). (a) A histogram of the lateral distance of tracer particles from the ridge axis when they enter the melting region. Bins are 2.5 km wide and are plotted according to their maximum distance from the axis. Black bars represent particles advected to the ridge axis; white bars represent particles frozen into the lithosphere. Total tracers per bin increase towards the ridge because faster matrix upwelling there introduces a larger flux of particles into the melting region, where they travel rapidly to the ridge axis. (b) Points represent transit time of fluid parcels from 54 km depth (just within the melting region) to the ridge axis as a function of their initial lateral distance from the axis. Dashed lines represent the transit time along fluid streamlines that start at the same depth.

 
Figure 5b compares the transit time along (tracer) path lines with the transit time along streamlines. If the flow field was in perfect steady state then the two would give the same result. Good agreement close to the ridge axis for K₀=10^–7 m² suggests that streamline transit times, which can be calculated instantaneously, are an acceptable substitute for expensive tracer transport calculations for simulations that approach a steady state. As expected, transit time decreases as the permeability constant K₀ is increased. For K₀ between 10^–6 and 10^–7 m², transit times for melts that originate beneath the ridge near the bottom of the melting region (54 km depth) range from 50 to 100 kyr. For comparison, the half-life of ²³⁰Th is about 75 kyr.

Estimation of the distance over which magma is focused to the ridge is possible using results in Fig. 5a. A similar result can be obtained by finding the width of the region that produces enough melt to balance the surface melt flux. In steady state one can calculate the balance of melt production and melt extraction at the ridge by integrating equation (1). Let us consider a rectangular region extending from the surface to a depth below the melting region and from the ridge axis to some distance, x, as shown in Fig. 4a. Using Gauss' theorem, the integral of (1) can be written

(31)

where is given by equation (23), represents the boundary of the region and is a unit normal vector pointing out of the region. This equation states that the net melting within the region must be balanced by magmatic flux into or out of the region.

The utility of equation (31) is clarified by expanding its two terms. The surface integral can be rewritten as the sum of the melt flux out of the domain at the ridge and the melt flux into the domain along the vertical boundary at x; there is no melt flux over the bottom boundary or the ridge-axis boundary. The volume integral of can be rewritten in terms of the integral of the melting (_mf > 0) and freezing (_fm < 0) rates. With these expansions, equation (30) becomes

(32)

where u and w are the horizontal and vertical components of the velocity of the fluid and z_D is the depth of the domain. Equation (32) is true for any value of x. At some distance , shown in Fig. 6, the flux out of at the ridge axis [term 2 in equation (32)] is balanced by melting (term 4) and the flux into through the vertical boundary at (term 1) is balanced by freezing (term 3). This distance is named the equivalent focusing distance (EFD), the distance within which the melt produced is equal to the melt delivered to the ridge. Because streamlines initiating farther from the ridge axis always travel through the compaction channel at shallower depth than those originating closer to the ridge axis (see Fig. 4), these distal melts are frozen into the lithosphere first. The EFD is therefore approximately equal to the maximum lateral distance over which magma is focused to the axis in simulations that have reached steady state.

View larger version (22K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 6. A plot of the terms in equation (32) as a function of the width x of the region . This result is from the simulation shown in Fig 4. The equivalent focusing distance is denoted . By definition, the cumulative melting curve is always positive and the cumulative freezing curve is always negative. The outflux curve is positive and it is constant for all x greater than the width of the ridge outflux region. The melting curve and the outflux curve must intersect because all the melt erupted at a ridge is produced in the melting region below it. Because equation (32) states that at any x the y-values of the four curves must sum to zero, the freezing curve must intersect the influx curve, at steady state, where the melting curve intersects the outflux curve. This point defines the equivalent focusing distance.

 

View larger version (13K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 7. Plots showing the variation of equivalent focusing distance, with model parameters. (a) EFD vs permeability constant K0 for three grid spacings, 0.5, 0.75 and 1 km. Evidently, for R = 100, a 1 km grid resolution is sufficient to resolve the integrated behaviour of the system. (b) EFD versus K0 for values of R varying from 12.5 to 200. At lower R the equivalent focusing distance is small, even for large permeability. This is due to magmatic trapping in mushy pools along the solidus, as shown in Fig. 4b. (Note the agreement between the EFD for R = 100, K0 = 10–7 m2 here and in the tracer focusing region in Fig. 5a). (c) EFD vs compaction length calculated with equation (25) at 30 km depth directly below the ridge axis. Each point represents the state of an independent simulation at the final output time of that simulation (i.e. at steady state, if applicable). The line is a least-squares fit to the data and has a y-intercept of 6 km, the prescribed axial half-width of the ridge. The simulations vary in K0 and R but all have a grid spacing of 0.75 km.

 
Figure 7b shows the EFD for an ensemble of simulations with different values of K₀ and the viscosity ratio, _R. For values of _R below 25, increased permeability cannot compensate for reductions in _R. Larger values of _R and K₀, however, yield a larger compaction length and more efficient melt focusing. Moreover, Fig. 7c shows that the EFD is correlated with the compaction length at 30 km depth below the ridge axis. Evaluation of the compaction length at this position avoids perturbations caused by the decompaction channel and gives a value representative of the behaviour of the simulation.

	DISCUSSION

TOP ABSTRACT INTRODUCTION MODEL FORMULATION SOLUTIONS TO THE EQUATIONS DISCUSSION SUMMARY AND CONCLUSIONS Appendix A: Conservation of... Appendix B: Non... Appendix C: The enthalpy... REFERENCES

 
An important check of a model of melt transport at a mid-ocean ridge is the predicted crustal thickness as a function of time. In this model, crustal thickness is computed as the integrated volumetric magma flux through the surface divided by the half-spreading rate times the ratio of crustal density to magmatic density (which is taken as unity here) and is shown in figure 8. Mid-ocean ridges with full spreading rates greater than 2 cm/yr have seismically measured crustal thickness of 5–7 km (Bown & White, 1994). In simulations, crustal thickness is sensitive to the melt production rate in the melting region and the efficiency of focusing of melt to the ridge. All simulations reported here use the same mantle potential temperature, composition and phase diagram, and hence the maximum degree of melting (20%) and melt production is roughly constant between them.

View larger version (37K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 8. Crustal thickness as a function of time. Each panel shows the evolution of crustal thickness from four simulations with K0 = 10–9, 10–8, 10–7 and 10–6 m2 and a grid spacing of 0.75 km. (a) R = 100, (b) R = 50, (c) R = 25, (d) R = 12:5. Simulations with R = 200 are not shown but display similar evolution to those with R = 100. The continuous grey lines for K0 = 10–7 m2 in (a) and (c) correspond to the simulations shown in Fig. 4.

 
Quasi-steady-state solutions produce crustal thickness that is either constant or oscillates around a constant value, as evident in Fig. 8. Oscillations may be due to low-amplitude, large-wavelength solitary waves in the melting region or the decompaction channel. Numerical simulations reported by Ghods & Arkani-Hamed (2000) also produced oscillations of crustal thickness, although of a much larger amplitude than those reported here. This may be due to their recipe for extracting melt at the ridge, which differs from the open top boundary in the present simulations.

The results in Fig. 8 show that predicted crustal thickness depends on the permeability and resistance to compaction of the mantle matrix. It is evident that smaller values of _R lead to oscillatory crustal thickness that is consistently below observed values. Both K₀ and _R exert an important control on melt transport within the decompaction layer that, in turn, controls the efficiency of magmatic focusing to the ridge. The streamlines in Fig. 4a and b demonstrate how a decrease in efficiency occurs in simulations. As noted by Spiegelman (1993c), melt is deflected into the compaction channel more efficiently if the width of the freezing interval around the solidus is small relative to the compaction length. Figure 9 shows that this is clearly the case when the permeability and bulk viscosity prefactors, K₀ and _R, are large.

View larger version (29K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 9. A map of the compaction length c from the same simulation as in Fig. 12a with K0 = 10–7 m2 and R = 100. It should be noted that although the bulk viscosity decreases in the sub-lithospheric decompaction layer, the effect of increased permeability is larger and leads to a factor of 5 increase in compaction length there. The black band around the melting region is where the porosity is very small, giving an extremely large bulk viscosity and hence an extremely large compaction length. Simpson (2008) has suggested that in the limit of vanishing porosity, compaction length should go to zero; in that work this was accomplished by regularizing the permeability and bulk viscosity functions.

 
The total focusing efficiency, averaged over a period from 250 kyr into each simulation until its end, is shown in Fig. 10 for a range of simulation parameters. These results can be compared directly with figure 6 of Sparks & Parmentier (1991). In general, total focusing efficiency is higher than their prediction for a given value of K₀/µ where the ratio of bulk to shear viscosity is greater than about 50. For the maximum degree of melting (20%) that is imposed in present simulations by the choice of mantle potential temperature, phase diagram parameters, spreading rate, etc., the total melt production can yield a crust of about 7 km thickness; hence an efficiency of greater than 70% is required to produce a crust thicker than 5 km. If, however, more melt were produced overall then the necessary efficiency would be smaller.

View larger version (20K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 10. Focusing efficiency and crustal production in simulations with various values of K0 and R. Each point represents the mean efficiency over a range of simulated time from 250 kyr after initiation of the simulation until the maximum time simulated (see Fig. 8). Compared with results reported by Sparks & Parmentier (1991) for K0/µ = 10–7 m2, the efficiency of focusing in this figure is higher for R = 100 and lower for R = 50 (see their figure 6, for a spreading half-rate of 3 cm/yr). The contribution of melts formed in the corner of the melting region not included within the domain of the present calculations is negligibly small. The R = 100 curve uses simulations at a grid resolution of 0.5 km; all others have a resolution of 0.75 km.

 
Figure 11 shows the composition of the evolving and steady-state crust for _R=100 and a range of K₀ in terms of the concentration of the high-temperature melting component. The variation is 5% around the mean, whereas crustal thickness varies by more than a factor of two. Higher permeability corresponds to a thicker crust with a larger fraction of the high-temperature-melting component. One-dimensional upwelling columns produce crust of a single composition, independent of permeability or other flow parameters (Ribe, 1985a). In 2D simulations, lateral melt focusing injects magma into the mantle column directly beneath the axis. This exotic magma equilibrates with the sub-ridge mantle by reactive melting. Greater focusing efficiency yields a larger flux of magma (and heat) through the mantle immediately beneath the ridge. Melting reactions increase the local degree of melting and deplete the solid of its low-T-melting component. Therefore, more melt focusing yields a more depleted bulk composition and a more depleted magma composition beneath the ridge axis at steady state. This is shown in Fig. 11a and explains the positive correlation in Fig. 11b between crustal thickness (magmatic flux) and depletion of that crust in the low-T-melting component.

View larger version (22K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 11. Crustal composition assuming perfect mixing of the magmatic flux through the top of the domain. Simulations have R = 100 and K0 ranging from 10–9 to 10–6 m2. (a) Secular trends in the concentration of the high-temperature melting component in the crust. (b) Steady-state crustal composition vs steady-state crustal thickness for the four simulations from (a).

 
Observations of uranium-series disequilibrium in young lavas provide another constraint on magmatic transport processes. If fractionation of uranium and thorium takes place mainly within the garnet stability field at 60 km depth and magmatic transport is by diffuse porous flow, then preservation of a ²³⁰Th excess at the surface requires mantle porosities below 1% and transport times of the order of a few half-lives (Spiegelman & Elliott, 1993). The present model does not include explicit calculations of radiogenic transport. It does, however, achieve ambient porosities and magma transit times that are small enough to be consistent with observed ²³⁰Th disequilibrium. Figure 12 shows mantle residual porosity predicted by simulations with different values of K₀ and _R. Consistent with 1D column model results presented above, permeability is a strong control on ambient porosity. Resistance to compaction provided by the bulk viscosity does not exert an important influence on residual porosity because in the melting region, where the melting rate varies slowly with respect to the compaction length, the zero-compaction length approximation (neglecting compaction stresses) is valid.

View larger version (21K): [in this window] [in a new window] [Download PowerPoint slide]   Fig. 12. Fine lines represent the porosity at 30 km depth and 10 km distance from the ridge axis as a function of K0 for different values of R. The bold continuous line is the porosity at 30 km depth in 1D simulations (see Fig 2). The bold dashed line is the curve = K0–1/3/(3 x 104).

 
²³¹Pa and ²²⁶Ra, with half-lives of only about 32 500 and 1600 years respectively, may represent a tighter constraint on magmatic transport if excesses are generated within the melting region and not the magma chamber (Stracke et al., 2003). According to models by Stracke et al. (2006), disequilibrium transport of magma with a velocity between 2 and 100 m/yr is required to match data from Iceland. Stracke et al. (2006) also argued for mixing of rapidly transported melts from beneath the ridge with melts transported over a greater transit time, perhaps those originating at a larger lateral distance from the ridge. Melt focusing along the base of the lithosphere might provide the mechanism for such mixing. Porosity in the decompaction channel, however, is predicted to be larger than that in the ambient melting region, and this will affect the degree of secular disequilibrium of the melts that pass through it. In particular, if the porosity of the decompaction channel is significantly larger than the distribution coefficients of the radiogenic nucleides then decay back to secular equilibrium will occur. Future work will incorporate direct calculations of U-series fractionation, transport and disequilibrium.

	SUMMARY AND CONCLUSIONS

TOP ABSTRACT INTRODUCTION MODEL FORMULATION SOLUTIONS TO THE EQUATIONS DISCUSSION SUMMARY AND CONCLUSIONS Appendix A: Conservation of... Appendix B: Non... Appendix C: The enthalpy... REFERENCES

 
The simulations described here combine the fluid-mechanical theory of McKenzie (1984) for mass and momentum conservation with the Enthalpy Method for modelling energy transport and the thermodynamics of melting. A pressure-dependent phase diagram describing a two-component system with full solid solution is employed as a representation of thermodynamic equilibrium. The equations are solved numerically in one and two dimensions and the results are analysed in terms of the efficiency of melt focusing.

Melt focusing at mid-ocean ridges may have important 3D characteristics, especially near the ends of offset segments of ridge. Carbotte et al. (2004) discovered a global pattern of asymmetry in axial depth across ridge transform faults and showed that it is correlated with the direction of ridge migration over the mantle in the hotspot reference frame. These observations motivated a model of asymmetric upwelling, melt production and 3D magmatic focusing by Katz et al. (2004). In that work, magmatic flow was not calculated explicitly but was parameterized according to inferences about the behaviour of magma beneath a migrating ridge with a transform fault (Magde & Sparks, 1997). A 3D extension of the current simulations would provide a basis for evaluating these inferences and exploring along-axis variations in melt supply. Such a simulation, however, would require a substantial increase in the number and/or speed of processors to maintain current simulation run-times.

The choice of constitutive equations has an important influence on the behaviour of the system. Simulations presented here use a constant shear viscosity and a bulk viscosity that varies with the inverse porosity. The later has a singularity when the porosity is zero, which is handled here such that the equations reduce to incompressible Stokes' flow in that limit. This is different from the approach of Ghods & Arkani-Hamed (2001), which used a non-singular (but temperature-dependent) relation for the bulk viscosity. The difference suggests a possible explanation for the relatively low focusing efficiency predicted by Ghods & Arkani-Hamed (2000): in their simulations, magma was able to migrate into the sub-solidus region above the decompaction layer and solidify, instead of being deflected by a gradient in compaction pressure. The temperature dependence of the bulk viscosity might be expected to enhance magmatic focusing. However, this is true only if a sharp change in matrix viscosity occurs at the temperature of the mantle solidus—not the case for experimentally constrained rheological parameters. Hence the incorporation of a temperature-dependent bulk viscosity into our model would probably not change the results significantly.

Temperature and porosity dependence of the shear viscosity might have interesting consequences for the flow field of the matrix phase. In particular, it will be instructive to study the effect of porosity-weakening viscosity in the context of the stress and deformation field of a ridge. There has been much speculation (e.g. Rabinowicz & Vignersse, 2004; Katz et al., 2006; Holtzman & Kohlstedt, 2007) on the relevance of a mechanical porosity-banding instability for the melting region beneath ridges. Although this instability is clearly active in experimental deformation of partially molten rocks (e.g. Holtzman et al., 2003), experiments show it emerges only at strains larger than unity. Flow within the region of partial melting beneath ridges may not achieve such large strains, or it may achieve large enough strains but over a time-scale that is long relative to the time-scale for annealing to textural equilibrium. Furthermore, although buoyancy-driven segregation of melt does not occur in experiments because of their small size and rapid deformation rate, it is clearly active in the mantle. Butler (2009) has shown that buoyancy will modify the signature of a porosity-banding instability, although it should not erase it altogether. An extension of the models described here will provide a basis for addressing these questions.

Assuming that the neglected dependences of shear and bulk viscosity are not leading-order controls on melt focusing, and hence that the current models capture the physics of melt focusing at ridges, we can draw conclusions about the magnitude of the parameters K₀ and _R. If the melting rates calculated here are within 10% of natural values such that the total melt production rate is approximately correct, then it is reasonable to rule out any combination of K₀ and _R that yield a crustal thickness smaller than 5 km or a focusing efficiency less than about 70%. From Fig. 10 we can see that this can be achieved for K₀ as small as 10^–8 m² if _R 100. If we take account of the observed ²³⁰Th excesses and assume that the current models capture the physics of magma migration (i.e. no channel or dike flow), then K₀10^–7 m² is required to give magma transport times that would preserve disequilibrium formed deep in the melting column. In this case, _R could be as low as about 25 and the bulk viscosity prefactor _R0 would be in good (although perhaps fortuitous) agreement with previous estimates (Batchelor, 1967; Bercovici et al., 2001; Hewitt and Fowler, 2008; Simpson, 2008).

Although porosity-induced matrix buoyancy is included in the full governing equations [ in equation (19)], it has been excluded from the equations that are solved numerically. Many researchers, beginning with Buck & Su (1989) and Scott & Stevenson (1989), have considered the possibility of active, porosity-driven upwelling beneath ridges. Melting and magmatic flow within an actively convecting sub-ridge mantle were modeled by Scott & Stevenson (1989) and Spiegelman (1996); the pattern of convection is similar to that modeled by Rabinowicz et al. (1984). These models prescribe isoviscous conditions throughout the domain, which may be a significant deficiency. In the natural ridge system the cold, rigid lithosphere above and the viscosity reduction caused by ambient porosity in the melting region could have important effects. Porosity-induced buoyancy, along with temperature and porosity-dependent viscosity can be included in simulations that extend those presented here.

Reactive channelization of magmatic flow is a chemical instability that is thought to occur in the melting region beneath ridges (Kelemen et al., 1995). High-flux magma channels have implications for uranium-series disequilibrium, as well as trace element distribution and variability in erupted lavas (Spiegelman & Kelemen, 2003). Past models have considered reactive channelization in the context of a static or uniformly upwelling mantle without internal deformation. Combining tectonic-scale models, such as those described here, with calculations of channelized melt transport would illuminate the behaviour of a channelized system in a deforming mantle with a cold, impermeable lid. Calculations of geochemical transport layered on top of such simulations would provide an input for investigations of magma mixing and fractionation in magma chambers beneath volcanoes (e.g. Maclennan, 2009).

Some workers have interpreted geological and geochemical evidence to suggest that melt migration is extremely rapid. For example, Maclennan et al. (2002) estimated magmatic ascent rates of >50 m/yr based on eruption rates in Iceland after the the last glacial period. It is unlikely that porous magmatic flow, even if it is channelized, could reproduce such rapid melt migration. If these estimates prove correct, an alternative fluid-mechanical model of melt migration involving flow in cracks and dikes may be required.

Cracks and dikes are clearly responsible for melt transport at shallow depths across the lithosphere. Off-axis magmatism may tap pools of magma trapped along the solidus boundary, as in Fig. 4b. Whether over-pressures in such pools are sufficient to initiate hydrofracture is a question for future work. Moreover, it is interesting to note that for often-cited estimates of the compaction length of 10–1000 m, melt focusing is predicted to be inefficient and melt pooling on the solidus is expected. If melting beneath ridges is due to passive upwelling and if melt focusing occurs according to the physics of porous flow and compaction, then the compaction length remains a key parameter in explaining the high efficiencies that are expected.

Melting in subduction zones is more complex than beneath ridges because of the effects of water. Subduction-related magmatism has therefore received less attention. Reactive flow may be a useful framework for considering magma genesis in arcs, as it is for ridges (Grove et al. 2006). Along these lines, simulations by Katz (2005) of infiltration of reactive, aqueous fluid into the mantle wedge predict that channelization of fluid or melt occurs above the slab where temperature increases upward. These simulations consider a static mantle, however, and thus neglect advection of porosity and solid depletion through the wedge. A model of melt transport in a deforming mantle wedge was described by Cagnioncle et al. (2007). This work invoked a melt-focusing mechanism based on that of Sparks & Parmentier (1991). However, by neglecting compaction stresses and freezing, the model of Cagnioncle et al. (2007) excluded the possibility of actually resolving such behaviour. A subduction-zone model implemented using the approach described here, augmented by the inclusion of water as a thermodynamic component, would provide a means for investigating reactive melting as well as magmatic focusing in arcs.

The number of possible extensions to the work described here indicates the utility of the Enthalpy Method. It provides a clean, transparent approach to incorporating thermochemical complexity into models of magma dynamics. Such complexity is required to address physical phenomena such as magmatic focusing in a tectonic context. This power comes at a cost of making the questionable assumption of thermodynamic equilibrium between magma and the mantle matrix. Coupling this approach with disequilibrium geochemical transport, however, may provide an effective tool for making geochemical predictions that are testable against measurements of lava chemistry.

	Appendix A: Conservation of bulk enthalpy and composition

TOP ABSTRACT INTRODUCTION MODEL FORMULATION SOLUTIONS TO THE EQUATIONS DISCUSSION SUMMARY AND CONCLUSIONS Appendix A: Conservation of... Appendix B: Non... Appendix C: The enthalpy... REFERENCES

 
The derivation begins by considering an arbitrary Eulerian volume V bounded by a surface V containing both fluid and matrix material. The change in time of the energy contained in this volume is related to the flux across its boundary. Internal sources of heat, radioactive decay, and viscous dissipation are not included in this derivation. Rewriting equation (7),

(A1)

where we have also neglected the contribution of kinetic energy because, for the mantle, it is very small relative to the other terms. Applying Gauss' theorem and allowing the volume to shrink in size until it is small compared with the length-scale of variation of any phase averaged property but large relative to the grain size of the mantle, we can write a differential form of equation (A1) as

(A2)

Using the sum of conservation of mass equations (1) and (2) and expanding the term describing the divergence of the enthalpy flux in (A2) gives

(A3)

This equation can also be derived by recasting equation (A39) of Mckenzie (1984) as an equation for conservation of enthalpy and neglecting viscous dissipation and internal heating (J. Rudge, personal communication).

Assuming that the phase densities _f and _m are constant, the sum of equations (1) and (2) can be written

(A4)

where = _m – _f. Furthermore, for changes between equilibrium states, the total differential of enthalpy can be expressed in terms of total differentials of temperature and pressure as (e.g. Denbigh, 1981)

(A5)

where is the coefficient of thermal expansion. Substituting (A4) and (A5) into (A3) and rearranging gives

(A6)

where L=h_f–h_m is the latent heat per unit mass. The specific heat, thermal diffusivity and thermal expansivity have been assumed equal between phases.

Applying the extended Boussinesq approximation, _m =_f = , and assuming a static pressure gradient, P=g,

(A7)

which states that changes in volumetric enthalpy are due to advection of latent heat, advection of sensible heat, adiabatic pressure changes and diffusion of sensible heat.

Equation (A7) can be simplified by substitution of the potential temperature, which implicitly accounts for adiabatic changes. Potential temperature is defined as

(A8)

Using equation (A8) and taking the z-coordinate to be depth into the Earth, equation (A7) becomes

(A9)

In this equation, terms proportional to (g/c_P) << 1 have been neglected.

Mass conservation for one thermochemical component can be written for the same Eulerian volume as considered above. If the concentration of this component is C_f in the fluid phase and C_m in the matrix phase,

(A10)

Applying the same reasoning as for the equation (7) and taking _f =_m =, this equation can be rewritten in differential form as

(A11)

where C=C_f + (1–) C_m.

	Appendix B: Non-dimensionalization

TOP ABSTRACT INTRODUCTION MODEL FORMULATION SOLUTIONS TO THE EQUATIONS DISCUSSION SUMMARY AND CONCLUSIONS Appendix A: Conservation of... Appendix B: Non... Appendix C: The enthalpy... REFERENCES

 
Following Katz et al. (2007) we introduce a decomposition of the fluid pressure into three parts

(A12)

where gz is the lithostatic pressure, is the compaction pressure and P is the remaining, ‘dynamic’ pressure. Moreover, we introduce the following dimensionless variables:

The velocity scale w₀ is given by

(A13)

and the length scale, , is

(A14)

Using equations (3), (4) and (A12), substituting nondimensional variables, and dropping primes gives the equations governing matrix shear and compaction:

(A15)

(A16)

(A17)

Here, is the normalized gravity vector and = – 2/3 is the compaction viscosity. The nondimensional fluid velocity is

(A18)

To nondimensionalize the conservation of energy equation we introduce the nondimensional temperatures

(A19)

where T₀ and T₁ are the minimum and maximum melting temperatures over all compositions at atmospheric pressure. The nondimensional bulk enthalpy is defined as H = c_P TH' with T = T₁–T₀. Using nondimensional variables in equations (8) and (9) and dropping primes gives

(A20)

(A21)

Here, _t is a partial derivative with respect to time, is the adiabatic parameter, is the Stefan number, Pe_T = w₀/ is the thermal Peclet number, and is the compositional Peclet number. Equation (12), which describes the bulk enthalpy, becomes

(A22)

where *=T₀/T.

	Appendix C: The enthalpy method for systems with more than two components

TOP ABSTRACT INTRODUCTION MODEL FORMULATION SOLUTIONS TO THE EQUATIONS DISCUSSION SUMMARY AND CONCLUSIONS Appendix A: Conservation of... Appendix B: Non... Appendix C: The enthalpy... REFERENCES

 
The Enthalpy Method can be generalized to N > 2 components in two phases (fluid and matrix). To do so requires mathematical descriptions of the solidus and liquidus surfaces in N dimensions. We seek a closure condition for porosity and temperature as functions of bulk enthalpy and composition. In this case, however, instead of solving for the concentration of a single component in each phase we must solve for all the N concentrations. It is convenient to define N-component vectors C_f, C_m, and C for the fluid, matrix, and bulk concentrations respectively. These vectors have the property that

(A23)

where a superscript i represents the ith component of a vector. This property allows us to trivially obtain the Nth component of any vector if we know the other N – 1 components. As in the case of a two-component system where one PDE for bulk composition is solved, for an N-component system N – 1 partial differential equations are used, each identical to equation (9) but each for a different component. The equation that describes the conservation of energy is unchanged from equation (8).

The problem is then to find closure conditions for 2N+2 unknowns: , T, C_f, and C_m as functions of H and C. To do so, in addition to (A23), we have the equations

(A24)

(A25)

(A26)

and equation (12), which defines bulk enthalpy as the sum of sensible and latent heat. Here f_S and f_L are N-component vector functions that give the matrix and fluid phase compositions on the solidus and liquidus as a function of pressure, temperature and bulk composition—they define the phase diagram.

Combining equations (12) and (A24)–(A26) with the constraint on bulk composition from equation (A23) gives an expression for porosity:

(A27)

A value of from equation (A27) is substituted into equation (12) to obtain the temperature, which is then be used in equations (A25) and (A26) to calculate the phase compositions. These values of , T, C_f, and C_m are used in the discretized partial differential equations for bulk enthalpy and composition to construct a residual vector suitable for use in Newton's; method (see Katz et al., 2007).

The Enthalpy Method may be simplified by prescribing a phase diagram in terms of (H,P,C) instead of (T,P,C). To my knowledge, neither case has been implemented for more than two thermochemical components and both raise the difficulty of reasonably parameterizing the surfaces described by equations (A25) and (A26).

	ACKNOWLEDGEMENTS

 
Guidance in derivation of the enthalpy equation by M. G. Worster is gratefully acknowledged. Discussions with J. Maclennan and B. Wood, and comments by J. Rudge and J. Neufeld helped to clarify the manuscript. An initial review by M. Spiegelman helped to prepare the manuscript for submission. Detailed reviews by P. Asimow, O. rámek and M. Rabinowicz were useful in revising it. Further thanks are due to Asimow for follow-up discussions and Rabinowicz for his even-handedness. This work was performed with financial support from the US National Science Foundation's; International Research Fellowship Program and an Academic Fellowship from the Research Councils UK.

	FOOTNOTES

 
Present address: Department of Earth Science, University of Oxford, Parks Road, Oxford OX1 3PR, UK

---

*Corresponding author. E-mail: richard.katz{at}earth.ox.ac.uk

	REFERENCES

TOP ABSTRACT INTRODUCTION MODEL FORMULATION SOLUTIONS TO THE EQUATIONS DISCUSSION SUMMARY AND CONCLUSIONS Appendix A: Conservation of... Appendix B: Non... Appendix C: The enthalpy... REFERENCES

 
Aharonov E, Whitehead JA, Kelemen PB, Spiegelman M. Channeling instability of upwelling melt in the mantle. Journal of Geophysical Research—Solid Earth (1995) 100:20433–20450.[CrossRef]

Ahern JL, Turcotte DL. Magma migration beneath an ocean ridge. Earth and Planetary Science Letters (1979) 45:115–122.[CrossRef][Web of Science]

Albers M. A local mesh refinement multigrid method for 3-d convection problems with strongly variable viscosity. Journal of Computational Physics (2000) 160(1):126–150. doi: 10.1006/jcph.2000.6438.[CrossRef][Web of Science]

Alexiades V, Solomon AD. Mathematical Modeling of Melting and Freezing Processes. (1993) New York: Hemisphere.

Asimow PD, Stolper EM. Steady-state mantle–melt interactions in one dimension: I. Equilibrium transport and melt focusing. Journal of Petrology (1999) 40(3):475–494.[CrossRef][Web of Science]

Asimow PD, Hirschmann MM, Stolper EM. An analysis of variations in isentropic melt productivity. Philosophical Transactions of the Royal Society of London, series A (1997) 355(1723):255–281.[CrossRef]

Balay S, Buschelman K, Gropp WD, Kaushik D, Knepley M, McInnes LC, Smith BF, Zhang H. PETSc Web page (2001) http://www.mcs.and.gov/petsc.

Balay S, Buschelman K, Gropp WD, Kaushik D, Knepley M, McInnes LC, Smith BF, Zhang H. PETSc users manual. In: Argonne National Laboratory, Technical Report ANL-95/11—Revision 2.1.5 (2004).

Barcilon V, overa O. olitary waves in magma dynamics. Journal of Fluid Mechanics (1989) 204:121–133.[CrossRef][Web of Science]

Batchelor GK. An Introduction to Fluid Mechanics (1967) Cambridge: Cambridge University Press.

Bear J. Dynamics of Fluids in Porous Media (1972) Amsterdam: Elsevier.

Bercovici D, Ricard Y. Energetics of a two-phase model of lithospheric damage, shear localization and plate-boundary formation. Geophysical Journal International (2003) 152:581–596.[CrossRef][Web of Science]

Bercovici D, Ricard Y, Schubert G. A two-phase model for compaction and damage 1. General theory. Journal of Geophysical Research—Solid Earth (2001) 106:8887–8906.[CrossRef]

Bouhifd MA, Besson P, Courtial P, Gerardin C, Navrotsky A, Richet P. Thermochemistry and melting properties of basalt. Contributions to Mineralogy and Petrology (2007) 153:689–698. doi:10.1007/s00410-006-0170-8.[CrossRef][Web of Science]

Bown JW, White RS. Variation with spreading rate of oceanic crustal thickness and geochemistry. Earth and Planetary Science Letters (1994) 121:435–449.[CrossRef][Web of Science]

Buck WR, Su WS. Focused mantle upwelling below mid-ocean ridges due to feedback between viscosity and melting. Geophysical Research Letters (1989) 16(7):641–644.[Web of Science]

Butler S. The effects of buoyancy on shear-induced melt bands in a compacting porous medium. In: Physics of the Earth & Planetary Interiors (2009) (Submitted).

Cagnioncle A.-M, Parmentier EM, Elkins-Tanton LT. Effect of solid flow above a subducting slab on water distribution and melting at convergent plate boundaries. Journal of Geophysical Research—Solid Earth (2007) 112. doi: 10.1029/2007JB004934.

Carbotte SM, Small C, Donnelly K. The influence of ridge migration on the magmatic segmentation of mid-ocean ridges. Nature (2004) 429:743–746.[CrossRef][Web of Science][Medline]

Cheadle MJ, Elliott MT, McKenzie D. Percolation threshold and permeability of crystallizing igneous rocks: The importance of textural equilibrium. Geology (2004) 32:757–760.[Abstract/Free Full Text]

Choblet G, Parmentier EM. Mantle upwelling and melting beneath slow spreading centers: effects of variable rheology and melt productivity. Earth and Planetary Science Letters (2001) 184:589–604.[CrossRef][Web of Science]

Cooper RF. Differential stress-induced melt migration—an experimental approach. Journal of Geophysical Research—Solid Earth (1990) 95(B5):6979–6992.[CrossRef]

Cordery MJ, Morgan JP. Melting and mantle flow beneath a midocean spreading center. Earth and Planetary Science Letters (1992) 111(2–4):493–516.[CrossRef][Web of Science]

Cordery MJ, Morgan JP. Convection and melting at midocean ridges. Journal of Geophysical Research—Solid Earth (1993) 98(B11):19477–19503.[CrossRef]

Daines MJ, Kohlstedt DL. Influence of deformation on melt topology in peridotites. Journal of Geophysical Research-Solid Earth (1997) 102:10257–10271.[CrossRef]

Denbigh KG. The Principles of Chemical Equilibrium. (1981) 4th edn. Cambridge: Cambridge University Press.

Elliott T, Spiegelman M. Melt migration in oceanic crustal production: a U-series perspective. In: The Crust, Treatise on Geochemistry—Holland HD, Rudnick RL, Turekian KK, eds. (2003) Vol. 3. Oxford: Elsevier–Pergamon. 465–510.

Faul UH. Melt retention & segregation beneath mid-ocean ridges. Nature (2001) 410:920–923.[CrossRef]

Faul UH. Permeability of partially molten upper mantle rocks from experiments and percolation theory. Journal of Geophysical Research—Solid Earth (1997) 102:10299–10311.[CrossRef]

Faul UH, Toomey DR, Waff HS. Intergranular basaltic melt is distributed in thin, elongated inclusions. Geophysical Research Letters (1994) 21:29–32.[CrossRef][Web of Science]

Forsyth DW, Scheirer DS, Webb SC, Dorman LM, Orcutt JA, Harding AJ, Blackman DK, Morgan JP, Detrick RS, Shen Y, Wolfe CJ, Canales JP, Toomey DR, Sheehan AF, Solomon SC, Wilcock WSD. Imaging the deep seismic structure beneath a mid-ocean ridge: The MELT experiment. Science (1998) 280:1215–1218.[CrossRef][Web of Science][Medline]

Fowler A. A mathematical model of magma transport in the asthenosphere. Geophysical and Astrophysical Fluid Dynamics (1985) 33:63–96.[CrossRef]

Ghiorso M. Algorithms for the estimation of phase stability in heterogeneous thermodynamic systems. Geochimica et Cosmochimica Acta (1994) 58(24):5489–5501.[CrossRef][Web of Science]

Ghiorso MS, Sack RO. Chemical mass-transfer in magmatic processes. 4. a revised and internally consistent thermodynamic model for the interpolation & extrapolation of liquid–solid equilibria in magmatic systems at elevated temperature and pressure. Contributions to Mineralogy and Petrology (1995) 119:197–212.[Web of Science]

Ghods A, Arkani-Hamed J. Melt migration beneath mid-ocean ridges. Geophysical Journal International (2000) 140:687–697.[CrossRef][Web of Science]

Grove TL, Chatterjee N, Parman SW, Medard E. The influence of H₂O on mantle wedge melting. Earth and Planetary Science Letters (2006) 249:74–89. doi:10.1016/j.epsl.2006.06.043.[CrossRef][Web of Science]

Hall CE, Parmentier EM. Spontaneous melt localization in a deforming solid with viscosity variations due to water weakening. Geophysical Research Letters (2000) 27(1):9–12.[CrossRef][Web of Science]

Hewitt IJ, Fowler AC. Partial melting in an upwelling mantle column. In: Philosophical Transactions of the Royal Society of London, Series A (2008) doi:10.1098/rspa. 2008.0045.

Hirth G, Kohlstedt D. Rheology of the upper mantle and the mantle wedge: A view from the experimentalists. In: The Subduction Factory.—Eiler J, ed. (2003) Geophysical Monograph, American Geophysical Union, pp. 83–105.

Holtzman BK, Kohlstedt DL. Stress-driven melt segregation and strain partitioning in partially molten rocks: Effects of stress and strain. Journal of Petrology (2007) 48:2379–2406. doi:10.1093/petrology/egm065.[Abstract/Free Full Text]

Holtzman BK, Groebner NJ, Zimmerman ME, Ginsberg SB, Kohlstedt DL. Stress-driven melt segregation in partially molten rocks. Geochemistry, Geophysics, Geosystems (2003) 4. article number 8607.

Jull M, Kelemen PB, Sims K. Consequences of diffuse and channelled porous melt migration on uranium series disequilibria. Geochimica et Cosmochimica Acta (2002) 66:4133–4148.[CrossRef][Web of Science]

Karato S, Wu P. Rheology of the upper mantle—a synthesis. Science (1993) 260:771–778.[Abstract/Free Full Text]

Katz RF. The deep roots of volcanoes: models of magma dynamics with applications to subduction zones. (2005) New york: Columbia University. PhD thesis.

Katz RF, Worster MG. Simulation of directional solidification, thermochemical convection, and chimney formation in a Hele–Shaw cell. Journal of Computational Physics (2008) doi:10.1016/j.jcp.2008.06.039.

Katz RF, Spiegelman M, Carbotte SM. Ridge migration, asthenospheric flow and the origin of magmatic segmentation in the global mid-ocean ridge system. Geophysical Research Letters (2004) 31:L15605.[CrossRef]

Katz RF, Spiegelman M, Holtzman B. The dynamics of melt & shear localization in partially molten aggregates. Nature (2006) 442:676–679. doi:10.1038/nature05039.[CrossRef][Web of Science][Medline]

Katz RF, Knepley MG, Smith B, Spiegelman M, Coon ET. Numerical simulation of geodynamic processes with the Portable Extensible Toolkit for Scientific Computation. Physics of the Earth and Planetary Interiors (2007) 163:52–68. doi:10.1016/j.pepi.2007.04.016.[CrossRef][Web of Science]

Kelemen P, Parmentier M, Rilling J, Mehl L, Hacker B. Thermal convection of the mantle wedge. Geochimica Cosmochimica Acta (2002) 66(15A):A389–A389.

Kelemen PB, Shimizu N, Salters VJM. Extraction of mid-ocean-ridge basalt from the upwelling mantle by focused flow of melt in dunite channels. Nature (1995) 375(6534):747–753.[CrossRef]

Kelemen PB, Hirth G, Shimizu N, Spiegelman M, Dick HJB. A review of melt migration processes in the adiabatically upwelling mantle beneath oceanic spreading ridges. Philosophical Transactions of the Royal Society of London, Series A (1997) 355(1723):283–318.[CrossRef]

Langmuir CH, Klein E, Plank T. Petrological systematics of mid-oceanic ridge basalts: constraints on melt generation beneath ocean ridges. In: Mantle Flow and Melt Generation at Midocean Ridges, Geophysical: Monograph,—Phipps Morgan J, Blackman DK, Sinton JM, eds. (1992) American Geophysical Union 71. 183–280.

Maclennan J. Concurrent mixing and cooling of melts under Iceland. In: Journal of Petrology (2009) (submitted).

Maclennan J, Jull M, McKenzie D, Slater L, Gronvold K. The link between volcanism and deglaciation in Iceland. Geochemistry, Geophysics, Geosystems (2002) 3:1062. doi:10.1029/2001GC000282.[CrossRef]

Magde LS, Sparks DW. Three-dimensional mantle upwelling, melt generation, and melt migration beneath segment slow spreading ridges. Journal of Geophysical Research—Solid Earth (1997) 102:20571–20583.[CrossRef]

McKenzie D. The generation and compaction of partially molten rock. Journal of Petrology (1984) 25(3):713–765.[Abstract/Free Full Text]

McKenzie D. ²³⁰Th–²³⁸U disequilibrium and the melting processes beneath ridge axes. Earth and Planetary Science Letters (1985) 72(2–3):149–157. doi:10.1016/0012821X(85)90001-9.[CrossRef][Web of Science]

Phipps Morgan J. Melt migration beneath mid-ocean spreading centers. Geophysical Research Letters (1987) 14(12):1238–1241.[Web of Science]

Oertling AB, Watts RG. Growth of and brine drainage from NaCl–H₂O freezing: A simulation of young sea ice. Journal of Geophysical Research—Oceans (2004) 109. doi:10.1029/2001JC001109.

Rabinowicz M, Ceuleneer G. The effect of sloped isotherms on melt migration in the shallow mantle: a physical and numerical model based on observations in the Oman ophiolite. Earth and Planetary Science Letters (2005) 229:231–246.[CrossRef][Web of Science]

Rabinowicz M, Vigneresse JL. Melt segregation under compaction and shear channeling: Application to granitic magma segregation in a continental crust. Journal of Geophysical Research—Solid Earth (2004) 109:B04407. doi:10.1029/2002JB002372.[CrossRef]

Rabinowicz M, Nicolas A, Vigneresse JL. A rolling-mill effect in the asthenosphere beneath oceanic spreading centers. Earth and Planetary Science Letters (1984) 67:97–108.[CrossRef][Web of Science]

Ramberg H. Temperature changes associated with adiabatic decompression in geological processes. Nature (1971) 234:539–540.

Ribe NM. The generation and composition of partial melts in the Earth's mantle. Earth and Planetary Science Letters (1985a) 73:361–376.[CrossRef][Web of Science]

Ribe NM. The deformation and compaction of partial molten zones. Geophysical Journal of the Royal Astronomical Society (1985b) 83:487–501.[Web of Science]

Richardson CN. Melt flow in a variable viscosity matrix. Geophysical Research Letters (1998) 25(7):1099–1102.[CrossRef][Web of Science]

Saal AE, van Orman JA. The Ra-226 enrichment in oceanic basalts: Evidence for melt–cumulate diffusive interaction processes within the oceanic lithosphere. Geochemistry, Geophysics, Geosystems (2004) 5:Q02008. doi:10.1029/2003GC.[CrossRef]

Schmelling H. Partial melting & melt segregation in a convecting mantle. In: Physics and Chemistry of Partially Molten Rocks—Bagdassarov N, Laporte D, Thompson AB, eds. (2000) Dordrecht: Kluwer Academic Publisher. 141–178.

Schubert G, Turcotte DL, Olsen P. Mantle Convection in the Earth and Planets (2001) Cambridge: Cambridge University Press.

Scott DR, Stevenson DJ. Magma ascent by porous flow. Journal of Geophysical Research—Solid Earth (1986) 91:9283–9296.[CrossRef]

Scott DR, Stevenson DJ. A self-consistent model of melting, magma migration and buoyancy-driven circulation beneath mid-ocean ridges. Journal of Geophysical Research—Solid Earth (1989) 94:2973–2988.[CrossRef]

Simpson G. The Mathematics of magma migration (2008) . PhD thesis, Columbia University, New York.

Sims KWW, Goldstein SJ, Blichert-Toft J, Perfit MR, Kelemen P, Fornari DJ, Michael P, Murrell MT, Hart SR, DePaolo DJ, Layne G, Ball L, Jull M, Bender J. Chemical and isotopic constraints on the generation and transport of magma beneath the East Pacific Rise. Geochimica et Cosmochimica Acta (2002) 66:3481–3504.[CrossRef][Web of Science]

Sparks DW, Parmentier EM. Melt extraction from the mantle beneath spreading centers. Earth and Planetary Science Letters (1991) 105(4):368–377.[CrossRef][Web of Science]

Spiegelman M. Linear analysis of melt band formation by simple shear. Geochemistry, Geophysics, Geosystems (2003) 4. doi:10.1029/2002GC000499.

Spiegelman M. Flow in deformable porous-media. part 1 Simple analysis. Journal of Fluid Mechanics (1993a) 247:17–38.[CrossRef][Web of Science]

Spiegelman M. Flow in deformable porous media. part 2. Numerical analysis—The relationship between shock waves and solitary waves. Journal of Fluid Mechanics (1993b) 247:39–63.[CrossRef][Web of Science]

Spiegelman M. Physics of melt extraction: theory, implications, and applications. Philosophical Transactions of the Royal Society of London, Series A (1993c) 342:23–41.

Spiegelman M. Geochemical consequences of melt transport in 2-d: The sensitivity of trace elements to mantle dynamics. Earth and Planetary Science Letters (1996) 139:115–132.[CrossRef][Web of Science]

Spiegelman M, Elliott T. Consequences of melt transport for uranium series disequilibrium in young lavas. Earth and Planetary Science Letters (1993) 118:1–20.[CrossRef][Web of Science]

Spiegelman M, Kelemen PB. Extreme chemical variability as a consequence of channelized melt transport. Geochemistry, Geophysics, Geosystems (2003) 4. doi:10.1029/2002GC000336.

Spiegelman M, McKenzie D. Simple 2-D models for melt extraction at mid-ocean ridges and island arcs. Earth and Planetary Science Letters (1987) 83:137–152.[CrossRef][Web of Science]

Spiegelman M, Kelemen PB, Aharonov E. Causes and consequences of flow organization during melt transport: the reaction infiltration instability in compactible media. Journal of Geophysical Research—Solid Earth (2001) 106(B2):2061–2077.[CrossRef]

Stevenson DJ. Spontaneous small-scale melt segregation in partial melts undergoing deformation. Geophysical Research Letters (1989) 16:1067–1070.[Web of Science]

Stevenson DJ, Scott DR. Mechanics of fluid-rock systems. Annual Reviews of Fluid Mechanics (1991) 23:305–339.[CrossRef]

Stracke A, Zindler A, Salters JM, McKenzie D, Grönvold K. The dynamics of melting beneath Theistareykir, northern Iceland. Geochemistry, Geophysics, Geosystems (2003) 4(10)). doi:10.1029/2002GC000347.

Stracke A, Bourdon B, McKenzie D. Melt extraction in the Earth's mantle: Constraints from U–Th–Pa–Ra studies in oceanic basalts. Earth and Planetary Science Letters (2006) 244:97–112. doi:10.1016/j.epsl.2006.01.057.[CrossRef][Web of Science]

Tackley PJ. Mantle convection and plate tectonics: Toward an integrated physical and chemical theory. Science (2000) 288:2002–2007.[Abstract/Free Full Text]

rámek O, Ricard Y, Bercovici D. Simultaneous melting and compaction in deformable two-phase media. Geophysical Journal International (2007) 168:964–982. doi:10.1111/j.1365-246X.2006.03269.x.[CrossRef][Web of Science]

Wark DA, Watson EB. Grain-scale permeabilities of texturally equilibrated, monomineralic rocks. Earth and Planetary Science Letters (1998) 164:591–605.[CrossRef][Web of Science]

Wark DA, Williams CA, Watson EB, Price JD. Reassessment of pore shapes in microstructurally equilibrated rocks, with implications for permeability of the upper mantle. Journal of Geophysical Research—Solid Earth (2003) 108. doi:10.1029/2001JB001575.

Zhu WL, Hirth G. A network model for permeability in partially molten rocks. Earth and Planetary Science Letters (2003) 212:407–416.[CrossRef][Web of Science]

CiteULike    Connotea    Del.icio.us    What's this?

This article has been cited by other articles:

				Y. Liang and E. M. Parmentier A Two-Porosity Double Lithology Model for Partial Melting, Melt Transport and Melt-rock Reaction in the Mantle: Mass Conservation Equations and Trace Element Transport J. Petrology, January 7, 2010; (2010) egp086v1. [Abstract] [Full Text] [PDF]

				D. S. H. King, M. E. Zimmerman, and D. L. Kohlstedt Stress-driven Melt Segregation in Partially Molten Olivine-rich Rocks Deformed in Torsion J. Petrology, January 1, 2010; 51(1-2): 21 - 42. [Abstract] [Full Text] [PDF]

This Article Abstract FREE Full Text (PDF) All Versions of this Article: 49/12/2099    most recent egn058v1 Alert me when this article is cited Alert me if a correction is posted Services Email this article to a friend Similar articles in this journal Alert me to new issues of the journal Add to My Personal Archive Download to citation manager Request Permissions Google Scholar Articles by Katz, R. F. Search for Related Content GeoRef GeoRef Citation Social Bookmarking          What's this?

Online ISSN 1460-2415 - Print ISSN 0022-3530

Copyright © 2010 Oxford University Press

Oxford Journals Oxford University Press

• Site Map
• Privacy Policy
• Frequently Asked Questions

Other Oxford University Press sites: Oxford University Press Oxford Journals China Oxford Journals Japan American National Biography Booksellers' Information Service Children's Fiction and Poetry Children's Reference Corporate & Special Sales Dictionaries Dictionary of National Biography Digital Reference English Language Teaching Higher Education Textbooks Humanities International Education Unit Law Medicine Music Online Products Oxford English Dictionary Reference Rights and Permissions Science School Books Social Sciences Very Short Introductions World's Classics