# Simultaneous melting and compaction in deformable two phase media

Ondrej Sramek, Yanick Ricard, David Bercovici

► To cite this version:

Ondrej Sramek, Yanick Ricard, David Bercovici. Simultaneous melting and compaction in deformable two phase media. Geophysical Journal International, 2007, 168 (3), pp.964-982. 10.1111/j.1365-246X.2006.03269.x. ensl-00139360

HAL Id: ensl-00139360

https://ens-lyon.hal.science/ensl-00139360

Submitted on 17 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Simultaneous melting and compaction in deformable two-phase media

Ondřej Šrámek,<fnref n="1" />,<fnref n="2" /> Yanick Ricard<fnref n="3" />,<fnref n="1" /> and David Bercovici<fnref n="1" />

Accepted 2006 October 22. Received 2006 October 16; in original form 2005 July 8

## SUMMARY

Melt generation and extraction are typically modelled using the two-phase equations developed by McKenzie or Scott and Stevenson. Various approximations are often made to simplify the problem which may lead to some unphysical results (e.g. thermodynamically inconsistent conditions of melting and unrealistic porosity profiles). We discuss a generalized version of the set of equations introduced by Bercovici et al. that allows for mass transfer between the two phases in a single component system and consider a self-consistent set of equations. In our description the two phases, solid and melt, are submitted to individual pressure fields whose difference is related to the surface tension at the interfaces, changes in porosity and the melting rate. A kinetic relation for the melting rate arises from the second law of thermodynamics. The condition of chemical equilibrium corresponds to the usual univariant equality of the chemical potentials of each phase when the matrix and melt are motionless. In the most general form, phase equilibrium is influenced by both the Gibbs–Thomson effect that arises naturally from thermodynamic considerations on surface tension and by the viscous deformation of the phases. We apply these new equations to a steady state problem of pressure release melting in a univariant system. We treat melting and compaction simultaneously and observe several new effects including multiple domains near the onset of melting that correspond to various force balances. A consequence of matrix compaction and melt expulsion (or matrix dilation and melt accumulation) is a pressure difference between melt and solid that facilitates (inhibits) melting. For parameters corresponding to mid-oceanic ridge magmatism, compaction permits melting to start as much as ~2 km below the standard solidus. Numerical solutions are necessary to determine the magnitude of the melt zone shift. Numerical results support the boundary layer solutions obtained analytically and suggest that in most of the melting zone the movement of melt and matrix should be close to the Darcy equilibrium where the buoyancy of melt is balanced by the viscous drag between the phases. The Darcy equilibrium follows an initial stage where the matrix viscous stresses balance Darcy drag. In all situations the steady state porosity profile remains a monotonic function of depth. The existence of a compaction layer following a melting zone where the porosity is maximum as described in various earlier publications has never been found.

**Key words:** melt generation, mid-ocean ridges, phase transitions, thermodynamics.

## 1 INTRODUCTION

Melting and melt migration are important processes within the Earth to transport heat and mass and to drive chemical differentiation. Below mid-ocean ridges upwelling material undergoes pressure release melting and new oceanic lithosphere is created (e.g. McKenzie & Bickle 1988). Partial melting processes are responsible for continental crust formation (e.g. White & McKenzie 1989). There are suggestions that dehydration-induced partial melting occurs in the transition zone (Bercovici & Karato 2003; Ohtani et al. 2004). Partial melting may also occur in the region above the core–mantle boundary (Williams & Garnero 1996). A rigorous theoretical framework is, therefore, necessary to account for melting in geodynamic modelling of various processes in the Earth.

Description of melting and subsequent melt migration requires a continuum model with at least two distinct phases. Since the pioneering work of Frank (1968) and Sleep (1974) the dynamics of two-phase flows has been extensively discussed in the geophysical literature (e.g. Ahern

<footnote n="1">¹Department of Geology and Geophysics, Yale University, PO Box 208109, New Haven, CT 06520-8109, USA. E-mail: ondrej.sramek@yale.edu</footnote>

<footnote n="2">²Laboratoire de Sciences de la Terre, CNRS UMR 5570, École normale supérieure de Lyon, 46 allée d'Italie, 69364 Lyon cedex 07, France.</footnote>

<footnote n="3">³Laboratoire de Sciences de la Terre, CNRS UMR 5570, Université de Lyon 1, 43 bd du 11 Novembre, bat. géode, 69622 Villeurbanne cedex, France.</footnote>
& Turcotte 1979; McKenzie 1984; Ribe 1985a,b; Scott & Stevenson 1986; Rabinowicz et al. 1987; Spiegelman 1993a,b). However, in previous papers (Bercovici et al. 2001a; Bercovici & Ricard 2003), we have shown that the self-consistent introduction of surface tension and surface free energy between two phases imposes additional rigor to the theoretical modelling. With surface tension, the two phases are clearly submitted to two different pressure fields. The partitioning of the surface energy and interface forces between the two phases also has to be defined. This led us to propose a new set of equations that is also applicable to non-equilibrium situations and provides a model of damage (Bercovici et al. 2001b; Bercovici & Ricard 2003; Ricard & Bercovici 2003). Close to thermodynamic equilibrium, neglecting surface tension, assuming that $\mu_m$, the viscosity of the solid matrix phase, is infinitely larger than $\mu_f$, the viscosity of the fluid phase, and assuming small porosity, our equations become equivalent to those proposed by McKenzie (1984), when the bulk viscosity, $\zeta$, is replaced by the porosity dependent factor proportional to $\mu_m/\phi$. Such a porosity-dependent bulk viscosity was also suggested by others (Nye 1953; Fowler 1985; Scott & Stevenson 1986; Schmeling 2000); it hinders the squeezing of the matrix and expulsion of melt at low porosity and favours it at large porosity where bulk and shear viscosity become comparable.

In Bercovici et al. (2001a) and Bercovici & Ricard (2003) however, we were not concerned with the possible transfer of matter between phases (i.e. with melting of the matrix or solidification of the fluid melt). While some previous studies (e.g. Ribe 1985a; Schmeling 2000) discussed the effect of a prescribed melting rate on matrix compaction, to our knowledge, no previous study considered the full feedback between viscous two-phase deformation and the thermodynamic conditions for melting/solidification. The usual simplifying assumptions (ascribing identical pressure to both phases or maintaining the matrix pressure to the lithostatic pressure, neglecting the interfacial effects, decoupling melting from viscous compaction, using empirical laws to relate degree of melting to temperature and pressure) can only be justified by comparison with results based on a more general description. In this study, we offer such a general framework. After formulating the theory for a single component system in Section 2, with the mass transfer between phases explicitly taken into account, we apply our new set of equations to a simple 1-D steady-state model of pressure-release melting in Section 3. Concluding remarks follow in Section 4.

## 2 GENERAL THEORY

We describe a two-phase continuum consisting of a solid matrix and a fluid (melt) phase. We assume that the individual phases have constant densities, $\rho_f$ in the fluid phase and $\rho_m$ in the solid, and constant Newtonian viscosities $\mu_f$ and $\mu_m$. Despite this notation that uses indices $f$ and $m$, however, in its most general form the theory is symmetrical with respect to swapping the phases. For any quantity $q$ with values $q_f$ and $q_m$ in the fluid and matrix phase, we define average and difference quantities by $\bar{q} = \phi q_f + (1-\phi)q_m$ and $\Delta q = q_m - q_f$.

Equations presented below result from averaging the ‘true’ or microscopic quantities of the mixture assumed isotropic over a control volume. This control volume is sufficiently small to be regarded as a ‘point’ in the continuum mechanical sense but also sufficiently large to contain many pores or grains (see Bear 1988; Bercovici et al. 2001a).

### 2.1 Mass conservation

The derivation of the two mass conservation equations is fairly standard (see McKenzie 1984; Bercovici et al. 2001a) and leads to

$$ \frac{\partial \phi}{\partial t} + \nabla \cdot [\phi \mathbf{v}_f] = \frac{\Delta \Gamma}{\rho_f}, \qquad (1) $$

$$ -\frac{\partial \phi}{\partial t} + \nabla \cdot [(1-\phi)\mathbf{v}_m] = -\frac{\Delta \Gamma}{\rho_m}, \qquad (2) $$

where $\phi$ is porosity or the volume fraction of the fluid phase, the fluid and matrix phases have volume averaged velocities $\mathbf{v}_f$ and $\mathbf{v}_m$ and $\Delta\Gamma$ is the rate of mass exchange per unit volume (kg m$^{-3}$ s$^{-1}$) between the phases. A positive $\Delta\Gamma$ corresponds to melting, a negative $\Delta\Gamma$ to solidification.

### 2.2 Momentum conservation

The momentum equations remain identical to those without phase change, at least in the approximation of infinite Prandtl number. Strictly speaking, new source terms of momentum could arise from the phase change. When matrix material melts, it undergoes a velocity jump from $\mathbf{v}_m$ to $\mathbf{v}_f$ which imposes an impulse discontinuity (Drew 1983). These terms are however of the order of the kinetic energy which is neglected in this study, as are the terms related to acceleration.

As discussed in Bercovici et al. (2001a), the total momentum conservation is

$$ -\nabla \bar{P} + \nabla \cdot \left( \bar{\mathbf{\Gamma}} - \bar{\rho} g \hat{\mathbf{z}} \right) + \nabla(\sigma\alpha) = 0. \qquad (3) $$

In this equation, $\bar{P}$, $\bar{\mathbf{\Gamma}}$ and $\bar{\rho}$ are the average pressure, viscous stress and density, $\hat{\mathbf{z}}$ is the unit vector in vertical direction, $g$ is the magnitude of gravitational acceleration, and $\sigma$ is the coefficient of surface tension. The quantity $\alpha$ is the melt-solid interfacial area per unit volume of the mixture; it has units of m$^{-1}$ and basically represents the inverse of the average size of grains and pores in the mixture. Here we do not consider grain-grain contacts; these are subject of a study by Hier-Majumder et al. (2006).
The total momentum equation is the sum of equations for fluid and matrix which, individually, are (Bercovici & Ricard 2003)

$$-\phi[\nabla P_f + \rho_f g \hat{z}] + \nabla \cdot [\phi \underline{\tau}_f] + c \Delta \mathbf{v} + \omega[\Delta P \nabla \phi + \nabla(\sigma \alpha)] = 0, \quad (4)$$

$$-(1 - \phi)[\nabla P_m + \rho_m g \hat{z}] + \nabla \cdot [(1 - \phi)\underline{\tau}_m] - c \Delta \mathbf{v} + (1 - \omega)[\Delta P \nabla \phi + \nabla(\sigma \alpha)] = 0, \quad (5)$$

where $\Delta \mathbf{v} = \mathbf{v}_m - \mathbf{v}_f$ is the difference in velocities, and likewise $\Delta P = P_m - P_f$ is the difference in pressures between the phases. The interaction coefficient $c$ is related to permeability which is itself a function of porosity. A symmetrical form compatible with the usual Darcy term is (see Bercovici et al. 2001a)

$$c = \frac{\mu_f \mu_m}{k_0 [\mu_f (1 - \phi)^{n-2} + \mu_m \phi^{n-2}]}, \quad (6)$$

where the permeability of the form $k_0\phi^n$ was used (the exponent $n$ is usually found to be about 2–3). In the following we assume $n = 2$ and also $\mu_m \gg \mu_f$. The interaction coefficient then becomes $c = \mu_f/k_0$ (independent of porosity). The parameter $\omega$ is discussed in Section 2.2.1

The relationship between stress and velocities does not include an explicit bulk viscosity term (Bercovici et al. 2001a), and for each phase $j$ the deviatoric stress is simply

$$\underline{\tau}_j = \mu_j \left( \nabla \mathbf{v}_j + [\nabla \mathbf{v}_j]^T - \frac{2}{3} \nabla \cdot \mathbf{v}_j \mathbf{I} \right), \quad (7)$$

where $j$ stands for $f$ or $m$, $[\cdot]^T$ denotes tensor transpose and $\mathbf{I}$ is the identity tensor.

Instead of using the two force balance eqs (4) and (5), it is often easier to use the total momentum eq. (3) and the ‘action–reaction’ equation that can be deduced by combining the two former equations,

$$-\phi(1-\phi)[\nabla \Delta P + \Delta \rho g \hat{z}] + \nabla \cdot [\phi(1-\phi)\Delta \underline{\tau}] - \bar{\underline{\tau}} \cdot \nabla \phi - c \Delta \mathbf{v} + (\phi - \omega)[\Delta P \nabla \phi + \nabla(\sigma \alpha)] = 0. \quad (8)$$

