Testing the effects of basic numerical implementations of water migration on models of subduction dynamics

Subduction of oceanic lithosphere brings water into the Earth’s upper mantle. Previous numerical studies have shown how slab dehydration and mantle hydration can impact the dynamics of a subduction system by allowing a more vigorous mantle flow and promoting localisation of deformation in the lithosphere and mantle. The depths at which dehydration reactions occur in the hydrated portions of the slab are well constrained in these models by thermodynamic calculations. However, computational models use different numerical schemes to simulate the migration of free water. We aim to show the influence of the numerical scheme of free water migration on the dynamics of the upper mantle and more specifically the mantle wedge. We investigate the following three simple migration schemes with a finite-element model: (1) element-wise vertical migration of free water, occurring independent of the flow of the solid phase; (2) an imposed vertical free water velocity; and (3) a Darcy velocity, where the free water velocity is a function of the pressure gradient caused by the difference in density between water and the surrounding rocks. In addition, the flow of the solid material field also moves the free water in the imposed vertical velocity and Darcy schemes. We first test the influence of the water migration scheme using a simple model that simulates the sinking of a cold, hydrated cylinder into a dry, warm mantle. We find that the free water migration scheme has only a limited impact on the water distribution after 1 Myr in these models. We next investigate slab dehydration and mantle hydration with a thermomechanical subduction model that includes brittle behaviour and viscous waterdependent creep flow laws. Our models demonstrate that the bound water distribution is not greatly influenced by the water migration scheme whereas the free water distribution is. We find that a bound water-dependent creep flow law results in a broader area of hydration in the mantle wedge, which feeds back to the dynamics of the system by the associated weakening. This finding underlines the importance of using dynamic time evolution models to investigate the effects of (de)hydration. We also show that hydrated material can be transported down to the base of the upper mantle at 670 km. Although (de)hydration processes influence subduction dynamics, we find that the exact numerical implementation of free water migration is not important in the basic schemes we investigated. A simple implementation of water migration could be sufficient for a first-order impression of the effects of water for studies that focus on large-scale features of subduction dynamics.

Abstract.Subduction of oceanic lithosphere brings water into the Earth's upper mantle.Previous numerical studies have shown how slab dehydration and mantle hydration can impact the dynamics of a subduction system by allowing a more vigorous mantle flow and promoting localisation of deformation in the lithosphere and mantle.The depths at which dehydration reactions occur in the hydrated portions of the slab are well constrained in these models by thermodynamic calculations.However, computational models use different numerical schemes to simulate the migration of free water.We aim to show the influence of the numerical scheme of free water migration on the dynamics of the upper mantle and more specifically the mantle wedge.We investigate the following three simple migration schemes with a finite-element model: (1) element-wise vertical migration of free water, occurring independent of the flow of the solid phase; (2) an imposed vertical free water velocity; and (3) a Darcy velocity, where the free water velocity is a function of the pressure gradient caused by the difference in density between water and the surrounding rocks.In addition, the flow of the solid material field also moves the free water in the imposed vertical velocity and Darcy schemes.We first test the influence of the water migration scheme using a simple model that simulates the sinking of a cold, hydrated cylinder into a dry, warm mantle.We find that the free water migration scheme has only a limited impact on the water distribution after 1 Myr in these models.We next investigate slab dehydration and mantle hydration with a thermomechanical subduction model that includes brittle behaviour and viscous waterdependent creep flow laws.Our models demonstrate that the bound water distribution is not greatly influenced by the water migration scheme whereas the free water distribution is.
We find that a bound water-dependent creep flow law results in a broader area of hydration in the mantle wedge, which feeds back to the dynamics of the system by the associated weakening.This finding underlines the importance of using dynamic time evolution models to investigate the effects of (de)hydration.We also show that hydrated material can be transported down to the base of the upper mantle at 670 km.Although (de)hydration processes influence subduction dynamics, we find that the exact numerical implementation of free water migration is not important in the basic schemes we investigated.A simple implementation of water migration could be sufficient for a first-order impression of the effects of water for studies that focus on large-scale features of subduction dynamics.

