## Abstract

Decompression melting of hot upwelling rock in the mantle creates a region of partial melt comprising a porous solid matrix through which magma rises buoyantly. Magma transport and the compensating matrix deformation are commonly described by two-phase compaction models, but melt production is less often incorporated. Melting is driven by the necessity to maintain thermodynamic equilibrium between mineral grains in the partial melt; the position and amount of partial melting that occur are thus thermodynamically determined. We present a consistent model for the ascent of a one-dimensional column of rock and provide solutions that reveal where and how much partial melting occurs, the positions of the boundaries of the partial melt being determined by conserving energy across them. Thermodynamic equilibrium of the boundary between partial melt and the solid lithosphere requires a boundary condition on the effective pressure (solid pressure minus melt pressure), which suggests that large effective stresses, and hence fracture, are likely to occur near the base of the lithosphere. Matrix compaction, melt separation and temperature in the partially molten region are all dependent on the effective pressure, a fact that can lead to interesting oscillatory boundary-layer structures.

## 1. Introduction

Beneath mid-ocean ridges and in isolated mantle plumes, upwelling mantle rock undergoes partial melting. The production of melt and its subsequent migration are responsible for the creation of new oceanic lithosphere and for the occurrence of submarine and subaerial volcanism. Understanding the physical processes involved in this generation and movement of magma is important in helping to understand the thermal and chemical properties of the oceanic crust, and in describing how magma can come to be emplaced in magma chambers and feed volcanic eruptions. This paper aims to provide a consistent mathematical model for these governing physical processes.

Partial melting occurs as a result of decompression melting—hot crystalline rock ascending adiabatically from lower in the mantle finds itself above the pressure-dependent melting temperature and begins to melt (figure 1). As the rock continues to rise, it forms a two-phase mixture of solid and melt, which undergoes continual melting until, near the surface, the temperature drops and the rock solidifies. Thermodynamically, the melt forms preferentially at the intersections of individual mineral grains that can be expected to form an interconnected network (McKenzie 1984). The melt is therefore able to move through the porous solid matrix and, being less dense, will rise buoyantly.

The whole situation can be described by the general theories of two-phase flows with mass exchange between the phases (Drew 1983; Bercovici *et al*. 2001), and many authors have proposed equations to model partially molten material (Turcotte & Ahern 1978; Ahern & Turcotte 1979; McKenzie 1984; Fowler 1985; Ribe 1985; Scott & Stevenson 1986; Spiegelman 1993; Bercovici *et al*. 2001). These differ in some specifics, but have the same general form, and have been put to a variety of uses. Surprisingly little attention, however, has been given to posing and solving a full model for the partial melting process, which must include, for example, not only the governing equations but also consistent boundary conditions.

The first attempts to model the process were by Turcotte & Ahern (1978), who considered Darcy flow through a deforming solid matrix with the melting prescribed by an equilibrium relationship between temperature, pressure and melt fraction. The principal assumption they make is that the pressure in melt and solid are the same, stating that viscous deformation of the grains will readily occur over short length scales (less than 100 m) to quickly equalize pressures. Effectively this allows the matrix to freely compact as the melt migrates upwards through it.

