Devolatilization of Subducting Slabs, Part II: Volatile Fluxes and Storage

TL;DR

This study uses a reactive transport model to analyze how water and carbon are released and stored during subduction, highlighting the effects of fluid flow, lithological interfaces, and sediment removal on volatile fluxes.

Contribution

It introduces a two-dimensional reactive transport model to quantify volatile fluxes and storage during subduction, considering fluid flow, lithology, and sediment removal effects.

Findings

Up-slab fluid flow peaks at subarc depths.

Hydrated slab mantle enhances decarbonation and carbonation.

Sediment removal increases surface volatile fluxes.

Abstract

Subduction is a crucial part of the long-term water and carbon cycling between Earth's exosphere and interior. However, there is broad disagreement over how much water and carbon is liberated from subducting slabs to the mantle wedge and transported to island-arc volcanoes. In the companion paper Part I, we parameterize the metamorphic reactions involving H$_2$O and CO$_2$ for representative subducting lithologies. On this basis, a two-dimensional reactive transport model is constructed in this Part II. We assess the various controlling factors of CO$_2$ and H$_2$O release from subducting slabs. Model results show that up-slab fluid flow directions produce a flux peak of CO$_2$ and H$_2$O at subarc depths. Moreover, infiltration of H$_2$O-rich fluids sourced from hydrated slab mantle enhances decarbonation or carbonation at lithological interfaces, increases slab surface fluxes, and redistributes CO$_2$ from basalt and gabbro layers to the overlying sedimentary layer. As a result, removal of the cap sediments (by diapirism or off-scraping) leads to elevated slab surface CO$_2$ and H$_2$O fluxes. The modelled subduction efficiency (the percentage of initially subducted volatiles retained until $\sim$200 km deep) of H$_2$O and CO$_2$ is increased by open-system effects due to fractionation within the interior of lithological layers.

Tables

Table 1. Table 1: Notations of Symbols

Symbol	Meaning	Unit
$D$	diffusion coefficient of ${\mathrm{CO}}_2$ or $H_{2}O$	$\sim$${10}^{-8}$ $m/s^2$
$P$	pressure	GPa
$T$	temperature	∘C
$H$	slab thickness	m
$i$	$H_{2}O$ or ${\mathrm{CO}}_2$
$t$	time	s
$x$	coordinate along slab	m
$z$	coordinate normal to slab	m
$\theta$	uniform flow angle
$\varphi$	porosity
${𝒗}_f$	fluid velocity	$m/s$
$v_f$	magnitude of fluid velocity	$m/s$
${𝒗}_s$	solid velocity	$m/s$
$v_s$	subduction rate	$m/s$
${\rho }_f$	fluid phase density	$\mathrm{kg}/m^3$
${\rho }_s$	solid phase density	$\mathrm{kg}/m^3$
$c_f^i$	mass fraction of $i$ in fluid phase
$c_s^i$	mass fraction of $i$ in solid phase
$\overline{c^i}$	bulk content of $i$ per unit rock volume	$\mathrm{kg}/m^3$
$\overline{c_{th}^i}$	bulk mass fraction of $i$
$\Gamma$	reaction rate	$\mathrm{kg}/m^3/s$
${\Gamma }_i$	reaction rate for $i$	$\mathrm{kg}/m^3/s$

Equations

$$\frac{\partial \varphi }{\partial t}+\nabla \cdot ({\mathbf{v}}_f\varphi )=\frac{\Gamma }{{\rho }_f},$$

$$\frac{\partial (1-\varphi )}{\partial t}+\nabla \cdot [{\mathbf{v}}_s(1-\varphi )]=-\frac{\Gamma }{{\rho }_s},$$

$$\frac{\partial ({\rho }_f\varphi c_f^i)}{\partial t}+\nabla \cdot ({\rho }_f{\mathbf{v}}_f\varphi c_f^i)={\Gamma }_i,$$

$$\frac{\partial [{\rho }_s(1-\varphi )c_s^i]}{\partial t}+\nabla \cdot [{\rho }_s{\mathbf{v}}_s(1-\varphi )c_s^i]=-{\Gamma }_i,$$

$$-({\rho }_s-{\rho }_f)\frac{\partial \varphi }{\partial t}+{\rho }_f\nabla \cdot ({\mathbf{v}}_f\varphi )={\rho }_s{\mathbf{v}}_s\cdot \nabla \varphi ,$$

$$\frac{\partial [{\rho }_f\varphi c_f^i+{\rho }_s(1-\varphi )c_s^i]}{\partial t}+{\rho }_f\nabla \cdot [\varphi {\mathbf{v}}_f{c}_f^i]+{\rho }_s{\mathbf{v}}_s\cdot \nabla [(1-\varphi )c_s^i]=0.$$

$$\begin{array}{rl}{\mathbf{v}}_s& =(v_s,0),\\ {\mathbf{v}}_f& =v_f(\cos \theta ,\sin \theta ),\end{array}$$

$$-({\rho }_s-{\rho }_f)\frac{\partial \varphi }{\partial t}+{\rho }_f\left[\cos \theta \frac{\partial (v_f\varphi )}{\partial x}+\sin \theta \frac{\partial (v_f\varphi )}{\partial z}\right]-{\rho }_s{v}_s\frac{\partial \varphi }{\partial x}=0,$$

$$\frac{\partial \overline{c^i}}{\partial t}+{\rho }_f\cos \theta \frac{\partial (v_f\varphi c_f^i)}{\partial x}+{\rho }_f\sin \theta \frac{\partial (v_f\varphi c_f^i)}{\partial z}+{\rho }_s{v}_s\frac{\partial }{\partial x}\left[(1-\varphi )c_s^i\right]=0,$$

$$\overline{c_{th}^i}=\frac{\overline{c^i}}{{\rho }_f\varphi +{\rho }_s(1-\varphi )},$$

$$\left(\varphi ,c_s^{H_{2}O},c_f^{H_{2}O},c_s^{C{O}_2},c_f^{C{O}_2}\right)=\Phi \left(P,T,\overline{c_{th}^{H_{2}O}},\overline{c_{th}^{C{O}_2}}\right).$$

$${\rho }_f\nabla \cdot ({\mathbf{v}}_f\varphi )={\rho }_s{\mathbf{v}}_s\cdot \nabla \varphi .$$

$$F={\int }_{\mathrm{path}\mathrm{start}}^{\mathrm{path}\mathrm{end}}\frac{{\rho }_s}{{\rho }_f}{v}_s\frac{\partial \varphi }{\partial x}\mathrm{d}l.$$

Full text

I Introduction

Earth distinguishes itself from other solar-system planets through its habitability that is maintained by its surface climate. Over geological time, water and carbon modulate the climate through geochemical cycles between Earth’s exosphere and interior [1]. Subduction is a tectonic process that brings altered, near-surface rock into the deep Earth and therefore participates in the long-term geochemical cycles. However, the flux of carbon accompanying subduction into deep Earth is still actively debated. Dasgupta and Hirschmann [2] and Hirschmann [3] argue that subducting slabs don’t experience significant degassing or partial melting and so sequester their carbon into the deep Earth over geological history. On the other hand, Kelemen and Manning [4] contend that most carbon is liberated from slabs and migrates into the subduction-zone mantle lithosphere, rather than being recycled into the deep mantle.

Detailed field and modelling studies also give disparate views on the fate of subducting carbon. Kerrick and Connolly [5, 6, 7] used thermodynamics to construct petrological phase diagrams for representative lithologies in subduction zones (i.e., hydrothermally altered slab mantle, metabasalts, metasediments). They concluded that all representative lithologies bring a significant amount of $\mathrm{CO_{2}}$ into the deep mantle, except that clay-rich slab sediments undergo complete decarbonation at forearc depths along hot subduction geotherms. Subsequently, Gorman et al. [8] evaluated the open-system effects induced by aqueous fluid infiltration on the thermodynamic modelling of subduction-zone dehydration and decarbonation. They concluded that $\mathrm{CO_{2}}$ liberation is still limited, in spite of $\mathrm{H_{2}O}$-rich fluid infiltration. Field studies on subducted carbonates, however, suggest considerable carbon release by carbonate dissolution [9, 10]; Piccoli et al. [11] show that the dissolved carbon can be re-precipitated within the slab or proximal mantle wedge. Further thermodynamic models considering aqueous ionic species also suggest that carbon release from subducting slabs is significant because dissolution of carbon in the form of organic ions can enhance carbon removal from rocks [12]. In particular, Connolly and Galvez [13] show that consideration of non-molecular species can roughly double the carbon solubility in fluids co-existing with sediments subducted along cool geotherms.