### 2.2.1 Interfacial quantities

The interface between the phases exerts a surface force on the two-phase mixture. The coefficient $\omega$, $0 < \omega < 1$, controls the partitioning of the interface surface force between the two phases (Bercovici & Ricard 2003) and represents the fraction of the volume-averaged surface force exerted on the fluid phase. The exact value of $\omega$ is related to the microscopic behaviour of the two phases (molecular bond strengths and thickness of the interfacial region or ‘selvedge’ layer) and represents the extent to which the interface is embedded in one phase more than the other. The only general physical constraints that we have is that $\omega$ must be zero when the fluid phase disappears (when $\phi = 0$) and when the fluid phase becomes unable to sustain stresses (when $\mu_f = 0$). A symmetrical form like

$$\omega = \frac{\phi \mu_f}{\phi \mu_f + (1 - \phi) \mu_m} \quad (9)$$

would satisfy these conditions.

Given the surface force partitioning, the velocity that controls the surface work is consequently found to be the $\omega$-weighted sum of the phase velocities

$$\mathbf{v}_{\omega} = \omega \mathbf{v}_{f} + (1 - \omega) \mathbf{v}_{m}, \quad (10)$$

which reflects the extent to which the interface is advected by each of the phases (Bercovici & Ricard 2003).

On the microscopic level, a phase change contributes to the interface motion with the velocity of the propagating phase front. However, this phase front velocity has opposite directions at opposite sides of a grain/pore. The velocity of the interfaces due to a phase change averaged over a control volume is proportional to the integral $\int \Delta\Gamma \hat{n} dS$, where the integration is over all the interfacial surfaces with normal $\hat{n}$ (pointing say, from the matrix phase to the fluid phase). Assuming that the two-phase mixture is isotropic then the interfaces are randomly oriented, in which case the integral $\int \Delta\Gamma \hat{n} dS$ is zero if $\Delta\Gamma$ is uniform or slowly varying relative to fluctuations in $\hat{n}$. The average interface velocity, therefore, remains $\mathbf{v}_{\omega}$ with or without melting, that is, is independent of the instantaneous melting rate.

Bercovici & Ricard (2003) show that the pressure $P_i$ on the interface between the two fluids is (see also Drew 1983)

$$P_i = (1 - \omega)P_f + \omega P_m. \quad (11)$$

Note that the weighting factors $\omega$ and $1 - \omega$ appear in a reverse order in (11) compared to (10). The behaviour of $P_i$ can be understood by rewriting the pressure jump $\Delta P$ as the difference between the pressure jump across the matrix selvedge layer and that across the fluid selvedge layer, $\Delta P = (P_m - P_i) - (P_f - P_i)$. In the case where the surface effects are totally embedded in the matrix ($\mu_f = 0$ and $\omega = 0$), the entire pressure drop across the interface should be accounted for by the matrix contribution, and therefore, controlled by the first term $P_m - P_i$ containing the matrix pressure. One thus gets $P_i = P_f$. Later (in Section 2.3) we will find it necessary to introduce an effective interface density weighted in the same way as the interface pressure $P_i$

$$\rho_i = (1 - \omega)\rho_f + \omega\rho_m. \quad (12)$$
## 2.3 Energy conservation

No terms related to phase change appear in the creeping momentum equations, and only appear in a trivial way in the two mass conservation equations. In contrast, the energy equation is modified in a substantial way and is now discussed in detail.

Each phase has a specific (per mass) internal energy, $\varepsilon_f$ and $\varepsilon_m$, respectively. In our model where individual phases are incompressible, the internal energy density of each phase is uniquely determined by the temperature. We assume the same temperature $T$ in both phases (local thermal equilibrium). We can write $d\varepsilon_j = C_jdT$, where $C_j$ is the heat capacity of phase $j$. To account for its thermodynamic presence, the interface is treated as an independent phase. The mechanical and thermodynamic state of the interface is described by its surface energy $\xi_i$, entropy $s_i$ (both per unit area), and surface tension coefficient $\sigma$; these quantities are (for chemically homogeneous systems) related according to (Bailyn 1994; Bercovici et al. 2001a)

$$s_i = -\frac{d\sigma}{dT}, \quad \xi_i = \sigma - T\frac{d\sigma}{dT}. \qquad (13)$$

The interfacial surface energy is transported by the moving continuum with the average interface velocity $v_\omega$ defined by eq. (10).

The conservation of energy is expressed by the following equation where the left-hand side represents the temporal change of energy content in a fixed control volume and the right-hand side represents the different contributions to this change, namely internal heat sources, $Q$, loss of energy due to diffusion, $\nabla \cdot q$, advection of energy, and rate of work of both surface and body forces.

$$\begin{aligned} \frac{\partial}{\partial t} [\phi \rho_f \varepsilon_f + (1-\phi) \rho_m \varepsilon_m + \xi_i \alpha] \\ = Q - \nabla \cdot q - \nabla \cdot [\phi \rho_f \varepsilon_f v_f + (1-\phi) \rho_m \varepsilon_m v_m + \xi_i \alpha v_\omega] \\ + \nabla \cdot [-\phi P_f v_f - (1-\phi) P_m v_m + \phi v_f \cdot \tau_f + (1-\phi) v_m \cdot \tau_m + \sigma \alpha v_\omega] \\ - \phi v_f \cdot \rho_f g \hat{z} - (1-\phi) v_m \cdot \rho_m g \hat{z}. \end{aligned} \qquad (14)$$

The last equation is manipulated in the standard way using the mass and momentum equations. We also use expressions for $d\varepsilon_j$ and relationship between interfacial energy and entropy, (13). After some algebra we get

$$\begin{aligned} \phi \rho_f C_f \frac{D_f T}{Dt} + (1-\phi) \rho_m C_m \frac{D_m T}{Dt} - T \frac{D_\omega}{Dt} \left(\alpha \frac{d\sigma}{dT}\right) - T \alpha \frac{d\sigma}{dT} \nabla \cdot v_\omega \\ = Q - \nabla \cdot q + \Psi - \left(\Delta P + \sigma \frac{d\alpha}{d\phi}\right) \frac{D_\omega \phi}{Dt} + \left[\Delta \varepsilon + \frac{P_m}{\rho_m} - \frac{P_f}{\rho_f}\right] \Delta \Gamma, \end{aligned} \qquad (15)$$

where $\Delta\varepsilon = \varepsilon_m - \varepsilon_f$ is the difference between the specific internal energies of the phases, and $\Psi$ is the rate of deformational work,

$$\Psi = c \Delta v^2 + \phi \nabla v_f : \tau_f + (1-\phi) \nabla v_m : \tau_m. \qquad (16)$$

The fundamental derivatives are defined by

$$\frac{D_j}{Dt} = \frac{\partial}{\partial t} + v_j \cdot \nabla, \qquad (17)$$

where $j$ is to be substituted with the appropriate subscript $f$, $m$ or $\omega$.

Eq. (15) is arranged in such a way as to group the terms related to temporal entropy variations on the left-hand side ($-d\sigma/dT$ is the interfacial entropy). The right-hand side then contains terms related to entropy sources and fluxes. In addition to the usual heat production, diffusion and deformational work, two other expressions appear with clear physical meanings. One of them includes $\Delta P + \sigma d\alpha/d\phi$. As the quantity $d\alpha/d\phi$ is the sum of the average curvatures of grain/pores, $2/(R)$ (here curvature is defined positive when concave to the fluid; see Bercovici et al. 2001a), the Laplace's condition on the pressure drop across a curved interface is recovered when the pressure difference between the two phases is $\Delta P = -\sigma d\phi/d\alpha$ (Landau & Lifshitz 1959). The Laplace's equilibrium condition can only be satisfied if the phases are stationary. Thus $\Delta P + \sigma d\phi/d\alpha$ is the out-of-equilibrium pressure difference, that is due to the viscous deformation of the phases. The other term contains the difference in the specific enthalpies $\Delta h = h_m - h_f$, where the enthalpy of phase $j$ is defined by $h_j = \varepsilon_j + P_j/\rho_j$. The energy eq. (15) corresponds to the one found in Bercovici et al. (2001a, eq. 59). It contains one additional source term, $\Delta h \Delta \Gamma$, which arises from the phase change.

The mass conservation eqs (1) and (2), the momentum eqs (4) and (5), augmented by the viscous stress relations (7), and the energy eq. (15) need to be supplemented with two more relations. In particular, we need to determine the pressure difference $\Delta P$ between the phases and the melting rate $\Delta \Gamma$. To this end, we examine entropy production in the system and employ the method of non-equilibrium thermodynamics (e.g. de Groot & Mazur 1984) to arrive at the two requisite equations, similar to the approach that was presented in Bercovici et al. (2001a). In the case with no phase change ($\Delta \Gamma$ identically zero) this non-equilibrium thermodynamics approach led to a relationship between the interphase pressure drop $\Delta P$ and viscous compaction (Bercovici et al. 2001a).

The entropy conservation writes in the most general case (de Groot & Mazur 1984)

$$\frac{\partial}{\partial t} \left[ \phi \rho_f s_f + (1-\phi) \rho_m s_m - \frac{d\sigma}{dT} \right] = -\nabla \cdot \left[ \phi \rho_f s_f v_f + (1-\phi) \rho_m s_m v_m - \frac{d\sigma}{dT} \alpha v_\omega \right] - \nabla \cdot \mathbf{J} + S, \qquad (18)$$

where $\mathbf{J}$ represents the microscopic non-convective entropy flux and $S$ is the internal entropy production. According to the second law of thermodynamics, $S$ cannot be negative. To identify the entropy sources and the microscopic flux in our particular case, we compare the energy and the entropy eqs (15) and (18) taking into account that, for each incompressible phase, $ds_j = C_jdT/T = d\varepsilon_j/T$. After some algebra, one
gets

$$ \mathbf{J} = \frac{\mathbf{q}}{T}, \qquad (19) $$

$$ TS = Q - \frac{1}{T} \mathbf{q} \cdot \nabla T + \Psi - \left( \Delta P + \sigma \frac{d\alpha}{d\phi} \right) \frac{D_{\omega}\phi}{Dt} + \Delta\Gamma\Delta\mu, \qquad (20) $$

where we have introduced the difference in chemical potential between the two phases

$$ \Delta\mu = \Delta\epsilon + \frac{P_m}{\rho_m} - \frac{P_f}{\rho_f} - T\Delta S, \qquad (21) $$

where $\Delta S = S_m - S_f$ is the change in specific entropies. In addition to the usual entropy sources related to heat production, diffusion, and dissipation $\Psi$ (see eq. 16), we find two additional entropy sources: one is related to the out-of-equilibrium pressure difference and goes as $-(\Delta P + \sigma d\alpha/d\phi)$ multiplied by the temporal change in porosity $D_{\omega}\phi/Dt$; the other arises from the phase change rate $\Delta\Gamma$, which is multiplied by the difference of chemical potential of the bulk phases $\Delta\mu$.

It is tempting to recognize two thermodynamic scalar forces and their conjugate fluxes in the last two terms of the entropy production (20). However, with this choice, the individual forces and fluxes are not independent; the pressure jump affects the chemical potential difference, and the porosity change is affected by the melting rate. Therefore, we reorganize the entropy production term $S$ in a somewhat different form. We write

$$ TS = Q - \frac{1}{T} \mathbf{q} \cdot \nabla T + \Psi - \left( \Delta P + \sigma \frac{d\alpha}{d\phi} \right) \left( \frac{D_{\omega}\phi}{Dt} - \frac{\rho_i}{\rho_f \rho_m} \Delta\Gamma \right) + \Delta\Gamma \left[ \Delta\mu - \frac{\rho_i}{\rho_f \rho_m} \left( \Delta P + \sigma \frac{d\alpha}{d\phi} \right) \right]. \qquad (22) $$

In this expression, the pressure jump in excess of Laplace's condition appears as conjugate to the porosity change in excess of melting, which can be expressed, using (1), (2), (10) and (12), as

$$ \frac{D_{\omega}\phi}{Dt} - \frac{\rho_i}{\rho_f \rho_m} \Delta\Gamma = (1-\omega)(1-\phi)\nabla \cdot \mathbf{v}_m - \omega\phi\nabla \cdot \mathbf{v}_f. \qquad (23) $$

Using (11), (12) and (21), we can write

