# Slab dehydration and fluid migration at the base of the upper mantle: implications for deep earthquake mechanisms

Guillaume Richard, Marc Monnereau, Michel Rabinowicz

► To cite this version:

Guillaume Richard, Marc Monnereau, Michel Rabinowicz. Slab dehydration and fluid migration at the base of the upper mantle: implications for deep earthquake mechanisms. Geophysical Journal International, 2007, 168 (3), pp.1291-1304. 10.1111/j.1365-246X.2006.03244.x. hal-02566101

HAL Id: hal-02566101

https://hal.science/hal-02566101

Submitted on 8 Nov 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
# Slab dehydration and fluid migration at the base of the upper mantle: implications for deep earthquake mechanisms

Guillaume Richard,<fnref n="1" /> Marc Monnereau<fnref n="2" /> and Michel Rabinowicz<fnref n="2" />

Accepted 2006 October 3. Received 2006 October 2; in original form 2005 July 2

## SUMMARY

Water enters the Earth's mantle via subduction of oceanic lithosphere and sediments. A lot of this water immediately returns to the atmosphere through arc volcanism, but part, retained in Dense Hydrous Magnesium Silicates (DHMSs) and Nominally Anhydrous Minerals (NAMs) like olivine, is expected to be subducted as deep as the bottom of the upper mantle (660 km depth). Then, due to its low solubility in lower mantle minerals, water is likely to be released as a hydrated fluid during the spinel–post-spinel phase change. The dynamics of this fluid phase is investigated through a 1-D model of compaction, in which a source term has been introduced to take the fluid precipitation into account.

The competition between the advective transport by the descending slab and the buoyant rise of the fluid results in three distinct situations, depending on the properties of the solid and the fluid phases. Low matrix permeability and high fluid viscosity inhibit compaction and favour the entrainment of fluid towards the deep mantle. In this case, the entire slab water content would enter the lower mantle and would be mixed at large scale. However, realistic values of the fluid viscosity and matrix permeability make this possibility unlikely. When effective, compaction results in an accumulation of fluid at and below the phase boundary. Then, depending on the value of the matrix viscosity, the situation evolves differently. Above $10^{20}$ Pa s, accumulation of fluid extends below the phase boundary and the pressure difference between the fluid and the matrix increases continuously, exceeding the yield strength of rocks. As a result, cracks would form and evolve towards the formation of dykes. In case of very low mantle viscosity, possibly due to strong grain size reduction during phase change, compaction becomes very efficient and the fluid remains confined within the phase change horizon, without increasing pressure. In the long term, this last situation appears unstable and would also evolve towards the formation of dykes.

Thus, it is expected that water returns to the upper mantle by dykes propagating in the direction of the maximum compressive stress. Since maximum compressive stress appears to follow the dip of the slab below 410 km depth, we predict the formation of dykes extending from the 660 km phase change to 410 km depth. The possible existence of such dykes in slabs offers the necessary conditions for strong double couple component earthquake in the deepest part of the upper mantle.

**Key words:** deep focus earthquake, fluids in rocks, mantle convection, melt segregation, percolation, phase transitions, rock fracture, seismicity, subduction, transition zone.

## 1 INTRODUCTION

For more than a decade, the presence of water in the Earth's mantle has been thought to play a key role in geodynamic processes (Ahrens 1989; Bell & Rossman 1992; Thompson 1992). Its role in the processes related to the origin of the continents has been recognized for a long time (Campbell & Taylor 1983). Moreover, the dramatic effects of water on various rock properties suggest a control on

other geophysical processes and over regions wider than the direct vicinity of slabs. In particular, if the water injected in the mantle at trenches were not released by the dehydration accountable for the arc volcanism, the remaining part would be carried and dispersed at greater depths. The amount of water recycled into the deep mantle remains a matter of debate, and depends essentially on the geotherm followed by the slab and the stability of Dense Hydrous Magnesium Silicates (DHMS) (Luth 1995; Angel et al. 2001; Komabayashi et al.

<footnote n="1">¹Department of Geology and Geophysics, Yale University, 210 Whitney Avenue, New Haven, Connecticut 06520, USA. E-mail: guillaume.richard@yale.edu</footnote>

<footnote n="2">²UMR5562, Observatoire Midi-Pyrénées, 14 Av. E. Belin, 31400 Toulouse, France</footnote>
2005). As a consequence, the properties of hydrated regions will be strongly different from dehydrated ones. For instance, the low velocity zone below the lithosphere, initially interpreted as a partially melted zone, may be viewed as a hydrated zone (Karato & Jung 1998). Further, some hotspots might be 'wet spots' (Bonatti 1990), and more recently, it has been proposed that the high depletion in trace elements of the upper mantle would result from dehydration melting above the 410 km discontinuity (Bercovici & Karato 2003).

These hypotheses, although attractive, remain based on speculations concerning the water distribution inside the mantle. Unfortunately, the high diffusivity of hydrogen ions in minerals (just one order of magnitude lower than the thermal diffusivity) makes extraction of information on the water content from deep mantle regions extremely difficult. Thus, assumptions on the water distribution mainly derive from our knowledge on the solubility of water in the various minerals composing the mantle rocks. It is largely accepted that the water solubility (2–3 wt per cent) in the transition zone minerals (410–660 km depth) is higher than in the upper part of the mantle (10–410 km depth) (Koga et al. 2003). The water solubility in the lower mantle is still controversial (Litasov et al. 2003), but a consensus is arising around the weak storage capacity of deep minerals, that is, two orders of magnitude lower than in the transition zone (see Bolfan-Casanova 2005, for review). These values indicate that the transition zone is an important potential reservoir for water on Earth. Most models where water plays an important role assume this reservoir relatively full, that is, that water content reflects the solubility. This assumption is justified only at equilibrium. Despite the strong solubility drop across the spinel–post-spinel phase change, the diffusion of water is not high enough to overcome the convective transport and to hold back the water in the transition zone (Richard et al. 2002). As a consequence, a probable exsolution of water from minerals occurs when they enter the lower mantle. Incidentally, it also means that if some hydrous phases are still able to hold water in the lower mantle (in a very cold slab for example) this water will be well-mixed by convection into the whole mantle. The possible fluid precipitation at 660 km depth is of the first importance for the water content of the transition zone, and could also concerns the origin of deep Earthquakes. Based on depth distribution of seismicity, garnet or DHMS dehydration reactions have been invoked to explain earthquakes at intermediate depth in subducting slabs (Estabrook 2004; Omori et al. 2004).

In the context of subduction, the amount of water carried down to the transition zone is unclear, it mainly depends on the degree of hydration of the two peridotite layers that are the subducting lithospheric mantle and the downdragged mantle at the base of the mantle wedge. The maximum concentration is controlled by the phase A (Angel et al. 2001; Komabayashi et al. 2004) and the 10 Angstrom phase (Fumagalli & Poli 2005). For the sake of simplicity, we have chosen a value of 1 wt per cent that is in the upper range of models estimates (Schmidt & Poli 1998; Fumagalli & Poli 2005). This remains below the average saturation of the transition zone minerals (wadsleyite, ringwoodite) of around 3 wt per cent. Crossing the spinel–post-spinel phase change, the hydrated minerals undergo a drastic reduction of their water storage capacity (to $10^{-2}$ wt per cent) and much of the water they contain cannot be accomodated by the new mineral structure (Bolfan-Casanova et al. 2003; Komabayashi et al. 2004). Then, two opposite behaviours, still debated, may occur. Either water remains in non-connected fluid inclusions, or it is exsolved from minerals at grain boundaries, creating an interconnected network. In the former case, water will be driven down into the lower-mantle and eventually homogenized within the whole mantle. In the latter, the water concentrated in the fluid phase may

percolate upward. Indeed, the percolation of fluid in an intergranular network is credible, even at very low concentrations (Mibe et al. 1998). For instance, at low pressure (1 GPa), Minarik & Watson (1995) have shown that in a dunite with a 1 mm grain size, a melt network can remain interconnected for concentrations as low as 0.01 vol per cent.

The goal of this paper is to investigate the mechanism and path of transport of a hydrous fluid phase released at the post-spinel mineral transition. For this purpose, we have developed a two-phase flow model of fluid percolation within a permeable and plastic rocky matrix, in which water precipitation is introduced through a source term. In spite of its simplicity, this model provides insights on the behaviour of the precipitated fluid phase. In particular, it predicts, for a large range of parameters an accumulation of the fluid below the phase change, going with an increase of the fluid overpressure. Hydro-fracturing and dyke formation are the expected consequences of such a case. Thus, a significant part of the discussion is devoted to the probable pathway of the dykes discharging the hydrous fluid accumulated below the phase change. This speculated mode of transport has consequences for the large scale distribution of water in the upper mantle and also for conditions favouring deep earthquakes.