Purely thermodynamic models of the fate of subducted carbon are zero-dimensional in that the system is assumed to be closed, with no directional mass transfer. In contrast, the one-dimensional, open-system model by Gorman et al. [8] treated $\mathrm{H_{2}O}$ & $\mathrm{CO_{2}}$ allowing for vertical fluid migration. Two-dimensional geodynamic models of porous fluid migration indicates substantial fluid migration nearly parallel to subducting slabs [14]. This focused flow is caused by the formation of a high-permeability channel in the dewatering layer and a compaction pressure gradient that helps contain liquids in the slab. Given that flow directions within the slab are uncertain and reactive flow is path-dependent, it is important to assess the effect of fluid flow direction on the fluxes of $\mathrm{H_{2}O}$ & $\mathrm{CO_{2}}$ out of subducting slabs.

In the computational treatment of open-system behaviors, an extra challenge exists when incorporating metamorphic reactions involving volatiles into fluid flow modelling. Since fluid movement constantly changes the bulk composition of each subdomain within the modelled slab, the computational cost may be prohibitively high if a traditional phase diagram calculation is applied to each subdomain repeatedly throughout the model evolution. In the companion paper Part I, we have parameterized the coupled dehydration and decarbonation processes for representative subducting lithologies (i.e., sediments, MORB, gabbro, and peridotite). This light-weight thermodynamic module focuses on the behaviors of $\mathrm{H_{2}O}$ & $\mathrm{CO_{2}}$ and can readily capture the fractionation and infiltration effects typical in open systems. Thus it forms the basis of the efficient reactive flow model in the current study. We note here that the liquid volatile phase in the model is a molecular fluid residing in the H2O–CO2 binary, so it excludes other carbon species from consideration. The limitation of this assumption is discussed in section V and in the companion Part I.

In this paper, we provide a two dimensional model of reactive fluid flow in subducting slabs. Since the dynamics of within-slab flow remains highly uncertain [e.g., 15, 14, 16, 17], we prescribe the flow direction in our model and investigate model behavior as a function of this parameter, rather than solving equations for momentum conservation. The model incorporates open-system equilibrium thermodynamics and enables us to assess the factors controlling slab dehydration and decarbonation. In particular, we find that nearly up-slab fluid flow produces a peak in the volatile flux at subarc depths and might thus be relevant for arc magmatism. We find that a sedimentary layer can act as a cap that absorbs $\mathrm{CO_{2}}$ released from underlying slab lithologies, if the sediments are not removed during subduction. Furthermore, slab lithospheric mantle, if extensively serpentinized, can cause significant $\mathrm{H_{2}O}$ & $\mathrm{CO_{2}}$ fluxes at subarc depths in warm subduction settings. In all the cases we explored, 20–90 wt% H2O and 80–100 wt% CO2 in the slabs subduct beyond a mantle depth of $\sim$200 km.

In the following, we start with a description of our model setup in section II, followed by the numerics of model solution in section III. Model results with closed-system behavior are presented in section IV.1, which is included to facilitate comparison of our model with previous work assuming constant bulk compositions. We present results on a reference, open-system model in section IV.2, both for comparison with the closed-system model and to serve as a reference case to investigate parameter sensitivity. The subsequent models presented explore the variability of H2O and CO2 fluxes on slab surface in response to changes in fluid flow direction (section IV.3), slab age (section IV.4), extent of slab mantle serpentinization (section IV.5), and removal of slab surface sediments (section IV.6). For all the open-system models we run, the efficiency of H2O and CO2 subduction into the deep mantle is summarized in section IV.7, which is followed by a discussion (section V) of the model limitations.

II Model Setup

The geometry of the model is illustrated in Figure 1; we focus only on the subducting slab. Cold temperatures in the slab create strong resistance to viscous (de-)compaction [14] and hence the slab experiences negligible isotropic deformation over the model times in this study ($\sim$6 Ma). We therefore treat the slab as a rigid plate in our model. Slab deformation has also been modelled using visco-elastic [17] and visco-elasto-plastic rheologies [15], yielding diverse fluid flow patterns within slab. By treating the slab as a rigid plate, the flow direction can be considered as a free parameter that is to be explored (see below), and we can thus circumvent the uncertainties of slab deformation and coupled fluid flow that are difficult to constrain.

Under the rigid plate assumption, the solid velocities ($\boldsymbol{v}_{s}$) in the slab are uniformly of magnitude equal to the subduction rate and can be prescribed as a model parameter. If a slab age and dip are further chosen, the steady-state temperature ($T$) and pressure ($P$) of the slab can be calculated using canonical thermo-mechanical models [e.g., 18]. Our model assumes such a steady-state $P$–$T$ structure of the slab and extracts it from the geodynamic model by England and Katz [19]. Figure 2 illustrates an example of this for a 10-Ma-old slab with a convergence rate of 5 cm yr*-1* and a dip of 45*∘*. This $P$–$T$ field is employed in the following sections except section IV.4 where different thermal structures, dependent on slab age, are tested. Note that this is a young, warm, and relatively fast-converging oceanic lithosphere. The relatively fast subduction speed implies a fast spreading rate at mid-ocean ridges, which correlates with small extent of upper mantle serpentinization [20]. The chosen H2O content of 1 wt% for the modelled slab mantle is thus consistent with this spreading rate.

Superimposed on this rigid and thermally steady-state slab is the reactive H2O & CO2 transport model. This transport model assumes local equilibrium between the solid and percolating liquid phases. The solid phase in the slab is lithologically layered, as shown in Fig. 1b, and consists of sediments, mid-ocean ridge basalts (MORB), gabbros, and peridotites from top to bottom. In the companion paper Part I, we parameterize the equilibrium partitioning of H2O & CO2 between the liquid and solid for each of the four representative lithologies; this parameterization serves as a thermodynamic module coupled to the fluid flow model detailed below.

Following previous studies [22, 23, 24, 25], mass conservation for liquid and solid phases in porous media are expressed as:

[TABLE]

[TABLE]

where densities of both phases are assumed to be constant, and the meaning of symbols is listed in Table 1. In addition, the conservation of volatile species $\mathrm{H_{2}O}$ & $\mathrm{CO_{2}}$ in both phases are:

[TABLE]

[TABLE]

where $i$ represents either $\mathrm{H_{2}O}$ or $\mathrm{CO_{2}}$. Chemical diffusion is neglected here because the Péclet number ($v_{f}H/D$) is $\sim$250 if slab thickness ($H$) of 8 km, diffusion coefficient ($D$) of \sim$$10^{-8} $\mathrm{m^{2}\,s^{-1}}$, porosity ($\phi$) of \sim$$10^{-2} (e.g., Fig. 4g), and a very conservative flux ($v_{f}\phi$) estimate of $\sim$0.1 m kyr*-1* (e.g., Fig. 4h) are used. Adding equation (1) to (2) and (3) to (4) leads to:

[TABLE]

[TABLE]

Note that, in deriving equations (5) and (6), the zero-compaction-rate assumption ($\nabla\cdot\boldsymbol{v}_{s}=0$) is used.

In the coordinate system where the $x$-axis is parallel to and the $z$-axis is normal to the slab extension, if the angle of flow direction ($\theta$) is uniform across the slab (Fig. 1b), then $\boldsymbol{v}_{s}$ and $\boldsymbol{v}_{f}$ in 2D can be written as:

[TABLE]