$$ \Delta\mu - \frac{\rho_i}{\rho_f \rho_m} \left( \Delta P + \sigma \frac{d\alpha}{d\phi} \right) = \Delta\epsilon - T\Delta S - P_i \frac{\Delta\rho}{\rho_f \rho_m} - \sigma \frac{d\alpha}{d\phi} \frac{\rho_i}{\rho_f \rho_m}. \qquad (24) $$

The expressions (23) and (24) are chosen as the thermodynamic forces. Following the standard procedure of non-equilibrium thermodynamics (e.g. de Groot & Mazur 1984), we assume that there is a coupled linear relationship between the two thermodynamic scalar fluxes, and the two thermodynamic forces, to write

$$ - \left( \Delta P + \sigma \frac{d\alpha}{d\phi} \right) = L_{11} \left( \frac{D_{\omega}\phi}{Dt} - \frac{\rho_i}{\rho_f \rho_m} \Delta\Gamma \right) + L_{12} \left( \Delta\epsilon - T\Delta S - P_i \frac{\Delta\rho}{\rho_f \rho_m} - \sigma \frac{d\alpha}{d\phi} \frac{\rho_i}{\rho_f \rho_m} \right), \qquad (25) $$

$$ \Delta\Gamma = L_{21} \left( \frac{D_{\omega}\phi}{Dt} - \frac{\rho_i}{\rho_f \rho_m} \Delta\Gamma \right) + L_{22} \left( \Delta\epsilon - T\Delta S - P_i \frac{\Delta\rho}{\rho_f \rho_m} - \sigma \frac{d\alpha}{d\phi} \frac{\rho_i}{\rho_f \rho_m} \right). \qquad (26) $$

The phenomenological coefficients $L_{ij}$ form a 2 by 2 matrix $L$. The off-diagonal coefficients of this matrix can be constrained using the Onsager's theorem (de Groot & Mazur 1984). The two forces (23) and (24) in our problem are odd and even functions of velocities, respectively; then following Onsager's findings we have $L_{21} = -L_{12}$, that is, the matrix $L$ is antisymmetric. To insure the positivity of the entropy production, the symmetric part of the matrix $L$ must be positive definite; this constraint gives us $L_{11} > 0$ and $L_{22} > 0$.

The existence of a non-zero coupling phenomenological coefficients like $L_{12}$ is possible in irreversible thermodynamic but is not always the case. In our context, a thought experiment suggests that the two fluxes decouple ($L_{12} = L_{21} = 0$). Let us consider a situation of homogeneous isotropic melting without surface tension where the melt has such a low viscosity that it cannot sustain viscous stresses and cannot interact with the solid by Darcy terms. For such an inviscid melt $\omega = 0$ according to (9). In this case, as the melt can escape instantaneously, the matrix should not dilate (i.e. using 23, the first term on the right-hand side of 25 is zero), and thus the two pressures should be the same (therefore, the left-hand side should disappear). The only solution for the eq. (25) to hold is to have $L_{12} = 0$ as the term containing $\Delta\mu$ can be arbitrarily imposed.

Another requirement is that (25) should correspond to what was found in Bercovici et al. (2001a) and Bercovici & Ricard (2003) when $\Delta\Gamma = 0$. Simple micromechanical models (e.g. Nye 1953) allow us to evaluate $L_{11}$,

$$ L_{11} = K_0 \frac{\mu_f + \mu_m}{\phi(1 - \phi)}. \qquad (27) $$

The dimensionless constant $K_0$ accounts for grain/pore geometry and is of $O(1)$.

Using these constraints we can, therefore, restate the two phenomenological equations as

$$ \Delta P + \sigma \frac{d\alpha}{d\phi} = -K_0 \frac{\mu_f + \mu_m}{\phi(1-\phi)} \left( \frac{D_{\omega}\phi}{Dt} - \frac{\rho_i}{\rho_f \rho_m} \Delta\Gamma \right), \qquad (28) $$

$$ \Delta\Gamma = L_{22} \left( \Delta\epsilon - T\Delta S - P_i \frac{\Delta\rho}{\rho_f \rho_m} - \sigma \frac{d\alpha}{d\phi} \frac{\rho_i}{\rho_f \rho_m} \right). \qquad (29) $$
The first equation establishes a general relationship controlling the pressure drop between phases. The left-hand side represents the out-of-equilibrium pressure difference due to deformation of the two-phase mush. This deformation is represented on the right-hand side by the temporal change in porosity less the contribution from the phase transition. This condition is in agreement with that found in previous papers without phases changes (Bercovici et al. 2001a; Ricard et al. 2001; Bercovici & Ricard 2003). In the former, it was argued that some deformational work could also affect this pressure drop. This hypothesis led to a damage theory developed in Bercovici et al. (2001b), Ricard & Bercovici (2003) and Bercovici & Ricard (2005). Here we assume that the system remains close enough to mechanical equilibrium that no damage occurs and we use (28). The generalization would anyway be straightforward.

The second equation gives the melting rate in terms of departure from thermodynamic equilibrium between the phases, represented by the parentheses on the right-hand side. The linear relationship assumes that the departures from equilibrium are small. The positive coefficient $L_{22}$ is related to the usual rate constant $k_r$ of kinetic theory (see e.g. Lasaga 1998).

## 2.4 The equilibrium condition