Introduction
Dehydration of subducting lithosphere and the related hydration of the mantle wedge can influence the dynamics of subduction, as water has a weakening effect on viscous and brittle rheologies (e.g.Sibson et al., 1975;Peacock, 1987;Hirschmann, 2006;Connolly, 2005;Gerya et al., 2008).The amount of fluids carried by subducting oceanic lithosphere is debated, but it is thought that water content can reach up to 3 wt% at the surface, decreasing with depth (Rüpke et al., 2004).Water is acquired through near-surface hydration, which is aided by flexure-related extensional fractures, and by hydrothermal activity with circulation of hot water and vapour in the upper section of the oceanic crust (Staudigel, 2003;Rüpke et al., 2004;Faccenda et al., 2008).
Published by Copernicus Publications on behalf of the European Geosciences Union.
The subducting crust carries water in two phases: (1) free fluids that are contained in the porosity of the rock and can percolate along grain boundaries (Stern, 2002;Bercovici and Karato, 2003;Wark et al., 2003;Rüpke et al., 2004;Cheadle et al., 2004) and (2) mineralogically bound fluids in the form of hydroxyl complexes (OH) (Schmidt and Poli, 1998;Hirschmann, 2006).Once a slab starts to subduct, it undergoes dehydration processes due to the increase in pressure and temperature.Most of the water contained in the porosity of the sediments is expelled near the trench through compaction and is not transported into the mantle.This is known as fore arc volatile discharge (Stern, 2002;Rüpke et al., 2004).Hydrated minerals include, among others, amphiboles, chlorite and serpentine.It has been well documented that mineralogically bound water is released when hydrated minerals undergo certain phase transitions (Schmidt and Poli, 1998;Iwamori, 1998;Kerrick and Connolly, 2001;Hacker et al., 2003;Ohtani et al., 2004;Rüpke et al., 2004;Syracus and Abers, 2006;Hirschmann, 2006).At the same time, experimentally determined phase diagrams suggest that mineralogically bound water can be transported to the base of the upper mantle, or perhaps even greater depths (Iwamori, 1998;Schmidt and Poli, 1998;Stern, 2002;Ohtani et al., 2004).
Dehydration processes can influence subduction in multiple ways.For example, the depth at which major dehydration occurs determines the location of volcanic arcs, which are located ca. 110-120 km above the surface of subducted slabs (England et al., 2004;Syracus and Abers, 2006).It is the melting of mantle wedge materials that is thought to lead to arc volcanism.But water released from the subducting slab decreases the pressure and temperature at which melting occurs, thus enhancing mantle wedge melting and causing volcanism.Mantle material that is hydrated by water released from the slab may also form weak, positively buoyant "wet plumes" that rise upwards and efficiently hydrate the mantle wedge (Billen and Gurnis, 2001;Billen, 2008;Gorczyk et al., 2007;Richard and Iwamori, 2010).These fluids can then cause a more vigorous flow in the mantle wedge.Arcay et al. (2005) showed that mantle wedge hydration can result in thermal erosion and softening of the overriding lithosphere.
Once subduction has started, (de)hydration processes may further influence the evolution of subduction by enforcing an asymmetrical geometry of subduction zones, causing subduction to be one-sided (Gerya et al., 2008).Dehydration processes can in addition aid the exhumation of high-and ultrahigh-pressure metamorphic rocks by creating a wide and weak subduction channel through which rocks are exhumed (Gerya et al., 2002).Dehydration of the subducting slab may in turn increase slab strength, but this effect may be overwhelmed by the strong impact of water on the mantle wedge.
The models mentioned above use similar methods to determine the conditions of pressure and temperature at which dehydration processes occur.These are usually based on thermodynamic calculations (de Capitani and Brown, 1987;Holland and Powell, 1998;Powell et al., 1998;Connolly, 2005) or high-pressure experiments (Schmidt and Poli, 1998;Ohtani et al., 2004;Komabayashi et al., 2005;Iwamori, 2007), and the locations of the dehydration fronts during subduction do not greatly vary between models.
Water migration can be described by two-phase flow conservation equations (Spiegelman , 1993a, b;Bercovici et al., 2001).However, these are not routinely used in numerical simulations of subduction zone dynamics at the scale of the upper mantle as it adds a fairly complex set of equations to an already highly non-linear model.Usually simplifications are therefore made to migrate water in large-scale subduction models.
Bound water is advected along the solid flow field.Bound water can in addition be affected by diffusion mechanisms which allow water to migrate through the solid-phase field as a function of the chemical and temperature gradients (Richard et al., 2006).Because we focus on first-order behaviour of free water migration on subduction dynamics, we keep our models relatively simple and neglect the effect of bound water diffusion.Free water can migrate in the interconnected pore space of the solid phase (Stern, 2002;Wark et al., 2003;Cheadle et al., 2004;Rüpke et al., 2004), create its own hydrated channels (Katz et al., 2006), or be absorbed by non-saturated rocks of the mantle wedge (Iwamori, 1998) to be potentially transported with the mantle flow into the lower mantle (Bercovici and Karato, 2003;Iwamori, 2007;Richard and Bercovici, 2009;Fujita and Ogawa, 2013).Free water is also advected by the solid flow field through which it migrates, and this can result in cases where part of the free water migrates up through the mantle wedge while the rest is carried with the solid flow and subducted into the mantle (Cagnioncle et al., 2007).
Numerical studies of hydration of the mantle by slab dehydration have used different simplified numerical approximations for the migration of free water in the mantle: I Free water migrates vertically in the upper mantle and is not coupled to the solid-phase flow in the mantle wedge (Arcay et al., 2005).
II The migration of free water is vertical but coupled to the mantle flow.The effective migration path of water is therefore no longer purely vertical, but can include a horizontal component.This method has been implemented as an imposed vertical velocity added to the velocity of the solid-phase flow (Gorczyk et al., 2007) or as a dehydration front with an imposed horizontal and vertical velocity (Gerya et al., 2002).
III Free water migrates as a Darcy flow, following the density gradient between the solid phase and the fluid phase in the mantle wedge (Cagnioncle et al., 2007).Darcy flow changes the migration paths of fluids which are now no longer necessarily vertical (Cagnioncle et al., 2007).Also here the solid flow phase may advect free water in addition to the Darcy mechanism.
Studies that use the above water migration schemes show differences in the spatial distribution of hydrated material in the mantle wedge as subduction evolves (Arcay et al., 2005;Gerya et al., 2002;Gorczyk et al., 2007;Cagnioncle et al., 2007).However, as the numerical setup of the subduction models also differs between these studies, it is difficult to evaluate the possible effects of the numerical implementation of water migration.So far, the influence of the basic numerical implementation of water migration on the dynamics of a subduction model has not been investigated.
We aim to investigate the effects of the three numerical water migration schemes described above (schemes I, II and III) on the dynamics of a subducting slab and its overlying mantle wedge.These models are kept simple, allowing us to focus on the first-order effects of dehydration and water migration.Therefore our models do not include melting, shear heating or adiabatic heating.We first illustrate the effects of dehydration and water migration for a simple model of a cold and hydrated cylinder sinking in a warm mantle.Our second series of models examines the effects of (de)hydration and water migration on a thermo-mechanical subduction model at the scale of the upper mantle.

Thermo-mechanical equations
We solve the equations for conservation of mass (assuming incompressibility) (Eq.1), momentum (Eq.2) and energy (Eq.3): v is the velocity vector, ρ density, t time, P pressure (mean stress), σ the deviatoric stress tensor, g gravitational acceleration (g x = 0 and g y = −9.81m s −2 ), C p specific heat, T temperature, k thermal conductivity, and H radioactive heat production per unit volume.In the subduction models, the Boussinesq approximation is assumed, i.e. ∂ρ ∂t = 0 but ρ = ρ 0 (1 − α(T − T 0 )), where ρ 0 is the reference density at T = T 0 and α is the volumetric thermal expansion coefficient.
Materials are either viscous or brittle.Our viscous rheologies are linear or pressure-and temperature-dependent:

;
(4) 2 ), A a material constant, n the power law stress exponent, d grain size, p grain size exponent, C OH water content in ppm, r the water content exponent, Q activation energy, V activation volume and R the molar gas constant.df and ds refer to deformation by diffusion creep (p > 0 and n = 1) and dislocation creep (p = 0 and n > 1) respectively.Diffusion and dislocation creep are assumed to act in parallel in all materials, resulting in a composite viscosity (η comp ) (Karato and Li, 1992;van den Berg et al., 1993): In our models, only bound water influences the viscosity and we only consider the impact of water on sub-crustal materials.A water content of 0.4 wt % results in a viscosity decrease by ca. 2 orders of magnitude when using the dislocation or diffusion creep flow law for wet olivine from Hirth and Kohlstedt (2003).This can result in viscosities that are lower than the minimum viscosity of 10 18 Pa s which is imposed in our models.We therefore assume that viscosity no longer decreases further once water content exceeds 0.4 wt %.Brittle behaviour in the subduction models follows a Drüker-Prager criterion (Handin, 1969;Jaeger and Cook, 1976;Twiss and Moores, 1992): σ e is the effective deviatoric shear stress (σ e = ( 1 2 σ ij σ ij ) 1 2 ), φ the angle of internal friction, and C the cohesion.φ undergoes a linear decrease with total effective plastic strain (measured as the square root of the second invariant of the strain tensor) to simulate strain weakening.Such strain weakening is thought to result from a reduction in fault rock grain size (Handy et al., 2007), mineral transformations (White and Knipe, 1978;Tingle et al., 1993) or the development of foliation or high fluid pressures (Hubbert and Rubey, 1959;Sibson, 1977).The effective viscosity for plastic flow is (Lemiale et al., 2008) In our thermo-mechanical subduction models, we use a minimum viscosity value of 10 18 Pa s and a maximum cut-off of 10 24 Pa s.These values ensure efficient convergence of our mechanical solution, while allowing for viscosity contrasts of 6 orders of magnitude.The minimum viscosity of 10 18 Pa s is low enough to capture almost all viscosity variations in our model.The cut-off value is only reached very locally in the mantle wedge.We solve the thermal and mechanical equations with a 2-D version of SULEC, which is an arbitrary Lagrangian-Eulerian (ALE) finite-element code (Buiter and Ellis, 2012).The mesh consists of quadrilateral elements which have linear continuous velocity and constant discontinuous pressure fields.Particles are used to track the Figure 1.Water content as a function of pressure and temperature calculated using Perple_X for (A) bulk oceanic crust, and (B) serpentinised harzburgite lithologies (Chemia et al., 2010).Bulk compositions are given in Table 1.
material field (through a material identifier) and properties such as particle strain, stress and water content.The particles are advected with the solid flow field using a second-order Runge-Kutta scheme.We use harmonic viscosity averaging and arithmetic density averaging schemes from particles to elements (Schmeling et al., 2008).The subduction models have a free surface and we use the surface stabilisation algorithm of Kaus et al. (2010) and Quinquis et al. (2011).

Calculation of water content
Our subduction model includes three lithologies: a bulk oceanic crust (BOC), serpentinised harzburgite (SHB) for the lithospheric mantle (Chemia et al., 2010), and pyrolite for the sub-SHB mantle (Schmidt and Poli, 1998).Water content is determined in wt% as a function of pressure, temperature and bulk composition (i.e. the average chemical composition of each lithology, Table 1).The water contents of BOC and SHB are calculated using Perple_X (Connolly, 2005) by Chemia et al. (2010) (Fig. 1).Perple_X is a thermodynamic code that minimises the Gibbs free energy of a chemical system to determine the stability fields of the phases which compose the mineral assemblages.Once these stability fields are calculated, it is possible to determine the maximum allowed water content of each phase as a function of pressure and temperature.The thermodynamic calculations of Perple_X are valid up to pressures of 7 GPa and temperatures of 1300 • C and therefore do not cover upper-mantle conditions.We base our water contents in the upper mantle on Schmidt and Poli (1998), who experimentally determined water storage capacity for pyrolite up to 8 GPa and 1100 • C. We extrapolate these data to 11 GPa and 1400 • C by linearly continuing the Clapeyron slopes of the stability fields, similar to Arcay et al. (2005) (Fig. 2).Serpentinisation of the mantle wedge can locally increase the water storage capacity up to 7 wt% H 2 O (Iwamori, 1998;Rüpke et al., 2004;Connolly, 2005;Férot and Bolfan-Casanova, 2012).We furthermore assume that the water storage capacity in the sublithospheric mantle (i.e. the pyrolytic material) does not converge to zero as shown by the phase diagrams determined by Schmidt and Poli (1998) because this would not allow the transport of bound water into the upper mantle.We impose the minimum storage capacity of the mantle to be 0.2 wt%, following Bercovici and Karato (2003) (Fig. 2).