where $v_{s}$ is the subduction rate, $v_{f}$ is the magnitude of fluid velocity, and $\theta$ represents flow direction, a model parameter that will be explored in section IV. Substituting equation (7) into (5) and (6) yields:

[TABLE]

[TABLE]

where $\overline{c^{i}}$ is volatile mass per unit rock volume and can be converted to mass fraction via:

[TABLE]

which is used as input for the thermodynamic module from the companion paper Part I (the subscript “$th$” indicates input variables for the thermodynamic module):

[TABLE]

Equations (8)–(11) constitute the governing equations in this study.

The unknowns in the governing equations are the magnitude of fluid velocity ($v_{f}$), porosity ($\phi$), bulk compositions ($\overline{c^{i}}$), and $\mathrm{H_{2}O}$ & $\mathrm{CO_{2}}$ content in both phases ($c^{i}_{f}$, $c^{i}_{s}$)—8 unknowns in total. The thermodynamic equation (11) provides solutions to 5 unknowns, so it can be counted as 5 equations. Including equations (8)–(9), the total number of equations is 8, so the set of governing equations is closed. Note that the closedness of governing equations is premised on the prescription of flow angle ($\theta$) and steady-state $P$–$T$ structure of the slab. To be consistent with the steady-state $P$–$T$ field, we solve the governing equations (8)–(11) until a steady state is reached (after $\sim$6 Ma model time for a 5 cm yr*-1* convergence rate).

The model yields the magnitude distribution of fluid velocity with pre-defined uniform flow direction by mass conservation, without resorting to the momentum conservation equation (i.e., Darcy’s law). With the $P$–$T$ path of a subducting slab pre-determined, thermodynamics of local equilibrium dictates the amount of volatiles liberated from or absorbed to every rock parcel during each timestep of the model. Under steady state, with flow direction ($\theta$) prescribed, mass conservation dictates that the magnitude of fluid flux is simply the integration of volatiles released or absorbed (the term $v_{s}(\partial\phi/\partial x)$ in equation (8)) along flow trajectories. Readers are referred to Appendix A for further details.

We note that avoiding the momentum conservation equation by prescribing the flow direction ($\theta$) has positive and negative consequences. The downside is the lack of emergent, dynamic flows; instead, fluid flows are determined by choice of flow direction and conservation of mass along the selected flow paths. As briefly discussed in the introduction, the upside is that the large uncertainties associated with the diversity of possible flow patterns yielded by various dynamic considerations is avoided. For example, two-phase dynamic model using viscous rheology suggests a high-permeability flow channel at the interface between slab and mantle wedge [14], whereas those using visco-elastic rheology show the development of porosity waves within the slab [17]. On smaller spatial scales (e.g., centimeter to meter), Plümper et al. [16] and Malvoisin et al. [26] show that flow channels and porosity waves can develop as well. Besides, visco-elasto-plastic models suggest that faults formed during slab bending exert a strong control on the fluid flow directions and overall flow pattern within the slab [15]. Indeed, there is various field evidence attesting to the complexity of fluid flows in subducting slabs [e.g., 27, 28, 29, 30, 31, 32, 33, 11, 34]. Designating flow direction ($\theta$) as a free model parameter, however, enables us to approximate the overall trend of fluid flows and explore it by setting $\theta$ to different values. Nevertheless, ignoring the flow dynamics does not revert the current model to the conventional open-system models. For example, the fraction of exsolved volatiles that is moved out and reacted upwards is determined ad hoc [e.g., 8], but this fraction is constrained by mass conservation in the current model such that porosity is maintained at a finite steady level.

In summary, the assumptions regarding the thermo-mechanical aspects of the slab in our model are: (i) rigid plate with a prescribed subduction rate $v_{s}$; (ii) constant densities for liquid ($\rho_{f}$) and solid ($\rho_{s}$) phases; (iii) steady-state $P$–$T$ structure of the slab pre-calculated from canonical thermo-mechanical models. The assumptions in the reactive flow model are: (i) local equilibrium between liquid and solid phases; (ii) no diffusion or dispersion of H2O and CO2; (iii) flow direction $\theta$ is prescribed and uniform across the slab; (iv) the reactive flow is solved until a steady state is reached.

III Numerical Method

III.1 Initial and Boundary Conditions

We undertake timestepping for equations (8)–(11) to achieve the steady-state solution. The initial conditions are: bulk volatile content $\overline{c^{i}_{th}}$ are uniformly set to those of the incoming rocks on the left boundary; porosity ($\phi$) and volatile content in both phases ($c^{i}_{f}$, $c^{i}_{s}$) are determined via the thermodynamic module (eq. (11)); initial fluid speed ($v_{f}$) is zero across the modelled slab domain.

The incoming slab is assumed to be layered in lithology (Fig. 1b) and thus has bulk volatile content: $c^{\mathrm{CO_{2}}}_{th}=3.01$ wt% and $c^{\mathrm{H_{2}O}}_{th}=7.29$ wt% [35] for a sedimentary layer from 0 km to 0.5 km; $c^{\mathrm{CO_{2}}}_{th}=2.95$ wt% and $c^{\mathrm{H_{2}O}}_{th}=2.68$ wt% [7] for a MORB layer from 0.5 km to 2.5 km; $c^{\mathrm{CO_{2}}}_{th}=2.84$ wt% and $c^{\mathrm{H_{2}O}}_{th}=2.58$ wt% for a gabbroic layer from 2.5 km to 6.5 km; and $c^{\mathrm{CO_{2}}}_{th}=0.02$ wt% and $c^{\mathrm{H_{2}O}}_{th}=1$ wt% [36] for a slab-mantle layer from 6.5 km to 8 km. Note that $1$ wt% H2O content in serpentinized peridotite roughly corresponds to a serpentinization degree of $\sim$8% [37]. The H2O and CO2 content in the altered gabbro is adopted from typical metabasalt values [7], and re-normalized using the non-volatile composition for gabbro in Hacker [38]. This very likely represents an upper limit because hydrothermal alteration decreases with depth [4]. In fact, gabbros altered in mid-ocean ridges are estimated to contain 0.2–1.3 wt% H2O [39, 40], but all these estimates exclude the contribution from volatile addition in the outer rise near subduction [41]. If hydrothermal changes at both mid-ocean ridges and outer rises are taken into account, the high gabbro volatile content adopted here could represent an upper limit.

As for boundary conditions, since the governing equation for flow speed ($v_{f}$) is first order in space, only one boundary condition along the flow paths is needed. For flow paths originating from the slab base, no-flux ($v_{f}=0$) boundary condition is prescribed (Fig. 1b). However, for the ghost points above the upper boundary in the numerical grid, Neumann boundary conditions ${\partial v_{f}}/{\partial z}$ and ${\partial\overline{c^{i}_{th}}}/{\partial z}$ are adopted for flow speed and bulk volatile content. The left and right boundary conditions depend on flow direction ($\theta$). When $\theta\leq 90^{\circ}$, the flow paths originate from the left and basal boundaries, so we let $v_{f}=0$ and \overline{c^{i}_{th}}=\overline{c^{i}_{th}}\Big{|}_{\mathrm{incoming}} on the left, and ${\partial v_{f}}/{\partial x}=0$ and ${\partial\overline{c^{i}_{th}}}/{\partial x}=0$ on the right. When $\theta>90^{\circ}$, the flow paths originate from the right and basal boundaries, so we let ${\partial v_{f}}/{\partial x}=0$ and \overline{c^{i}_{th}}=\overline{c^{i}_{th}}\Big{|}_{\mathrm{incoming}} on the left, and $v_{f}=0$ and ${\partial\overline{c^{i}_{th}}}/{\partial x}=0$ on the right.

III.2 Solution Procedure