In the case of mechanical equilibrium (i.e. when Laplace's condition holds), (29) recovers the usual proportionality of the melting rate to the difference in chemical potentials of single phases $\Delta\mu$ (see also 24). In general the non-equilibrium pressure drop affects the melting rate. We thus define a new chemical equilibrium state from (29) with $\Delta\Gamma/L_{22} \to 0$, which requires

$$ \Delta\epsilon - T\Delta S - P_i \frac{\Delta\rho}{\rho_f \rho_m} - \sigma \frac{d\alpha}{d\phi} \frac{\rho_i}{\rho_f \rho_m} = 0. \quad (30) $$

This equilibrium is controlled by the interface properties, interface pressure $P_i$ and interface density $\rho_i$ (see 11 and 12).

In the following we assume that the melting occurs at equilibrium and use (30) rather than the kinetic relation (29) (this is equivalent to using (29) with a very large kinetic coefficient $L_{22}$). We also make the assumption that the heat capacities are constant; in this case, for each phase j one has

$$ \epsilon_j = \epsilon_j^0 + C_j(T - T_0), \quad (31) $$

$$ s_j = s_j^0 + C_j \ln \left( \frac{T}{T_0} \right), \quad (32) $$

where $T_0$, $\epsilon_j^0$ and $s_j^0$ are reference values. We assume that the two heat capacities are equal, $C_f = C_m = C$, and the equilibrium condition (30) thus simply writes

$$ (T - T_0)\Delta S + P_i \frac{\Delta\rho}{\rho_f \rho_m} + \sigma \frac{d\alpha}{d\phi} \frac{\rho_i}{\rho_f \rho_m} = 0, \quad (33) $$

where $\Delta S = \Delta s^0$ and $T_0 = \Delta\epsilon^0/\Delta s^0$; $T_0$ is, therefore, the melting temperature of a flat interface at zero pressure. In this equation one can introduce the classical Clapeyron slope $\gamma' = -\Delta\rho/(\rho_f\rho_m\Delta S)$ (in K Pa$^{-1}$). We can also express the melting temperature in terms of the measurable average pressure $\bar{P}$:

$$ T = T_0 + \gamma' \bar{P} - \frac{\sigma}{\Delta S} \frac{\bar{P}}{\rho_f \rho_m} \frac{d\alpha}{d\phi} + \gamma' K_0 (1 - \omega - \phi) \frac{\mu_m + \mu_f}{\phi(1 - \phi)} \left( \frac{D_\omega \phi}{Dt} - \frac{\rho_i}{\rho_f \rho_m} \Delta\Gamma \right). \quad (34) $$

This equation shows that two new terms affect the Clapeyron static melting condition, $T = T_0 + \gamma'\bar{P}$, that would hold in a single phase modelling. The third term on the right side of (33), containing the surface tension $\sigma$, represents the Gibbs–Thomson effect. For grains of silicates in contact with the melt, the values of surface tension are around 1 J m$^{-2}$. With $\rho \sim 3000$ kg m$^{-3}$ and $\Delta S \sim -(300-400)$ J K$^{-1}$ kg$^{-1}$, this term is quite small for grain size of 10 $\mu$m–1 mm, spanning a range of about 0.001–0.1 K. In this calculation we took 1/grain size as the measure of characteristic interface curvature. In theory the Gibbs–Thomson effect could be significant, depending on the behaviour of $d\alpha/d\phi$, which is a more appropriate measure of the average curvature, that is, depending on the microscopic model of the interfaces. In the case where all interfaces are between matrix and magma (i.e. when grain–grain interfaces are not considered), the dependence of $\alpha$ on $\phi$ is such that $d\alpha/d\phi$ becomes very large as $\phi \to 0$, remains positive at small $\phi$ and changes sign at some $\phi_c$. The nucleation of a new phase (at $\phi = 0$) requires lower pressure or larger temperature than those of the Clapeyron curve to overcome the pressure jump associated with interface surface tension. Similar behaviour was found in the model that includes grain–grain interfaces in the case of large (> 60°) dihedral angle. If the dihedral angle is small (< 60°), the presence of grain–grain interfaces causes the sum of average curvatures $d\alpha/d\phi > 0$ to become negative and nucleation of melt occurs at higher pressure or lower temperature compared to Clapeyron conditions (Hier-Majumder et al. 2006). The question remains as to the nature of the first melt to appear. If the first liquid forms at some localized nucleation sites as tiny droplets, the pressure will be increased due to large positive interface curvature, and the Gibbs–Thomson effect will act against melting. This case could be, at least qualitatively, modelled using a simple interface model without grain–grain interfaces (discussed above); the model could be adjusted such that the change in curvature sign occurs at some very small porosity $\phi_c$. On the other hand, if the first melt behaves as a wetting grain boundary phase, then melting is probably facilitated by the Gibbs–Thomson effect.

The last term on the right side of (34) is the dynamic pressure perturbation caused by matrix dilation or compaction (that is not associated with melting of freezing). With $\mu_f \ll \mu_m$ and thus $\omega = 0$, the dynamic pressure perturbation simplifies to $K_0\mu_m(1 - \phi)/\phi \nabla \cdot \mathbf{v}_m$. In a
compacting melting matrix ($\nabla \cdot \mathbf{v}_m < 0$), the dynamic pressure is negative and melting thus occurs at larger mean pressure $\bar{P}$ or at lower temperature than in a static case. Reciprocally, in a dilating matrix ($\nabla \cdot \mathbf{v}_m > 0$) melting requires lower pressure or higher temperatures. When $\rho_f < \rho_m$, the bulk volume increases upon melting. The melt is extracted since it is lighter than the matrix and rises faster. Whether the matrix compacts or dilates depends on how efficient the melt extraction is relative to the bulk volume increase. We have thus identified two effects that change the Clapeyron melting condition based on the average pressure $\bar{P}$: Gibbs-Thomson effect and dynamic pressure perturbation. Both can influence the phase equilibrium in either direction. The dynamic pressure effect will be quantified in Section 3.3.

## 2.5 Simplifying assumptions and governing equations

The introduction of surface tension was crucial to derive our equations properly by forcing us to be explicit in the definitions of the various quantities, pressures, densities and velocities, pertaining to surfaces or phases. However, the surface tension also induces the well-known Gibbs-Thomson effect in which case melting may be totally forbidden by capillary forces unless a more sophisticated model of heterogeneous nucleation is considered. This physical difficulty could be removed by considering a more realistic case wherein the matrix itself contains both grain-fluid and grain-grain interfaces, each with different surface tension (see Hier-Majumder et al. 2006, for a complete discussion). We leave the discussion of the surface tension effects for future studies and assume in the rest of the paper that $\sigma = 0$.

We further assume that the matrix is much more viscous than the fluid phase ($\mu_f \ll \mu_m$) as typical for melting scenarios, which implies that $\tau_f = 0$, $\omega = 0$, $\rho_i = \rho_f$, $P_i = P_f$, and $\mathbf{v}_\omega = \mathbf{v}_m$. For convenience we summarize the governing equations with these assumptions ($\sigma = 0$, $\mu_f \ll \mu_m$).

By combining eqs (1) and (2) we obtain the mass conservation for the two phase material

$$ \nabla \cdot \bar{\rho} \mathbf{v} = 0, \qquad (35) $$

where $\bar{\rho}\mathbf{v} = \rho_f\phi\mathbf{v}_f + \rho_m(1-\phi)\mathbf{v}_m$. It must be supplemented by another mass conservation equation, for example, the matrix equation

$$ -\frac{\partial \phi}{\partial t} + \nabla \cdot [(1-\phi)\mathbf{v}_m] = -\frac{\Delta \Gamma}{\rho_m}. \qquad (36) $$

The equations of conservation of momentum for the fluid phase is

$$ -\phi \nabla [P_f + \rho_f g z] + c \Delta \mathbf{v} = 0, \qquad (37) $$

and the force-difference (or action-reaction) equation is

$$ -\nabla[(1-\phi)\Delta P] - (1-\phi)\Delta\rho g \hat{\mathbf{z}} + \nabla \cdot [(1-\phi)\tau_m] - \frac{c\Delta\mathbf{v}}{\phi} = 0. \qquad (38) $$

The deviatoric stress in the matrix is given by

$$ \tau_m = \mu_m \left( \nabla \mathbf{v}_m + [\nabla \mathbf{v}_m]^T - \frac{2}{3} \nabla \cdot \mathbf{v}_m \mathbf{I} \right). \qquad (39) $$

The pressure jump between phases becomes

$$ \Delta P = -K_0 \frac{\mu_m}{\phi} \nabla \cdot \mathbf{v}_m. \qquad (40) $$

Here $K_0\mu_m/\phi$ represents the bulk viscosity. The value of the dimensionless factor $K_0$ depends on the geometry of the two-phase mixture. In the case of non-interacting cylindrical inclusions Bercovici et al. (2001a) found $K_0 = 1$. Some more complicated functional dependence of bulk viscosity on porosity was suggested (see discussion in Schmeling 2000) that takes into account the interaction between individual pores/grains. Here we will assume the simple $1/\phi$ dependence as it captures the most important characteristic of the bulk viscosity in the limit of small porosity. For the sake of compact mathematical expressions here we adopt the value $K_0 = 4/3$.

In general, when the phases are not in thermodynamic equilibrium, the rate of melting is determined by (29), which, with (31) and (32), yields

$$ \Delta\Gamma = -L_{22} \left[ (T - T_0)\Delta S + P_f \frac{\Delta\rho}{\rho_f \rho_m} \right], \qquad (41) $$

where $P_f$ could be replaced by $\bar{P} + K_0\mu_m(1-\phi)\nabla \cdot \mathbf{v}_m/\phi$ to emphasize the effect of compaction on the melting condition. In the following we assume thermodynamic equilibrium, that is, an infinitely fast reaction rate ($L_{22} \to \infty$).

The energy equation is

$$ \rho_f \phi C \frac{D_f T}{D_t} + \rho_m (1-\phi) C \frac{D_m T}{D_t} - T \Delta S \Delta \Gamma = Q - \nabla \cdot \mathbf{q} + \frac{\Delta \Gamma^2}{L_{22}} + K_0 \mu_m \frac{1-\phi}{\phi} (\nabla \cdot \mathbf{v}_m)^2 + \Psi. \qquad (42) $$

On the right-hand side one recognizes the various dissipative terms, in particular, the chemical relaxation term $\Delta\Gamma^2/L_{22}$ (de Groot & Mazur 1984), and the bulk compression term $K_0\mu_m(1-\phi)/\phi(\nabla \cdot \mathbf{v}_m)^2$.

In the case of thermodynamic equilibrium ($L_{22} \to \infty$) the Gibbs equilibrium condition (41) is no longer an equation for the melting rate. It can be interpreted as an equation for the temperature while the energy eq. (42) controls the melting rate.
## 3 1-D STEADY-STATE MELTING

Various studies have dealt with 2-D simulations of melt migration with a simplified set of equations (Ribe 1988; Spiegelman 1993c; Choblet & Parmentier 2001). However, to understand the physics of two-phase medium near the onset of melting we restrict our study to a simple 1-D melting problem. Similar problems were discussed by McKenzie (1984), Ribe (1985a) and Turcotte & Phipps Morgan (1992). Although some of these studies described simultaneous melting and compaction, they did not account for the feedback between the viscous deformation and the thermodynamics of melting. McKenzie (1984) considered isentropic batch melting (i.e. no melt extraction) and prescribed a dependence of the degree of melting on temperature and pressure. Turcotte & Phipps Morgan (1992) and Ribe (1985a) used an a priori melting relation; Turcotte & Phipps Morgan (1992) actually decoupled the processes of melting and compaction. Only Fowler (1989) discussed the effect of compaction on the melting behaviour in the case of simultaneous melting and compaction. However, Fowler's basic equations differ from ours in that he assumed the matrix pressure is lithostatic. The Darcy interaction term present in his fluid momentum equation that goes as $v_f - \bar{v}$ does not have its equal and opposite counterpart in the matrix momentum equation. Therefore, Fowler's formulation violates Newton's 3rd law; in particular, when the two momentum equations are added to describe the force balance on the entire mixture, the Darcy drag term remains as a spurious body force (i.e. the drag of the fluid on the matrix must be equal and opposite to the drag of the matrix on the fluid).

In our model we consider the melting of a simple univariant system. In the melting region temperature follows the Clapeyron slope, which describes equilibrium between the phases. Melt extraction occurs simultaneously with melting, which leads to viscous deformation of the matrix and consequently to different pressure fields in the two phases.

We assume matrix material initially ascends with an upward velocity $V$ before melting starts. No heat is supplied to the ascending matrix. As the matrix is incompressible, before melting starts, there is no adiabatic gradient and the far-field temperature $T_0$ remains uniform. The melting starts when the ascending solid reaches the melting pressure. This point is considered as the origin of the coordinate $z$, which is positive upward. From this point on two phases coexist until all the matrix has melted. Above the melting zone, melt is transported isothermally again. As in the previous section we assume that the matrix is much more viscous than the fluid and we neglect surface tension. We only study the steady-state solution of this problem.

First, we write the equations in the dimensionless form. The natural length scale for this melting problem is the width of the melting zone. It can be estimated from a simplified form of the energy eq. (42) assuming that the drop in temperature is purely due to heat consumed by melting, that is, $\Delta\Gamma = \rho_m V C(dT/dz)/(T \Delta S)$, where the temperature gradient $dT/dz$ is the Clapeyron slope. For the purpose of scaling, $T$ can be replaced with $T_0$, the initial temperature, and the melting rate $\Delta\Gamma$ can be assumed constant; in this case, integration of (36) in 1-D steady state over the depth of the partial-melt layer $H$ leads to $\Delta\Gamma = \rho_m V/H$. The previous two expressions for $\Delta\Gamma$ combine to give $H = \rho_f T_0 \Delta S^2/(\Delta\rho g C)$. Velocity and temperature are scaled with their initial values $V$ and $T_0$. Pressure is scaled with $\rho_m g H$. Performing this scaling, the dimensionless thermal conductivity $k_T = \rho_f k'_T T_0/(\rho_m \Delta\rho g V H^2)$ (where $k'_T$ is dimensional conductivity), buoyancy velocity $V_B = \Delta\rho g/(cV)$, dimensionless compaction length $\delta = \sqrt{4\mu_m/(3cH^2)}$ (McKenzie 1984), and dimensionless Clapeyron slope $\gamma = \gamma' \rho_m g H/T_0$ appear. The dimensionless Clapeyron slope is also simply the ratio of entropy difference to heat capacity, $\gamma = -\Delta S/C$. At last we introduce $R$ as the density ratio $R = \rho_f/\rho_m$. We assume $\rho_f < \rho_m$ and $\Delta S < 0$, therefore, all the dimensionless parameters are positive. Note that in our description, the Darcy interaction coefficient $c$ is constant. More generally it is written as $c = \mu_f \phi^2/k(\phi)$ in the case of $\mu_m \gg \mu_f$ (see eq. 6) with $k(\phi)$ a porosity-dependent permeability. We assume $k(\phi) = k_0 \phi^2$ which leads to a constant $c = \mu_f/k_0$.

With these definitions, the 1-D steady-state governing equations are the dimensionless form of the action-reaction eq. (38) (with eqs 39 and 40),

$$ \delta^2 \phi \frac{d}{dz} \left( \frac{1-\phi^2}{\phi} \frac{dv_m}{dz} \right) - \Delta v = V_B \phi (1-\phi), \quad (43) $$

where $\Delta v = v_m - v_f$ can be expressed in terms of porosity and matrix velocity from integrated mass conservation eq. (35) (continuous mass flux implies $\overline{\rho v} = \text{const.} = \rho_m V$),

$$ \Delta v = \frac{[1 - (1 - R)\phi]v_m - 1}{R\phi}, \quad (44) $$

the equation for the fluid pressure (37),

$$ \frac{dP_f}{dz} = -R + \frac{1-R}{V_B} \frac{\Delta v}{\phi}, \quad (45) $$

the equilibrium condition (41) assuming $L_{22} \to \infty$,

$$ T = 1 + \gamma P_f, \quad (46) $$

the energy eq. (42) with equilibrium melting ($L_{22} \to \infty$),

$$ \frac{1}{\gamma^2} \frac{dT}{dz} + \frac{1}{\gamma} T \Delta\Gamma = k_T \frac{d^2 T}{dz^2} + \frac{R \Delta v^2}{V_B} + \frac{R \delta^2}{V_B} \frac{1-\phi^2}{\phi} \left( \frac{dv_m}{dz} \right)^2, \quad (47) $$

and the matrix mass conservation (36),

$$ \frac{d}{dz}[(1-\phi)v_m] + \Delta\Gamma = 0. \quad (48) $$
This set of equations must be supplemented by their boundary conditions, which, as usual in two phase flows, are not trivial. The porosity $\phi$ goes continuously from 0, for $z \le 0$, to 1, for $z \ge z_f$, where $z_f$ is the as yet unknown position of where melting is complete ($\phi = 1$). The matrix and melt fluxes, $(1 - \phi)v_m$ and $R\phi v_f$ (or $\rho_m(1 - \phi)v_m$ and $\rho_f \phi v_f$ in dimensional form), are also continuous; this implies that $v_m$ is continuous and equal to 1 (the upward velocity) at $z = 0$, while $v_f$ is continuous and equal to $R^{-1}$ (the upward velocity augmented by the matrix to melt density ratio) at $z = z_f$. However, $v_f$ is not necessarily continuous at $z = 0$ nor is $v_m$ at $z = z_f$. By definition $T = 1$ and $P_f = 0$ at $z = 0$. Eliminating $\Delta v$ between (43) and (45) we obtain

$$ \frac{dP_f}{dz} = -1 + (1 - R)\phi + \frac{(1 - R)\delta^2}{V_B} \left. \frac{d}{dz} \left( \frac{1 - \phi^2}{\phi} \frac{dv_m}{dz} \right) \right|_{z=0} \quad (49) $$

To assure $P_f = 0$ at $z = 0$, we introduce an integration constant

$$ z_0 = \frac{(1 - R)\delta^2}{V_B} \left. \frac{1}{\phi} \frac{dv_m}{dz} \right|_{z=0} \quad (50) $$

which has to be finite for the solution to be non-singular at $z = 0$. We thus have $dv_m/dz = 0$ at $z = 0$ since $\phi = 0$ there. Comparing (50) to the 1-D dimensionless form of (40), one finds that $z_0 = -\Delta P|_{z=0} = -P_m|_{z=0} = -\bar{P}|_{z=0}$, the second equality coming from $P_f|_{z=0} = 0$ and the third using also $\phi = 0$ at $z = 0$. The finite pressure difference at the onset of melting results from the viscous deformation of the matrix; it is the dynamic pressure perturbation discussed in Section 2.4. The consequence is that the melting zone is shifted by $z_0$ (upwards for positive $z_0$) relative to what suggests the simple Clapeyron curve based on the average pressure. Finally, eq. (43) at $\phi = 1$ implies that $2\delta^2(d\phi/dz)(dv_m/dz) + \Delta v = 0$ at $z = z_f$.