It is now realized that this compaction is in fact *due* to the difference in pressure between the phases, which must therefore be accounted for in any model of the process, as is generally the case for other two-phase flows and particularly, for example, in soil mechanics. Following the widely used form of the equations suggested by McKenzie (1984), this pressure difference is commonly parameterized in terms of a *bulk viscosity ζ*, related in some way to the normal shear viscosity *η*_{s}. The bulk viscosity essentially describes how easily a two-phase material can be squeezed to extract the melt and, in the context of melt generation, will depend strongly on the melt fraction *ϕ*. A number of arguments, based on microscopic models, suggest that *ζ*∼*η*_{s}/*ϕ* (Batchelor 1967; Fowler 1985; Sleep 1988; Bercovici *et al*. 2001).

Much of the work following McKenzie's original equations has neglected the melting process and concentrated on the melt migration and the matrix motion when melt is somehow already *in situ* (Spiegelman & McKenzie 1987; Spiegelman 1993; Spiegelman *et al*. 2001). Fully describing the onset of partial melting and placing this within the wider context of mantle circulation clearly requires the thermodynamics to play a major part, and the often neglected energy equation must therefore be included. The partial melt occupies only part of the ascending mantle and its position should be found by considering the thermodynamics; the boundaries are ‘free boundaries’, in the sense that their location is unknown *a priori* and must be found as part of the solution.

The free boundary nature of this problem was actually identified in the early work of Ahern & Turcotte (1979), who locate the depth at which melting begins by considering the surrounding temperature field and ensuring continuity of temperature gradients. Fowler (1989) describes the procedure to uncouple the free boundary location from the interesting dynamics, which all occur *within* the partial melt region, a method that we will employ later in this paper. Besides this work, there has been no comparable attempt to solve the partial melt problem in its proper context until the recent work by Šrámek *et al*. (2007), whose aims were similar to our own. Based on the two-phase formalism of Bercovici *et al*. (2001), they considered the total melting of an ascending column of rock, correctly noting the free boundary at which melting starts, although attempting to determine its position without suitable jump conditions.

In this paper we will review the equations describing the partial melt region including the temperature as a principal variable and discuss appropriate boundary conditions to apply to the partial melting that occurs beneath a mid-ocean ridge. We provide solutions for the case of a one-dimensional upwelling column of mantle such as might be appropriate directly beneath the ridge. The situation differs from that considered by Šrámek *et al*. (2007), in that we apply thermal boundary conditions at the Earth's surface that acts as a ‘lid’ on the partial melt below. This causes a solidified boundary layer (the lithosphere) to form and its position, just as the position of the onset of melting, must be determined by conserving energy across the boundary. The requisite dynamical boundary condition to apply to the partial melt at this upper boundary has previously been addressed only by Fowler (1989, 1990*a*) and is not entirely obvious; we show below that the requirement that the boundary itself is in thermodynamic equilibrium necessitates that the pressures in the solid and the melt be equal there. Given the complexity of the modelling, the emphasis here is on physics rather than chemistry and we treat the rock as a single thermodynamic component. We hope that the additional effects when the mantle is more realistically treated as multi-component can be incorporated at a later stage.

Given a prescribed mantle ascent rate (this we consider to be externally driven by the large-scale mantle circulation) and the temperature of the upwelling rock, our results allow for the full determination of where and how much partial melting occurs.

## 2. Mathematical model

The situation we consider is shown in figure 2; the upwelling of mantle rock causes a region *D* to become partially molten while the region *D*^{+}, which will be referred to as ‘subsolidus’, is solid. The boundary ∂*D* is unknown and our model therefore consists of equations to describe the dynamics within *D* and *D*^{+}, together with the boundary conditions across ∂*D*, at the fixed Earth surface and in the deep mantle.

### (a) Conservation equations for a compacting partially molten region

The equations we use to model the partially molten region are for Darcy flow through a deforming solid matrix (Fowler 1985). More complicated equations giving a systematic description of two-phase flow are possible in which many extra surface effects can be included (Drew 1983; Fowler 1985, 1990*a*; Bercovici *et al*. 2001; Šrámek *et al*. 2007); but once the usual simplifying assumptions are made concerning interactive drag coefficients and the partitioning of surface forces, these are often reduced to Darcy's law. Since such details are available elsewhere and we do not wish to confuse the issue by introducing more physics than necessary, we start from the outset by assuming that melt flow obeys Darcy's law.

#### (i) Mass and momentum conservation

Conservation of mass for each phase is expressed by(2.1)(2.2)in which ** u** and

**are the velocities of melt and solid, respectively;**

*V**ρ*

_{l}and

*ρ*

_{s}are the densities of melt and solid, respectively; and

*m*(kg m

^{−3}s

^{−1}) is the melt rate converting solid into melt.

We write Darcy's law in the form(2.3)where the permeability *k*_{ϕ} is(2.4)Here, *a* is a typical grain size; *b* is a tortuosity factor; and *η*_{l} is the viscosity of the melt.

Along with Darcy's law to describe the relative motion of the two phases, we require the conservation of total momentum, which we write as(2.5)in which we take the following constitutive laws for the stress tensor ** σ** in melt and solid, based on the assumption that the fluid supports negligible deviatoric stress:(2.6)(2.7)where

*p*

_{l}and

*p*

_{s}are the pressure in melt and solid, respectively, and

*η*

_{s}is the viscosity of the solid rock.

#### (ii) Compaction

As in any two-phase flow, we must close the problem by prescribing some relation between the pressures of the phases. This is commonly done using a bulk viscosity,(2.8)but this expression is more properly viewed as the *definition* of the bulk viscosity, which we should therefore derive. Expression (2.8) follows from a microscopic model of the deformation of individual ‘tubules’ in the porous matrix (Nye 1953; Sleep 1988); if we consider these as cylinders of radius *a*, the walls of the cylinder will enlarge by melting and close down by viscous creep at a rate , thus(2.9)Associating the area of these tubules with the porosity, we have(2.10)and now combining with (2.2) and ignoring the small term proportional to *ϕ* (this argument is only appropriate with the assumption of small porosity) gives (2.8) in the form(2.11)

Similar expressions can be derived from alternative microscopic models; Batchelor (1967) derived the same expression for the bulk viscosity by considering energy dissipation when a two-phase material is compressed, and Šrámek *et al*. (2007) also recovered (2.8) with a more in-depth discussion of the interface thermodynamics.

#### (iii) Energy conservation

In order to model the melting process and prescribe the melting term *m* in (2.1) and (2.2), we must consider energy conservation. Since the phase change is all important, we separate out the internal energy of the solid *e*_{s} and melt *e*_{l}, and conservation of these requires(2.12)where to avoid complications we assume that the thermal conductivity , specific heat capacity *c* and thermal diffusivity *κ* are the same in each phase. The source term *Ψ* is the work done by viscous forces.

Using the conservation of mass and momentum in the usual way and the definition of latent heat(2.13)(*v*=1/*ρ* is the specific volume), we can rewrite (2.12) as(2.14)Here *β* is the thermal expansion coefficient and *Φ* is the viscous dissipation, given with the rheologies above, by(2.15)

#### (iv) Local thermodynamic equilibrium

We make the assumption that the partially molten region is in local thermodynamic equilibrium—that is, we assume that the flow of heat across the solid–melt interface is virtually instantaneous and given the local stress field the interfacial temperature will quickly reach equilibrium. Our assumptions concerning the pressure difference mean that the stresses do not equilibrate so readily, so local thermodynamic equilibrium means that the temperature will depend on the interfacial stress which, since , is the liquid pressure *p*_{l} (Kamb 1961).

Thermodynamic equilibrium requires the continuity of free energy across the interface, thus we require(2.16)where is the difference in specific Gibbs free energy; is the specific enthalpy difference; and is the specific entropy difference. By considering variations in the pressure *p*_{l} and melting temperature *T*_{S} from a reference state, we have the equilibrium Clapeyron condition(2.17)in which *Γ* is the usual Clapeyron slope. The subscript ‘S’ here refers to the *solidus* temperature; more generally, this will also depend on the composition of the rock, but we consider only one component for simplicity, so that the solidus depends only on pressure.

### (b) Subsolidus rock

Outside the partial melt region *D*, the solidified rock obeys simpler and more standard equations. The continuity and momentum equations are(2.18)(2.19)(2.20)while the energy equation, including the same terms as (2.14), is(2.21)

### (c) Boundary conditions

Two types of boundary conditions are required to consider the problem set out in figure 2; prescribed conditions at the ‘fixed’ boundaries—the surface and the deep mantle—and conservation laws across the interface ∂*D* between subsolidus and partially molten rock. The former will simply prescribe the temperature at the surface and of the upwelling mantle, but we delay their discussion until the next section.

#### (i) Conservation laws

Across the boundary ∂*D*, we must apply the same conservation laws as in the partial melt and subsolidus regions. Working from the integral form of the conservation laws in (2.1), (2.2), (2.5) and (2.12), we derive (see Fowler 1985 for instance) the following ‘jump’ conditions between values on either side of the boundary:(2.22)(2.23)(2.24)Here ^{+} refers to subsolidus variables, while the subscripts *s* and *l* refer to values in the partial melt; ** n** is the outward pointing normal to the boundary; and

**is the velocity of the boundary. denotes the jump in the quantity from partial melt to subsolidus regions. Ignoring the small deviatoric stress in (2.23), these can be rearranged to give the following conditions:(2.25)(2.26)(2.27)We must also require the continuity of temperature across the boundary,(2.28)**

*v*#### (ii) Thermodynamic equilibrium of the boundaries

Conditions (2.25)–(2.28) tell us something about the velocity, temperature and pressures at ∂*D*, but are not sufficient to close the problem; condition (2.26) can provide only one condition on the pressures *p*_{s} and *p*_{l} and since these are different we still need another condition to relate them. We also need a condition on the melt fraction *ϕ*; it seems intuitive to assume that the melt fraction at the onset of melting is 0 so that the necessary condition is *ϕ*=0 on the lower part of ∂*D* (this may be defined specifically as that part of ∂*D* where ). We can in fact derive this condition, and the extra condition on the pressure, by the requirement that the boundary itself is in thermodynamic equilibrium; we have previously assumed a local thermodynamic equilibrium within the partial melt, where the temperature is thus related to the local stress conditions. At the boundary ∂*D*, we also suppose that macroscopic thermodynamic equilibrium must be achieved between the partially molten and subsolidus regions, by ensuring the continuity of the average Gibbs free energy across the boundary.

We can write the average free energy on either side of the boundary (see Fowler 1990*a*) as(2.29)In the reference state in which , so continuity requires(2.30)which is the local equilibrium condition (2.16). Perturbations *δp*_{s}, *δp*_{l}, *δp*^{+}, *δT* and *δϕ* to pressures, temperature and melt fraction must maintain this equilibrium and, noting that for local equilibrium, we therefore require(2.31)Since in the reference state *p*_{s}=*p*^{+}, we find the thermodynamic boundary condition *p*_{s}=*p*^{+}. From (2.26), this provides the extra condition on ∂*D*,(2.32)Note that this condition requires that *either* the melt fraction *ϕ* is zero *or* the pressure difference *p*_{s}−*p*_{l} is zero. In practice, we expect the former to apply to the lower part of ∂*D* where the melting begins and the latter to apply to the upper boundary. From (2.27), this implies a continuous temperature gradient at the lower boundary, but solidification and a jump in temperature gradient at the upper boundary.

We are forced to assume here that melt will solidify at the upper boundary of the partial melt region. That this is not necessarily the case is manifestly true; magma erupts from the surface and is known to be emplaced in magma chambers within the lithosphere. However, the processes that allow this to happen are not entirely clear; certainly, it seems that magmafracturing and transport up pre-existing conduits within the lithosphere play a major part, but how these are connected to the partial melt zone below is very much open to debate. Some form of localization of the melt flow, or tendency for the partial melt itself to undergo fracture, seem to be required.

In the absence of any definite idea of what precisely does happen at the top of the partial melt region, we look for the simplest thermodynamically viable solution, and therefore make the naive assumption that all the melt solidifies—solidification of this sort is referred to as *underplating* of the lithosphere.

### (d) Non-dimensionalization

We expect to prescribe the mantle velocity that results from larger scale mantle convection, and our solutions will focus only on the one-dimensional column of mantle directly beneath a ridge or plume axis. We label the surface of the Earth *y*=0, with *y* being the vertical coordinate. If no partial melting occurred, the mantle would rise steadily, the pressure profile would be lithostatic and the melting temperature would be given instead of (2.17), by the *lithostatic solidus*,(2.33)The depth *y*_{m} here is a reference depth at which the lithostatic pressure is *p*_{m} and the lithostatic solidus is *T*_{m}. The adiabatic geotherm that the ascending rock would follow would in some places be above the lithostatic solidus (figure 1), and this is (approximately) the region that must therefore undergo partial melting. We therefore choose the reference depth *y*_{m} to be the depth at which the lithostatic solidus and adiabatic geotherm first intersect—this should be near to (though as we shall see, not exactly) the depth at which melting starts. We take this depth to be known, which effectively means that we prescribe the temperature of the upwelling rock.

Somewhere close to *y*_{m} partial melting begins, and somewhere close to the surface *y*=0, the cold surface temperature is noted and the partial melting region ends. We label these positions *y*_{b} and *y*_{a}, respectively, as shown in figure 2. As discussed above, the exact positions of these boundaries must form part of the solution, being found from the conservation laws in §2*c*. However, the interesting dynamics are all contained *within* the partial melt region that therefore attracts the majority of our attention; for this reason, we define another vertical coordinate so that the partial melt occupies . By taking a good guess for this depth *l*, we can find the solution for the partial melt region, and then use the jump conditions at *y*_{b} and *y*_{a} to precisely locate these boundaries that then define the actual depth *l*. In this way, we can effectively separate out the solution within the partial melt region and the determination of the free boundaries. A similar procedure was outlined by Fowler (1989).

The relevant length scale is the depth *l* (*a priori* we have only a guess at exactly what it should be, but it can be revised later) to which all lengths are scaled. The temperature within the partial melt follows the solidus, which will be approximately given by (2.33); we thus define and scale .

We take the mantle velocity scale *V*_{0} to be prescribed—it can be inferred by measuring the rate of spreading at mid-ocean ridges—and the melting scale follows from balancing the Clapeyron slope with the latent heat consumption in (2.14). The melt velocity is determined principally by buoyancy in (2.3), and the melt fraction then follows from balancing divergence with melt production in (2.1). Thus, we scale(2.34)With the reference point *y*_{m}, *p*_{m} and *T*_{m} defined above, we scale the solid pressure . We expect the melt pressure *p*_{l} to follow a similar scale, but in fact it is the pressure difference or *effective* pressure(2.35)which is crucial to the compaction dynamics, and we therefore use it as a primary variable. The appropriate scale for this effective pressure comes from balancing terms in the compaction relation (2.11),(2.36)

These balances define the following scales:(2.37)

With these scales, we rewrite the equations for the partial melt region in the non-dimensional form(2.38)(2.39)(2.40)(2.41)The compaction relation becomes(2.42)and the energy equation is(2.43)The solidus (2.17), to which the temperature in the partial melt is confined, is(2.44)

Several new non-dimensional parameters have been defined as(2.45)*r* is the density ratio; *λ* is the ratio of adiabatic to solidus temperature gradients; *μ* is the ratio of temperature variation across the partial melt to absolute temperature; *Pe* is the usual Péclet number; and *ν* is the ratio of gravitational potential energy to latent heat; *St* is the usual Stefan number, the ratio of latent to sensible heat, which also turns out in this problem to give the ratio of solid mass flux to melt mass flux—further emphasizing that the whole system is thermodynamically driven; *ϵ* is the ratio of typical velocities of solid and melt; *δ* is the ratio of effective pressure gradients to buoyancy; and *δ*_{s} is the ratio of effective pressure gradients to the lithostatic pressure gradient.

Table 1 shows typical values of the various constants, which we use to estimate the size of the scales and non-dimensional parameters. Typical estimates of the rate of mid-ocean ridge spreading suggest that new crust is generated at a rate between 1 and 10 cm yr^{−1}. We take the mantle velocity scale consistent with this range. The depth *l* of the partial melting region is yet to be found exactly but will be given approximately by −*y*_{m}, for which we take the nominal value 50 km.

With the values given in table 1, we find(2.46)and the parameters are(2.47)

There is a certain degree of uncertainty and variability in many of the parameters used here. Properties of the mantle and melt such as density and viscosity depend considerably on composition; rhyolitic magma is many orders of magnitude more viscous than basaltic magma for example, but somewhere between 1 and 100 Pa s is probably an average value for basalt (Ahern & Turcotte 1979). Density and viscosity also vary with temperature (this is what causes mantle convection in the first place), but such effects will be ignored for the purposes of this study. We take the values in (2.47) to be representative of typical mantle conditions, but bear in mind that there may be significant variations from these. Many of the parameters are small, which will allow approximate analytic solutions to be found.

With the same scales and parameters, the boundary conditions on ∂*D* become(2.48)(2.49)(2.50)(2.51)and(2.52)

We can now also consider the boundary conditions for the whole problem. We assume that the limiting vertical velocity exists and is known. Then from the dimensionless version of (2.21), noting that *μ* is small, the adiabatic temperature variation of the deep mantle requires(2.53)The temperature *T*_{s} at the surface *y*=0 is known and therefore we have(2.54)where can be expected to be large and negative, being measured with respect to the reference melting temperature *T*_{m}.

To summarize our model, (2.38)–(2.44) provide the equations for the partial melt region *D*, dimensionless versions of (2.18)–(2.21) provide the equivalent equations for the subsolidus regions *D*^{+}, (2.48)–(2.52) provide the conditions across the boundary ∂*D* and (2.53) and (2.54) are the fixed thermal boundary conditions. The equations for the partial melt region are essentially the same as those proposed by Fowler (1985) and also by McKenzie (1984), when the porosity dependence of the bulk viscosity is realized. In comparing with McKenzie's original formulation, note that the *compaction length* is given here by(2.55)The equations also contain the same information as those used by Šrámek *et al*. (2007), but are written in terms of the effective pressure *N* rather than fluid pressure *p*_{l}. Our simplification of the equations will be slightly different from theirs because we make use of the fact that the melt fraction is everywhere small by neglecting terms of relative order *ϕ*_{0}.

The equations will be further simplified in our analysis of the next section by making the Boussinesq approximation *r*=1, by setting *λ* to zero (we do not expect the adiabatic effects to be important), and by neglecting the terms of order *ϵδ* in the momentum equations. This last approximation is appropriate because the bulk viscosity is considerably larger than shear viscosity when the melt fraction is small; in the momentum equations (2.40) and (2.41), this means matrix stress gradients are less important than effective pressure gradients by a factor *ϵ*. This highlights the dynamical importance of the effective pressure and justifies Fowler's earlier assumption (Fowler 1985, 1989, 1990*b*) that the solid pressure (from (2.41)) is approximately lithostatic (Šrámek *et al*. 2007).

## 3. Steady one-dimensional solutions

In one dimension, we have ** V**=

*W*

**and**

*k***=**

*u**w*

**. With the approximations**

*k**λ*=0,

*r*=1 and

*δ*

_{s}=0 and ignoring terms of order

*ϕ*

_{0}, the partial melt equations become(3.1)(3.2)(3.3)(3.4)(3.5)We have defined the new parameter(3.6)which we do not neglect, on the basis that this diffusive term may be important near the boundaries; with our typical values

*P*is

*O*(1). Equation (3.1) is Darcy's law, (3.2) and (3.3) are mass conservation equations, (3.4) is the compaction relation and (3.5) is the energy equation, in which the remaining physics are latent heat consumption, heat advection and heat conduction. Equations (3.1)–(3.5) are to be solved on the domain 0<

*z*<1, and we have the boundary conditions(3.7)(3.8)The remaining jump conditions will be used later to find the exact position of the boundaries. Note that (3.4) is elliptic for

*N*, and we might therefore expect a condition on

*N*at

*z*=0, but since it is degenerate as

*ϕ*→0, it seems conditions (3.7) and (3.8) are sufficient.

### (a) Boundary-layer solutions

From (3.3) and (3.7), we find immediately(3.9)and noting that and we are neglecting the terms of order *ϕ*_{0}, the remaining equations reduce to give(3.10)

(3.11)

#### (i) Outer solution

From (2.47), we have typical values *ϵ*∼0.03, *δ*∼0.03 and *P*∼0.1. The steady outer solution obtained by taking *ϵ,δ*→0 satisfies(3.12)It turns out that in order to find a sensible boundary-layer solution we should choose this outer solution to satisfy the condition *ϕ*=0 at *z*=0, although we will still require a boundary layer there since this suggests *N*→∞ as *z*→0. The outer solutions are therefore(3.13)

Šrámek *et al*. (2007) describe this as the ‘Darcy solution’ since the buoyancy force is balanced by viscous drag. The effective pressure is able to freely adjust to allow the necessary compaction of the matrix, which balances the melt divergence—the melt flow is essentially decoupled from the matrix deformation, and exactly the same behaviour is therefore found in Ahern & Turcotte's (1979) results. Near the boundaries, however, the effective pressure cannot adjust in the same way and the relation between melt and matrix flow is less straightforward.

#### (ii) Boundary layer at *z*=0

Near *z*=0, guided by the outer solution (3.13), we write(3.14)Then(3.15)(3.16)with the boundary and matching conditions,(3.17)Assuming and , we have the single ordinary differential equation for ,(3.18)It is clear that depending on the value of at , solutions to this equation may either blow up at finite (if is large enough) or decrease to 0 at finite (if is small enough). The behaviour varies monotonically with the value of between these extremes and there must be a unique value, , such that the solution matches up with the outer solution as . Numerical solution gives .

With this initial value for *N,* we no longer satisfy *ϕ*=0 at *z*=0, but have . As suggested by the first term in (3.15), we need to rescale *z* again to find the smaller inner boundary region in which *ϕ* goes to 0. Writing , the equations become(3.19)with as . The solution is straightforward,(3.20)

In earlier work on this problem, the boundary layer at the bottom of the melt region has been called the *compaction layer* (McKenzie 1984; Fowler 1985; Ribe 1985), a somewhat misleading name since compaction occurs everywhere within the region and is actually restricted in this bottom layer by the large bulk viscosity; the melt pressure here is almost hydrostatic and melt production balances the convective derivative of *ϕ* in (3.2). This corresponds to Šrámek *et al*.'s (2007) ‘viscogravitational equilibrium’.

#### (iii) Boundary layer at *z*=1

At *z*=1 there must be another boundary layer to satisfy *N*=0. We write(3.21)and defining(3.22)rewrite the equations as(3.23)(3.24)with matching conditions(3.25)

Now by taking *δ*→0, a first integral of (3.24) produces the coupled set of first-order ordinary differential equations,(3.26)for which we are required to find a trajectory joining *N*=0 to as *Z* goes from 0 to ∞.

Represented on a phase plane these equations are shown in figure 3, where we define the nullclines *K*=0 and *G*=0 and the line *V*=0 on which *ϕ*_{Z} becomes infinite. Provided , we find that the fixed point is a saddle point and there is one unique trajectory that reaches it from *N*=0, with *ϕ* and *N* varying monotonically. The implications when will be discussed in §4.

### (b) Numerical solutions

Our simplified equations (3.7)–(3.10) can also be solved numerically. To find the solution, we actually solve the time-dependent version of the equations, which we allow to evolve to a steady state. Equation (3.11) is treated as a quadrature for *N* given *ϕ* and this is used to step the *ϕ* solution forward in time using (3.10). We use a uniform grid on which (3.11) is discretized to second-order accuracy.

A solution found in this way is shown in figure 4. In this we can clearly see the general parabolic profile of the melt fraction and its inverse relation with the effective pressure (3.13), the boundary layer at *z*=1 of width and the inner boundary layer at *z*=0 of order . Less obvious is the outer part of this boundary layer of width , although the limiting values and can be seen.

### (c) Free boundary location

We have now found solutions for the melt fraction and effective pressure (and therefore also temperature, melt and matrix velocities) within the partially molten region *assuming* that we knew the size of this region. As has already been discussed, its position and size (indeed whether partial melting occurs at all) must be found by considering the temperature of the surrounding rock.

In the one-dimensional situation we consider, and with the same approximations as before, the equations (2.18)–(2.21) for the subsolidus region simply tell us that the velocity is constant , the pressure is lithostatic and the dimensionless temperature therefore satisfies the steady-state equation(3.27)This is to be solved on and , and the relevant boundary conditions are (2.53) and (2.54) along with the jump conditions (2.50) and (2.51) (with in the steady state),(3.28)(3.29)(3.30)(3.31)These six boundary conditions provide all the information we need to solve for the temperature profile within the two disjoint subsolidus regions, with the extra two conditions determining the position of the two boundaries *y*_{a} and *y*_{b} (non-dimensionally, , so locating *y*_{a} is equivalent to fixing the length scale *l*).

The solutions of (3.27) for the temperature in the subsolidus regions are exponentials; applying all six boundary conditions (3.28)–(3.31) results in the simultaneous equations for *y*_{a} and *y*_{b},(3.32)(3.33)The parameters and the scalings for the variables in (3.32) and (3.33) are themselves dependent on the unknown length scale *l*, and, since , these two equations may be rewritten as one nonlinear equation for *l*.

The procedure we adopt is as follows: we take a guess *l*^{*} at the depth of the partial melt region and define the non-dimensional parameters *δ*^{*}, *Pe*^{*} and *ϵ*^{*} and variable scales as in (2.37) using this length scale. With these values, we find the solutions *ϕ*^{*}(*z*) and *N*^{*}(*z*) numerically, as above. Then, writing *X*=*l*/*l*^{*}, so that, for example, *Pe*=*Pe*^{*}*X*, (3.32) and (3.33) combine to give an equation for *X*(3.34)in which all the ^{*} variables are known. This can be solved to find *X* and therefore the true depth scale *l*. Provided our original estimate of the length scale was good, the solutions within the partial melt still hold with the true non-dimensional solutions being given by and .