The nonlinear character of the governing equations (8)–(11) is evident by inspection. For example, equation (10) shows that bulk volatile content depends on fluid flow ($v_{f},\theta$), porosity ($\phi$), and volatile content in each phase ($c^{i}_{f}$, $c^{i}_{s}$), among which porosity and phase volatile content depend further on bulk compositions through thermodynamics expressed in equation (11). To solve this equation set, we employ PETSc [Portable, Extensible Toolkit for Scientific computing, 42], following the procedure illustrated in Katz et al. [43] to write down the residuals of governing equations (8)–(9). Every timestepping solve is handled by the nonlinear solver provided by PETSc, and the special treatment in our solution procedure is that we emplace the thermodynamic module (eq. (11)) within every evaluation of residuals at each Newton iteration.

At the sites where devolatilization onsets, there will be pulses of volatile production owing to the increase of porosity from ideally zero to some finite value determined by thermodynamics. These pulses of volatile production along the envelope of the devolatilizing region eventually contribute to integrated volatile fluxes, and cause small fluctuations in the computed fluxes. To focus on the general trend and overall pattern of volatile fluxes within and atop the slab, the results returned by the solution procedure above are post-processed through a smoothing step. Different length scales of smoothing have been experimented with and we choose 12 numerical grid points (corresponds to 3 km) as the spatial span for smoothing. Details on the smoothing are provided in Appendix A.

IV Results

Before presenting the results for the open-system reference model, we first provide the closed-system results (Fig. 3) with a two-fold goal: firstly to compare our closed-system results with the previous, purely thermodynamic model that assumes a closed system, and secondly for later comparison with the open-system reference model to elucidate the open-system effects.

IV.1 Closed System

The closed system is characterized by the fact that there are no mass fluxes into or out of the rocks of interest, leaving the bulk compositions unaltered. As a result, the bulk $\mathrm{CO_{2}}$ & $\mathrm{H_{2}O}$ content of the incoming rocks stay unchanged as subduction proceeds with a speed $v_{s}$; that is, these bulk values in the slab remain equal to those at the inflow boundary (left side in Fig. 1b). Such a closed-system model corresponds to that typically used in earlier studies, where a representative but fixed bulk composition is used for each rock type to construct thermodynamic phase diagrams over which various subduction geotherms are superimposed [44, 45, 39, 46, 38].

With the $P$–$T$ field of a hot slab (Fig. 2), Figure 3 shows our closed-system results. Left panels are results on CO2, and right panels on H2O. From top to bottom, the rows of panels show CO2 and H2O content in the solid rock phase, liquid volatile phase, and their loss relative to the starting values on the inputting boundary. Of course in a closed system, there should be no volatile loss. However, in closed-system devolatization models it is conventionally assumed that volatiles, once exsolved, are immediately extracted from the bulk system (without altering the bulk composition) [e.g., 44, 45, 39, 46, 38]. In this context, a general trend of CO2 and H2O loss from rocks can be seen with along-slab distance in Figure 3a–b, and more clearly in Figure 3e–f. However, there is a reversal of this trend near the slab distance $\sim$210 km in the top sedimentary layer (e.g., Fig. 3a). It can also be seen from Figure 3c–d that, within each lithological layer, the coexisting liquid phase is gradually more CO2-rich as subduction goes deeper, except where the trend reversal takes place.

Under the model assumption regarding extraction from a closed system, the volatile loss in Figure 3e–f is the change in volatile content of the solid phase from its initial value. Since the amount of liquid phase in equilibrium depends on $P$–$T$ conditions, the non-uniform $P$–$T$ field within the slab yields varying extent of volatile loss in the closed-system model. The trend reversal in Figure 3 can be understood by superimposing the top-layer geotherms over the phase diagram for closed systems [e.g., 46]. The $P$–$T$ conditions in the segment beyond $\sim$210 km in the top layer exceed the $P$–$T$ curve for the onset of devolatilization, but to a lesser extent than the $P$–$T$ conditions in the segment before $\sim$210 km, leading to a smaller degree of devolatlization. Similar trend reversal of closed-system volatile loss can also be observed for the warm subduction geotherm in Kerrick and Connolly [6, figure 3a].

Depths of volcanic arcs to slab surface range from $\sim$70 km to $\sim$170 km [21], which correspond to along-slab distances of $\sim$100 km to $\sim$250 km in Figure 3. Closed-system model predictions, shown in Figure 3a–b and e–f, indicate that slab sediments almost completely lose $\mathrm{CO_{2}}$ & $\mathrm{H_{2}O}$ at forearc depths; the basaltic and gabbroic layers can supply significant $\mathrm{H_{2}O}$ but only limited $\mathrm{CO_{2}}$ at subarc depths; the slab-mantle lithosphere can release almost all of its $\mathrm{H_{2}O}$ at subarc depths. All these results are consistent with the previous studies that assume closed systems and consider hot subduction geotherms [5, 6, 7].

To facilitate later comparison, it is worth noting the implications of the closed-system model for the open-system behaviors in the succeeding sections. Firstly, the basal slab-mantle layer in this model contains negligible CO2 and serves essentially as a water supplier. Because hydrated slab mantle has only a small set of hydrous minerals (i.e., talc, brucite, serpentine, chlorite), H2O loss in this layer is more abrupt and complete relative to in basaltic and sedimentary layers. A sudden onset of the supply of the basal H2O will, in later sections on open-system flow, lead to significant H2O infiltration. Secondly, at the peridotite–gabbro and gabbro–basalt interfaces (6.5 and 2.5 km deep into the slab, respectively), there is a sharp CO2 concentration gradient in the liquid phase (Fig. 3c–d). Fluid ascent in the open systems would inevitably cause H2O-rich and CO2-poor fluid infiltration that enhances decarbonation [8]. Thirdly, the CO2 concentration in the liquid phase decreases from the basaltic to the sedimentary layer, suggesting the potential for carbonation by fluid flow down the gradient of CO2 concentration.

IV.2 Open-System Reference Model

Taking all the parameters for the closed-system model above, the reference model additionally considers fluid flows with a uniform direction $\theta=90^{\circ}$ (Fig. 1b), namely, fluid flow is normal to the direction of slab motion ($\boldsymbol{v}_{s}$). The results of the steady-state solution to this open-system model are presented in Figure 4. The panel layout of Figure 4 is the same as that of Figure 3 except for the two additional panels in the bottom row. These show the porosity (Fig. 4g) and volatile flux (Fig. 4h) distribution within the slab. To elucidate and further highlight the open-system effects, Figure 5 maps the differences in CO2 and H2O loss (or gain) between the closed- and open-system models, that is, volatile loss in the open system minus that in the closed system. We next consider the reference model results in each layer of the slab, from bottom to top.

The basal slab-mantle layer does not experience fluid infiltration because the flows are upward and there is no flux at the bottom boundary. Figure 4c–d shows that the devolatilized liquid phase is almost pure H2O, given that the inputting peridotitic rock contains negligible bulk CO2 (0.02 wt%). The H2O loss is thus abrupt and complete (Fig. 4f), similar to that for the closed-system model.

On top of the base layer, the gabbro layer is fluxed by nearly pure H2O sourced from below (Fig. 4c–d). As a result, the infiltration and fractionation effects discussed in Part I come into play. The very H2O-rich and CO2-poor fluids will infiltrate the base of the gabbro layer, causing enhanced decarbonation and inhibited dehydration. The enhanced CO2 loss is evident in Figure 5a, but the expected reduction in H2O loss is not as clear in Figure 5b. This can be explained by the fractionation effect detailed in Part I. Because increased decarbonation reduces the ratio of CO2 over H2O in the bulk system of solid plus liquid, the onset temperature of devolatilization decreases at the base of the gabbro layer. Given that temperatures within the modelled slab are held fixed, such a decrease in the onset temperature means devolatilization of the bulk system becomes easier, so the resultant increase of devolatilization extent can offset the inhibitive effect on H2O loss induced by the infiltration. The elevated extent of devolatilization is also reflected by the higher porosity at the base of the gabbro layer in Figure 4g.

The infiltration effect is strongest at the peridotite–gabbro interface because the concentration gradients of $\mathrm{CO_{2}}$ and $\mathrm{H_{2}O}$ in the liquid phase are the steepest (Fig. 4c–d). Above the interface, the liquid phase compositions have been adjusted by the buffering reactions drastically occurring at the interface, so the infiltration effect weakens. In consequence, the bulk interior of the gabbro layer is dominated by the fractionation effect. At and after the onset of devolatilization, preferential $\mathrm{H_{2}O}$ loss causes increased difficulty in devolatilization and thus inhibits overall volatile loss in the bulk interior (Fig. 5).