In order to constrain the dimensionless parameters of the governing equations, we estimate values of various material parameters for typical conditions of oceanic spreading centres. The viscosity of the solid $\mu_m$ could be about $10^{18}–10^{19}$ Pa s. For the viscosity $\mu_f$ of the basaltic melt we adopt 10 Pa s. These values of viscosities are however not well constrained and they span a wide range in the literature. Moreover, a shift of several orders of magnitude in viscosity might be expected in the presence of water (e.g. Karato 1989; Mei & Kohlstedt 2000; Karato & Jung 2003). Here we will assume dry conditions. We choose $k_0 = 5 \times 10^{-10}$ m$^2$, which corresponds to grain size of the order of 1 mm. The value of the compaction length $\delta$ then spans a range 8–26 km. A characteristic upwelling velocity of the mantle below the melting region is few cm yr$^{-1}$, up to about 10 cm yr$^{-1}$. We take $\rho_m = 3200$ kg m$^{-3}$ and $\rho_f = 2800$ kg m$^{-3}$ for the densities of the solid and melt. The buoyancy velocity scale $V_B$, therefore, spans a range of 60–160, the lower value corresponding to high upwelling velocity. We note that these values of $\delta$ and $V_B$ were obtained for relatively low viscosity basaltic melt. These numbers would be smaller for more viscous siliceous melts. The usual values of heat capacity and entropy of fusion found in literature fall into $C = 1000–1300$ J K$^{-1}$ kg$^{-1}$, and $\Delta S = -(250–400)$ J kg$^{-1}$ K$^{-1}$ (see Kojitani & Akaogi 1995, and references therein). The Clapeyron slope of melting can be calculated from the given parameters. Alternatively, one can look at the solidus values inferred from experiments on peridotite samples. Compilation by Hirschmann (2000) of experimental data gives a depth-dependent solidus, where the slope changes from 133 K GPa$^{-1}$ at ambient pressure to 100 K GPa$^{-1}$ at 6.5 GPa (200 km depth). The parameters are listed in Table 1.

For calculations in this study we use the following preferred values: $\mu_m = 10^{18}$ Pa s, $\mu_f = 10$ Pa s, $V = 10$ cm yr$^{-1}$, $k_0 = 5 \times 10^{-10}$ m$^2$, $\Delta S = -340$ J kg$^{-1}$ K$^{-1}$, $C = 1200$ J K$^{-1}$ kg$^{-1}$, $k'_T = 3.7$ W m$^{-1}$ K$^{-1}$. We take the Clapeyron slope $\gamma' = 1.3 \times 10^{-7}$ K Pa$^{-1}$ (or 3.6 K km$^{-1}$). For initial matrix temperature $T_0 = 1673$ K we get a length scale $H = 115$ km. This choice of parameters leads to dimensionless numbers $\delta = 0.07$, $k'_T = 0.03$, $V_B = 60$, $R = 0.875$ and $\gamma = 0.28$.

<table><caption>**Table 1.** Table of parameters applicable to dry melting below mid-ocean spreading centres.</caption><thead><tr><th>Symbol</th><th>Description</th><th>Definition</th><th>Possible ranges</th><th>Preferred value</th><th>Units</th></tr></thead><tbody><tr><td>μ<sub>f</sub></td><td>Fluid shear viscosity</td><td></td><td></td><td>10</td><td>Pa s</td></tr><tr><td>μ<sub>m</sub></td><td>Matrix shear viscosity</td><td></td><td>10<sup>18</sup>–10<sup>19</sup></td><td>10<sup>18</sup></td><td>Pa s</td></tr><tr><td>k<sub>0</sub></td><td>Constant in permeability relationship</td><td></td><td>10<sup>-10</sup>–10<sup>-9</sup></td><td>5 × 10<sup>-10</sup></td><td>m<sup>2</sup></td></tr><tr><td>c</td><td>Darcy interaction coefficient</td><td>μ<sub>f</sub>/k<sub>0</sub></td><td>10<sup>10</sup>–10<sup>11</sup></td><td>2 × 10<sup>10</sup></td><td>Pa s m<sup>-2</sup></td></tr><tr><td>ρ<sub>f</sub></td><td>Fluid density</td><td></td><td></td><td>2800</td><td>kg m<sup>-3</sup></td></tr><tr><td>ρ<sub>m</sub></td><td>Matrix density</td><td></td><td></td><td>3200</td><td>kg m<sup>-3</sup></td></tr><tr><td>g</td><td>Gravitational acceleration</td><td></td><td></td><td>9.8</td><td>m s<sup>-2</sup></td></tr><tr><td>ΔS</td><td>Entropy of fusion</td><td>s<sub>m</sub> – s<sub>f</sub></td><td>(-250–400)</td><td>-340</td><td>J K<sup>-1</sup> kg<sup>-1</sup></td></tr><tr><td>C</td><td>Heat capacity</td><td></td><td>1000–1300</td><td>1200</td><td>J K<sup>-1</sup> kg<sup>-1</sup></td></tr><tr><td>T<sub>0</sub></td><td>Initial temperature of upwelling</td><td></td><td>1573–1673</td><td>1673</td><td>K</td></tr><tr><td>V</td><td>Initial velocity of upwelling</td><td></td><td>4–10</td><td>10</td><td>cm yr<sup>-1</sup></td></tr><tr><td>H</td><td>Length scale</td><td>ρ<sub>f</sub> T<sub>0</sub>ΔS<sup>2</sup>/(Δρg C)</td><td>60–150</td><td>115</td><td>km</td></tr><tr><td>k'<sub>T</sub></td><td>Thermal conductivity</td><td></td><td></td><td>3.7</td><td>W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>k<sub>T</sub></td><td>Dimensionless thermal conductivity</td><td>ρ<sub>f</sub>k'<sub>T</sub>T<sub>0</sub>/(ρ<sub>m</sub>ΔρgVH<sup>2</sup>)</td><td></td><td>0.03</td><td>—</td></tr><tr><td>γ'</td><td>Clapeyron slope</td><td>Δρ/(ρ<sub>f</sub>ρ<sub>m</sub>ΔS)</td><td>100–133</td><td>130</td><td>K GPa<sup>-1</sup></td></tr><tr><td>γ</td><td>Dimensionless Clapeyron slope</td><td>γ'ρ<sub>m</sub>gH/T<sub>0</sub> = -ΔS/C</td><td></td><td>0.28</td><td>—</td></tr><tr><td>δ'</td><td>Compaction length</td><td>√4μ<sub>m</sub>/(3c)</td><td>8–26</td><td>8</td><td>km</td></tr><tr><td>δ</td><td>Dimensionless compaction length</td><td>δ'/H</td><td>0.07–0.23</td><td>0.07</td><td>—</td></tr><tr><td>V<sub>B</sub></td><td>Buoyancy velocity scale</td><td>Δρg/(cV)</td><td>60–1000</td><td>60</td><td>—</td></tr><tr><td>R</td><td>Density ratio</td><td>ρ<sub>f</sub>/ρ<sub>m</sub></td><td></td><td>0.875</td><td>—</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Porosity profiles in the melting zone in the Darcy equilibrium approximation ($\delta = 0$). All are curves calculated with $\gamma = -\Delta S/C = 0.28$, $R = \rho_f/\rho_m = 0.875$. The values of dimensionless buoyancy velocity $V_B = \Delta\rho g/(cV)$ are indicated. For $V_B$ larger than ~5.9, porosity becomes a multivalued function in the Darcy approximation.</figcaption></figure>

## 3.1 Darcy equilibrium: the outer solution

In magmatic settings the compaction length is usually much smaller than the thickness of the melting zone. The dimensionless compaction length $\delta$ is a small parameter, and therefore, the force balance on the phases should be similar to that with a simpler state where $\delta = 0$. We call this state, where gravity is balanced by Darcy drag, the Darcy equilibrium. In Darcy equilibrium the force balance (43) becomes an algebraic relation where the velocity difference and the matrix and melt velocities are explicit functions of porosity, that is, combining (43) and (44) and using $\Delta v = v_m - v_f$ leads to

$$ \Delta v = -V_B \phi (1 - \phi) \qquad (51) $$

$$ v_m = \frac{1 - R V_B \phi^2 (1 - \phi)}{1 - (1 - R)\phi}, \qquad (52) $$

$$ v_f = \frac{1 + V_B \phi (1 - \phi)^2}{1 - (1 - R)\phi}. \qquad (53) $$

Numerical solutions for the Darcy equilibrium are easily obtained since four of the equations are algebraic in this case (46, 47 with thermal diffusion neglected, 52 and 53) and a single quadrature is needed in (45). Porosity as a function of vertical coordinate is shown in Fig. 1 for several values of the dimensionless buoyancy velocity $V_B$. When $V_B$ is large, (i.e. when $V_B \gtrsim 5.9$ for $R = 0.875$ and $\gamma = 0.28$), porosity becomes a multivalued function of $z$ and Darcy equilibrium cannot be a legitimate solution over the whole extent of the melting region. To understand the limit of the Darcy approximation near the end of the melting zone, a new quantity has to be introduced.

The degree of melting $X$ at a point $z$ is the mass fraction of matrix which has undergone melting all the way from 0 to $z$, which in dimensional quantities is

$$ X(z) = \frac{1}{\rho_m V} \int_0^z \Delta\Gamma(z') dz'. \qquad (54) $$

Contrary to the porosity, which is the instantaneous volume fraction of fluid, the degree of melting accounts for the melting history. Integrating the mass conservation eq. (36) with the appropriate boundary conditions and using eq. (35), the degree of melting can also be written in dimensionless quantities as

$$ X = R\phi v_f = 1 - (1 - \phi)v_m. \qquad (55) $$

The degree of melting is thus also the magma mass flux divided by the total mass flux (Ahern & Turcotte 1979; Ribe 1985b). If the densities are kept distinct but velocities are assumed the same as in the case of batch melting, then we have the relation $X = \phi\rho_f/\bar{\rho}$, which is what McKenzie (1984) and Iwamori et al. (1995) use. The batch melting approximation implies that $X \sim \phi$ as $\rho_m \sim \rho_f$. In the general case however $\phi$ and $X$ can be largely different as $v_f$ can be much larger than $v_m$. As long as melt is lighter than matrix and can escape easily, $X$ can be close to 1 even at small porosity. At the end of the melt zone, $z = z_f$, however, one must have $X = 1$, $\phi = 1$ and $v_f = 1/R$ simultaneously. The Darcy equilibrium solutions for large $V_B$ do not satisfy these boundary conditions; when $X$ reaches 1, $\phi < 1$ and $v_f \gg 1/R$. This implies that near the end of the melting zone $\phi$ has to increase and $v_f$ has to decrease rapidly relative to Darcy equilibrium values. The magnitude of Darcy drag $\Delta v$ thus decreases while the matrix viscous term $\nabla \cdot [(1 - \phi)\bar{\tau}_m]$, or at least its component $-\bar{\tau}_m : \nabla\phi$, increases. At the end of the melting zone, the viscous forces become important relative to Darcy drag even for small values of the compaction length $\delta$, and the Darcy approximation does not hold.
To understand the complete behaviour of the system, we must drop the assumption of Darcy equilibrium and take into account the viscous terms near the origin and near the end of the melting. We focus on the behaviour of our equations at low porosity as this case is more interesting geophysically. An analytical expression is useful to understand the limitations of the Darcy solution near $z = 0$. The matrix and fluid velocities (52) and (53) can be expressed at small $\phi$ as

$$v_m = 1 + (1 - R)\phi - [V_B R - (1 - R)^2]\phi^2 + O(\phi^3), \quad (56)$$

$$v_f = 1 + (V_B + 1 - R)\phi - [V_B(1 + R) - (1 - R)^2]\phi^2 + O(\phi^3). \quad (57)$$

According to (45) and (46) and using (51), we get $\gamma^{-1}dT/dz = -1 + (1 - R)\phi$. As $T = 1$ at $z = 0$, the initial melting rate deduced from (47) is $\Delta\Gamma \approx -\gamma^{-1}dT/dz = 1 - (1 - R)\phi$. Integration of the mass conservation (48) implies in turn that $1 - (1 - \phi)v_m = \int \Delta\Gamma dz'$. Combining these results and introducing the expression for $v_m$ from (56) we get, after some algebra,

$$\phi + (V_B + 1 - R)\phi^2 = \frac{z}{R} - \frac{1 - R}{R} \int_0^z \phi dz' + V_B \int_0^z \frac{\phi^2(1 - \phi)^2}{\gamma^{-1} - z + (1 - R) \int_0^{z'} \phi dz''} dz'. \quad (58)$$

In the typical case, $1 - R \ll 1 \ll V_B$. In the vicinity of $z = 0$ we neglect the two integral terms on the right of (58) of order $(1 - R)\phi z$ and $V_B\phi^2 z$; we also neglect $1 - R$ with respect to $V_B$. We get the z-dependence of porosity

$$\frac{z}{R} \approx \phi + V_B \phi^2. \quad (59)$$