## 2 MODEL

The model consists of a vertical slab crossing the spinel–post-spinel phase transition with a vertical velocity $v$. The phase change is assumed to occur progressively over a thickness $l$, with a constant rate equal to $v/l$. As the ringwoodite transforms into perovskite, water is exsolved, so that within the phase change region three phases are coexisting: ringwoodite, perovskite (plus magnesiowustite) and hydrous fluid. For simplicity, the four solid phases will be considered as a single matrix. Two opposing processes are at work in this system. On the one hand, the fluid and the matrix are advected down into the lower mantle by the descending slab; on the other hand the fluid percolates upward because of its lower density. To describe the evolution of this system, we will first develop a set of equations derived from Bercovici & Ricard (2003), in which the water exsolution will be introduced as a source term. All the notations used are summarized in Table 1.

### 2.1 Poro-viscous equations

#### 2.1.1 Mass conservation

In spite of the fact the model will be used to investigate the conditions of hydro-fracturing induced by fluid migration, the matrix and fluid will be considered as incompressible. Indeed, the effects of compressibility are only substantial when the dynamic pressure (the part resulting from the fluid/matrix movements and related to the fluid buoyancy $\delta\rho g$) is of the same order of magnitude as the fluid (or matrix) bulk modulus $\kappa$. A measure of this is given by $\delta\rho gL/\kappa$, called the Deborah number, where $L$ is a characteristic length of the fluid accumulation, (see eq. 21 for a definition of $L$). The bulk modulus of rocks or fluid at that depth are around $10^{11}$ Pa, leading to Deborah number less that $10^{-2}$ and to negligible effects on the compaction process (Vasilyev et al. 1998). Under these assumptions (incompressibility and consequently constant density of both phases), the equations of mass conservation gives:

$$ \frac{\partial \phi}{\partial t} + \nabla \cdot [\phi \mathbf{v}_f] = \frac{\tilde{\Gamma}}{\rho_f} \qquad (1) $$
<table><caption>**Table 1.** Notations used in the paper.</caption><thead><tr><th>Symbols</th><th>Meaning</th><th>Range of value</th><th>Unit</th></tr></thead><tbody><tr><td>φ</td><td>Porosity</td><td></td><td></td></tr><tr><td>φ<sub>0</sub></td><td>Reference porosity</td><td>0.01</td><td></td></tr><tr><td>v</td><td>Slab velocity</td><td>1</td><td>cm yr<sup>-1</sup></td></tr><tr><td>v<sub>f</sub>, v<sub>f</sub></td><td>Fluid velocity</td><td></td><td>m s<sup>-1</sup></td></tr><tr><td>v<sub>m</sub>, v<sub>m</sub></td><td>Matrix velocity</td><td></td><td>m s<sup>-1</sup></td></tr><tr><td>P<sub>f</sub></td><td>Fluid pressure</td><td></td><td>Pa</td></tr><tr><td>P<sub>m</sub></td><td>Matrix pressure</td><td></td><td>Pa</td></tr><tr><td>s, s</td><td>Separation velocity</td><td></td><td></td></tr><tr><td>c, c</td><td>Mean flow velocity</td><td></td><td></td></tr><tr><td>δP</td><td>Effective pressure</td><td></td><td></td></tr><tr><td>ρ<sub>f</sub></td><td>Fluid density</td><td>2 × 10<sup>3</sup></td><td>kg m<sup>-3</sup></td></tr><tr><td>ρ<sub>m</sub></td><td>Matrix density</td><td>4 × 10<sup>3</sup></td><td>kg m<sup>-3</sup></td></tr><tr><td>δρ</td><td>Density contrast</td><td>2 × 10<sup>3</sup></td><td>kg m<sup>-3</sup></td></tr><tr><td>η<sub>f</sub></td><td>Fluid viscosity</td><td>0.1–10</td><td>Pa s</td></tr><tr><td>η<sub>m</sub></td><td>Matrix viscosity</td><td>10<sup>18</sup>–10<sup>24</sup></td><td>Pa s</td></tr><tr><td>g</td><td>The gravity</td><td>10</td><td>m s<sup>-2</sup></td></tr><tr><td>Γ</td><td>Mass of water exsolved</td><td></td><td>kg m<sup>-3</sup> s</td></tr><tr><td>k</td><td>Permeability</td><td></td><td>m<sup>2</sup></td></tr><tr><td>k(φ<sub>0</sub>)</td><td>Reference permeability</td><td></td><td>m<sup>2</sup></td></tr><tr><td>L</td><td>Compaction length</td><td></td><td>m</td></tr><tr><td>l</td><td>Phase change thickness and length scale</td><td>10<sup>4</sup></td><td>m</td></tr><tr><td>l/v</td><td>Timescale</td><td>10<sup>6</sup></td><td>yr</td></tr><tr><td>η<sub>m</sub> v/l</td><td>Pressure scale</td><td></td><td></td></tr></tbody></table>

and

$$ -\frac{\partial \phi}{\partial t} + \nabla \cdot [(1 - \phi)v_m] = -\frac{\tilde{\Gamma}}{\rho_m}, \qquad (2) $$

where $\phi$ is the porosity, that is, the proportion of the volume occupied by the fluid, $\rho_i$ and $v_i$ the density and the velocity of phase $i$, and subscripts $f$ and $m$ indicate the fluid and the matrix, respectively. $\tilde{\Gamma}$ represents the amount of fluid exsolved during the phase change:

$$ \tilde{\Gamma} = \phi_0 \rho_f \frac{v}{l} \tilde{\Gamma}, \qquad (3) $$

in which $\phi_0$ is the volume of exsolved fluid at phase change completion. The hydration state of the lithospheric part of the slab at the 660 km discontinuity is difficult to assess. Here, we assume that the dehydration and water transport processes undergone by the slab have homogenized its water content. Thus, $\tilde{\Gamma}$ is constant with time and can be represented as a box car function centred around the depth at which the phase change occurs. This depth is taken as the origin of the vertical axis:

$$ \Gamma(z) = \begin{cases} 1 & \text{for } |z| < l/2, \\ 0 & \text{otherwise.} \end{cases} \qquad (4) $$

## 2.1.2 Momentum conservation

In a geological context, where the fluid is assumed to be much less viscous than the matrix, Bercovici & Ricard (2003) proposed the following equations for the momentum conservation in the fluid and the matrix:

$$ 0 = \phi \rho_f \mathbf{g} - \phi \nabla P_f - \eta_f \frac{\phi^2}{k(\phi)} (\mathbf{v}_f - \mathbf{v}_m) \qquad (5) $$

and

$$ 0 = (1 - \phi)\rho_m \mathbf{g} - (1 - \phi)\nabla P_m + \nabla \cdot [(1 - \phi)\bar{\bar{\tau}}_m] \\ + \eta_f \frac{\phi^2}{k(\phi)} (\mathbf{v}_f - \mathbf{v}_m) + (P_m - P_f)\nabla\phi, \qquad (6) $$

where $P$ is the pressure, $\bar{\bar{\tau}}_m$ the viscous stress tensor of the matrix:

$$ \bar{\bar{\tau}}_m = \eta_m \left[ \nabla \mathbf{v}_m + \left[ \nabla \mathbf{v}_m \right]^t - \frac{2}{3} (\nabla \cdot \mathbf{v}_m) \bar{\bar{\mathbf{I}}} \right], \qquad (7) $$

where $\eta$ the viscosity and $k(\phi)$ the permeability. These equations are not symmetric because the interfacial forces and the viscous dissipation are neglected in the fluid equation (eq. 5). Apart from the classical Navier–Stokes terms holding for the body forces, the pressure gradient and the viscous stresses, the two phases flow description introduces a Darcy term, proportional to the fluid–matrix velocity difference, and a term proportional to the matrix–fluid pressure difference that accounts for the viscous dissipation due to the interface deformation. In the equations of Bercovici & Ricard (2003), there is also a term accounting for the surface tension forces that has been neglected here.

## 2.1.3 Porosity equation

To close the system, an additional equation is required. On the basis of energy conservation and entropy considerations, Bercovici & Ricard (2003) proposed to relate the variation of the porosity to the pressure difference between the matrix and the fluid:

$$ \frac{D\phi}{Dt} = -\phi(1 - \phi)K \frac{P_m - P_f}{\eta_m + \eta_f}, \qquad (8) $$

where $K$ is a dimensionless factor accounting for pore or grain geometry and is typically $O(1)$ (see Bercovici & Ricard 2003, for discussion). Consequently $K$ disappears in the following equations.

$$ \frac{D}{Dt} = \frac{\partial}{\partial t} + \mathbf{v}_m \cdot \nabla, \qquad (9) $$

is the material derivative. In this equation, we introduce a source term accounting for porosity creation resulting from water exsolution. This volume that occupied by the matrix that disappears $\tilde{\Gamma}/\rho_m$ or that occupied by the fluid that is exsolved $\tilde{\Gamma}/\rho_f$. Both correspond to the mass exchange $\tilde{\Gamma}$. Only the smaller of the two actually accounts for positive entropy source and thus should be included in the porosity equation (Sramek et al. 2006). This gives:

$$ \frac{D\phi}{Dt} = -\phi(1 - \phi) \frac{P_m - P_f}{\eta_m} + \frac{\tilde{\Gamma}}{\rho_m}, \qquad (10) $$

where the fluid viscosity, several orders of magnitude lower than that of the matrix, has been neglected.

## 2.1.4 Separation velocity and effective pressure

The introduction of the effective pressure $\delta P$, the separation velocity $\mathbf{s}$ and the mean flow $\mathbf{c}$:

$$ \begin{aligned} \delta P &= P_m - P_f \\ \mathbf{s} &= \phi(\mathbf{v}_f - \mathbf{v}_m) \\ \mathbf{c} &= \phi \mathbf{v}_f + (1 - \phi) \mathbf{v}_m, \end{aligned} \qquad (11) $$

allows us to rewrite this set of three scalar and two vector equations in a more simple form. The mass conservation equations (eqs 1 and 2) may be rearranged, so that their sum gives an equation for the mean flow $\mathbf{c}$:

$$ \nabla \cdot \mathbf{c} = \frac{\delta \rho}{\rho_m \rho_f} \tilde{\Gamma} \qquad (12) $$

and the combination $(1 - \phi) \times (1 - \phi) \times (2) - (10)$, an expression for $\mathbf{s}$ as a function of $\delta P$:

$$ \nabla \cdot \mathbf{s} = \phi \frac{\delta P}{\eta_m} + \frac{\delta \rho}{\rho_m \rho_f} \tilde{\Gamma}. \qquad (13) $$
The force balance is highlighted by combination of momentum equations, $\phi \times (6) - (1 - \phi) \times (5)$:

$$0 = \frac{\eta_f}{k(\phi)} \mathbf{s} - \nabla [(1-\phi)\delta P] + (1-\phi)\delta\rho\mathbf{g} + \nabla \cdot [(1-\phi)\bar{\tau}_m], \quad (14)$$

whereas their simple sum leads to the hydrostatic equilibrium:

$$0 = -\nabla[(1-\phi)P_m + \phi P_f] + [(1-\phi)\rho_m + \phi\rho_f]\mathbf{g} + \nabla \cdot [(1-\phi)\bar{\tau}_m], \quad (15)$$

in which the viscous stress of the matrix is negligible in comparison to the global pressure gradient. This set of equations, (10) and (12–15), are very useful for 1-D models because it can be reduced to a set of three scalar equations, the solving for the mean flow and the hydrostatic pressure becoming unnecessary.

## 2.2 Simplified dimensionless equations

In reality, the geometry of the system (i.e. a slab entering the lower mantle) is mainly 2-D. However, the dynamics of the water after its release within the spinel–post-spinel phase change will result in a competition between the upward fluid buoyancy and the downward advective transport by the slab. As a first step, this can be approximated by a 1-D vertical model.

In 1-D, the viscous stress tensor is proportional to the divergence of the matrix velocity. It can be evaluated by combination of the matrix mass conservation (2) with the porosity eq. (10):

$$(1 - \phi)\bar{\tau}_m = \frac{4}{3}\eta_m(1 - \phi)\frac{\partial v_m}{\partial z} = -\frac{4}{3}\phi(1 - \phi)\delta P. \quad (16)$$

The force balance eq. (14) gives:

$$0 = \frac{\eta_f}{k(\phi)} s - \frac{\partial}{\partial z} [(1-\phi)\delta P] + (1-\phi)\delta\rho g - \frac{4}{3} \frac{\partial}{\partial z} [\phi(1-\phi)\delta P], \quad (17)$$

where the vertical axis is oriented downward and $z$ is the depth. One can see that the first part of the viscous stress tensor (the fourth term of the right hand side of equation) is $\phi$ times the order of magnitude of the gradient pressure term, that is, the second term. Since the water content of mantle rocks is likely not to exceed one percent (Schmidt & Poli 1998), we can assume that the porosity will remain small, so that $\phi \ll 1$. This allows us to neglect this first part of the viscous stress tensor. Subsequently, we will assume that $1 - \phi \simeq 1$ which simplifies eq. (17) to:

$$0 = \frac{\eta_f}{k(\phi)} s - \frac{\partial}{\partial z} \delta P + \delta \rho g, \quad (18)$$

Thus, a differential equation for the separation velocity, $s$, is obtained. However, dimensional analysis reveals that $s$ is simply proportional to the balance between the pressure gradient and the buoyancy forces. The dimensionless form of the equations is obtained scaling spatial quantities with $l$, the time with $l/v$, the pressure with $\eta_m v/l$ and the porosity with $\phi_0$. This gives:

$$\frac{s}{\phi^n} = N_L^2 \frac{\partial}{\partial z} \delta P + N_V, \quad (19)$$

for the momentum equation, where $\phi^n$ is the dimensionless permeability of rocks, whose expression will be discussed in the next section, and where $N_V$ and $N_L$ are two dimensionless numbers.

$N_V$, denotes the ratio of the interstitial velocity of the fluid $V$ to the velocity of the descending slab $v$:

$$N_V = \frac{k(\phi_0)\delta\rho g}{v\phi_0\eta_f} = \frac{V}{v}. \quad (20)$$

It expresses the capacity of the fluid to overcome the slab drag and to accumulate at a fixed depth. The quantity $N_L$, represents the ratio of the characteristic size of a compaction wave $L$ to the thickness of the phase transition $l$:

$$N_L = \frac{1}{l} \sqrt{\frac{k(\phi_0)\eta_m}{\eta_f\phi_0}} = \frac{L}{l}. \quad (21)$$

Its meaning is less obvious than that of $N_V$, and will be discussed in the light of the experiment results. Finally, the momentum equation can be written as a simple Darcy-like equation:

$$s = \phi^n \left[ N_L^2 \frac{\partial}{\partial z} \delta P - N_V \right]. \quad (22)$$

This last equation is used to eliminate $s$ in (13) and thus to obtain an equation for the effective pressure, $\delta P$:

$$\phi \delta P - N_L^2 \frac{\partial}{\partial z} \left[ \phi^n \frac{\partial \delta P}{\partial z} \right] = -N_V \frac{\partial \phi^n}{\partial z} - \frac{\delta \rho}{\rho_m} \Gamma. \quad (23)$$

In this equation, we note the presence of a source term proportional to the matrix–fluid density contrast and to the exsolution rate, a term also present in the mean flow eq. (12). It accounts for effects induced by the volume variation due to the fluid exsolution. As a matter of fact, in the Boussinesq approximation, the density variations are only considered as non zero in the buoyancy term. Thus, keeping this term would be only consistent with a compressible approach, where the elastic forces due to the volume variation related to the phase changes are taken into account. Indeed, we verified that this term induces a dramatic increase of the pressure difference between the fluid and the matrix ($-\delta P$). Thus, to stay in the Boussinesq approximation, we drop this term in both (12) and (23). This lead to a mean velocity $c$ that remains constant in the phase change region, and thus equal to $v$ everywhere.

Finally, the dynamics of the fluid exsolved during the phase change can be described in a 1-D framework by a set of two equations as function of the porosity $\phi$ and the effective pressure $\delta P$:

$$\frac{D\phi}{Dt} = -\phi\delta P + \frac{\rho_f}{\rho_m}\Gamma, \quad (24)$$

$$\phi \delta P - N_L^2 \frac{\partial}{\partial z} \left[ \phi^n \frac{\partial \delta P}{\partial z} \right] = -N_V \frac{\partial \phi^n}{\partial z}. \quad (25)$$