Further up, the H2O-rich and CO2-poor fluids (Fig. 4c–d) derived from the gabbro layer infiltrate the base of the basalt layer, similar to the scenario at the base of gabbro layer. Likewise, relative to the closed system, enhanced $\mathrm{CO_{2}}$ loss and approximately unmodified $\mathrm{H_{2}O}$ loss take place at the layer base (Fig. 5), and the overall $\mathrm{CO_{2}}$ and $\mathrm{H_{2}O}$ loss is inhibited in the layer interior. The porosity elevation immediately above the gabbro–basalt interface is even more discernible in this case (Fig. 4g). Moreover, elevated porosity levels reflect higher extent of devolatilization. According to our parameterization in the companion paper Part I, higher extent of devolatilization corresponds to higher CO2 and lower H2O concentrations in the coexisting liquid phase, which are also evident at $\sim$2.5 km depth in Figure 4c–d. Upward infiltration of these fluids, coupled with the fractionation effect, gives rise to the diminished volatile loss immediately above the region of enhanced decarbonation (blue stripes $\sim$2 km deep in Fig. 5). Note also that in Figure 5 the degree of decarbonation is higher at the basalt base than at the gabbro base, which is consistent with the higher integrated fluid fluxes at the shallow lithological interface (Fig. 4h).

Figure 4c–d shows that the fluid infiltration at the base of the topmost sedimentary layer is down the CO2 and up the H2O concentration gradients. As is inferred from the closed-system model, the infiltration will result in carbonation reactions in the sediments immediately above the lithological contact, which is demonstrated by the negative values of CO2 loss (i.e., CO2 gain by the sediments) in Figure 4e. Figure 5a compares the CO2 loss between the open- and closed-system models; carbon sequestration in the top sedimentary layer is more evident. In Figure 5b, it can be seen that H2O loss in the sedimentary layer diminishes before $\sim$250 km along the slab coordinate. This is due to the fractionation effect discussed in Part I: carbonation increases the ratio of bulk CO2 over H2O, which further raises the onset temperature of devolatilization and in turn inhibits dehydration. As for the enhanced H2O loss after $\sim$250 km in Figure 5b, it actually stems from the reversal of trend in the H2O loss calculated for the closed-system model (e.g., Fig. 3f).

A common feature above the gabbro–basalt and basalt–sediment interfaces in Figure 5a is the formation of paired stripes in red and blue, but the vertically reversed color order for the stripe pair in the sedimentary layer indicates that decarbonation follows carbonation along infiltration. This is because the former stripe pair is caused by infiltration down a $\mathrm{CO_{2}}$ gradient, whereas the latter is by that up the gradient. Another noteworthy feature in Figure 5 is that the starting positions of the infiltration-induced stripes in the basalt and gabbro layers coincide horizontally with the onset of dehydration of the slab upper mantle. This coincidence reflects the potential for fluids from hydrated lithospheric mantle to remobilize slab H2O and CO2. In fact, the enhanced infiltration into the top sedimentary layer causes considerable carbonation that substantially reduces the porosities and fluid fluxes immediately after $\sim$160 km along the slab (Fig. 4g–i).

Close examination of Figure 4h suggests that the volatile fluxes fluctuate laterally before the onset of slab mantle dehydration at $\sim$160 km. This is because the onset of devolatilization at different slab depths occurs at different horizontal positions that are laterally separated from one another. Since the onset of devolatilization corresponds to a pulse of volatile production, after being vertically integrated into volatile fluxes, these pulses cause lateral fluctuation in volatile flux distribution. Appendix A discusses this issue further.

IV.3 Effect of Flow Direction

Different flow trajectories pass through areas with different $P$–$T$ conditions that are associated with different volatile partition coefficients. Hence volatile loss and gain during reactive flow are path-dependent, and we assess the impact of fluid flow direction ($\theta$) on CO2 and H2O fluxes emanating from slab surface. Figure 6 shows the slab surface CO2 (left panels) and H2O (right panels) fluxes with flow angles of $22.5^{\circ}$, $45^{\circ}$, $90^{\circ}$, $135^{\circ}$, and $157.5^{\circ}$. The case of flow angle $\theta=90^{\circ}$ is the open-system reference model already presented in section IV.2. Other than flow angle, all the parameters for model runs in Figure 6 are identical to those in the reference model. The flow angles $22.5^{\circ}$, $45^{\circ}$, $90^{\circ}$, $135^{\circ}$, and $157.5^{\circ}$ are chosen to spread between $0^{\circ}$ and $180^{\circ}$, and thus representative of various flow scenarios; other choices of flow angles don’t affect our conclusion. We also note that the angle of $22.5^{\circ}$ corresponds to downward flow, which could occasionally occur in localized settings [e.g., 15]. The consideration of this flow scenario is primarily to explore as completely as possible the spectrum of within-slab flow patterns.

Figure 6 shows profiles of the slab surface volatile flux (i.e., $[v_{f}\phi c_{f}^{i}]_{z=0}\sin\theta$). All the profiles have a peak near $\sim$80 km along the slab, caused by the devolatilization of the top sedimentary layer. In addition, there is a dip in all the flux profiles around the starting position of slab-mantle devolatilization ($\sim$160 km). This sharp drop of surface flux stems from the boost of carbonation reactions by slab-mantle devolatilization, which reduces the top-layer porosity ($\phi$) and further the surface flux ($v_{f}\phi$). The shift in position of this dip is due to the change of flow direction that advances ($\theta>90^{\circ}$) or retreats ($\theta<90^{\circ}$) the carbonation reaction in the top sedimentary layer. A more striking feature in Figure 6 is a second flux peak that emerges at the average position of global arcs projected onto the slab surface (red triangle). As shown in Figure 6g–j, this peak appears only when flow angle $\theta>90^{\circ}$, that is, for fluid flows that are nearly slab-parallel and up-dip. If it is further assumed that magmas generated by flux melting traverse the mantle wedge vertically without deflection [e.g., 47], Figure 6 implies that up-slab fluid flow is responsible for the magma supply to volcanic arcs.

The emergence of the subarc flux peak in Figure 6 can be explained by Figure 7, which maps the difference in fluxes ($v_{f}\phi$) between models with flow angle $\theta\neq 90^{\circ}$ and the reference model with $\theta=90^{\circ}$. The flux increase in Figure 7 can be understood by inspecting equation (13): fluid fluxes accumulate along flow paths and longer paths tend to cause larger fluxes. Moreover, when $\theta>90^{\circ}$, additional flux increase is caused by flow trajectories tapping into the region of slab-mantle devolatilization (Fig. 7c–d). Conversely, when $\theta<90^{\circ}$, the flux diminution in Figure 7a–b is caused by flow trajectories passing through regions before the onset of slab mantle devolatilization. As flows with $\theta>90^{\circ}$ can transport more volatiles to slab surface, this explains the emergence of flux peak near $\sim$150 km in the flux profiles in Figure 6g–j.

IV.4 Effect of Slab Age

Published closed-system and 1D open-system studies [5, 6, 7, 8, 38, 46] show that high-temperature subduction geotherms promote slab devolatilization and thus increase slab-surface volatile fluxes. We test this hypothesis using our 2D open-system model. Results with flow angle $\theta=90^{\circ}$ and slab ages ranging from 10Ma to 60Ma are shown in Figure 8, where it can be seen that, for old slabs with cold geotherms, devolatilization is deferred to greater depths. For the coldest slab (60Ma old) in this study, there is almost no dehydration or decarbonation in the range of subarc depths. In general, relative to the 10-Ma slab in the reference model, old, cold slabs will release less volatiles around subarc depths, therefore promoting $\mathrm{CO_{2}}$ & $\mathrm{H_{2}O}$ recycling into deep mantle [48].