With the length scale known, the depth of the onset of melting now follows directly from (3.33), with the values of *N*_{z}(0) and *N*(0) coming from the partial melt solution. In fact, our analytic boundary-layer solution allows us to find *y*_{b} exactly; from (3.18), we had and . In terms of dimensional variables, (3.33) therefore becomes(3.35)We see that the partial melt region may begin at a shallower depth than that predicted by the lithostatic solidus due to heat conduction in the rock below, or at a greater depth due to the depression of the solidus temperature from its lithostatic value as a result of the effective pressure. In fact, these two effects may have counterbalancing effects—in terms of our dimensionless parameters, *y*_{b}<*y*_{m} if , and with the values in (2.37) this is the case.

The full solutions for the temperature of the ascending mantle column and the behaviour of the partial melt region are shown in figure 5. The predicted thickness of the solidified lithosphere is approximately 2 km, considerably less than the estimated thickness beneath mid-ocean ridges. This results from the one-dimensional assumption that the rock continues to ascend all the way to the surface, whereas in reality the cold rigid lithosphere will move sideways, driven by the convective motion of the overlying plate.

The model equations in §2 apply equally well to this more realistic two-dimensional situation, but the solutions are rather more complicated. The varying mantle velocity alters the stress gradients driving melt flow, and the boundaries *y*_{a} and *y*_{b} will depend on the lateral coordinate. The same technique to determine the location of the boundaries will therefore not work and an alternative method of solution needs to be found.