This system has been discretized over 30 000 equidistant grid points and solved with a numerical scheme that is fourth order in space and second order in time (Barcilon & Richter 1986). Boundary conditions are periodic and the initial porosity is equal to zero everywhere. The equations are solved in a Lagrangian system of coordinates that moves with the mean flow $c$, that is, the slab velocity $v$.

## 3 PARAMETERS AND DIMENSIONLESS NUMBERS

The estimation of the magnitude of the dimensionless numbers $N_L$ and $N_V$ is not easy, because they depend on physical properties that remain poorly constrained or simply unknown under the pressure of the spinel–post-spinel phase transition. This is the case for the permeability and the viscosity of the fluid whose existence and, a
fortiori, composition, still remain speculative. These parameters, as well as the lower mantle viscosity will serve as the free parame- ters of the problem, whereas other parameters can be constrained more precisely. Thus, the magnitude of N<sub>V</sub> is principally related to the estimation of the ratio k(φ<sub>0</sub>)/η<sub>f</sub> that controls the velocity of fluid inside the porous matrix, N<sub>L</sub> depending on both k(φ<sub>0</sub>)/η<sub>f</sub> and η<sub>m</sub>.

The dependence of matrix permeability on porosity is commonly
described by a power law (McKenzie 1984): k(φ) = d²/Cφⁿ, in
which d is the rock grain size, C a constant and n an exponent
that depends on the porous network geometry. For connected chan-
nels with cylinder shape, it is proposed that n = 2, but with values
for the constant C varying strongly among different authors. For
porosity of 1 per cent in a milimetric grain rock, the predicted per-
meability ranges from 1.5 × 10⁻¹⁵ m² with C = 64 000 (Maaloe &
Scheie 1982) to 3 × 10⁻¹³ m² with C = 3000 (McKenzie 1989).
Contrasting with these values, Faul (1997) proposed a much lower
estimation of the permeability of rocks at small percentage of fluid.
In laboratory experiments where the fluid fraction remains below
2 per cent, he observed that only 10 per cent of the fluid is con-
nected in tubules network, the most part being in intergranular films
and disc-shaped inclusions. In such cases, he suggests a power-law
exponent at least equal to 3, leading to permeability of ~10⁻¹⁷ m²
for porosity of 1 per cent. Note that uncertainties also come from
probable effects of deformation on permeability (Holtzman *et al.*
2003). With regard to the water concentration in deep minerals
(~1 per cent), we have chosen to set n = 3 with a possible per-
meability range of [10⁻¹⁷, 10⁻¹⁴] m².

The viscosity of the fluid expected at the spinel–post-spinel phase transition should be bounded by those of magmas known at the surface. Andesitic melts correspond to the upper limit, around 1000 Pa s. Intermediate values (~10 Pa s) are found for basaltic melts at subsurface pressure and the lowest ones (<1 Pa s) are those of kimberlitic melt, also known for their high content in water and carbon dioxide. This last argument could support the hypothesis of a low viscosity, however the composition of such a fluid could be silicate-rich (Bureau & Keppler 1999), and therefore, highly viscous.

In summary, the ratio k(φ₀)/ηf could span seven orders of mag-
nitude [10⁻²⁰, 10⁻¹³]. The lower mantle viscosity ηm is no more
constrained, varying from very low values (10¹⁸ Pa s), possible just
in the vicinity of the phase transition where a drastic grain size re-
duction may induce a super-plasticity effect (Solomatov 2001), to
high values (10²⁵ Pa s) expected in a cold slab. In light of such un-
certainties, a very careful discussion of the remaining parameters
would appear unnecessary:

(i) The thickness of the spinel–post-spinel transition zone can be reasonably set around 10 km. This value is a good compromise between the seismological (20–30 km) (Petersen *et al.* 1993) and mineralogical (<4 km) (Ito & Takahashi 1989) estimations.

(ii) Fluid dynamics shows that speed of subducting plate suffers a slowdown by a factor proportional to the square-root of the viscosity contrast when entering the lower mantle. A generally accepted value for this contrast is around 50, implying a decent velocity up to 2 cm yr⁻¹ in the lower mantle. For the sake of simplicity, the value of v will be set to 1 cm yr⁻¹.

(iii) The fluid density is not constrained at all, but we estimate that a density contrast with the matrix around 2000 kg m⁻³ is reasonable. Note that a fluid denser than the matrix δρ < 0, like the one proposed at the core–mantle boundary (Ohtani & Maeda 2001) will sink to the base of the mantle.

(iv) The reference porosity φ₀ is equal to the assumed average
water content of rocks entering the lower mantle, that is, 1 per cent.
Other initial slab concentrations can be easily tested remembering
that N_L (eq. 21) is proportional to φ₀² and N_V (eq. 20) to φ₀

With this set of parameters, the unit of porosity is 1 per cent, of
time 10⁶ yr, of length 10 km and of velocity 1 cm yr⁻¹.

4 NUMERICAL EXPERIMENTS

We present the results for intervals of N<sub>V</sub> and N<sub>L</sub> relevant to fluid migration, that is, [0.01, 100] for N<sub>V</sub> and [0.03, 3] for N<sub>L</sub>. These parameter ranges have been investigated through 221 experiments regularly distributed on a logarithmic scale (13 values for N<sub>V</sub> and 17 for N<sub>L</sub>). These experiments permit us to distinguish between different regimes.

4.1 Fluid flow dominated by the matrix drag

First, we look at experiments where the scale of the interstitial fluid velocity (k(φ₀)δρg/φ₀ηf) is very low compared to the subducting plate velocity v, that is, N<sub>V</sub> ≤ 0.1. In these cases, the fluid is swept along by the matrix, as shown in Fig. 1. During the exsolution (inside the shaded box on the top-left panel), the porosity increases at a constant rate, as would be the case in the absence of gravity. It exceeds the value expected at the phase change completion (φ = 1), to reach 1.155. This results from a residual upward migration of the fluid (v<sub>m</sub> − v<sub>f</sub> ~ .15), which remains too small to overcome the descending matrix drag. Observing the time evolution of the porosity, the velocity difference, and the effective pressure, we see that, at a given depth below the source zone, these quantities reach a maximum and then remain constant so that the region of maximum porosity propagates downward. In fact, this implies that 100 per cent of the fluid is advected downward, towards the deep mantle. Such a situation occurs as long as N<sub>V</sub> < 0.15. It is worth noting that this threshold does not depend on N<sub>L</sub>. It is also evident that similar mode of water transport will be seen if the fluid density is higher than the matrix density (δρ < 0).

4.2 Fluid flow dominated by percolation

For higher N<sub>V</sub> and any N<sub>L</sub> values, the migration of the fluid over-
comes the drag due to the matrix descending flow, and accumulates
at the phase change. In Fig. 2, where N<sub>V</sub> = 1 and N<sub>L</sub> = 10<sup>-1/2</sup>,
the porosity develops in a single stationary wave, anchored inside
the source region. By increasing N<sub>V</sub> or decreasing N<sub>L</sub>, the fluid
accumulation shrinks inside the source region and splits in several
stationary waves as it is shown in Fig. 3.

Fig. 4 is a summary of all the experiments that we ran. It repre-
sents the growth rate of the porosity, dφ/dt, when the porosity has
reached a non-dimensional value of 5, except in cases corresponding
to N<sub>V</sub> = 0.1, where the growth rate tends to zero and the porosity
never exceeds 1.15. In this figure, dashed lines draw a dimensioned
framework, obtained by fixing the value of some parameters. Set-
ting l = 10 km, δρ = 2000 kg m<sup>-3</sup> and v = 1 cm yr<sup>-1</sup>, allows us
to plot the matrix viscosity as a function of N<sub>L</sub> and N<sub>V</sub>: η<sub>m</sub> =
N<sub>L</sub><sup>2</sup>/N<sub>V</sub> δρgl<sup>2</sup>/v. Then, the choice of the reference porosity, φ<sub>0</sub> =
1 per cent, fixes the proportionality coefficient between K(φ)/η<sub>f</sub>
and N<sub>V</sub>.