Flow angles in Figure 8 are all $90^{\circ}$ and, according to our assessment of the directional effect of flows, it can be inferred that up-slab flows ($\theta>90^{\circ}$) will produce flux peaks. However, the positions of these peaks are controlled by the position of the onset of slab-mantle devolatlization. Given that the onset of slab-mantle devolatlization is also deferred to greater depths in old slabs (Fig. 8), up-slab flows would not affect the subduction efficiency of H2O and CO2 (see section IV.7). In other words, cold slabs promote H2O and CO2 recycling even if within-slab flows are upward along the slab.

IV.5 Effect of Serpentinization

It is commonly proposed that hydration of slab lithospheric mantle takes place in fast-spreading oceans at the outer rise near the trench [49]; subsequent dehydration leads to slab embrittlement, which has been invoked for interpreting subduction zone seismicity [50, 41, 51, 52]. Observational seismic studies infer the extent of serpentinization to be 5%–31% [53, 54, 55], which corresponds to $\mathrm{H_{2}O}$ content of 0.6–3.5 wt% for hydrated slab mantle. Given this uncertainty of $\mathrm{H_{2}O}$ content and its control on the within-slab flux pattern (e.g., Fig. 7), we test how this factor affects slab volatile storage and surface fluxes by varying the $\mathrm{H_{2}O}$ content in the incoming slab mantle lithosphere from 0.5 wt% ($\sim$4% serpentinization) to 10 wt% ($\sim$76% serpentinization). Panels in Figure 9 show the resultant variation of $\mathrm{CO_{2}}$ and $\mathrm{H_{2}O}$ fluxes along the slab surface, where 1 wt% H2O is the value used for the reference model. Firstly, the position of onset for lithosphere devolatilization moves from depth of $\sim$125 km to $\sim$80 km as the hydration extent of incoming slab increases. This is because increasing the bulk H2O content while keeping bulk CO2 constant (negligibly small 0.02 wt%) decreases the onset temperature of devolatilization. Secondly, the magnitudes of slab surface fluxes are substantially increased by the elevated basal $\mathrm{H_{2}O}$ content. Figure 9e–h show that serpentinite dehydration can dominate other controls of slab surface $\mathrm{CO_{2}}$ and $\mathrm{H_{2}O}$ fluxes and thus can play a crucial role in arc magmatism and subduction-zone volatile recycling [56].

In Figure 10, the differences in CO2 and H2O loss relative to that of the reference model are mapped for different hydration states of the slab mantle. For the smaller extent of hydration (Fig. 10a–b), the weakened infiltration leads to general reduction in reaction progress, as demonstrated by the decrease in both the CO2 gain at the sites of carbonation (red stripes in Fig. 10a), and the CO2 loss at the sites of decarbonation (blue stripe/regions in Fig. 10a). In particular, there is a stripe pair above the basalt–gabbro interface in Figure 10a, but vertically in reverse order to that in Figure 5a. This is still attributable to the diminished infiltration flux as explained above. By contrast, for higher extent of slab mantle hydration (Fig. 10c–f), stripe pairs similar to those in Figure 5a also appear in Figure 10b & e, but are pushed to upper levels due to the increased infiltration flux. Moreover, for very strong infiltration as in Figure 10e, decarbonation is so strong that even the carbonation within the basalt layer is offset (white stripe at $\sim$2 km in Fig. 10e). If inspected closely, a red–blue stripe pair appears at the peridotite-gabbro interface at $\sim$6.5 km in Figure 10e, which is not the case for the reference model in Figure 5a, and the reason is simply that the stronger infiltration makes the stripe pair discernible.

Quantitatively, for the intensively infiltrated basalt layer in Figure 10c & e, $\mathrm{CO_{2}}$ loss can be enhanced to as high as 100%, whereas for the even more intensively infiltrated sedimentary layer, the rise in $\mathrm{CO_{2}}$ gain can be as high as 300%. In other words, enhanced infiltration extracts an extra amount of $\mathrm{CO_{2}}$ from the gabbro and basalt layers along flow pathways, but much of the extracted $\mathrm{CO_{2}}$ is merely redistributed into the top sedimentary layer. By comparing Figure 9c and e, or 9c and g, it can be inferred that not all of the infiltration-induced CO2 transfer is deposited in the sedimentary layer, a fraction of it escapes the slab top and contributes to the surface flux increase.

$\mathrm{H_{2}O}$ loss or gain differs from that of $\mathrm{CO_{2}}$. For the basalt and gabbro layers, both are infiltrated by H2O-rich fluids that cause decarbonation above the layer contacts. As illustrated in the companion paper Part I, fractionation caused by decarbonation leads to a larger extent of overall devolatilization. The enhanced H2O loss at depths $\sim$2.5 km and $\sim$6.5 km in Figure 10d & f is a result of this. For the top sedimentary layer, however, it is infiltrated by flow down a CO2 concentration gradient, so the fractionation caused by carbonation leads to inhibited H2O loss, as reflected in Figure 10d. Different from that in Figure 10d, the top sedimentary layer in Figure 10f experiences more H2O loss. This is because the stronger infiltration down the CO2 gradient also induces dehydration that counteracts the fractionation effect that is dominant for the case in Figure 10d.

It is worth noting that field studies document the occurrence of a metasedimentary layer immediately on top of serpentinized slab mantle, where the sediments firstly experience decarbonation/hydration and secondly carbonation/partial dehydration [57, 11]. The flow in the outcrops are inferred to be confined into up-slab fluid channels. This field-documented two-stage flow event can be roughly elucidated by the current model. Taking Figure 10e for illustration, if there are channels tapping fluids from dehydrating mantle lithosphere in an up-slab direction (e.g., Fig. 7d), then the subducting sediments will be infiltrated by very H2O-rich fluids, leading to the first stage decarbonation/hydration. At greater depths, when the tapped fluids scavenge CO2 by passing through gabbro and basalt layers, carbonation reaction can take place in the sedimentary layer, accounting for the second-stage flow documented in the field.

IV.6 Effect of Sediment Removal

For all of our models so far, a sedimentary layer of 500 m is assumed to cover the slab top and not detach during subduction. However, the sediments are subject to off-scraping in the accretionary prism of subduction zones [e.g., 58]. In addition, the slab surface sediments that survive off-scraping at trenches are further subject to removal by diapirism. Previous studies [59, 60] showed that a sedimentary layer $>$100 m thick and $\sim$200 $\mathrm{kg\,m^{-3}}$ lighter than surrounding peridotites may form rising diapirs that traverse the mantle wedge. Kelemen and Manning [4] propose that diapirs entraining carbonates are an important avenue for carbon release to mantle lithosphere. We therefore assess the impact of sediment removal on slab volatile storage and emission based on a simplified model setup: (i) the topmost 500-m-thick sedimentary layer in the reference model is replaced by basalt of the same composition as the underlying layer, mimicking the scenario of complete removal of sediments; (ii) this replacement basaltic layer is excluded during calculations of the resultant volatile fluxes and subduction efficiency (see next section). (iii) fluid flow angles are varied among 22.5*∘, 45∘, 90∘, 135∘, 157.5∘* to examine the resultant slab surface fluxes.

For this three-layer model with a 90*∘* flow angle, Figure 11 maps the volatile loss within the slab. Regarding CO2 loss (Fig. 11a), except for the thin sheet of CO2 gain above the basalt–gabbro interface, the entire plate loses CO2, and the loss can reach a maximum of $\sim$70% at the infiltration site from gabbro to basalt. The thin sliver of up to 25% CO2 gain is due to the infiltration-related carbonation, as previously explained. In contrast, H2O loss in Figure 11b is similar to that of the reference model. This is because the sedimentary layer, even if it survives the removal, should dehydrate and not restrict $\mathrm{H_{2}O}$ transport from beneath.