Water migration schemes
Mineralogically bound water is advected along the solidphase flow field.In our models, each particle therefore carries not only a material identifier (which determines its material properties) but also the amount of bound and free water it contains.Free water migrates following one of these imposed migration schemes: (I) one element vertically up per time step (also referred to as "elemental"), (II) imposed vertical velocity, or (III) Darcy flow velocity.All water migration schemes follow three steps: (1) determine the maximum water storage capacity of each particle and the amounts of free and bound water; (2) determine the migration path for the free water, if present and; (3) distribute the free water along the migration path.For every particle we calculate the maximum bound water that the particle can contain using a standard bilinear interpolation in pressure and temperature of wt% H 2 O on gridded versions of Figs.capacity, the particle is undersaturated and no free water is produced.If the mineralogically bound water of the particle exceeds the water storage capacity, it is oversaturated in water and dehydration occurs.The amount of free water is the difference between the mineralogically bound water of the particle and the maximum allowed water.The free water migrates through the model following one of the schemes we are investigating.The first migration scheme (scheme I) assumes that free water moves vertically upwards owing to its negative buoyancy, with one element per time step ( t), and is not affected by the solid-phase flow (Arcay et al., 2005).This implies that the water migration velocity is purely vertical and is imposed as the local vertical grid size divided by the time step.A model using a variable grid size would not have a constant free water velocity, and therefore this should be avoided.If free water is present, undersaturated particles in the current element are hydrated first.If free water remains after this first step, it migrates to the element above; there it hydrates the undersaturated particles of that element from the bottom up.If all particles are saturated and free water is still present, it is evenly distributed over all particles of the element, waiting for the next time step for further upward migration.
The second migration scheme (scheme II) imposes the velocity for free water (v f,x and v f,y ) (Gorczyk et al., 2007).This method reduces the grid dependence of the migration scheme, though does not eliminate it completely as the migration path itself is grid-dependent.When a particle is over- saturated, it releases water distributed along the path defined by v f × t (Fig. 3).The horizontal and vertical components of the solid flow phase are added to the respective components of the free water velocity.Migration of free water is therefore no longer necessarily vertical.As in the case of the migration scheme I, the first step is to hydrate undersaturated particles in the current element.If free water is still present after this step, the remaining free water migrates to the next element, saturating the undersaturated particles from the bottom up, and so on.If all particles along the migration path of this time step are saturated (i.e.all particles in elements 1 to 6 of Fig. 3), the remaining water is distributed evenly over all particles of the last element.The motivation for element-wise distribution is that water migration paths are likely irregular and a linear path for free water would be unlikely (Rüpke et al., 2004;Katz et al., 2006).We only show examples with an imposed vertical velocity and v f,x = 0.
The third migration scheme (scheme III) follows a simplification of Darcy flow where the fluid follows the pressure gradient caused by the difference in density between the fluid and the solid it is percolating through (Turcotte and Schubert, 2002): where q is the Darcy velocity, ρ s and ρ f the density of the solid and fluid respectively, g the gravitational acceleration, η f the viscosity of the fluid, and k the permeability.The permeability follows the empirical definition of Wark et al.
where is the volume fraction of fluid and d the grain size (same as in Eq. 4).The fluid velocity is the sum of the Darcy velocity relative to the volume fraction of fluid and the solid velocity: where v f is the fluid velocity and v s the solid-phase velocity.
The water migration and distribution then follow scheme II.We assume, however, that water migration in scheme III is vertical.This assumption (also previously made by (Cagnioncle et al., 2007)) is reasonable, as the horizontal pressure gradient in our models in the mantle wedge is much smaller than the vertical pressure gradient (Fig. 4).
, the volume fraction of fluid, is determined from the initial water content and the grids for pressure, temperature and wt% H 2 O.However, this assumes that all the free water that is present in interconnected channels is used to calculate the fluid velocity (Eqs.8-10).This would result in unnaturally high fluid velocities.We therefore introduce an efficiency factor, ω, that corresponds to the percentage of interconnected channels of the network through which water can migrate.The effective permeability in Eq. 8 then becomes κ e = ωκ.This reduces the effective fluid velocity as water can only migrate through the interconnected network.
The Darcy water velocity is calculated for every particle of the model.The water velocity is not constant throughout the model, and areas with higher water content have higher water migration velocities.
We use the following output values to quantify the influence of the water migration schemes on the water distribution and the dynamic evolution of the models: (1) the water distribution is described by tracking the depths of the top-most, lower-most, and left-and rightmost particles of hydrated material.The top hydrated particle gives insight into the rate of advancement of the hydration front, which is not necessarily the same as the imposed water migration velocity.For example, assuming a high water migration velocity and a low amount of free water, free water will first saturate undersaturated materials, resulting in an effectively lower velocity of the hydration front.(2) The root-mean-square water contents (OH rms ) for slab or cylinder and the mantle provide information on the rate of dehydration of the slab or cylinder and the rate of hydration of the mantle.This gives the average water content over various domains in the model (e.g. in the mantle or in the slab) where OH b is the bound water content, OH f is the free water content, and A is the area of computation.
(3) For sinking cylinder models in which viscosity depends on water content, the bottom-most particle of the subducting cylinder is also tracked.This shows the influence of water content on the dynamic evolution of the model.
3 The effects of (de)hydration on a simple model of a sinking wet cylinder

Simple model of a sinking cylinder
We first investigate the effects of (de)hydration and water migration with a simple model which simulates the subduction of a detached piece of lithosphere by the sinking of a cold, hydrated cylinder into a warm, dry mantle.These experiments are based on simple Stokes flow and use linearly viscous rheologies.We solve for the advection and conduction of temperature in addition to the mechanical flow, but the temperature does not play a role in the mechanical flow as viscosities are linear viscous.The model domain is 300 km × 300 km and has a uniform Eulerian resolution of 1 km×1 km elements.The initial particle density is 25 particles per element.In this series of experiments no particle injection or deletion scheme is used.The cold cylinder has a radius of 20 km and is centred on the coordinates x = 150 km, y = −130 km (Fig. 5).A strong viscosity contrast between the cylinder and the mantle avoids deformation of the cylinder (Table 2).The mechanical boundary conditions are free-slip on all sides (i.e. the velocity component parallel to the boundary is free, whereas the velocity component perpendicular to the boundary is zero).A 20 Myr old lithosphere, 58 km thick, overlies the mantle.The initial thermal condition of this oceanic lithosphere is determined from the plate-cooling model (e.g.Turcotte and Schubert, 2002) for a surface temperature of 0 • C, a temperature of 1300 • C at 58 km, and a thermal diffusivity of 10 −6 m 2 s −1 .The surface temperature is held at 0 • C and the bottom temperature (at y = −300 km) at 1360.5 • C throughout model evolution, while the lateral sides are insulated (zero heat flux).A high conductivity (k = 105 W m −1 K −1 ) is used in the mantle to enforce a mantle adiabat of 0.25 • C km −1 (Pysklywec and Beaumont, 2004).The initial temperature of the cylinder is 400 • C, which gradually increases as the cylinder warms up.The mantle is of pyrolytic composition, while the lithosphere and the cylinder are composed of SHB.The material and thermal properties are in Table 2.The initial water content of the hydrated cylinder is imposed at 0.2 wt %.
This model setup is run using the three different migration schemes (schemes I, II and III; see Sect.2.3).First, we use an elemental vertical migration of free water.The water migration velocity is the vertical element size divided by the time step: 20 cm yr −1 .Second, we use migration scheme II, imposing vertical free water velocities of 10, 20 and 60 cm yr −1 .These velocities are in line with free water velocities reached in previous models (Gorczyk et al., 2007;Gerya et al., 2002).As we use a uniform grid resolution, the 20 cm yr −1 model should be similar to the scheme I.However, the schemes are not identical because the vertically imposed velocity scheme assumes that free water is also displaced by the solid flow phase field.Finally, we use migration scheme III, where the vertical water velocity is calculated from the Darcy equation (Eqs.8-10).