Also shown in figure 6 is the predicted region of partial melt when the parameters are changed; as the mantle ascent rate *W*_{0} is reduced the size of the partial melt region decreases until, if *W*_{0} is too small, there is no partial melting at all. The temperature of the ascending rock, which determines where it intersects the solidus, also influences the size of the partially molten region.

## 4. Discussion

### (a) Thermodynamic boundary layer

The thermodynamic condition *N*=0 at the upper boundary of the partial melt required a boundary layer in which the pressure has to adjust rapidly. This was described by the phase plane in (3.26) and figure 4, which admits a unique solution to match with the outer region provided . With our chosen parameters in table 1, this is indeed the case, but it would not require a large change for it not to be so.

If, conversely, , then the fixed point in the phase plane , to which we must match, becomes a stable node or spiral. This change is associated with the lines and crossing each other so that, when , the trajectory coming from *N*=0 must cross the line *V*=0. This is only possible by passing through the point *G*=*V*=0, which looks like a degenerate saddle point in the phase plane; thus although there are many trajectories that spiral into the fixed point, there is still a unique one that passes through this degenerate saddle. This suggests that a unique boundary-layer solution should exist in which both *ϕ* and *N* show oscillatory decay towards the outer solution.

This somewhat intriguing behaviour can be related to the characteristics of the equations; including the time dependence in the boundary-layer analysis gives(4.1)There is one real characteristic with speed *V*/*ϵ*. When , the solution has *V*<0 corresponding to the characteristics entering the boundary layer from below, whereas with the characteristics must change direction within the boundary layer and point outwards to the outer solution.

Such solutions are found numerically when we increase the value of the parameter *P*; figure 7 shows the steady state for *P*=2. As *K*=0 and *V*=0 move progressively further apart on the phase plane, the spiralling becomes more pronounced, and there is the possibility that the effective pressure becomes negative; this is shown in figure 7 for *P*=10.

An exploration of parameter values does not show any obvious criterion that produces these oscillatory solutions, and there appears to be no clear physical meaning attached to the condition , other than a certain combination of physical properties; being a scaled inverse Péclet number, larger values of *P* imply greater relative importance of heat conduction compared with matrix advection. The very large oscillations in which *N* approaches zero seem to be unlikely, but are within the bounds of possibility, particularly if the ascent rate is small (less than 1 cm yr^{−1}) and the melt is not too viscous.

The oscillations have their origin in the unusual thermodynamics; the effective pressure controls the temperature, the matrix compaction *and* the melt movement. These last two dynamical roles for the effective pressure require that near the boundary of the partial melt there must be a very sudden change in pressure, but the fact that this produces large temperature gradients which heat conduction would attempt to smooth out seems to force the solutions for *N* (and therefore also *ϕ*) to have the unusual oscillatory profile.

### (b) Fracture initiation

We have assumed in this study that melt remains as a distributed porous flow and solidifies at the base of the lithosphere. Continued transport into the lithosphere requires localized flow through a conduit or fracture and will need a sufficiently large melt flux to avoid solidification. Magmafracturing (Spence *et al*. 1987; Lister & Kerr 1991; Roper & Lister 2005), aided by the supply of large quantities of melt from below, may be possible if substantial flow localization occurs in the partial melt, and several mechanical and chemical instabilities have been suggested to do this (Spiegelman & McKenzie 1987; Stevenson 1989; Aharonov *et al*. 1995; Spiegelman *et al*. 2001; Katz *et al*. 2006). Alternatively, the fracture could occur *within* the partially molten region and draw in melt from the surroundings to establish an interconnected network of veins or dykes (Nicolas 1986; Sleep 1988; Ito & Martel 2002). In either case, something must cause these fractures to be initiated.

The Griffith criterion suggests that such fracture should occur when the stored elastic energy becomes larger than the surface energy that is created when individual grains are pulled apart. The stored energy has size approximately and the surface energy is approximately , where *a* is the typical grain size, *γ* the surface energy, *E* Young's modulus and *σ* is the applied stress. Fracture therefore occurs if(4.2)with typical values of *Σ* of the order of 1 MPa (Sleep 1988; Fowler 1990*b*). If *τ*_{1} is the largest principal deviatoric stress, the largest effective stress is , giving the condition(4.3)

This suggests that the matrix will fracture if the effective pressure becomes small enough. At a spreading mid-ocean ridge, a rough estimate of deviatoric stress is using our previously given values, but this does not take account of the temperature-dependent viscosity; as the rock cools towards the surface, the viscosity will increase significantly from the 10^{19} Pa used above and stresses of several MPa or more are probable there. In this case *τ*_{1}>*Σ*, and since the upper boundary condition on the partial melt requires *N* to reduce towards zero there, the conditions for fracture initiation may well be met. The possibility for oscillatory solutions to the boundary layer in §4*a* also offer the intriguing possibility that the fracture criterion may be met lower within the partial melt; this would enable the newly initiated fracture to grow by drawing in melt from the surroundings.

Fractures would occur normal to the largest extensional deviatoric stress *τ*_{1}; beneath spreading mid-ocean ridges this is horizontal, so fractures will be vertical and allow for efficient transport of melt.

## 5. Conclusions

We have proposed model equations for an upwelling mantle region, comprising both partially molten and solid rock, with consistent conditions at the interfaces between these. The equations for the partial melt are essentially the same as those of McKenzie (1984) and Šrámek *et al*. (2007), the differences arising from their simplification and use. Our work follows the related study of soil mechanics in using the effective pressure as a primary variable on account of its fundamental control on the compaction of the matrix. The chief difference between our analysis and other authors' is then to make use of the fact that the melt fraction is small, and the bulk viscosity is larger than the intrinsic shear viscosity.

There has been little previous discussion in the literature of the appropriate boundary conditions for the partial melt region, although they must form an integral part of the problem; thermodynamic equilibrium of the partial melt boundaries requires that either the melt fraction *ϕ* or effective pressure *N* must be zero there. The partial melt dynamics can be reduced to an unusual degenerate pair of equations for *ϕ* and *N*, (3.10) and (3.11), for which these conditions appear to determine a unique steady state.

The solutions for one-dimensional upwelling are shown in figures 5 and 6 and agree in many respects with other published results: the parabolic melt profile (Ahern & Turcotte 1979; Šrámek *et al*. 2007); the near-linear decrease in matrix velocity; the viscous boundary layer near the onset of melting (Fowler 1990*b*; Šrámek *et al*. 2007); and the presence of a temperature precursor beneath the partial melt (Ahern & Turcotte 1979). Future work will hope to extend these solutions to two dimensions, with lateral movement of both matrix and melt, but solving the full problem in this case is by no means trivial.

The boundary-layer treatment of temperature and pressure gradients at the bottom of the partial melt allow an expression (3.35) to be derived for the depth of onset of melting. Šrámek *et al*. (2007) pointed out that including the density difference *r*≠1 will cause some extra problems near this boundary; there will be a small region in which the lower density melt requires an infinite pressure gradient to drive it through the matrix. The same behaviour occurs in our solutions if the density difference is included, and it seems that some form of disequilibrium melting may be required. With realistic parameters however the region over which these issues arise is very small and would, for instance, make no discernible difference to the graph of figure 4.

The boundary layer required at the top of the partial melt suggests an interesting interplay between the thermal and dynamical roles of the effective pressure, which can lead to an oscillatory profile (figure 7). A steady oscillatory structure occurs when the relative effects of heat conduction to matrix advection within this layer are large enough. Very significant oscillations in which the effective pressure approaches zero may be possible if the mantle ascent rate is less than approximately 1 cm yr^{−1}. It is not clear how such effects will manifest themselves in the two-dimensional case or when more thermodynamic components are considered. We cannot therefore conclude that such structures will necessarily exist, but rather that the role of heat conduction, when the solidus is pressure dependent, may have some unexpected effects.

The boundary condition requiring *N* to go to 0 at the top of the partial melt will necessarily produce conditions that make fracture according to (4.3) a distinct possibility, and oscillations of the pressure may also cause this initiation criterion to be met deeper in the partial melt. The requirement that melt pressures must adjust to maintain thermodynamic equilibrium may therefore be the cause of fracture and consequent continued transport into the lithosphere.

## Acknowledgments

A.C.F. acknowledges the support of the Mathematics Applications Consortium for Science and Industry (www.macsi.ul.ie) funded by the Science Foundation Ireland mathematics initiative grant 06/MI/005. I.J.H. acknowledges the award of an EPSRC studentship. We thank the reviewers S. Hier-Majumder and R. F. Katz for their helpful comments on the manuscript.

## Footnotes

- Received January 31, 2008.
- Accepted April 10, 2008.

- © 2008 The Royal Society