Figure 12 provides the slab-surface CO2 and H2O flux profiles for the three-layer model with the various flow angles ($\theta$). It shows that removal of the capping sediments leads to an overall increase of the CO2 fluxes, and additionally to the emergence of a broad flux peak at the subarc depth range. In particular, when fluid flow is nearly up-slab (Fig. 12k & l), the position of the flux peak coincides with the average position of global arc volcanoes projected onto slab surface [21]. It is also interesting to compare between the flux profiles of the four-layer (Fig. 12a) and three-layer models (Fig. 12g), both with a flow angle of 90*∘*: the carbonation reaction in the sedimentary layer clearly suppresses the CO2 flux peak that would appear without it. Up-slab flows ($\theta>90^{\circ}$) can therefore avoid the more reactive sediments at greater depths, and eventually contribute to slab surface $\mathrm{CO_{2}}$ fluxes. On the other hand, removal of the capping sediments does not increase the slab surface H2O fluxes, except that it removes the flux dip at $\sim$160 km in the four-layer reference model (Fig. 12b). It is clear from our analysis of the open-system reference model (section IV.2) that the flux dip stems from porosity reduction due to enhanced carbonation in the sedimentary layer, thus the dip is eliminated after the removal of the sedimentary layer.

The simplified three-layer model represents an end-member case where the subducted slab is devoid of sediments. By comparing the total CO2 and H2O amount in the three-layered slab with those in the four-layered slab in the reference model, we find that sediment removal causes $\sim$11 % CO2 loss and $\sim$9 % H2O loss from the slab. It is noteworthy that these estimates of CO2 and H2O loss via sediment removal do not account for the thermal effect resulted from the direct contact between the basalt layer and mantle wedge. Further heating of the basalt layer by mantle wedge can increase the CO2 and H2O loss from the slab, and thus their fluxes along the slab surface.

IV.7 Subduction Efficiency of H2O and CO2

The subduction efficiency of H2O (and CO2) in this study is defined as the quotient of its total mass residing in the solid rock phase beyond $\sim$200 km deep over that residing in the incoming slab. To demonstrate how the various factors in the preceding sections affect the subduction efficiency, we plot them in Figure 13.

Volatile redistribution among different layers shown in Figure 10 implies that the subduction efficiency for individual layers depends on flow direction; however, Figure 13a and b suggest that flow direction exerts almost no influence on the subduction efficiency for the entire slab. Comparison between the open-system results in Figure 13a and b indicates that the subduction efficiency of CO2 is reduced by $\sim$7%, due to the absence of a sedimentary layer where carbonation reaction can occur. Figure 13c illustrates that slab age strongly influences volatile subduction efficiency—the older the slab, the higher the efficiency. In particular, the 60-Ma-old cold slab is nearly 100% and 93% efficient in recycling CO2 and H2O beyond mantle depth of $\sim$200 km. In addition, Figure 13d shows that the H2O content of hydrated slab mantle exerts a much stronger control on the efficiency of volatile subduction. For extensively hydrated slab mantle (e.g., 10 wt% H2O or $\sim$76 % serpentinization), the recycling efficiency of CO2 can be reduced to $\sim$80% and that of H2O to as low as $\sim$21%.

Given that reactive flows are path-dependent, the finding that the subduction efficiency in Figure 13a and b is insensitive to flow direction is surprising. However, it can be understood from the two assumptions in the model. The first is the local equilibrium assumption, which dictates that, as long as the $P$–$T$ structure of the slab is constant, steady-state volatile production within the entire slab is determined by $\rho_{s}\,v_{s}\partial\phi/\partial z$ (see Appendix A). Since porosity $\phi$ further depends on bulk composition, which is affected by flow directions, in theory the volatile production ($\rho_{s}\,v_{s}\partial\phi/\partial z$) should be path-dependent and thereby give rise to path-dependence of subduction efficiency. However, changes in bulk H2O and CO2 content induced by varying flow direction are small and spatially limited such that the path-dependence of subduction efficiency becomes indiscernible in Figure 13a and b. Thus, the flow angle ($\theta$) shifts the volatile flux peaks, but not the overall amount of devolatilization. The second relevant model assumption is the uniform flow direction across slab, which makes strongly focused flows unlikely and hence prevents substantial alterations of bulk volatile content.

The subduction efficiency of H2O has also been estimated by previous studies that don’t consider CO2. For example, van Keken et al. [46] estimate an average subduction efficiency of 32% for H2O in the present-day subducting slabs. In their model, the basal slab-mantle layer is 4 km thick and contains 2 wt% H2O. In terms of the total amount of H2O infiltration, this is roughly equivalent to a H2O content of 5.3 wt% in the 1.5 km thick basal layer in our model. As shown in Figure 13d, our model would yield a H2O subduction efficiency of $\sim$35%. Nonetheless, the average slab age in the model by van Keken et al. [46] is older than 10 Ma that is adopted in the calculation for Figure 13d, implying that our model would yield a H2O subduction efficiency higher than 35% for the present-day global subduction zones. The higher subduction efficiency of H2O is expected because our model accounts for open-system effects. As demonstrated in Figure 13, the subduction efficiency of H2O is generally higher in open-system models than in closed-system ones. This is mainly due to the fractionation-induced inhibition of H2O loss in the bulk interior of subducted lithological layers (see section IV.2). Only in the case of strong infiltration in Figure 13d does the H2O (and CO2) subduction efficiency become lower in open systems than in closed systems. Similar to H2O, CO2 subduction efficiency is also elevated by open-system fractionation. Moreover, comparison between Figure 13a and b suggests that, for closed systems, removal of top sedimentary layer exerts almost no effect on the CO2 subduction efficiency; however, for open systems, the CO2 subduction efficiency is reduced by $\sim$7% due to the sediment removal, highlighting again the importance of open-system effects.

V Discussion

By treating the subducting slab as a non-deformable plate and assuming local equilibrium and uniform flow direction, our model solves for the steady-state reactive flows of H2O and CO2 and quantifies their fluxes from and storage within the slab. The merit of the model is its capability to couple open-system thermodynamics with fluid flow in a more realistic and consistent way. However, in achieving this capability, the choices made in model setup also impose limitations.

The first limitation lies in neglecting partial melting. Various studies have focused on the effects of $\mathrm{H_{2}O}$ on partial melting, and suggest that melting primarily occurs in the overlying mantle wedge due to solidus depression caused by slab-derived $\mathrm{H_{2}O}$ [61, 47]. Field studies in volcanic arc settings, however, discovered adakites that are interpreted to be evidence of slab melting in hot geotherms [62, 63]. When $\mathrm{CO_{2}}$ is taken into account, experimental studies suggest that even in hot subduction zones, modern slab geotherms are still below the solidi of carbonated basalts and peridotites, and that only carbonated sediments are promising for carbonatite liberation [1, and references therein]. On the other hand, Poli [64] shows that partial melting of hydrous carbonated gabbros can produce carbonatitic liquids in warm subduction zones, as long as the gabbros are Ca-rich, a result likely from extensive hydrothermal alteration. In light of these experimental findings, neglecting the partial melting of slab-mantle peridotite in our model is a reasonable choice, but ignoring the melting of mafic and sedimentary layers could overestimate the subduction efficiency of CO2, especially when migration of carbonatitic melt can remove considerable CO2 from the slab. Moreover, absence of partial melts in the current model indicates that the thermal feedbacks from associated latent heat and melt transport is unaccounted for [e.g., 65]. The influences of perturbation to the slab thermal structure caused by the loss of partial melts and/or sediments (section IV.6) need to be further assessed in future research.

The second limitation stems from the binary H2O–CO2 molecular fluids assumed for the liquid phase in our model. Variations of pressure, temperature and particularly oxygen fugacity during subduction can give rise to the dominance of methane (CH4) or precipitation of graphite [e.g., 33, 66]. Moreover, ionic carbon species (e.g., HCO${}_{3}^{-}$, CO${}_{3}^{2-}$) are inferred to be common in subduction zone settings [e.g. 9, 10], although detailed thermodynamic modelling suggests that molecular carbon species can dominate over ionic ones in fluids derived from subducted basalts [67] and sediments [13]. The latter study also shows that potassium in ionic form can be depleted from the slab-top sediments, destablizing hydrous phases like micas and thus facilitating volatile release. In addition to potassium, the activity of dissolved silica and aluminum also affects the stability of silicates that participate in carbonate dissolution reactions [68]. Therefore, omission of the direct contribution from carbon species other than CO2 and indirect contribution from dissolved non-volatile elements (e.g., Si and Al) might underestimate the water and carbon fluxes from the slab, and in turn overestimate their recycling efficiency.