Sinking cylinder results
We first examine models in which water content does not influence the linear viscosity.The thermal and mechanical evolution of these models is therefore identical.These models are used to test cylinder dehydration and mantle hydration for the three free water migration schemes.
As only the water migration schemes are changed, the evolution of the cylinders dehydration is identical in all cases.Dehydration of the cylinder occurs due to the increase in temperature of the cylinder from the outer rim inwards (Fig. 6a).Dehydration processes are assumed to be nearly instantaneous, whereas water migration velocities are more than an order of magnitude larger than the flow velocities of the mantle.Therefore, the influence of pressure on dehydration is negligible as the cylinder does not sink as quickly as it dehydrates.The hydrated cylinder is initially undersaturated as the initial water content is 0.2 wt% (but it could contain up to 6.8 wt% at its initial pressure and temperature).The interior of the cylinder is therefore hydrated by water from the lower rim, which undergoes dehydration.This explains the differences observed in cylinder OH rms during the first 1 Myr of model evolution (Fig. 7).Once the centre of the subducting cylinder becomes saturated, dehydration in all models converges (Fig. 7c, f, i).
Hydration of the mantle initiates at the lateral sides of the cylinder, as the interior is being hydrated (Fig. 6b).This results in a horned-shape area of hydrated mantle, which progresses inwards as the cylinder dehydrates.The spatial distribution of bound water is affected by the flow of the mantle.This is visible in the area of hydrated mantle above the cylinder: a minimum in width of the hydrated domain occurs at a depth of ca.125 km (Fig. 6a, last stage).This thinning is accentuated as the cylinder sinks towards the bottom of the model domain.The spatial distribution of free water is not affected by the mantle flow as water migration velocities are much higher than the mantle flow velocities.
The final bound water distribution in the model is independent of the water migration scheme, due to the limited amount of initial water, but water migration schemes do influence the evolution of the water distribution in the mantle (Fig. 7d).Due to the limited amount of water that mantle material can absorb, the mantle is rapidly saturated by water released from the dehydrating cylinder.Therefore the faster the vertical migration of water, the faster the mantle hydrates (Fig. 7d).Due to the much higher water saturation values in the lithosphere, the rate of advancement of the hydration front in the overlying lithosphere is greatly reduced and differences in migration schemes are negligible there.Increasing the efficiency factor in the Darcy flow models (ω, in Sect.2.3) by an order of magnitude can locally increase the water velocity by an order of magnitude.However, this has a limited influence on the distribution of water as the average water velocity in the Darcy model stays close to 10 cm yr −1 (Fig. 7g).
To simulate the effect of water on viscosity in these linear viscous models, a linear decrease in the lithospheric and mantle viscosity of 2 orders of magnitude is used over a water content of 0 to 0.4 wt, %.This is of the same order of magnitude as obtained by including 0.4 wt % in the Hirth and Kohlstedt (2003) dislocation or diffusion creep flow law for wet olivine.Including the influence of water on mantle viscosity does not change the overall water distribution of the models (Fig. 7).It does, however, have a small effect on the mechanical evolution of the model (Fig. 8).The lower viscosity values of the hydrated mantle above the cylinder increase the mantle flow velocity in this region, which results in an increase in the sinking velocity of the cylinder.However, this effect is limited as it corresponds to a difference in velocity of 0.1 cm yr −1 .This is because hydration only occurs directly above the sinking cylinder.The sinking velocity of the cylinder is less sensitive to lowering the viscosity above the cylinder than it would be to changing the mantle viscosity below the cylinder.i) Figure 7. Evolution with time of the depth of the top-most particle of hydrated material (top row), mantle OH rms (middle row) and cylinder OH rms (bottom row) for our models of a wet sinking cylinder.The first column show results for scheme I (water migration velocity of 20 cm yr −1 ), the second for scheme II with different imposed velocities (v f,y = 10, 20 and 60 cm yr −1 ) and the third for scheme III, in which the water migration velocity is calculated using the Darcy law with efficiency factor ω = 0.01, 0.1 or 0.3 (Sect.2.3).Models labelled "visc" have a decrease in viscosity with increasing water content.

Sinking cylinder discussion
Our models of cylinder dehydration and mantle hydration show that different numerical water migration schemes do not result in large differences in the distribution of bound water in the mantle.The sinking cylinder induces a vertical solid flow field in the area of mantle hydration.The free water migration velocities are also vertical.The combined fluid migration velocity is therefore vertical and the difference between the schemes lies mainly in the rate of free water migration.The water migration velocities for our three migration schemes are within the same order of magnitude, resulting in a similar bound and free water distribution.However, changing the magnitude of the water migration veloc-ity does initially increase the hydration rate of the mantle and the hydration front rises at different velocities (Fig. 7d).This is a transient phenomenon as the models converge after ca.2.4 Myr.The low water absorption capability of the mantle material results in rapid effective migration velocities of bound water in the mantle, confirming the results of Arcay et al. (2005).Models that simulate a decrease in mantle viscosity with water content show an increase in mantle velocities and a slightly faster sinking of the cylinder.Our simple models of a sinking wet cylinder in a dry mantle indicate that the three schemes we investigated for numerical implementation of water migration lead to very similar results.Within these schemes, the exact implementation might be secondary to the first-order effect that water could have This figure is for models that follow free water migration scheme I (Fig. 7a-c).
on the system.We will test this in the next section with thermo-mechanical models that also include horizontal mantle (wedge) flow components.
Figure 6 shows a one-element-thick ring of hydrated mantle around the subducting cylinder, on the lateral and bottom sides.This is an artefact of the free water migration scheme.Once a particle dehydrates, undersaturated particles of the element are hydrated first.When the contact between cylinder and mantle lies within that element, mantle particles can absorb the released water, resulting in the hydrated ring around the cylinder.