The shaded area marks the zone where dispersion of porosity
waves occurs, as depicted in Fig. 3. This dispersion is favoured
at high N<sub>V</sub> and small N<sub>L</sub>. As expected, the porosity waves be-
come thinner as N<sub>L</sub> (the compaction length) decreases. This effect
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Evolution of the porosity $\phi$, the matrix-fluid velocity difference and the effective pressure for $N_v = 0.1$ and $N_L = 1$. The shaded area on the left-hand panel corresponds to the region where the fluid is exsolved. It extends from $z = -1/2$ to $1/2$. On the three top panels, the various curves show 5 different stages of the evolution, regularly sampled. The bottom panel shows the evolution of the maximum of the porosity (dashed line) and the effective pressure (plain line).</figcaption></figure>

is shown on the left-hand panel of Fig. 5, where porosity profiles corresponding to cases run at constant $N_v(=10^{3/2})$ and $N_L$ varying from $10^{-3/4}$ to $10^{1/2}$ are plotted. With this $N_v$ value, the width of the porosity wave is roughly equal to one compaction length. We verified that all these profiles superpose perfectly when plotted with a scale length equal to the compaction length, indicating that the width is strictly proportional to the compaction length. Besides, the proportionality coefficient depends on $N_v$. The central panel of Fig. 5 shows that at constant $N_L$, the porosity wave becomes thinner as $N_v$ increases. Furthermore, the width can be smaller than one compaction length, as in case of $N_v = 100$ and $N_L = 0.18$ where the wave size reaches roughly half of the compaction length.

On these various porosity profiles, we observe that the splitting occurs as soon as the wave size reduces to the source extent, or rather that it never occurs outside the source region. In fact, the compaction wave onsets when the fluid encounters an obstruction, that is, a porosity decrease along its migration (Spiegelman 1993). Here, the source term is such an obstruction due to the increase of the porosity with depth. This feature is only present where $\Gamma$ is non zero. Below this region, the porosity decreases with depth because

of the upward migration of the fluid, precluding the development of waves. When the porosity accumulation splits into 2 (or n) waves, the corresponding wavelength is $1/2$ (or $1/n$), that is, smaller than 1. As a consequence, the first term in eq. (25) appears negligible compared to the second one, the diffusion term. Neglecting this first term, the solutions of (25) becomes only dependent on the ratio $N_v/N_L^2$. This is highlighted on the right-hand panel of Fig. 5 showing four porosity profiles in which $N_v/N_L^2$ remains constant from one case to the other. As a consequence, for $N_v$ significantly higher than the threshold that delimits the onset of compaction, the evolution of the porosity seems only dependent on the $N_v/N_L^2$ ratio, that is, on the matrix viscosity. We may notice, in Fig. 4, that, for high $N_v$ values, the growth rate of the porosity wave remains constant along the dash lines corresponding to a constant matrix viscosity.

## 4.3 Compaction versus hydraulic fracturing

From these various experiments, we can discern three situations :

(i) The fluid is swept away in the deep mantle. This occurs in case of highly viscous fluid or very low permeability of the matrix
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Same as Fig. 1 but for $N_V = 1$. and $N_L = 10^{-1/2}$.</figcaption></figure>

(i.e. if $k(\phi_0)/\eta_f < 10^{-17} \text{ m}^2/\text{Pa s}$ or $N_V \le 0.1$) and whatever the viscosity of the matrix. Such a situation will be reinforced in case of fast downwellings.

(ii) Compaction is active and the fluid accumulates below the phase change, but a significant amount escapes and is dragged downward. This occurs for $k(\phi_0)/\eta_f > 10^{-17} \text{ m}^2/\text{Pa s}$ and $\eta_m > 10^{20} \text{ Pa s}$.

(iii) When the matrix viscosity is very low, $\eta_m < 10^{20} \text{ Pa s}$, compaction becomes strongly active and the fluid remains confined within the horizon of the phase transformation where it is produced. This last situation corresponds to the shaded area in Fig. 4, the ($N_L$, $N_V$) values for which the porosity waves split.

Compaction seems to be capable of retaining water in the upper mantle in the third situation only. In the other two, a moderate to significant amount of water is carried into the deep mantle. However, beyond a concentration of a few percent, the fluid may induce various changes in material properties with strong impacts on its migration. For instance, a possible softening of the matrix (Hirth & Kohlstedt 1995) would lead to a drastic enhancement of compaction, pushing the system from the second (ii) to the third regime (iii). Another major change in the evolution of the fluid concentration may come from the increase of the effective pressure with time. Hydraulic fracturing and dyke generation will be triggered as soon

as the difference between the fluid and the solid pressure will exceed the yield strength of rocks. This threshold is estimated to be around 50 MPa (Spera 1987; Rubin 1995). It will be reached at the top of the impregnation where the fluid pressure is always maximum, as shown in Figs 1–3. Generation of cracks would evolves towards the formation of dykes able to drain the fluid back to the upper mantle.

To determine in which case this ending is probable, we will first examine the magnitude of the fluid excess pressure, $P_f - P_m$, as a function of $N_L$ and $N_V$. As shown on the bottom panel of Figs 1 and 3, the maximum fluid excess pressure is an increasing function of time that tends to a limit. However, in most of cases, such as in Fig. 2, this limit is not reached before $\phi = 5$, value at which the experiments were stopped. Fig. 6(a) is a plot of the maximum fluid pressure reached at the end of the experiments. This value has been multiplied by $N_L^2/N_V$, in order to be easily dimensioned by a scale based on the buoyancy of the fluid $\delta\rho g l$, instead of $\eta_m v/l$ used previously. For instance, the 50 MPa threshold corresponds to the 0.25 curve for $l = 10$ km, or to the 1.0 curve for the same density contrast, but with a phase change thickness reduced by a factor of 4 ($l = 2.5$ km). Thus, a pair of values for $l$ and $\delta\rho$ fixes the threshold position in a ($N_L$, $N_V$) diagram. Further, the choice of the downward velocity $v = 1 \text{ cm yr}^{-1}$ and the reference porosity $\phi_0 = 1$ per cent allows us to draw a dimensioned framework as in Fig. 4. Figs 6(b)
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Same as Fig. 1 but for $N_v = 10$ and $N_L = 0.1$.</figcaption></figure>

and (c) represent the yield strength threshold (thick line) and the time necessary to reach it in such a dimensioned framework, where the only change concerns the phase transition thickness, from $l = 10$ km to 2.5 km, respectively.

First, when $l = 10$ km, which is an intermediate value between seismological and mineralogical estimations of the phase change thickness, fracturing happens in most of cases corresponding to situations (i) or (ii). Thus, the fluid either concentrates in several waves within the phase change horizon (the shaded area on the figure), or compaction trigger cracks after a few million years to fractions of million years. A reduction by a factor of 4 of the phase change thickness, that is, $l = 2.5$ km, shifts the threshold curve to $N_L$ values four times greater, leaving a gap between fracturing and splitting for cases corresponding to situation (ii), and thus opens some possibility of water leakage towards the lower mantle. However, whether the spinel–post-spinel transition zone thickness is $l = 10$ km or 2.5 km, fracturing is triggered for a matrix viscosity greater than $10^{20}$ Pa s. As a result, the yield strength threshold corresponds to a viscosity threshold. A comparison of these last two figures shows clearly that the threshold curve (thick line), the time curves and viscosity curves (dashed lines) drawn for $l = 10$ km are similar to $l = 2.5$ km, but shifted towards $N_L$ values exactly four times greater, indicating

that both plots would be the same if plotted as a function of the compaction length instead of $N_L$.

## 4.4 About the shape of impregnations

As the model is 1-D, no inference can be made on the evolution of the fluid concentration in the horizontal direction and its incidence on the pressure predicted here. Wiggins & Spiegelman (1995) showed that an homogeneous mush slowly develops spherical impregnation with a radius of several compaction lengths $L$. However, in case of obstructions to fluid percolation, that is, when the fluid concentration sharply drops upwards, instability growth in the vertical direction becomes much faster than in the horizontal one, so that sills develop instead of spheres. Clearly, in our model, the zero porosity above the phase change constitutes a total obstruction to the upward fluid migration.