The porosity has, therefore, a parabolic form and varies as $\phi \approx (\sqrt{1 + 4V_Bz/R} - 1)/(2V_B)$. The initial slope at $z = 0$ is equal to $1/R$. Practically for all simulations that we performed, $R$ gives a rough estimate of the total thickness of the two-phase melting zone. The parabolic form of the porosity profile is a consequence of $n = 2$ in the permeability-porosity law (see text after eq. 6). Schmeling (2006) finds identical $z(\phi)$ behaviour to our eq. (59) when he restricts his melt extraction model for plumes to a single component case.

The matrix velocity increases initially with $\phi$ increasing (see 56). This effect is due to the dilation of the matrix by the introduction of less dense melt that cannot easily escape. As soon as the porosity reaches $(1 - R)/(2RV_B)$ (see 56), the draining of the fluid towards the surface allows the matrix to compact and $v_m$ starts decreasing. With parameters appropriate for a ridge, the matrix would inflate until $\phi = 0.001$ which would occur at

$$z_D = \frac{1 - R}{2V_B} \quad (60)$$

after the first 120 m of melting.

Near $z = 0$, however, the Darcy equilibrium cannot be sustained. Indeed, the viscous compaction term that was neglected goes as $\delta^2\phi^{-1}dv_m/dz$ and is thus predicted by the Darcy approximation to vary as $\delta^2(1 - R)/z$; therefore, this term will dominate the Darcy term $\Delta v$ ($\Delta v \approx -V_Bz/R$) at small enough $z$.

## 3.2 Initiation of melting

How melting initiates is not obvious in the framework of the Darcy approximation. The presence of $\phi$ in the denominators of (43), (45) and (47) introduces a singularity at the onset of melting. The first drop of liquid requires more volume than the equivalent mass of solid it replaces. This raises the pressure in the fluid because the necessary dilation of the matrix is inhibited by the large equivalent bulk viscosity $\mu_m/\phi$. This excess pressure tends in turn to bring the thermodynamic conditions back to freezing. From the mass continuity eq. (35) at $z = 0$ where $\phi = 0$, $v_m = 1$ and $dv_m/dz = 0$, one has $[Rv_f - 1]d\phi/dz = 0$. One possibility to satisfy this equation would be $d\phi/dz = 0$ at $z = 0$. However, then from (48) the melting rate $\Delta\Gamma$ would also be zero, that is, at the onset of melting no melting would occur. Therefore, one necessarily has $v_f = 1/R$; unlike in the Darcy approximation there is a finite difference in velocities $\Delta v = -(1 - R)/R$ at the onset of melting. From (45) the gradient in $P_f$ is infinite at $z = 0$ and through the equilibrium melting condition (46), this singularity translates into a singularity of the temperature gradient that should not exist in the presence of thermal diffusion.

In fact, eqs (43)–(48) do not necessarily even have a solution when the starting condition is a pure matrix with zero porosity. There are various possibilities to circumvent this difficulty. One could consider that a finite porosity is present in the solid even before melting starts. This removes the difficulty but in a somewhat unphysical way since it assumes melt is already present before melting starts. One could account for the fluid phase compressibility (Connolly & Podladchikov 1998; Richard et al. 2006) since liquids are indeed more compressible than solids (Stolper et al. 1981). At high pressures this effect might even favour melt initiation by a release of pressure accompanying melting. Another possibility is to remove the discontinuity in temperature gradient at the onset of melting by taking into account the kinetics of the fusion reaction and replace (46) by (41). This would require the knowledge of the reaction rate factor $L_{22}$ which is likely to be a complex function of $\phi$, $T$ and $P_f$. We will further constrain $L_{22}$ in Section 3.2.3. A simpler approach is to accept the singularity of the thermal gradient but to neglect the thermal diffusion in the energy conservation (47). This is the approach that we take in the following of this paper.

### 3.2.1 The squirting solution

The finite value of $\Delta v$ at the onset of melting induces a Darcy friction that can only be balanced by the compactive viscous stress since the buoyancy forces which go as $\Delta\rho\phi$, are zero at $z = 0$ (see 43). We call this equilibrium between the viscous compaction term and the Darcy
drag the 'squirting' solution because the fluid pressure gradient has a singularity at $z = 0$, that is, by (45) and using $\Delta v = -(1 - R)/R$, goes as

$$ \frac{dP_f}{dz} \approx -\frac{(1-R)^2}{RV_B\phi}. \qquad (61) $$

According to the Clapeyron relation (46), the gradient in $P_f$ is also the gradient in temperature $dT/dz$, which, is therefore, also singular at $z = 0$. Furthermore, by matrix mass conservation (48), with the condition $dv_m/dz|_{z=0} = 0$, $\Delta\Gamma$ near $z = 0$ obeys

$$ \frac{d\phi}{dz} \approx \Delta\Gamma. \qquad (62) $$

In the energy eq. (47), assuming the absence of thermal conduction, the singularity in $dT/dz$ can only be balanced by an infinite rate of melting (all the other terms are negligible or zero at $z = 0$)

$$ \Delta\Gamma \approx -\frac{1}{\gamma} \frac{d \ln T}{dz} \approx \frac{(1-R)^2}{RV_B\phi}, \qquad (63) $$

since $T = 1$ at $z = 0$. From (62) and (63) one obtains

$$ \phi \approx (1-R) \sqrt{\frac{2z}{RV_B}}. \qquad (64) $$

With the initial behaviour of $\phi$ determined, we search for a solution of (43)–(48) in the form of a series in $\sqrt{z}$. After a bit of algebra we obtain the following terms of the expansions,

$$ \phi = (1-R) \sqrt{\frac{2z}{RV_B}} + \frac{2}{3} \left[ R + \frac{2\gamma(1-R)^2}{RV_B} \right] z + O(z^{3/2}), \qquad (65) $$

$$ P_f = -(1-R) \sqrt{\frac{2z}{RV_B}} - \frac{2}{3} \left[ R - \frac{\gamma(1-R)^2}{RV_B} \right] z + O(z^{3/2}), \qquad (66) $$

$$ v_m = 1 + \frac{2 z_0}{3 \delta^2} \sqrt{\frac{2 V_B}{R}} z^{3/2} + O(z^2). \qquad (67) $$

These analytical solutions will be compared to the numerical solutions in the Section 3.3. The expansion for $v_m$ remains a function of $z_0$ that cannot be determined by an analytical study at $z = 0$ but is controlled by the boundary condition at the end of the melting zone. Getting an analytical estimate for $z_0$ is feasible but cumbersome. It consists of using the boundary condition $2\delta^2(d\phi/dz)(dv_m/dz) + \Delta v = 0$ at $z = z_f$, obtained from (43) with the use of $\phi = 1$, with the expansions previously found for $\phi$ and $v_m$, the end of the melting zone $z_f$ being obtained by solving $\phi(z_f) = 1$. The sign of $z_0$ is, however, easy to estimate from a simple physical consideration. At low buoyancy velocity $V_B$, the melt is somewhat locked in the matrix, and the volume expansion associated with first melting induces matrix dilation at the base of the melting zone; therefore, according to (50), $z_0$ will be positive. At high buoyancy velocity $V_B$, the easy extraction of melt allows the matrix to readily compact upon first melting, and $z_0$ is negative. This change in sign of $z_0$ will be demonstrated later in Fig. 7. During the initial squirting equilibrium the porosity increases and for large $V_B$ the gravity term $V_B\phi(1 - \phi)$ can overcome the Darcy term $\Delta v = -(1 - R)/R$. This may lead to a new force balance.

### 3.2.2 Viscogravitational equilibrium

The potential of balance between viscous compaction and buoyancy is called viscogravitational equilibrium. When the Darcy term $\Delta v$ is negligible the pressure is simply (see 45 and 46)

$$ P_f = (T - 1)\gamma \approx -Rz. \qquad (68) $$

From the energy eq. (47) the melting rate is simply

$$ \Delta\Gamma \approx \frac{R}{1 - R\gamma z}. \qquad (69) $$

Knowing the melting rate, the other variables are readily computed; one obtains from (48) and (50)

$$ \phi \approx Rz, \qquad (70) $$

$$ v_m \approx 1 + \frac{z_0 V_B R}{2\delta^2(1-R)} z^2, \qquad (71) $$

where the expansion of both $\phi$ and $v_m$ was truncated after the first non-zero term in $z$. Again, as in the squirting approximation, the value of $z_0$ can only be obtained by solving the equations in the whole integration domain. However, the sign of $z_0$ can be determined: Darcy drag is negligible in this approximation, therefore, the easy extraction of melt will allow the matrix to compact and $z_0$ will be negative.
### 3.2.3 Ordering of boundary regions

Up to now we have defined three possible force equilibria: Darcy, squirting and viscogravitational. In each case one force term was neglected: the viscous compaction term, the buoyancy term, $V_B\phi(1 - \phi)$, and the Darcy drag term, $\Delta v$, respectively. From the first order solutions given in each domain (56, 57 and 59 for Darcy equilibrium; 65–67 for squirting approximation; 68–71 for viscogravitational balance), these neglected terms amount to $\delta^2(1 - R)/(Rz)$ (viscous compaction term), $(1 - R)\sqrt{2V_Bz}/R$ (buoyancy force), and $-(1 - R)/R$ (Darcy drag), respectively. We compare the magnitude of the neglected terms between pairs of equilibrium approximations and introduce three coordinates

$$z_1 = \frac{1}{2RV_B} \qquad \text{squirt. } (z < z_1) \text{ vs. viscograv. } (z > z_1), \qquad (72)$$

$$z_2 = \delta^2 \qquad \text{viscograv. } (z < z_2) \text{ vs. Darcy } (z > z_2), \qquad (73)$$

$$z_3 = \sqrt[3]{\frac{\delta^4}{2RV_B}} \qquad \text{squirting } (z < z_3) \text{ vs. Darcy } (z > z_3). \qquad (74)$$

The squirting solution is thus valid from $z = 0$ to min ($z_1, z_3$), above which is possibly the viscogravitational equilibrium up to $z_2$. These inner solutions then give place to the outer Darcy equilibrium. As $z_1z_2^2 = z_3^3$, the three domains have a common point where $2\delta^2RV_B = 1$.

This study of the squirting solution nearest the origin can be used to illustrate the problems when the thermal diffusivity is not neglected. The derivative of the fluid pressure remains always singular (proportional to $1/\phi$, see 61) although the presence of the thermal diffusivity in the energy eq. (47) forbids the singularity of $dT/dz$. There is, therefore, no possibility to insure the thermodynamic equilibrium. We could instead replace the thermodynamic equilibrium (46) by the kinetic law (41), but this puts some constraints on the functional form of the kinetic coefficient $L_{22}$. For example, if $\phi$ varies like $z^n$, the leading term of $P_f$ will be in $z^{n-1}$ (see 61) and that of $\Delta\Gamma$ in $z^{n-1}$ (see 62). To balance these terms into the kinetic equation we need to choose $L_{22} \propto z^{2n-2} \propto \phi^{2-2/n}$. To have $P_f(0) = 0$ we need $n < 1$, and therefore, $L_{22}$ should depend on a negative power of $\phi$. The kinetic rate should, therefore, be infinite at $z = 0$ for a solution with diffusion to exist. Rather than exploring this situation, and in the absence of a clear experimental or theoretical form for the kinetic rate, we simply assume that the very small parameter $k_T$ is zero. In fact a layer dominated by thermal conduction should exist close to $z = 0$. Its thickness