Subduction model setup
We investigate the effects of slab dehydration, water migration, and mantle wedge hydration using a model of a 70 Myr old oceanic lithosphere subducting under a 40 Myr oceanic lithosphere.The model without water follows Quinquis et al. (2013) and has been tested with a number of different numerical codes.
The oceanic plates are composed of two layers: (1) a 7 and 8 km crustal layer for the overriding and subducting lithospheres, respectively, composed of bulk oceanic crust (BOC) (Chemia et al., 2010), and (2) a 32 km thick serpentinised harzburgite (SHB) layer (Chemia et al., 2010).The rheological and thermal parameters are given in Tables 3 and 4. We assume that the upper 16 km of the oceanic plates is hydrated to a certain degree through fractures in the oceanic crust.The upper kilometre of BOC is fully hydrated, resulting in an initial water content of 2.68 wt% H 2 O. Few faults exceed 1 km depth, resulting in an undersaturation of the remaining BOC.The initial water content of the BOC below 1 km depth is set at 1.5 wt% H 2 O. Finally, between 7 (or 8) and 15 km (or 16 km) depth, the overriding and subducting SHB is also undersaturated at an initial water content of 2 wt% H 2 O.These values follow Rüpke et al. (2004) and Faccenda et al. (2012).
The model domain is 3000 km wide and 670 km deep (Fig. 9) and has the highest horizontal and vertical Eulerian resolution at the trench (1 km per element).The total number of elements is 473 × 269 (horizontal × vertical), and 16 particles per element are used initially.Due to the variable grid resolution, this series of models requires injection and deletion of particles to maintain elemental particle density between 12 and 36.The particle deletion scheme helps keep the code memory requirements reasonable.To maintain the overall water content of the model, particles are injected dry, whereas the bound or free water content of a deleted particle is distributed evenly over all other particles of the element.Subduction is initialised by a "weak seed" located at the interplate boundary.The weak seed simulates a preexisting shear zone in the oceanic lithosphere separating the overriding and subducting plates.It is 14 km thick (in the direction perpendicular to the dip angle), extends to a depth of 82 km and has a 35 • dip angle.The top mechanical boundary is a true free surface (both v x and v y are free), whereas the bottom and left boundaries of the model are free-slip.Balanced material in-and outflow is defined on the right boundary of the model domain (the boundary parallel component is again free).To avoid strong shearing at the transition between in-and outflow, a linear velocity gradient from in-to outflow is defined over a 20 km depth interval.The inflow velocity is 5 cm yr −1 over the thermal thickness of the 70 Myr old lithosphere.The outflow velocity is imposed from a depth of −130 km to the bottom of the model domain at −670 km.
The initial thermal conditions of the 40 Myr (from x = 0 to x = −1500 km) and 70 Myr (from x = 1500 to x = 3000 km) old oceanic lithospheres are determined from the plate-cooling model (Turcotte and Schubert, 2002) for a surface temperature of 0 • C, a mantle temperature of 1300 • C at 82 and 110 km depth respectively, and a thermal diffusivity of 10 −6 m 2 s −1 .The initial step in temperatures at x = 1500 km is rapidly diffused.During model evolution, the surface temperature is held at 0 • C and the bottom temperature (at y = −670 km) at 1440 • C, while the lateral sides are insulated (zero heat flux).A high conductivity (k = 183.33W m −1 K −1 ) is defined for the mantle to enforce the mantle adiabat of 0.25 • C km −1 (Pysklywec and Beaumont, 2004).
The subduction models are run using the three water migration schemes.Because vertical grid size varies (coarsening downwards from 1 to 7 km per element), water migration velocities are no longer constant in the elemental scheme and vary between 10 and 70 cm yr −1 .However, the bulk of the dehydration processes occurs in the mantle wedge, and there our model has a constant mesh resolution and therefore a a Angle of internal friction (φ) softens from first to second value over an effective strain interval of 0.5 to 1.5.
constant water migration velocity of 10 cm yr −1 .All schemes are used to investigate the evolution of models with or without the effect of bound water content on viscosity (Eq.4).The sinking cylinder models showed that an efficiency factor ω of 0.1 in scheme III (Sect.2.3) provides realistic water migration velocities.We therefore also use an efficiency factor ω of 0.1 in scheme III in the subduction models.For scheme II, two water velocities are investigated: v f,y = 5 cm yr −1 and 10 cm yr −1 .

Subduction model results
Subduction is initiated at the weak seed by pushing the 70 Myr oceanic plate inwards, resulting in ca.15 km of advance of the interplate contact until brittle failure helps localise deformation at the trench.The slab subducts at a fairly steep angle and with a more-or-less constant sinking velocity.
The hydrated portion of the slab is limited to the top 16 km, and most of the slab is dry and therefore stiff.As we discuss below, the evolution of this stiff subducting slab is affected by (de)hydration processes, but not significantly.
Main dehydration occurs at two locations in the slab: at ca. 150 km depth, where dehydration occurs at the phase tran-sition of blueschists to eclogite, and 210 km depth, where dehydration of chlorite occurs (Figs. 10 and 11).The mantle wedge is hydrated by these dehydration reactions.Hydrated mantle wedge material is entrained by the downwards flow above the subducting slab.This can bring bound water down to the transition zone.Figure 10 shows that the horizontal width of hydrated mantle above the slab decreases with depth.In the wedge, i.e, above 200 km depth, the mantle can be hydrated up to a distance of 100 km away from the surface of the slab.As the slab deepens, the distance of hydrated mantle from the slab surface decreases, from 50 km at ca. 200 km depth to less than 10 km at the slab tip (Fig. 10).Including the effects of bound water on viscosity increases the width of the hydrated mantle, BOC and SHB regions (Fig. 12e, g, h).The difference is, however, limited as the width of the hydrated regions increases by ca.40 km towards the overriding plate.
Due to inflow of hydrated oceanic lithosphere during model evolution, the OH rms of the subducting slab increases as the model evolves (Fig. 12j).After 3 Myr, the OH rms of the slab decreases slightly.This represents the onset of dehydration and is synchronous to the increase in OH rms of the mantle (Fig. 12i).The numerical water migration schemes cause small differences in the distribution of bound water in the mantle.The differences in mantle OH rms between the three water migration schemes is small initially, but it increases to ca. 1×10 −5 after 6 Myr and then remains constant (Fig. 12i).
The difference corresponds to variations in the lateral distribution of hydrated mantle material (Figs. 10 and 12).Including the effects of water on the viscosity does not change the OH rms of the slab (Fig. 12j), but it does increase the OH rms of hydrated mantle (Fig. 12i).The effects on the distribution of free water are more substantial (Fig. 11).The distribution of free water for the elemental and Darcy schemes is similar, but the free water domain is somewhat broader, at ca. 40 km (Fig. 11) in the Darcy models.This is caused by the horizontal component of the mantle flow that is added to the free water velocity.Larger variations occur in the free water distribution for scheme II and can result in locally large quantities of free water of up to 4 wt% H 2 O (Fig. 11).However, in our models, free water has no effect on rheology and therefore no influence on the dynamics of the system.Introducing the effect of water content on viscosity does not have a strong impact on the large-scale mechanical evolution of the model.This is because the evolution of the subducting slab is controlled by its stiffness.Due to the relatively little amount of water present in the slab (which is initially hydrated up to 16 km depth), dehydration processes will not greatly affect evolution of the already stiff slab.We do find that water weakening of viscosity increases flow in the mantle wedge and slightly reduces the curvature of the slab (Figs. 14 and 15).Corner flow is more pronounced in the models that have water-dependent viscosity, as a large part of the mantle above the slab is hydrated and thus weakened, promoting stronger mantle flow.