Connolly & Podladchikov (1998) studied conditions enhancing lateral instabilities. In case of purely viscous rheology, spherical impregnations elongate in the vertical direction when viscosity increases progressively with depth and flatten in the reverse situation. With more complex rheologies, they found that a yield stress to plastic deformation or/and a transition from plastic to elastic
<figure><img src="image_5.png" /><figcaption>**Figure 4.** *dφ/dt* when *φ* = 5 as a function of *N*<sub>L</sub> and *N*<sub>V</sub>. Isolines have been extrapolated from a compilation of 221 runs. Three different regimes are distinguishable. From and below *N*<sub>V</sub> = 0.1 the fluid is dragged downward by the slab and the porosity never exceeds 1.15. Above this threshold, the fluid accumulates in a stationary wave below the source term. In the shaded area, the fluid accumulation is restricted to the zone where the phase transition occurs and splits in several units. Note that far from the threshold *N*<sub>V</sub> = 0.1, the growth rate of the porosity waves seems only dependent on the matrix viscosity. Stars refer to cases shown in Figs 1–3.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Porosity profiles in various cases. In the left-hand panel, cases are run with the same *N*<sub>V</sub> and the width of the porosity accumulation is proportional to *N*<sub>L</sub>. In the central panel, *N*<sub>L</sub> remains constant, but the wavelength of the porosity decreases as *N*<sub>V</sub> increases. In the right-hand panel, the porosity profiles, for which the ratio *N*<sub>V</sub>/*N*<sub>L</sub><sup>2</sup> is constant, look similar. It is important to notice that the porosity splits in several waves when it develops inside the source region (the shaded area on the three panels).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** (a) The non-dimensional fluid excess pressure using $\delta\rho g l$ as pressure scale (i.e. $(P_f - P_m)N_L^2/N_V$), (b) and (c) time in Myr required to reach the yield strength, that is, $P_f - P_m > 50$ MPa, computed with two different values of the phase change thickness: (b) $l = 10$ km and (c) $l = 2.5$ km. In the bottom two panels, the thick line marks the position of the yield strength, with cases in which this threshold is never reached on the left side. Alike in Fig. 4, the shaded area corresponds to cases of porosity wave splitting, in which the fluid remains confined in the phase change horizon.</figcaption></figure>
rheology generates tube-like vertical impregnations, separated by a compaction length *L*. These tubes, equivalent to dykes, develop in absence of hydro-fracturing. The authors suggest that this situation applies to the transition between mantle and lithosphere, when the rheology changes from plastic to elastic behaviour, or when there is a stress threshold for the plastic deformation. We find difficult to assume such rheologies for the horizon of spinel-perovskite phase change, where the grain size reduction suggests the existence of a possible superplasticity.

Another limit of the 1-D formulation is that the phase change is assumed to be horizontal. In fact, because of the temperature gradient across the slab, the spinel-perovskite transition should be deeper in the coldest part of the slab, such that the zone of water precipitation may be tilted by few tens of degrees. In case of initially tilted impregnations, Rabinowicz & Ceuleneer (2005) showed that the fluid percolates upslope gathering all the eventual heterogeneities into a single sheet parallel to the slope. Even if it is difficult to transpose this result to the present context, we may suspect that our model underestimates the pressure that can be reached in case of tilted geometry.

# 5 CONJECTURES ON DYKE SWARMS ORIGINATING FROM THE SPINEL-PEROVSKITE BOUNDARY

The issue of the lower mantle water content depends on competition of three processes: the advective transport, driving the hydrous fluid exsolved at the post-spinel phase change down into the lower mantle, the percolation or compaction, returning this fluid up to the transition zone, and the hydraulic fracturing triggered by the increase of fluid excess pressure over the rock yield strength.

## 5.1 Conditions for dyke generation

At small values of $N_L$ and $N_V$, compaction is not able to retain the water swept by the mantle flow. Thus, upon entering the lower mantle, the water is possibly homogeneously distributed by mixing. (Richard *et al.* 2002). Though this situation cannot be excluded, it would occur only in case of very small permeability or large fluid viscosity and large plate velocity.

In most cases, the evolution of the fluid migration depends on the magnitude of the effective pressure, which is mainly related to the matrix viscosity. For a characteristic mantle flow velocity of 1 cm yr⁻¹ and a yield strength of 50 MPa, rock damage is avoided for solid viscosity below $10^{20}$ Pa s. Various dynamic models suggest a viscosity around $10^{21}$ Pa s in the transition zone and around $10^{23}$ Pa s in the upper part of the deep mantle (Kaufmann & Lambeck 2002). Nevertheless, some of these models agree on the possible existence of a low viscosity zone located in between (Ito & Sato 1991). It might result from superplasticity of perovskite (Karato *et al.* 1995), that is, a rheological weakening associated to grain size reduction during the post-spinel phase change. In that case, the fluid accumulate within the phase change horizon and can end up in a complete separation from the solid in a few million years (Rabinowicz *et al.* 2001). This situation appears unstable because of the slope of the phase change inside a slab. A more probable outcome is the growth of pure fluid lenses, at the top of which the fluid pressure will eventually reach the yield strength, inducing hydro-fracturing and emptying through dykes.

In case of matrix viscosity higher than $10^{20}$ Pa s, as expected if no superplasticity affects the rock rheology, fractures most probably

develop in the matrix during dehydration. Then, interstitial fluid concentrates in the early-formed dykes. A dyke made by grouped cracks is likely to drain most of the fluid in an area of thickness *L* (Sleep 1988), decreasing the fluid excess pressure. Thereafter, the fluid pressure rises again and a new cycle starts (Sleep 1995). The recurrence of fracturing is dominated by the time necessary to reach the yield strength. As shown in Figs 6(b) and (c), it ranges from a few millions years to hundred thousands years and could be shorter in case of high permeability and high viscosity of the matrix. This time indicates the maximum depth at which the water may be driven into the lower mantle before to be draining in dykes. This does not exceed a few kilometres below the phase change.

## 5.2 About the dykes' trajectories

To summarize, except in the case where the fluid is driven towards the deep mantle because of its high viscosity, compaction would eventually trigger dykes. Percolation within dykes is the most efficient migration path for a fluid and this process should be preponderant in terms of water distribution. As a consequence, an important issue is the direction of dykes. Usually, when cracks propagate, they follow the direction of the maximum compressive stress ($\sigma_1$) (Rubin 1995). The stress field in a slab crossing the transition zone depends on numerous parameters (temperature, phase transition, rheology, etc.) and is complicated. Nevertheless, deep earthquakes data (Isacks & Molnar 1971) and most of the studies agree upon the presence of a dominant region of downdip compression around the entry of the lower mantle (660 km depth) and a region of downdip tension above (410–500 km). More precisely, Yoshioka *et al.* (1997) locate the tilt of the $\sigma_1$ direction at the wadsleyite–ringwoodite phase transition (~600 km) and in a more realistic geometry (slab dip angle of 60°) Bina *et al.* (2001) predict it around the olivine–wadsleyite transition (~410 km). Based on these models we have drawn a simplified sketchmap of the the direction of dykes in a slab (see Fig. 7). (in other models (Goto *et al.* 1985), the compression region can be narrower and localized at the surface of the slab but the trend remains the same). The rotation of $\sigma_1$ is expected within the phase change horizon. In all our compaction experiments, the maximum fluid pressure is reached well inside this horizon. More precisely, it

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Expected direction of the dykes formed by slab dehydration at 660 km depth. The $\sigma_1$ directions are simplified from Bina *et al.* (2001). Fluid precipitated at the spinel–post-spinel phase change is likely to accumulate into these dykes and to be transported up to transition zone. Eventually, after having settled in the slab neighbourhood and diffused into the nearby mantle the water would be reintroduced into the large scale convective mantle circulation.</figcaption></figure>
is found at mid-depth for cases close to the compaction/fracturing limit (e.g. case of Fig. 2) and at the very top for cases with a short recurrence time. According to this geometry, the hydrated fluid is likely to be drained inside the transition zone following the dykes opened in the cold centre of the slab. First, water will be absorbed by the rocks of the dyke walls, up to saturation. Then, the following dykes will be channelled in this saturated zone and progressively will extend it upwards. The progression of these dykes would stop as soon as $\sigma_1$ becomes horizontal. In the case shown in Fig. 7, this happens at the top of the transition zone (410 km). The two major consequence of such a process is first to return the water into the transition zone and then to build a weak plane that can play an important role in the onset of deep earthquakes.

# 6 CONCLUSION

## 6.1 Implication for generation of deep Earthquakes

It has long been recognized that the occurrence of deep earthquakes is problematic since brittle fracture and frictional sliding on fault surfaces should be prohibited by confining pressure at great depth. As a matter of fact, these earthquakes show negligible isotropic components, contrary to the effect expected from the phase change volumetric variation (Wiens 2001). Various processes have been proposed to account for great depth double couple sources, indicative of shear movement along a plane: lubrication induced by frictional melting (Kanamori et al. 1998), runaway thermal creep along weakness plane caused by the small grain size of a slab interior after the olivine to $\beta$-spinel transformation zone (Karato et al. 2001), or by rapid anticrack growth (Green & Burnley 1989). According to the anticrak theory, earthquakes cannot occur below the spinel–post-spinel phase change because an endothermic transformation is unable to sustain the thermal runaway required for the faulting instability. Dehydration embrittlement is another mechanism that can operate at any pressure, given a suitable source of fluid (Meade & Jeanloz 1991). In fact, dehydration has long been recognized as a probable cause of intermediate depth earthquakes (75–300 km), but has been dismissed for the deepest ones, because beyond 300 km depth the residual water is mainly dissolved in minerals and thus has never been considered as available for faulting mechanism (Green & Houston 1995).

A recent semi-quantitative study describes the stability fields of hydrated minerals down to the upper-lower mantle (30 GPa) (Komabayashi et al. 2004), in which the authors propose the precipitation of an hydrous fluid phase, when hydrogen rich materials enter the lower mantle. If this assumption is correct, the resulting behaviour of the fluid, predicted here, can explain some observations concerning deep earthquakes. The weak planes, necessary to account for the strong double couple sources, naturally result from the migration of fluid released by the dehydration of slabs at the spinel-post-spinel phase transition. The expelled hydrous fluid accumulates in dykes, that converge towards the coldest part of the slab and then propagate upward along it. Downwards, their extension is restricted by the incapacity of the slab to drag along fluid deeper than a few kilometres below the phase change, limiting to this depth the conditions for earthquake triggering. Upwards, following the direction of maximum compressive stress $\sigma_1$, dykes can propagate far away and return water (or hydrous fluid) up to the base of the upper mantle. As shown on Fig. 7, $\sigma_1$ in a subducting slab is vertical between 660 and 410 km, becoming horizontal above. Seismicity of slabs also reveals consistently downdip compressional solutions for

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Number of earthquakes per year as a function of depth after Frohlich (1989). The figure has been redone with events referenced in the Harvard CMT Catalogue, from 1977 January to 2003 December.</figcaption></figure>

earthquakes deeper than 400 km (Isacks & Molnar 1971), whereas it is downdip tensional for shallower seismicity (Green & Houston 1995). In spite of there being no dehydration reaction below 300 km other than the phase change marking the entrance into the lower mantle, the propagation of hydrous dykes offers the conditions for earthquakes triggering in between. Thus, the bimodal aspect of the distribution of earthquakes with depth, shown by Frohlich (1989) (Fig. 8) and the apparent link between dehydration process and earthquake phenomenons reported by Omori et al. (2004) may find here a simple explanation: the seismicity minimum marking the progressive decrease with depth of dehydration reaction and the extend limit of dykes coming from the upper-lower mantle phase transition. It is worth noting that in this frame, the model proposed by Yoshioka et al. (1997) where the $\sigma_1$ direction tilt at 500 km depth gives a perfect explanation of the seismicity minimum.

## 6.2 Implications for large scale mantle water concentration

The geochemical implications of the proposed deep dehydration are not certain. Indeed, the incompatible elements will preferentially leave the matrix for the fluid phase. However, by a depth of 660 km, the slab has already been affected by various dehydration reactions removing much of the water it contained at the beginning of subduction. On the other hand, the solubility of the trace elements depends on the temperature and the chemical nature of the fluid, so that this last dehydration could have a non negligible impact on the formation of geochemical reservoirs.

The role of melting (dehydration) on the geochemical evolution of the mantle has been considered recently in a different context. Bercovici & Karato (2003) propose a water filter acting at the top of the transition zone. Their model, that can explain some features of the geochemical nature of the mantle, partially lies on the current assumption that the water solubility variation reflects the average water distribution in the mantle: they assume a concentration higher
in the transition zone than in the lower and upper mantle. However, the point of the water distribution issue is the existence of a mechanism able to concentrate then to keep the water inside the transition zone. We have shown that diffusion of hydrogen was inefficient in that role (Richard *et al.* 2002). The conclusion of the present study is that most of the water subducted towards the lower mantle will return upwards and only a minor part corresponding to the solubility of hydrogen in deep minerals (around 100 ppm) will actually go into the lower mantle. It is probable that the path of the fluid ends when the direction of the maximum compressive stress becomes horizontal. This probably happens in the transition zone. Giving the 1-D assumption we have made, it is risky to speculate further, but we can say that if the lateral transport of water by percolation is weak, the fluid is somehow trapped in centre part of the slab where dykes formed and diffuses in the nearby region. In the case of a stagnant slab, this type of water rich region could be related to the low P-wave velocity region recently described by Obayashi *et al.* (2006). Water enriched material could be eventually mixed into the whole upper mantle and not preferentially into the transition zone.

It is very difficult to predict all the consequences of the possible precipitation of an hydrous fluid phases at the ringwoodite to perovskite plus magnesiowustite phase transition. Thus, further speculation concerning the global repartition of water in the mantle and its consequences appears difficult, until laboratory experiments have clarify the existence, and composition of the precipitated fluid and the chemistry of the dehydration reaction.

## ACKNOWLEDGMENTS

This paper benefited from thorough and helpful reviews by Haro Schmeling, Lars Ruepke and an anonymous reviewer. Many thanks also to O. Sramek and S. Hier-Majumder for fruitful discussions and G. Leahy for his comments.

## REFERENCES

* Ahrens, T.J., 1989. Water storage in the mantle, *Nature*, **342**, 122-123.

* Angel, R.J., Frost, D.J., Ross, N.L. & Hemley, R., 2001. Stabilities and equations of state of dense hydrous magnesium silicates, *Phys. Earth planet. Int.*, **127**, 181-196.

* Barcilon, V. & Richter, F.M., 1986. Non-linear waves in compacting media, *J. Fluid. Mech.*, **164**, 429-448.

* Bell, D.R. & Rossman, G.R., 1992. Water in the Earth's mantle: the role of Nominally Anhydrous Minerals, *Science*, **255**, 1391-1397.

* Bercovici, D. & Karato, S.-I., 2003. Whole-mantle convection and the transition-zone water filter, *Nature*, **425**, 39-43.

* 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**, 581-596.

* Bina, C.R., Stein, S., Marton, F.C. & Van Ark, E.M., 2001. Implications of slab mineralogy for subduction dynamics, *Phys. Earth planet. Int.*, **127**, 51-66.

* Bolfan-Casanova, N., 2005. Water in the Earth's mantle, *Min. Mag.*, **69**(3), 229-257.

* Bolfan-Casanova, N., Keppler, H. & Rubie, D.C., 2003. Water partitioning at 660 km depth and evidence for very low water solubility in magnesium silicate perovskite, *Geophys. Res. Lett.*, **30**, doi:10.1029/2003GL017182.

* Bonatti, E., 1990. Not so hot 'Hot Spots' in the oceanic mantle, *Science*, **250**, 107-111.

* Bureau, H., Keppler, H., 1999. Complete miscibility between silicate melts and hydrous fluids in the upper mantle: experimental evidence and geochemical implications, *Earth planet. Sci. Lett.*, **165**, 187-196.

* Campbell, I.H. & Taylor, S.R., 1983. No water, no granites—no oceans, no continents. *Geophys. Res. Lett.*, **10**, 1061-1064.

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

* Estabrook, C.H., 2004. Seismic constraints on mechanisms of deep earthquake rupture, *J. geophys. Res.*, **109**, doi:10.1029/2003JB002449.

* Faul, U.H., 1997. Permeability of partially molten upper mantle rocks from experiments and percolation theory, *J. geophys. Res.*, **102**, 10299-10311.

* Frohlich, C., 1989. The nature of deep-focus earthquakes, *Annu. Rev. Earth planet. Sci.*, **17**, 227-254.

* Fumagalli, P. & Poli, S., 2005. Experimentally determined phase relations in hydrous peridotites to 6.5 GPa and their consequences on the dynamics of subduction zones, *J. Petrol.*, **46**, 555-578.

* Goto, K., Hamaguchi, H., Suzuki, Z., 1985. Earthquake generating stresses in a descending slab, *Tectonophysics*, **112**, 111-128.

* Green, H. W. II. & Burnley, P.C., 1989. A new, self-organizing, mechanism for deep-focus earthquakes, *Nature*, **341**, 733-737.

* Green, H.W. II. & Houston, H., 1995. The mechanics of deep earthquakes, *Annu. Rev. Earth planet. Sci.*, **23**, 169-213.

* Hirth, G. & Kohlstedt, D.L., 1995. Experimental constraints on the dynamics of the partially molten upper mantle: deformation in the diffusion creep regime, *J. geophys. Res.*, **100**(B2), 1981-2001.

* Holtzman, B.K., Groebner, N.J., Zimmerman, M.E., Ginsberg, G.B. & Kohlstedt, D.L., 2003. Stress-driven melt segregation in partially molten rocks, *Geochem. Geophys. Geosyst.*, **4**, 5.

* Huang, W.C. & Okal, E.A., 1998. Centroid-moment-tensor solutions for deep earthquakes predating the digital era: discussion and inferences, *Phys. Earth planet. Int.*, **106**, 191-218.

* Isacks, B. & Molnar, P., 1971. Distribution of stresses in the descending lithosphere from a global survey of focal-mechanism solutions of mantle earthquakes, *Rev. Geophys. Space Phys.*, **9**(1), 103-174.

* Ito, E. & Sato, H., 1991. Aseismicity in the lower mantle by superplasticity of the descending slab, *Nature*, **351**, 140-141.

* Ito, E. & Takahashi, E., 1989. Postspinel transformation in the system Mg₂SiO₄-Fe₂SiO₄ and some geophysical implications, *J. geophys. Res.*, **94**, 10637-10646.

* Kanamori, H., Anderson, D.L. & Heaton, T.H., 1998. Frictional melting during faulting, *Science*, **279**, 839.

* Karato, S., Riedel, M.R. & Yuen, D.A., 2001. Rheological structure and deformation of subducted slabs, *Phys. Earth planet. Inter*, **127**, 83-108.

* Karato, S.-I. & Jung, H., 1998. Water, partial melting and the origin of the seismic low velocity and high attenuation zone in the upper mantle. *Earth planet. Sci. Lett.*, **157**, 193-207.

* Karato, S.-I., Zhang, S. & Wenk, H.-R., 1995. Superplasticity in Earth's lower mantle: evidence from seismic anisotropy and rock physics, *Science*, **270**, 458-461.

* Kaufmann, G. & Lambeck, K., 2002. Glacial isostatic adjustment and the radial viscosity profile from inverse modeling. *J. geophys. Res.*, **107**(B11), 2280.

* Koga, K., Hauri, E., Hirschmann, M. & Bell, D., 2005. Hydrogen concentration analysis using SIMS and FTIR: comparison and calibration for nominally anhydrous minerals, *Geochem. Geophys. Geosys.*, **4**, 1019, doi:10.1019/2002GC00378.

* Komabayashi, T., Omori, S. & Maruyama, S., 2004. Petrogenetic grid in the system MgO-SiO₂-H₂O up to 30 GPa, 1600 C: applications to hydrous peridotite subducting into the Earth's deep interior. *J. geophys. Res.*, **109**(B03206), doi:10.1029/2003JB002651.

* Komabayashi, T., Hirose, K., Funakoshi, K. & Takafuji, N., 2005. Stability of phase A in antigorite serpentine) composition determined by in situ X-ray pressure observations. *Phys. Earth planet. Int.*, **151**(3-4), 276-289.

* Litasov, K., Ohtani, E., Langenhorst, F., Yurimoto, H., Kubo, T. & Kondo, T., 2003. Water solubility in Mg-perovskites and water storage capacity in the lower mantle, *Earth planet. Sci. Lett.*, **211**, 189-203.

* Luth, R.W., 1995. Is phase a relevant to the Earth's mantle? *Geochim. Cosmochim. Act.*, **59**, 679-682.

* Maaloe, S. & Scheie, A., 1982. The permeability controlled accumulation of primary magma, *Contrib. Min. Petrol.*, **81**, 350-357.
Meade, C. & Jeanloz, R., 1991. Deep focus earthquakes and recycling of water in the Earth's mantle, *Science*, **252**, 68–72.

Minarik, G.M. & Watson, E.B., 1995. Interconnectivity of carbonate melt at low melt fraction, *Earth planet. Sci. Lett.*, **133**, 423–437.

McKenzie, D., 1984. The generation and compaction of partially molten rocks, *J. Petrol.*, **25**, 713–765.

McKenzie, D., 1989. Some remarks on the movement of small melt fraction in the mantle, *Earth planet. Sci. Lett.*, **95**, 53–72.

Mibe, K., Fujii, T., Yasuda, A., 1998. Connectivity of aqueous fluid in the Earth's upper mantle, *Geophys. Res. Lett.*, **25**, 1233–1236.

Obayashi, M., Sugioka, H., Yoshimitsu, J. & Fukao Y., 2006. High temperature anomalies oceanward of subducting slabs at the 410-km discontinuity, *Earth planet. Sci. Lett.*, **243**, 149–158.

Ohtani, E. & Maeda, M., 2001. Density of basaltic melt at high pressure and stability of the melt at the base of the lower mantle, *Earth planet. Sci. Lett.*, **193**, 69–75.

Omori, S., Komabayashi, T. & Maruyama, S., 2004. Dehydration and earthquakes in the subducting slab: empirical link in intermediate and deep seismic zones, *Phys. Earth planet. Int.*, **146**, 297–311.

Petersen, N., Vinnik, L., Kosarev, G., Kind, R., Oreshin, S. & Stammler, K., 1993. Sharpness of the mantle discontinuities, *Geophys. Res. Lett.*, **20**, 859–862.

Rabinowicz, M. & Ceuleneer, G., 2005. The effect of slope isotherms on melt migration in the shallow mantle: a physical and numerical model based on observations in Oman ophiolites, *Earth planet. Sci. Lett.*, **229**, 231–246.

Rabinowicz, M., Genthon, P., Ceuleneer, G. & Hillairet, M., 2001. Compaction in a mantle mush with high melt concentrations and the generation of magma chambers, *Earth planet. Sci. Lett.*, **188**, 313–328.

Richard, G., Monnereau, M. & Ingrin, J., 2002. Is the transition zone an empty water reservoir? Inferences from numerical model of mantle dynamics, *Earth planet. Sci. Lett.*, **205**, 37–51.

Rubin, A.M., 1995. Propagation of magma-filled cracks, *Annu. Rev. Earth planet. Sci.*, **23**, 287–336.

Ruepke, L., Phipps Morgan, J., Hort, M. & Connolly, J., 2004. Serpentine and the subduction zone water cycle, *Earth planet. Sci. Lett.*, **223**, 17–34.

Schmidt, M. & Poli, S., 1998. Experimentally based water budgets for dehydrating slab and consequences for arc magma generation, *Earth planet. Sci. Lett.*, **163**, 361–379.

Sleep, N., 1988. Tapping of melt by veins and dykes, *J. geophys. Res.*, **93**, 10 255–10 272.

Sleep, N., 1995. Ductile creep, compaction, and rate and state dependent friction within major fault zones, *J. geophys. Res.*, **100**, 13 065–13 080.

Solomatov, V.S., 2001. Grain size-dependent viscosity convection and the thermal evolution of the Earth. *Earth planet. Sci. Lett.*, **191**, 203–212.

Spera, F.J., 1987. Dynamics of translithospheric migration of metasomatic fluid and alkaline magma, in *Mantle metasomatism*, pp. 1–19, eds Menzies, M. & Hawkesworth, C.J., Academics Press, London.

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

Sramek, O., Ricard, Y. & Bercovici, D., 2006. Discussion of two-phase equations with melting, *Geophys. J. Int.*, in press.

Thompson, A.B., 1992. Water in the Earth's mantle, *Nature*, **358**, 295–302.

Vasilyev, O.V., Podladchikov, Y.Y. & Yuen, D.A., 1998. Modeling of compaction driven flow in poro-viscoelastic medium using adaptive wavelet collocation method, *Geophys. Res. Lett.*, **25**(17), 3239–3242.

Wiens, D.A., 2001. Seismological constraints on the mechanism of deep earthquakes: temperature dependence of deep earthquake source properties, *Phys. Earth planet. Int.*, **127**, 145–163.

Wiggins, C. & Spiegelman, M., 1995. Magma migration and magmatic solitary waves in 3D, *Geophys. Res. Lett.*, **22**(10), 1289–1292.

Yoshioka, S., Daessler, R. & Yuen, D.A., 1997. Stress fields associated with metastable phase transitions in descending slabs and deep-focus earthquakes, *Phys. Earth planet. Int.*, **104**, 345–361.