The third limitation arises from the assumption of uniform flow directions. As explained in section II on the model setup, choosing the flow direction $\theta$ as a free parameter enables the model to explore the variation of general patterns of within-slab flows predicted by various dynamic models; however, it suffers from losing the details of flow dynamics within the slab. In particular, the strong flow focusing, widely documented in the field [e.g., 27, 31, 28, 29, 11, 34], and predicted by dynamic models [e.g., 15, 26, 16, 14], cannot be simulated in the current model. As demonstrated by Wada et al. [69], heterogeneous hydration associated with local flow focusing can promote H2O release and reduce its subduction efficiency. Therefore, further work remains to incorporate dynamics into the open-system flow model to evaluate how flow focusing can affect the slab-surface fluxes of CO2 and H2O, and their recycling efficiency by subduction.

The fourth limitation relates to the assumption of no compaction in the slab. To assess the effect of compaction, we compare the magnitude of fluid fluxes in our reference model with that in the compaction model by Wilson et al. [14]. Our slab model (section IV.2) initially contains $\sim$2.60 wt% H2O on average, yielding a maximum surface flux of $\sim$0.1 m kyr*-1* (Fig. 6f). In contrast, the compaction model initially contains $\sim$1.85 wt% H2O on average, yielding a maximum surface flux of $\sim$0.8 m kyr*-1*. Taken at face value, the compaction model generates a flux that is eight times larger than that from our non-compaction model. However, as shown in Figure 5b, the H2O flux drop in our reference model is caused not only by neglecting compaction, but also by the open-system fractionation, which can account for up to $\sim$70% of the flux drop. If this open-system effect were considered by Wilson et al. [14], the magnitude of H2O fluxes in their model would become comparable to that in our reference model, indicating that the influence of compaction is insignificant with respect to the magnitude of slab surface fluxes. Nevertheless, as stressed above, compaction may be important if it is a control on fluid flow pathways within the slab.

In interpreting our model results on the slab surface fluxes, an additional assumption is made that magmas produced by flux melting of the mantle wedge traverse the wedge vertically without any lateral deflection to feed arc volcanoes. It follows that, if the average position of global arc volcanoes is projected downwards onto the slab surface, it should correspond to a peak in slab surface flux. It is under this assumption that the up-slab flows in Figure 6 are necessary to produce the melt feeding arc volcanism. However, two-phase geodynamic models of the mantle wedge suggest that the melt and fluid transport through the wedge is more complex than simply vertical, and depends on many factors including solid rheology and grain size in the wedge [14, 70, 71]. In addition, the coincidence of flux peaks and the projected volcanic arc position in Figure 6 and 12 is premised on the hot thermal structure of the 10-Ma-old slab. As demonstrated in Figure 8 and 9, old slab shifts the flux peak to greater depths, and higher hydration extent of the slab base shifts the peak to shallower depths. Since slower convergence rate leads to a warmer subducting slab [72], it can be expected that the position of the flux peak also varies with convergence rate. How the shift in flux-peak position influences the dynamics of melt or volatile transport in the wedge, and eventually the genesis of arc magmas, remains to be explored by future models that include both the dynamics and thermodynamics in the wedge.

VI Conclusion

In summary, we assess the controlling factors for $\mathrm{CO_{2}}$ and $\mathrm{H_{2}O}$ release from subducting slabs. We find that up-slab flows are necessary to form a flux peak at subarc depths, which could be responsible for the magmatic genesis of arc volcanoes if the flux-induced melts traverse the mantle wedge vertically. Such a result aligns with the findings of dynamic studies by Wilson et al. [14] where a flow channel is formed to induce up-slab flows. Sufficiently hydrated slab mantle can relax the above requirement for flow directions, but when the slab mantle is additionally carbonated, the onset position of its devolatilization becomes deeper, still making up-slab flows important for supplying maximum volatile release at subarc depths. Moreover, H2O infiltration from dehydrating slab mantle can mobilize CO2 in the overlying gabbroic and basaltic layers, but only redistribute the carbon to the slab surface sediments where significant carbonation takes place. The subducting sediments, if passing the accretionary prism and surviving diapiric removal and partial melting, can be an important barrier to CO2 liberation from slabs. Conversely, loss or partial melting of the sediments can be crucial avenues for CO2 transfer from subducting slab to mantle wedge. The subduction efficiency of H2O is $\sim$20%–90% in our model, generally higher than that in previous closed-system models that consider only H2O. The CO2 subduction efficiency is $\sim$80%–100%. The higher volatile subduction efficiency compared to closed-system models stems from the open-system fractionation effect that inhibits overall devolatilization in the interior of lithological layers. The absolute values are likely overestimates because the effects of flow focusing and fluid species other than H2O and CO2 are ignored in the model, which should be considered in future studies.

Appendix A Appendix

At steady state, the fluid fluxes can be calculated according to equation (5) that further reduces to:

[TABLE]

Figure 14 shows part of the numerical grid established on the slab coordinate system (Fig. 1b), where the fluid fluxes $F=v_{f}\phi$ can then be expressed through integration along flow path ($l$):

[TABLE]

For the reference model with flow direction $\theta=90^{\circ}$, flow paths are perpendicular to the slab, so the integration is along the $z$-axis direction in Figure 14. At any specific depth (dashed lines in Fig. 14), the onset of devolatilization leads to an abrupt elevation of porosity ($\phi$) and thus a pulse in the $\partial\phi/\partial x$ profile along the iso-depth traverse. Figure 15 shows the $\partial\phi/\partial x$ and flux profiles along various iso-depth traverses within the modelled slab, and Figure 15a illustrates that the pulses in $\partial\phi/\partial x$ appear at different horizontal positions in the slab. According to equation (12), when flow trajectories pass the locations of the pulses, the calculated fluid fluxes will inherit these pulses. For the reference model where flow paths are normal to $x$-axis and the pulses are horizontally apart from one another (Fig. 14), the inherited pulses lead to the fluctuation in the fluid flux distribution (Fig. 15b–c and Fig. 4g). Beyond the slab distance $\sim$160 km where the basal slab mantle dehydrates, the envelope of the onset of devolatilization (red solid line in Fig. 14) is surpassed, so there are no pulses in $\partial\phi/\partial x$ at any depth. The flux fluctuations disappear accordingly as in Figure 15b–c and Figure 4h.

The analysis above indicates that refining numerical grid will produce more but weaker fluctuations in slab surface fluxes before $\sim$160 km, but the trend and main features of flux profiles stay unchanged. To focus on the general trend and main features, we smooth the flux profiles by averaging over a span of neighboring grid points. As illustrated in Figure 15b–c, smoothing with a span of 12 grid points (3 km) well preserves the flux trend and features, so it is applied to all the slab surface flux profiles in this study.

Acknowledgements.

We thank the three anonymous reviewers for their comments that helped us improve this manuscript. We thank the Isaac Newton Institute for Mathematical Sciences for holding the Melt in the Mantle program sponsored by EPSRC Grant Number EP/K032208/1. Support from Deep Carbon Observatory funded by the Sloan Foundation is acknowledged. M.T. acknowledges the Royal Society Newton International Fellowship (NF150745). D.R.J acknowledges research funding through the NERC Consortium grant NE/M000427/1, NERC Standard grant NE/I026995/1, and the Leverhulme Trust. This project has also received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement n*∘* 772255). This contribution is about numerical modelling, so it does not depend on experimental or field data. Relevant data and equations for reproducing the model results are already contained in the text. Nonetheless, we provide an example code on the usage of the thermodynamic parameterization: https://bitbucket.org/meng_tian/example_code_thermo_module/src/master/.