Subduction model discussion
Our models show a similar slab dehydration evolution to Rüpke et al. (2004), Arcay et al. (2005) andCagnioncle et al. (2007).This is because a similar method for determining the locations of dehydration reactions are used in these experiments.The depths at which dehydration reactions occur are slightly different between these studies because of differences in the thermal structure of the slab in the models.
We find that the overall dynamics of our subduction model are not strongly influenced by the viscosity decrease of mantle materials due to increase in water content.We suggest that this may be caused by the subducting slab being fairly stiff.Our slab is largely dry and the flow law of Hirth and Kohlstedt (2003) results in an average viscosity of ca. 5 × 10 23 Pa s, which is up to 5 orders of magnitude above the viscosity of the mantle.The evolution of this strong slab is not greatly influenced by a further viscosity increase caused by dehydration processes.Our models do show an increase in the corner flow (Fig. 14) but not on the same scale as that observed in the models of Arcay et al. (2005).The slab in the models of Arcay et al. (2005) dips at a shallower angle, which could focus corner flow and cause a larger effect of water on the flow field.
The numerical implementation of water migration has a significant effect on the distribution of free water in the mantle wedge (Fig. 11).In our models, this effect is not visible in the overall evolution of the model because free water does not affect the rheology of the mantle materials.However, the distribution of free water in the mantle wedge could influence the dynamics of subduction by changing the pore pressure, thereby changing the stress, and the viscosity.Similarly, the pore pressure effect of free water could reduce the plastic yield stress, thereby reducing effective viscosities in the slab and brittle parts of the mantle wedge.During fluid flow, compaction and dilation of the solid matrix may occur, related to pore pressure effects.This could locally change the pressure field and thus in return effect water migration paths, changing the hydration patterns in the solid matrix.
To obtain a first-order assessment of the potential impact of free water on the rheology of our models, we calculated a scheme II water migration model with v f,y = 5 cm yr −1 in which the free water content was added to the bound water content in the calculation of the creep viscosity (Eq.4).As bound water in the model can locally already reach values of 2-3 wt % (Fig. 11), we increased the maximum value that is allowed for C OH in the viscosity calculation from 0.4 to 1 wt %.In the initial stages of subduction, where only small amounts of free water are present, the difference in viscosity between the bound water viscosity model and the bound and free water viscosity model is small.However, as the model evolves and the amount of free water increases, the viscosity field changes (Fig. 16).Low viscosity values are localised in the mantle wedge.The amount of hydrated mantle material that reaches the bottom of the model domain is greatly reduced.This is because mantle flow is localised higher up in the mantle wedge.We emphasise that this model is a prelim-inary result, as the effect of free water (as opposed to bound water) on creep flow laws is not established.
In our simplified models of water migration we did not include melts.Melts are sinks that are thought to absorb most of the excess free water.This would then decrease the potential effects of free water on pore pressures.However, the resulting melt would build up pore pressure instead of a pore pressure increase that would have been caused by the now dissolved free water.Therefore, a potential weakening effect of free water on viscosity could remain.Water also decreases the temperature at which melting can occur, encouraging melting, and because melts have a low viscosity, this could impact subduction dynamics.A logical next step would therefore be to include the effects of melts in our models of subduction with (de)hydration processes.
Including a decrease of viscosity with water content in the models influences the bound water distribution in the mantle close to the surface of the subducting slab (Fig. 14).Due to the limited water absorption capabilities of the mantle and the weak mantle wedge, the flow in the mantle wedge is increased, causing a further increase in the area over which the mantle is hydrated and weak.This suggests that subduction and mantle wedge studies that investigate (de)hydration pro-cesses should preferably include the dynamic effects of water on viscosity and thus mantle flow during model evolution.
We assumed that the mantle material in our models can contain up to 0.2 wt % water at upper-mantle pressures and temperatures (Bercovici and Karato, 2003).We find that hydrated mantle material up to 0.2 wt % is entrained by the flow caused by the subducting slab down to the bottom of the model domain (Figs. 10 and 13).This therefore supports the initial assumptions used in the studies of slab dehydration at the transition zone (i.e. between 670 and 410 km) of Richard et al. (2006Richard et al. ( , 2007)).It also agrees with experimentally determined phase diagrams that suggest water could be present in the transition zone or deeper (Ohtani et al., 2004;Komabayashi et al., 2005).The amount of water reaching the transition zone could be greatly increased by including other chemical reactions in the mantle wedge, such as serpentinisation (Iwamori, 1998).
We show that including the weakening effect of bound water on viscosity increases the amount of water brought down to the bottom of the model domain (Fig. 14).This is because the viscosity reduction causes a stronger corner flow that en-trains more hydrated mantle material in the downward flow above the subducting slab.(de)hydration influences our models but that the exact manner of water migration is not that important.We suggest that studies of especially large-scale subduction dynamics may use a simple implementation of free water migration to capture the first-order effects of (de)hydration, as, for example, an imposed velocity or simple Darcy flow water migration scheme.Elemental water migration (simply moving free water up one element per time step) can be used as long as grid resolution is constant and it is realised that grid size determines the water migration velocity.
We find that the different water migration schemes influence the distribution of free water in subduction models.Elemental water migration results in a localised distribution of free water, while an imposed water migration velocity and Darcy flow result in a broader distribution of free water in the mantle wedge.This effect is caused by the free water being moved by the mantle flow in the latter migration schemes.Free water could affect pore pressure and thus material strength, but melting could decrease the amount of free water, and therefore the effect of free water on a subduction system requires future study.Including the effects of bound water content on viscosity does not strongly impact the overall evolution of the subducting slab as it is controlled by the slab stiffness.We do find that the decrease in mantle viscosity with increasing water content as slab dehydration continues causes a more vigorous corner flow.If we assume that the upper mantle contains an average of 0.2 wt % water, this allows saturated hydrated mantle material to be transported down to the base of the upper mantle, supporting previous assumptions of hydrated material residing in the mantle transition zone.
Figure 16.Viscosity fields in the mantle wedge for subduction models calculated using the imposed free water migration scheme II with v f,y = 5 cm yr −1 .(A) Only the effects of bound water are taken into account in the viscosity calculations (using a maximum cut-off for C OH of 0.4 wt % in Eq. 4).This is the same model as Fig. 10b.(B) The effects of both bound and free water are taken into account in the viscosity calculations (using a maximum cut-off for C OH of 1 wt % in Eq. 4).Joya Tetreault for many numerical discussions.This manuscript benefitted from the constructive reviews by Manuele Faccenda, Tobias Keller and Guillaume Richard.