$$z_T = \frac{k'_T}{\rho_m C V} \qquad (75)$$

should be of order 300 m.

We summarize the various domains near $z = 0$ in Fig. 2. The first domain of thickness $k'_T/(\rho_m CV)$ corresponds to a diffusive layer. It has not been explicitly studied as it requires one to take into account the kinetics of melting which would introduce various additional unknowns and difficulties. The second domain corresponds to the domain of the squirting equilibrium. It is replaced either by the Darcy equilibrium at low compaction lengths or by the viscogravitational equilibrium at larger compaction lengths. Near the end of the melting zone, where $z > R$ and $\phi \lesssim 1$, the viscous force becomes again important and limits the validity of the Darcy approximation. This diagram has been computed for $V_B = 60$ and changes with $V_B$, however its topology remains similar. For smaller values of $V_B$, one can distinguish a dilational and a compacting regime in the Darcy domain ($z_D$, which was defined in Section 3.1, becomes larger than the thermal diffusive layer thickness $z_T$). In the situation of the Earth the force equilibrium changes directly from squirting regime to Darcy equilibrium even for the least viscous magmas.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Graphical representation of the various regions along the melting zone (dimensionless $z$ coordinate plotted on the vertical axis in logarithmic scale, melting starts at $z = 0$) as a function of the dimensionless compaction length $\delta = \sqrt{4\mu_m/(3cH^2)}$ (on horizontal axis with logarithmic scaling). Plotted for buoyancy velocity $V_B = \Delta\rho g/(cV) = 60$. The compaction length appropriate for melting under a ridge is indicated by an arrow.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Calculated porosity profile in the melting zone (thick solid line) compared to the Darcy approximation (dotted line) and the initial squirting equilibrium approximation (dashed line). The details of the curves near the origin are zoomed in the embedded panel. All curves are calculated with $\gamma = 0.28$, $R = 0.875$ and $V_B = 10$, $\delta = 0.26$. As expected the viscous equilibrium prevails in a very thin boundary layer near the origin.</figcaption></figure>

### 3.3 Numerical solutions

The full system of eqs (43)–(48) can be solved numerically. We choose a domain of integration from $z = 0$ to a sufficiently large $z_{max}$ (remember that the width of the melting zone is yet unknown). From a guessed porosity we obtain the matrix velocity using a tridiagonal solver of eq. (43) with boundary conditions $v_m = 1$ at $z = 0$ and $dv_m/dz = 0$ at $z = z_{max}$. Fluid pressure is then calculated from (45) and the temperature from (46). We use the inner squirting solution, (66), to estimate $P_f$ for the first two grid points to avoid the numerical indeterminacy of the integral of $\Delta v/\phi$ near $z = 0$. The energy eq. (47) is used to calculate the melting rate. Finally, porosity is updated from eq. (48) (we actually use the time dependent form of this equation, i.e. 36). We iterate this scheme to obtain the steady-state solution.

This numerical solution can be compared to the analytical solutions (65) and (59). For this case we choose a rather low buoyancy velocity $V_B = 10$ and rather large compaction length $\delta = 0.26$ in order to emphasize the thickness of the boundary region near $z = 0$ where the viscous term is dominant. The results are depicted in Fig. 3 where the complete numerical solution (solid line) is compared to the Darcy solution (dotted line) and the squirting solution (dashed line). This figure confirms the accuracy of the analytical solutions and shows that the Darcy approximation remains very close to the real solution over the 2/3 of the melting zone. The blow up panel shows the initial behaviour of the numerical solution and confirms the validity of the squirting approximation in the vicinity of $z = 0$.

Porosities for three different values of the compaction length $\delta$ are shown in Fig. 4. When $\delta$ is small compared to the width of the melting zone, the porosity follows the $\delta = 0$ Darcy equilibrium profile (gravity is balanced by Darcy resistance). At the end of the melting zone there is a narrow region where the contribution of viscous forces becomes important again relative to Darcy drag. In this region porosity increases rapidly with $z$ until it reaches 1 where melting is complete. The overall thickness of the two-phase melting zone increases with $\delta$.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Porosity profiles in the melting zone for different values of the compaction length $\delta$. All curves are calculated with $\gamma = 0.28$, $R = 0.875$ and $V_B = 100$. Values of the compaction length are 2.6 (300 km; dashed line), 0.87 (100 km; dot-dashed line) and 0.26 (30 km; solid line). Porosity in Darcy equilibrium ($\delta = 0$) is shown as thin dotted line.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Degree of melting in the melting zone for different values of the compaction length $\delta$. Values of parameters and line textures correspond to those in Fig. 4.</figcaption></figure>

This behaviour is due to the contribution to the fluid pressure gradient by the viscous deformation (the last term in eq. 49 that is multiplied with $\delta^2$). In the limit of very large compaction length the force balance in the melting zone is basically between the viscous compaction and buoyancy, and the difference in velocities remains small between the two phases (numerical results for velocities are shown in Fig. 6). Then the pressure gradient in the melt is $dP_f/dz \approx -R$ (i.e. the fluid pressure is controlled by the magma density; see 68, also 49). In the limit of very small compaction length one can see from eq. (49) that $dP_f/dz \approx -1$ (i.e. the fluid pressure is controlled by the matrix density, since porosity also remains small in most of the melt zone; see Fig. 4). When $\delta$ is comparable to (or exceeds) the size of the melting zone, the deviation from the Darcy equilibrium is more pronounced. This case however is far from geophysical interest: in the ridge case when the buoyancy velocity is large and the Darcy equilibrium is close to the real force balance, the porosity profile is essentially identical to the Darcy solution, from which it deviates only at the very end of the melting zone. Taking the limit $V_B \gg 1$ at the end of the melting zone ($z \approx R$) in the approximate Darcy equilibrium solution (59) suggests that the porosity jumps abruptly from a value close to $\phi \approx R^{-1}V_B^{-1/2}$ (or about 0.1 for the characteristic values of $R \approx 0.9$ and $V_B \approx 100$) to $\phi = 1$. This means that even in the case of extensive melting (degree of melting of several tens per cent) the porosity remains quite small ($<10$ per cent).

Curves in Fig. 4 depict porosity in the case of complete melting of a univariant material. In real geophysical situations the matrix is a multicomponent material which partially melts up to some limited degree of melting. In the case of pressure release melting under mid-ocean ridges the degree of melting typically does not exceed 20–30 per cent (Turcotte & Phipps Morgan 1992). In the numerical simulations the degree of melting increases in an approximately linear fashion throughout the melting zone as is seen in Fig. 5. This behaviour is very different from the parabolic evolution of $\phi$. In particular, when $V_B$ is very large, $\phi$ is much smaller than $X$ everywhere except at the end of the melting zone. Using the values of parameters suggested in the previous section, porosity is only ~0.05 when $X = 0.2$. Although the porosity profile is very dependent on the compaction length, the melting rate is not.

Fig. 6 shows matrix (left panel) and fluid (right panel) velocities. The matrix enters the melting zone at dimensionless velocity $v_m = 1$ and decreases regularly during the compaction process. The first drop of melt is already moving at velocity $v_f = 1/R$ as already seen in Section 3.2. This velocity is also the melt velocity at the end of the melting zone. For the figure we used compaction lengths larger than

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Matrix (left panel) and fluid (right panel) velocities in the melting zone. Values of parameters and line textures correspond to those in Fig. 4.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Shift of the melting $z_0$ zone as a function of the compaction length $\delta$ for the following values of the buoyancy velocity: $V_B = 100$ (solid line, diamonds), $V_B = 10$ (dashed line, circles), $V_B = 1$ (dot-dashed line, triangles), $V_B = 0.1$ (dotted line, squares). Both dimensionless values (bottom and left axes) and values in km (top and right axes) are shown. Positive (negative) $z_0$ means shallower (deeper) melting. Calculated with $\gamma = 0.28$, $R = 0.875$.</figcaption></figure>

appropriate for a ridge. The case with $\delta = 0.07$ (8 km) would however be very similar to the case $\delta = 0.26$ but the final viscous boundary layer would be somewhat thinner. The maximum fluid velocity is about 10 times the velocity of the upwelling mantle, that is of order of $1 \text{ m yr}^{-1}$. The matrix velocity curves are terminated at the end of the melting zone; no matrix is present above the partial-melt zone, and matrix velocity, is therefore, undefined here.

In the present theory the two phases are submitted to different pressure fields. The fluid pressure $P_f$, the matrix pressure $P_m$ and the average pressure $\bar{P}$ are different. The pressure difference $\Delta P$ is a dynamic effect that comes from changes in porosity in excess to those due to melting as discussed in Section 2.3. In the limit where the fluid phase is much less viscous than the matrix, the pressure difference is directly related to deformation of the matrix. It is positive when matrix compacts and negative when matrix dilates (eq. 40). The observer at the surface of the Earth does not know the variation of the fluid pressure with depth. He may believe to have a better understanding of the variations with depth of the matrix pressure, $P_m = P_f + \Delta P$, or the average pressure, $\bar{P} = P_f + (1-\phi)\Delta P$ and may expect melting to start at the melting conditions $P_m = 0$ or $\bar{P} = 0$ when $T = 1$. However, at the onset of melting where $P_f = 0$, $P_m = \bar{P} = \Delta P = -z_0$ (see discussion around eq. 50). As a consequence the melting does not start where expected but the onset of melting is shifted by the distance $z_0$. When $z_0$ is negative, the fluid pressure is lower than the average pressure which favours melting at greater depths. On the contrary a positive $z_0$ implies higher pressure in the melt, and therefore, the melting condition is only reached at shallower depth.

Fig. 7 shows $z_0$ as a function of the compaction length $\delta$ for several values of the dimensionless buoyancy velocity $V_B$. No simple functional dependence of $z_0$ on $\delta$ of $V_B$ is found. As seen from Fig. 7, $z_0$ can be positive or negative depending on the value of $V_B$, the change of sign occurring near $V_B \sim 1$. When the buoyancy velocity is small, the melt extraction from the matrix is difficult. The only way to compensate for the bulk volume increase upon melting is that both matrix and fluid velocities increase. Pressure in the fluid is in this case larger than in the matrix, and melting starts at a shallower depth ($z_0 > 0$). This corresponds to the case of squirting equilibrium. On the other hand when $V_B \gg 1$ (viscogravitational balance), the melt extraction is efficient and the matrix compacts while the fluid velocity increases. Pressure is in this case decreased in the fluid and the melting can start deeper ($z_0 < 0$). When $\delta$ decreases, the Darcy equilibrium prevails and as expected, $z_0$ tends to zero. An asymptotic behaviour is obtained at very large compaction lengths. For values appropriate for a ridge the melting is likely to occur at most a few kilometres deeper than what could be expected from the average pressure.

The sign of $z_0$ is indicative of the sign of the difference between pressure of the two phases $\Delta P$. In Fig. 8 we show the pressure difference in the melting region as a function of $z$ for three values of the buoyancy velocity $V_B$ and a fixed value of the compaction length $\delta = 0.26$. This value corresponding to 30 km is somewhat too large for a ridge but it is used to emphasize the difference between the various cases. In the extreme case of very large or very small buoyancy velocity the melt is either uniformly under-pressured or over-pressured with respect to the matrix. For intermediate values of $V_B$ the pressure difference can change sign twice. The pressure difference between matrix and fluid is about 5 MPa (50 bar) even far from the incipient melting for typical ridge parameters. The fact that the fluid is under-pressured should somewhat inhibit hydro-fracturing of the matrix by the rising magma. For moderate buoyancy velocities, however (dot-dashed line), there is significant tendency for the fluid pressure to overcome that of the matrix.

# 4 DISCUSSION AND CONCLUSION

In this paper we developed a theoretical model of two-phase flow in the presence of melting in a single component system. Our model properly accounts for the feedback between the viscous deformation of the phases and the thermodynamic conditions of melting. As in previous papers
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Pressure difference $\Delta P = P_m - P_f$ in the melting zone calculated with $\delta = 0.26$ (30 km) and the following values of the buoyancy velocity: $V_B = 100$ (solid line), $V_B = 1$ (dot-dashed line), $V_B = 0.1$ (dotted line). Both dimensionless values of pressure difference (bottom axis) and values in MPa (top axis) are shown. Calculated with $\gamma = 0.28$, $R = 0.875$.</figcaption></figure>

(Bercovici et al. 2001a; Ricard et al. 2001; Bercovici & Ricard 2003), we accounted for the presence of surface tension on the interface between the two phases which imposes a clear distinction between the properties of each phase and those of the interface (pressures, velocities and densities). The conditions of equilibrium between the two phases are naturally deduced from the second law of thermodynamics. The usual Clapeyron slope is affected by the presence of the surface tension (Gibbs–Thomson effects) and by the dynamic pressure difference between the phases. This pressure difference is proportional to the rate of matrix compaction and to the inverse of porosity; its role is somewhat equivalent to the bulk viscosity of other theories (McKenzie 1984). It is only in the case of motionless mixture without surface tension that the usual Clapeyron slope is recovered. In the motionless case with surface tension, the pressure difference between the phases verifies the Laplace’s condition but the phase equilibrium is affected by the Gibbs–Thomson effect. The model also allows for non-equilibrium situations, in which case a kinetic relation links the melting rate to the departure from equilibrium. The general non-equilibrium treatment will be useful in future extended models where the realistic chemical composition of Earth material will be taken into account.

We use our equations in the simple case of 1-D melting. The equations are very complex by the presence of various singularities near $\phi = 0$. Neither the length of the melting zone nor the position of melt initiation ($z_0$) are known *a priori*. The solutions are primarily controlled by two parameters, the compaction length $\delta$ and the buoyancy velocity $V_B$. We discuss the existence of various possible boundary regions near the onset of melting. The onset of melting corresponds to a balance between the viscous forces and the Darcy friction. This first boundary region is followed at large enough $\delta^2 V_B$ by a zone of balance between gravity forces and viscous matrix stresses. Then the Darcy balance (equilibrium between gravity force and Darcy friction) dominates in most of the melting region. The Darcy equilibrium can itself be divided into a domain where the fluid does not escape fast enough to avoid matrix dilation (caused by volume increase upon melting), and a domain where matrix compacts while the fluid is extracted. This already very complex structure should also include a thin thermal diffusion boundary zone (that was neglected in this study) which cannot be understood at exact thermodynamic equilibrium unless the kinetics of the melting is taken into account. The effect of dynamic pressure on the onset of melting can raise or depress the depth of melting by few kilometres with respect to what would be expected from the average pressure. With parameters appropriate for melting under an oceanic ridge ($\delta \sim 0.1$, $V_B \sim 100$) all these boundary zones occur in the first few kilometres of the two-phase region. The melting starts up to a few km deeper that what could be expected from matrix pressure only.

The 1-D model suggests that for $V_B \gtrsim 1$, $z_0$ is negative, and therefore, melting begins deeper than standard Clapeyron slope predicts; likewise for $V_B \lesssim 1$, $z_0$ is positive, and melting does not begin until lower pressure is reached. Parameters that determine the value of $V_B$ are the density difference between phases, average grain size, upwelling velocity, and melt viscosity, the latter parameter exhibiting the largest variation in Earth. Thus, for example, differences in melt viscosity at two otherwise similar geological locations would result in different depths of melting; a more mafic magma of lower viscosity at oceanic arcs would facilitate deeper melting compared to more felsic magma of higher viscosity at continental arcs. We stress that these predictions are based on a fairly simple model and that further refinement of the description is needed to generate more reliable predictions. In particular, the effect of surface tension, which was neglected in the 1-D model application, could either facilitate or inhibit melting, depending on the geometry of the first melt to appear; localized spheroidal pockets of first melt of positive average curvature would hamper melting by Gibbs–Thomson effect, while melting along grain boundaries resulting in melt network of negative curvature would favour melting. More experimental evidence is needed on the geometrical distribution of the first melt within the solid matrix. Moreover, the effect of viscous compaction/dilation on solid–melt phase equilibria of multicomponent Earth
material should be studied. In fact, the main asset of the present study, rather than offering accurate predictions of geological observables, is the development of a framework to consistently account for the interplay between two-phase viscous flow and the melting condition.

The results of our 1-D model cannot be directly compared to those obtained by other authors. Most previous studies had not attempted to solve the self-consistent set of equations but have employed various approximations that do not seem to be supported by our study. Various authors have imposed the melting rate as a simple function of depth (average pressure). This decouples the melting problem from the compaction problem and may generate a spurious layer of compaction where the porosity decreases (Turcotte & Phipps Morgan 1992). In our simulation the porosity is always increasing function of the vertical coordinate and the porosity remains much smaller than the degree of melting ($\phi < V_B^{-1/2}$). With ridge parameters this porosity corresponds to a threshold ~10 per cent after which the melt accumulates instantaneously at some depth and the porosity jumps to 1 above it. This contradicts the adiabatic batch melting model where porosity and degree of melting are closely related (McKenzie 1984). Our results are more similar to those of Fowler (1989), at least in their goals, although significant difference in the starting equations makes a detailed comparison difficult: however, we agree with him about the presence of significant pressure difference (5 MPa or 50 bar) between the two phases all along the melting zone.

The present theory offers a framework for treatment of physical situations where melting of viscously deformable media is concerned. Although our focus was on geological settings, the description is applicable to much wider spectrum of problems. Various extensions to what was discussed in the present study are possible, including melting of multivariant material, non-equilibrium time-dependent melting, and 2-D (3-D) flow modelling.

## ACKNOWLEDGMENTS

The authors benefited from reviews by Harro Schmeling and Michel Rabinowicz. Support was provided by the National Science Foundation (NSF, grant EAR-0537599), and the Centre National de la Recherche Scientifique (CNRS).

## REFERENCES

* Ahern, J.L. & Turcotte, D.L., 1979. Magma migration beneath an ocean ridge, *Earth planet. Sci. Lett.*, **45**(1), 115-122.

* Bailyn, M.A., 1994. *Survey of Thermodynamics*, AIP Press, Woodbury, NY.

* Bear, J., 1988. *Dynamics of Fluids in Porous Media*, Dover Publications, New York.

* Bercovici, D. & Karato, S., 2003. Whole mantle convection and the transition-zone water filter, *Nature*, **425**(6953), 39-44.

* Bercovici, D. & Ricard, Y., 2003. Energetics of a two-phase model of lithospheric damage, shear localization and plate boundary formation, *Geo-phys. J. Int.*, **152**(3), 581-596.

* Bercovici, D. & Ricard, Y., 2005. Tectonic plate generation and two-phase damage: Void growth versus grain size reduction, *J. geophys. Res.*, **110**(B3), B03401.

* Bercovici, D., Ricard, Y. & Schubert, G., 2001a. A two-phase model for compaction and damage, part 1: general theory, *J. geophys. Res.*, **106**(B5), 8887-8906.

* Bercovici, D., Ricard, Y. & Schubert, G., 2001b. A two-phase model for compaction and damage, part 3: application to shear localization and plate boundary formation, *J. geophys. Res.*, **106**(B5), 8925-8939.

* Choblet, G. & Parmentier, E.M., 2001. Mantle upwelling and melting eneath slow spreading centers: effects of variable rheology and melt productivity, *Earth planet. Sci. Lett.*, **184**(3-4), 589-604.

* Connolly, J.A.D. & Podladchikov, Y.Y., 1998. Compaction-driven fluid flow in viscoelastic rock, *Geodinamica Acta*, **11**(2-3), 55-84.

* de Groot, S.R. & Mazur, P., 1984. *Non-Equilibrium Thermodynamics*, Dover Publications, New York.

* Drew, D.A., 1983. Mathematical modelling of two-phase flow, *Ann. Rev. Fluid Mech.*, **15**, 261-291.

* Fowler, A.C., 1985. A mathematical model of magma transport in the asthenosphere, *Geophys. Astrophys. Fluid Dyn.*, **33**(1), 63-96.

* Fowler, A.C., 1989. Generation and creep of magma in the Earth, *SIAM J. Appl. Math.*, **49**(1), 231-245.

* Frank, F.C., 1968. Two-component flow model for convection in the Earth's upper mantle, *Nature*, **220**(5165), 350-352.

* Hier-Majumder, S., Ricard, Y. & Bercovici, D., 2006. Role of grain boundaries in magma migration and storage, *Earth planet. Sci. Lett.*, **248**(3-4), 735-749.

* Hirschmann, M.M., 2000. Mantle solidus: Experimental constraints and the effects of peridotite composition, *Geochem. Geophys. Geosyst.*, **1**(10).

* Iwamori, H., McKenzie, D. & Takahashi, E., 1995. Melt generation by isentropic mantle upwelling, *Earth planet. Sci. Lett.*, **134**(3-4), 253-266.

* Karato, S., 1989. Defects and plastic deformation in olivine, in *Rheology of Solids and of the Earth*, pp. 176-208, eds Karato, S. & Toriumi, M., Oxford University Press, New York.

* Karato, S. & Jung, H., 2003. Effects of pressure on high-temperature dislocation creep in olivine, *Philos. Mag.*, **83**(3), 401-414.

* Kojitani, H. & Akaogi, M., 1995. Measurement of heat of fusion of model basalt in the system diopside-forsterite-anorthite, *Geophys. Res. Lett.*, **22**(17), 2329-2332.

* Landau, L.D. & Lifshitz, E.M., 1959. Fluid mechanics, *Course of Theoretical Physics*, Vol. 6, Pergamon Press, New York.

* Lasaga, A.C., 1998. *Kinetic Theory in the Earth Sciences*, Princeton University Press, Princeton, New Jersey.

* McKenzie, D., 1984. The generation and compaction of partially molten rock, *J. Petrol.*, **25**(3), 713-765.

* McKenzie, D. & Bickle, M.J., 1988. The volume and composition of melt generated by extension of the lithosphere, *J. Petrol.*, **29**, 625-279.

* Mei, S. & Kohlstedt, D.L., 2000. Influence of water on plastic deformation of olivine aggregates, 2. dislocation creep regime, *J. geophys. Res.*, **105**(B9), 21471-21481.

* Nye, J.F., 1953. The flow law of ice from measurements in glacier tunnels, laboratory experiments and the Jungfraufirn borehole experiment, *Proc. R. Soc. London A*, **219**, 477-489.

* Ohtani, E., Litasov, K., Hosoya, T., Kubo, T. & Kondo, T., 2004. Water transport into the deep mantle and formation of a hydrous transition zone, *Phys. Earth planet. Int.*, **143**, 255-269.

* Rabinowicz, M., Ceuleneer, G. & Nicolas, A., 1987. Melt segregation and flow in mantle diapirs below spreading centers: evidence from the Oman ophiolite, *J. geophys. Res.*, **92**(B5), 3475-3486.

* Ribe, N.M., 1985a. The deformation and compaction of partial molten zones, *Geophys. J. R. astr. Soc.*, **83**(2), 487-501.

* Ribe, N.M., 1985b. The generation and composition of partial melts in the Earth's mantle, *Earth planet. Sci. Lett.*, **73**(2-3), 361-376.

* Ribe, N.M., 1988. On the dynamics of mid-ocean ridges, *J. geophys. Res.*, **93**(B1), 429-436.

* Ricard, Y. & Bercovici, D., 2003. Two-phase damage theory and crustal rock failure: the theoretical 'void' limit, and the prediction of experimental data, *Geophys. J. Int.*, **155**(3), 1057-1064.

* Ricard, Y., Bercovici, D. & Schubert, G., 2001. A two-phase model for compaction and damage, part 2: applications to compaction, deformation and
the role of interfacial surface tension, *J. geophys. Res.*, **106**(B5), 8907–8924.

Richard, G., Monnereau, M. & Rabinowicz, M., 2006. Slab dehydration and fluid migration at the base of the upper mantle: implication for deep earthquake mechanisms, *Geophys. J. Int.*, in press.

Schmeling, H., 2000. Partial melting and melt segregation in a convecting mantle, in *Physics and Chemistry of Partially Molten Rocks*, pp. 141–178, eds Bagdassarov, N., Laporte, D. & Thompson, A.B., Kluwer Academic Publisher, Dordrecht.

Schmeling, H., 2006. A model of episodic melt extraction for plumes, *J. geophys. Res.*, **111**, B03202.

Scott, D.R. & Stevenson, D.J., 1984. Magma solitons, *Geophys. Res. Lett.*, **11**(11), 1161–1164.

Scott, D.R. & Stevenson, D.J., 1986. Magma ascent by porous flow, *J. geo-phys. Res.*, **91**(B9), 9283–9296.

Sleep, N.H., 1974. Segregation of magma from a mostly crystalline mush, *Geol. Soc. Am. Bull.*, **85**(8), 1225–1232.

Spiegelman, M., 1993a. Flow in deformable porous media. part 1 Simple analysis, *J. Fluid Mech.*, **247**, 17–38.

Spiegelman, M., 1993b. Flow in deformable porous media. part 2: numerical analysis—the relationship between shock waves and solitary waves, *J. Fluid Mech.*, **247**, 39–63.

Spiegelman, M., 1993c. Physics of melt extraction: Theory, implications, applications, *Trans. R. Soc. Lond. A*, **342**(1663), 23–41.

Stolper, E., Walker, D., Hager, B.H. & Hays, J.F., 1981. Melt segregation from partially molten source regions: the importance of melt density and source region size, *J. geophys. Res.*, **86**(B7), 6261–6271.

Turcotte, D.L. & Phipps Morgan, J., 1992. The physics of magma migration and mantle flow beneath a mid-ocean ridge, in *Mantle Flow and Melt Generation at Mid-Ocean Ridges*, pp. 155–182, eds Phipps Morgan, J., Blackmann, D.K. & Simpson, J.M., American Geophysical Union, New York.

White, R. & McKenzie, D., 1989. Magmatism at rift zones—the generation of volcanic continental margins and flood basalts, *J. geophys. Res.*, **94**(B6), 7685–7729.

Williams, Q. & Garnero, E.J., 1996. Seismic evidence for partial melt at the base of Earth’s mantle, *Science*, **273**(5281), 1528–1530.