Figure 3 .
Figure3.Schematic migration of free water and its distribution along a prescribed path for our velocity controlled water migration schemes II and III (see text for further explanation).V f is water velocity.

Figure 4 .
Figure 4.A close-up of the mantle wedge region at 5 Myr in a subduction model that uses free water migration scheme I in which viscosity decreases with bound water content.Shown is the elemental bound water content with contours at every 2 MPa of the horizontal pressure difference, dP dx .Horizontal pressure variations are low in the region with highest bound water content.
Figure 5. (A) Model setup for a cold, hydrated cylinder sinking in a warm, dry mantle.All materials are linear viscous, and their rheological parameters are given in Table 2. (B) Initial geotherm for model in (A).

Figure 6 .
Figure 6.Elemental water content for the Stokes flow models using scheme I (20 cm yr −1 for a grid resolution of 1 km × 1 km) showing (A) bound water and (B) free water.
Scheme II

Figure 8 .
Figure8.Depth of the bottom of the subducting cylinder versus time for linear viscous models with or without including the effects of water on viscosity.This figure is for models that follow free water migration scheme I (Fig.7a-c).

Figure 9 .
Figure 9. Model setup for subduction of a 70 Myr old oceanic plate under a 40 Myr old oceanic plate.The top boundary is free.The bottom and left boundary are free-slip, while the right boundary condition includes material in-and outflow.SHB: serpentinised harzburgite; BOC: bulk oceanic crust.

Figure 10 .
Figure10.Bound water distribution for subduction models after 5, 7.5 and 10 Myr for the following free water migration schemes: (A) elemental water migration scheme I (velocity v f,y increases from 10 cm yr −1 just below the lithosphere to 70 cm yr −1 at the base of the model where the grid is coarsest), (B) imposed vertical velocity v f,y = 5 cm yr −1 (scheme II), (C) imposed vertical velocity v f,y = 10 cm yr −1 and (D) Darcy velocity with efficiency factor ω = 0.1 (scheme III).Bound water moves with the solid-phase flow.Viscosity of all materials changes with water content following the flow law of Eq. 4.

Figure 11 .
Figure 11.Free water distribution for subduction models after 5, 7.5 and 10 Myr for the following water migration schemes: (A) elemental water migration scheme I (velocity v f,y varies between 10 and 70 cm yr −1 ), (B) imposed vertical water migration scheme II with v f,y = 5 cm yr −1 , (C) Imposed vertical water migration scheme II with v f,y = 10 cm yr −1 , and (D) Darcy water migration scheme III with efficiency factor ω = 0.1.Viscosity of all materials changes with water content following flow law of Eq. 4.

Figure 12 .
Figure 12.Evolution of: (A) depth of the topmost particle of hydrated mantle, (B) depth of the lowermost particle of hydrated mantle, (C) depth of the lowermost particle of hydrated BOC, (D) depth of the lowermost particle of hydrated SHB for the subduction models, (E) horizontal coordinate of the leftmost particle of hydrated mantle, (F) horizontal coordinate of the rightmost particle of hydrated mantle, (G) horizontal coordinate of the leftmost particle of hydrated BOC, (H) horizontal coordinate of the leftmost particle of hydrated SHB, (I) mantle OH rms , and (J) slab OH rms ."visc" indicates models that have a decrease in viscosity with increasing water content.

Figure 13 .
Figure 13.Viscosity field for models that follow the Darcy free water migration scheme III with efficiency factor ω = 0.1 after 7.5 Myr for (A) viscosity not dependent on water content and (B) viscosity decrease with water content.

Figure 14 .
Figure 14.Bound water distribution for subduction models using Darcy free water migration scheme III (with efficiency factor ω = 0.1) at 5, 7.5 and 10 Myr.(A) Water content does not change viscosity; (B) water content decreases the viscosity (following flow law Eq.4).

Figure 15 .
Figure 15.Free water distribution for subduction models using the Darcy water migration scheme (with ω = 0.1) at 5, 7.5 and 10 Myr.(A) Water content does not change viscosity; (B) water content decreases the viscosity (following flow law Eq.4).

Table 2 .
Input parameters for the linear viscous sinking cylinder model.For all materials thermal expansivity α = 0 and heat production H = 0.

Table 3 .
Subduction model parameters.The viscosity ranges from a minimum of 10 18 to a maximum of 10 24 Pa s.Heat production H = 0.