Regimes of carbonate subduction as a function of lithospheric mantle hydration state

,


INTRODUCTION
When oceanic lithospheres sink into the mantle in subduction zones, a fraction of the chemical elements (volatiles, metals and metalloids) are progressively dissolved out of the slabs by fluids at increasingly higher pressure and temperature (P-T) conditions.The chemical budget of the downgoing plate, mantle wedge, continental crust, and atmosphere hinges on the temporal and spatial sequence of this mass transfer, but several aspects of this process remain obscure.The destabilization of hydrous minerals, such as lawsonite, amphibole, phengite, serpentine and chlorite in the downgoing slab (Cook-Kollars et al., 2014;Grove et al., 2006;Schmidt and Poli, 1998;Turcotte and Schubert, 2002) from shallow (<80 km, fore-arc) to sub-arcs depths (<150 km) produce the fluids that mediate subsolidus mass transfers.The fate of carbonate in particular, which is the main carbon carrier to subduction zones, depends on the availability of hydrous fluids (Ague and Nicolescu, 2014;Connolly and Galvez, 2018;Dasgupta and Hirschmann, 2010;Gorman et al., 2006;Johnston et al., 2011;Plank and Manning, 2019;Sverjensky et al., 2014b).These fluids may remobilize carbon and other chemical elements throughout the slab, or transport it out of the slab to the mantle wedge.
Informed by the above observational and experimental constraints, quantitative subduction zone chemical budgets have relied on coupled petrological-thermodynamical numerical models.Most models (e.g.Connolly, 2005;Gorman et al., 2006;Tian et al., 2019) have relied on molecular (usually binary) CO 2 -H 2 O fluid models for lack of better alternatives.Although these petrological-numerical models of subduction zone chemical dynamics provide an accurate estimate of C solubility over a wide range of P-T conditions, only CO 2 and H 2 O are treated as mobile elements.Therefore, silica and other rock components have been implicitly considered as immobile, or infinitely compatible with rocks.We refer to this important limitation as the "infinite compatibility" limit, and it has important consequences for the temporal evolution of slab mineralogy, and slab decarbonation in particular.
Previous studies of slab decarbonation have all monitored a class of devolatilization reactions (commonly called "decarbonation") represented, e.g. by the incongruent dissolution reaction below (Aranovich and Newton, 1999;Harker and Tuttle, 1956): SiO 2 (solid) + CaCO 3 (+H 2 O) = CaSiO 3 + CO 2 (+H 2 O) (Eq. 1) where the chemical formula represents thermodynamic components of the system (as opposed to physical entities), and (+H 2 O) expresses that the reaction requires water as a solvent and catalyst at typical subduction zone conditions.The silica component is typically contained in silicate minerals, but is always part of the rock, SiO 2 (solid).In this context, slab decarbonation may result from prograde evolution (i.e.rising P-T conditions change the thermodynamic activity of all rock components), or via transport of fluid CO 2 and H 2 O into or out of the system, modifying the thermodynamic activity of CO 2 , and driving equation (1) to the right or left (Ferry and Dipple, 1991;Rumble, 1982).However, subduction zone fluids also contain dissolved solutes such as Na, K, Si, Al, Ca, Mg and Fe (Manning, 1994;Spandler et al., 2007;Wohlers et al., 2011).Subduction fluids are therefore not strictly binary, nor purely molecular.More fundamentally, this means that changes in chemical potentials (μ) of the rock components SiO 2 , CaO, K 2 O and CO 2, that drive Eq. 1 may also result from non-prograde (i.e.fixed P-T conditions) transport of, e.g., dissolved Si, Ca and alkali in and out of the system by metasomatism.Hence, if rock-forming components other than CO 2 and H 2 O are at least partly mobile, Eq. (1) may also be written as: SiO 2 (aq) + CaCO 3 = CaSiO 3 + CO 2 (Eq.2) where SiO 2 (aq) denotes that the component SiO 2 may now be found in a mobile fluid phase as well.
Of course, any more complex decarbonation reaction involving Al, Ca, Mg reactants may be rearranged the same way (Galvez and Pubellier, 2019).
Decarbonation reactions driven by Si metasomatism, in particular (e.g.Eq. 2), have been reported on various occasions in natural settings (Ague, 2003;Bucher, 1999;Ferry et al., 2011).However, recognizing the role of metasomatic-decarbonation processes in the overall subduction zones chemical budget has long remained elusive due to the lack of computational tools to connect volatile and nonvolatile (i.e.alkali, alkali earth, transition metals and metalloids) behaviour at a more fundamental level.Recently, Sverjensky et al. (2014) called attention to the capacity of the HKF model (Helgeson et al., 1981) to handle electrolyte chemistry at elevated pressure conditions owing to the knowledge of water dielectric constant over a range of P-T conditions (Fernández et al., 1997(Fernández et al., , 1995;;Pan et al., 2013).Meanwhile, Galvez et al. (2015Galvez et al. ( , 2016) ) have shown that this electrostatic formalism is compatible with conventional petrological-thermodynamical modelling tools relying on Gibbs Energy Minimization of complex mineral/fluid systems (GEM, Perple_X (Connolly, 2005)).They showed that electrolytic fluid solutions are compatible with mixed-volatile fluids where H 2 O activity may significantly depart from unity (Galvez et al., 2016(Galvez et al., , 2015)), which is commonly observed at the elevated P and T conditions of the Earth interior (Aranovich and Newton, 1999;Connolly and Cesare, 1993).As a result, a new metric of acidity-basicity for high P fluids has been defined, ΔpH (Galvez et al., 2016(Galvez et al., , 2015) ) which accounts for changing dissociation constant, activity and dielectric properties of water in high P and T aqueous solutions (c.f.Method section for details).But the second and most important advantage of this hybrid approach, which we exploit here, is that it relaxes the long-standing assumption of "infinite compatibility" of non-volatile components in most previous petrologicalthermodynamic numerical models of subduction zones Applying this hybrid approach to the average GLOSS sediment composition (Plank and Langmuir, 1998), Connolly and Galvez (2018) recently suggested that the mobility of potassium in fluids released from a subducting sediment may lead to early sediment depotassification, and therefore dehydration.But this result may have been amplified by considering pure H 2 O entering the sediment (2018), which is unrealistic.This bias can be easily resolved with a systematic approach of sub-solidus slab mass transfers where all lithologies, sediment, mafic and ultramafic are partly soluble in the fluid phase.
Here we report the first a first attempt at a quantitative chemical budget of subduction zones built on a petrological-thermodynamical numerical approach where all rock components have a finite, experimentally constrained, solubility in the fluid phase.We examine the process of slab decarbonation specifically, qualitatively and quantitatively.For the sake of simplicity, we adopt a simple model that considers dissolution, fractionation and vertical transport of aqueous fluid and their dissolved solutes within the slabs from global subduction zones, and we explicitly consider a range of upper mantle hydration states (Schmidt and Poli, 1998;Van Keken et al., 2011).Our model relies on a set of geological, chemical, geodynamical and computational assumptions which are discussed and tested quantitatively when possible.Overall, our work supports the notion that the subduction efficiency, measuring the fraction of subducted mass flux returned to the mantle (beyond sub-arc depth), is not a fixed and uniform value of the global subduction zone as is commonly assumed, for example, in global biogeochemical models that consider subduction processes (Caldeira, 1991;Hoareau et al., 2015;Volk, 1989).It is a variable that differs for all elements and varies in both space and time, as conjectured by e.g.Galvez and Pubellier (2019).We identify three main regimes of carbonate cycling as a function of upper mantle hydration state and slab thermal structure.Our qualitative and quantitative results exploit the effect of an intraslab metasomatic redistribution of alkali, alkali earth, alumina and silica on carbonate subduction and mantle wedge refertilization, and we report first estimates of fluxes for those elements through and out of the slab, and we discuss their uncertainties.These insights inform specific regional targets for future research into the cycle of lithospheric carbon and other major rock-forming elements.

Chemical structure of the subducting slab: model setup
For each slab, the oceanic lithosphere is composed of four layers: a sedimentary layer underlain by a basaltic upper crust, a gabbroic lower crust, and a partially serpentinized upper mantle layers at the bottom (see Fig. 1).The composition of the sediment is from Plank (2013), including H 2 O and CO 2 , with a minor correction to a few slabs to ensure consistent Ca/C ratio where data of C wt% is lacking or particularly uncertain (see Fig. S1 with equation given in the caption).The total CO 2 trench input is ca.186 Mt CO 2 /yr.Results of a model with higher C contents, ca.267 Mt CO 2 /yr from Clift (2017), are also presented for comparison in the supplementary materials.The altered basalt (1.5 km thick) and gabbro (5 km thick) compositions are from Pearce (1976).It has been shown that seafloor alteration of the oceanic lithosphere enriches the top basaltic layer in potassium (Alt and Honnorez, 1984;Staudigel et al., 1989).To simulate this process, we added a 300 m thin layer containing 1.5 wt% K 2 O, underlain by 1.2 km of basalt with a concentration of 0.1 wt% (Supplementary Table S1).
Overall, this corresponds to a K 2 O content of the altered basalt layer of around 0.37 wt%, which is consistent with natural observations (Alt and Teagle, 2003).The altered upper mantle composition is from Hart and Zindler (1986), and its thickness and water content is a variable and discussed below.
The global input of CO 2 associated with the basalt + gabbro (igneous) layer is ca.86 Mt CO 2 /yr (23 Mt reflect the global estimate from e.g.Kelemen and Manning, (2015).The Makran subduction zone is not included in our simulations for a lack of thermal model (Syracuse et al., 2010).
The H 2 O content and distribution within the slab determines the abundance of solvent available for element transport during the subduction process, e.g.see reaction 1 and 2. We have adopted a range of hydration states for our slabs.Our trench sedimentary H 2 O composition is from Plank (2013) and the mass of subducted sediment per year is from Clift (2017), leading to a trench sedimentary H 2 O input ca.269 Mt H 2 O/yr (Table 1).However, we find that about 74 Mt H 2 O/yr is lost as a free-fluid phase at the first grid of our simulation.This reduces the effective sedimentary input flux to ca. 195 Mt H 2 O/yr (Table 2).When combined with the igneous input of 638 Mt H 2 O/yr (Jarrard, 2003), the overall sedimentary and igneous flux in our models is 833 Mt H 2 O/yr (Table 2).This is consistent with the general consensus on the amount of H 2 O structurally locked in sediment and altered igneous minerals (Hacker, 2008;Jarrard, 2003;Rüpke et al., 2004;Van Keken et al., 2011), i.e. ~720 Mt H 2 O/yr.Fuca plate beneath Cascadia (Clift, 2017).The integrated H 2 O flux in d at the slab surface is ca.28 Mt/yr down to 120 km.
The main uncertainty in trench water input is that associated with the serpentinized altered upper mantle (AUM) on a global scale.In detail, Schmidt and Poli (2003) propose an upper-estimate AUM input of ca.860 Mt H 2 O/yr, while the upper-estimate of Rüpke et al. (2004) 1 and 2), i.e. an average of 2 wt% H 2 O over a 4 km AUM layer using a linear variation scheme from 0 wt% (AUM bottom) to 4 wt% (AUM top).
For our low-end hydration state system, we noted the consensus between the low-end AUM hydration estimates of Schmidt and Poli (2003), 490 Mt/yr; Rüpke et al. (2004), 240 Mt/yr;and Van Keken et al., (2011) 1 and 2).Note that our global estimates may locally overestimate (Canales et al., 2017), or underestimate (Cai et al., 2018) the water input for some hot or cold slabs, respectively.This point has been tested on a case by case basis when possible.
As a first-order approximation, pressure within the subducting slab is lithostatic using an average rock density of 3300 kg/m 3 .Estimates of slab thermal structures have been obtained from both analytical or numerical methods in previous works (e.g.McKenzie, 1969;Rüpke et al., 2004;Chemia et al., 2015;England and Wilkins, 2004;Furukawa, 1993;Abers et al., 2006;Van Keken et al., 2011).
Here we use the numerical results obtained by Syracuse et al. (2010) using input data of plate convergence rate and slab geometry.The slab name, slab length and sediment thickness are taken from Clift (2017).
Most of our article refers to the extraction efficiency of a given oxide i, it is denoted η i .The extraction efficiency of carbon, in particular, is a a key parameter of long-term climate and biogeochemical models.This variable corresponds to the fraction of the trench input of component i released from the slab within a certain depth interval.For fore-arcs, η forearc is defined at depth <80 km, and for sub-arcs, η subarc is defined at depth <150 km.The beginning of the model is defined at 20 km depth.The extraction efficiency η , the opposite to the subduction efficiency σ=1-η, describes the fraction of the total input released to the mantle wedge before a given depth (Galvez and Pubellier, 2019;Volk, 1989).

Computational strategy
Simulating the processes of internal thermodynamic equilibration between fluids and solids by Gibbs energy minimization require that the system is closed, that is, that its composition is fixed (Connolly, 2009).However, the geological fluids we simulate here contain solvents (typically H 2 O and CO 2 ) and solutes (volatile and non-volatile) that may migrate and modify rock compositions along their path (open system).Open system processes that involve changes in environmental conditions (P, T) and bulk rock chemistry may still be addressed by discretization of the subducting slab (Connolly, 2005;Gorman et al., 2006).This is used here in a simple dissolution-fractionation-transport model.Our fluid dynamics is similar to the model of Gorman et al. (2006), but our thermodynamics contain essential improvements with respect to the nature and composition of the fluid phase.This makes our predictions of rock mineralogies and element fluxes more geologically accurate, and it distinguishes our model from previous works.
In our initial model setup, we discretize a column of rock into grids (resolution provided below).We assign the bulk composition to each grid according to its position within the four distinct horizontal layers (see supplementary materials for detailed composition of all the layers) composing the slab, and we proceed as follows.In the first step, we compute the amount of free molecular H 2 O and CO 2 solvent by internal equilibration of the bottom grid in the column of rock using the Perple_X software (Connolly, 2005).The solid solution models used for this Gibbs energy minimization are listed in supplementary materials.Melts are not considered.Chemical potentials of rock components are obtained from Perple_X, and the concentration of dissolved aqueous species are then computed following the procedures developed in Galvez et al. (2015).This back-calculation strategy, as described in Galvez et al. (2015) (see also Connolly and Galvez (2018)) implies that the fluid electrolytes are not part of the general Gibbs energy minimization algorithm.It is therefore particularly accurate in predicting phase equilibria and equilibrium fluid compositions (electrolytes, solvents) when the Gibbs energy of a fluid component (non-solvent) represent a negligible fraction of the Gibbs energy of that component in the system (i.e.rock dominated limit, Galvez et al. (2015)), which is the case most of the time in the present work.In particular, here we exploit the fact that cycles of fluid back-calculations conducted in the rock-dominated limit followed by fluid fractionation may lead to profound, but yet overlooked, changes in rock composition.Our speciation algorithm may lose accuracy in predicting the partitioning of a given element between fluid and rock when the fluiddominated condition is reached (Connolly and Galvez, 2018;Galvez et al., 2015).This may typically occur when a component of the fluid-rock system reaches vanishingly low concentrations in the rock compared to that dissolved in the fluid.In that case, if the amount of dissolved component is predicted to be higher in the fluid than the leftover in the rock, only the leftover of this specific component from the rock is dissolved in the fluid.This does ensure mass-conservation during the entire transport and subduction sequence, which is most critical.More importantly, while remaining computationally efficient, our strategy is able to provide enhanced geological consistency and accuracy of subduction zone processes because all chemical components have a finite mobility in the fluid phase as constrained by experimental data.The thermochemical database of mineral solutions is from Holland and Powell (1998).Aqueous solutes follow the revised HKF equations of state from the SUPCRT database (Helgeson et al., 1981;Johnson et al., 1992;Shock et al., 1992;Tanger and Helgeson, 1988) as revised by (Sverjensky et al., 2014a).This hybrid approach is only possible if the Gibbs energies of minerals, solvents and electrolytes at elevated P and T are computed with a consistent convention (Galvez et al., 2015).We use the Benson-Helgeson convention (Anderson, 2005;Benson, 1976;Helgeson et al., 1978), where the apparent molar Gibbs free energies of species (minerals, solvent, solutes) at elevated P and T are computed relative to the Gibbs energy of formation from their constituent elements at 298.15 K and 1 bar (Benson, 1976;Helgeson et al., 1978;Shock and Helgeson, 1990;Tanger and Helgeson, 1988).We use the extended Debye-Huckel activity model for electrolytes (Anderson and Crerar, 1993;Galvez et al., 2015), and the activity model of (Holland and Powell, 2003) for the mixed-volatile molecular H 2 O-CO 2 solvent.In the second step, the mass of solvent H 2 O, CO 2 and dissolved solutes released are fractionated from the bulk composition at local grids.Hence, the mass of all major components (NCKFMASH) in the evolving rock system is treated as a variable.
It reflects the spatial sequence of devolatilization reactions, and the process of internal equilibration within the fluid-rock system as the slab sinks deeper into the subduction zone.In the third step, the free aqueous solution produced, if at all, is transported to the grid immediately above, and the composition of the fluid is added to local bulk composition before recomputing the stable phase assemblage.Complete fluid expulsion assumes high permeability within the reactive zone in the slab (Ague, 2003).Our model implicitly assumes vertical fluid transport within the slab as driven by buoyancy force.Mottl et al.(2004) supported this model at Mariana fore-arc.The dissolutionfractionation-transport procedures are repeated until a depth of 220 km is reached.The grid resolution along in-slab direction is ~20 m within the thin sedimentary and basaltic layers and increases to ca. 200 m within the thicker gabbroic crust and altered upper mantle (AUM) layers.Subsequently, the entire column of rock is moved downward according to the thermal model considered.The descending rate of individual subduction slab is used here to calculate the total mass of rock column being transported via the slab per year (Syracuse et al., 2010).We have chosen an increment of ca. 1 km/step to move down the slab from 20 km depth until 220 km depth, as preliminary tests have shown that this discretization optimizes computation time while warranting sufficient resolution to capture the mass transport process.
Temperature at each grid point is interpolated based on the thermal model from Syracuse et al. (2010).Hotter alternatives, which may be relevant (Penniston-dorland et al., 2015), are addressed in the discussion.Pressure is calculated based on the depth assuming a lithostatic pressure gradient in the slab.Therefore, no specific pressure-temperature step applies due to the heterogeneity in temperature field and slab-surface geometry of each subduction slab.

Devolatilization pattern of Cascadia vs magnetotelluric data
Figure 1c and 1d shows a match between fluid release events -associated with the carbonate sediment at ca. 40 km, basalt layer at 50~80 km, and with the serpentinized mantle at 90-100 km, respectively-and the depth of low resistivity area for the Cascadia subduction zone (McGary et al., 2014).Our model also predicts the formation of a region of high H 2 O flux spreading over several tens of km depth (~50 km horizontally shown in Fig. 1c), linked to upper mantle dehydration, and this feature, too, matches an extensive zone of low resistivity located between 80 and 110 km depth.
Finally, a minor fluid release event within the sedimentary layer at ca. 40 km is also qualitatively consistent with a weak resistivity response observed in this depth interval.This match is taken as a qualitative validation of our subduction model setup.

Interslab variability of element release out of slab: H 2 O
The location of slabs dehydration events varies dramatically across different subduction zones (Fig 2 ).
For example, the main dehydration event of the Nicaragua subduction zone occurs along a rather broad depth range between 125 and 160 km.The serpentine/chlorite dehydration reactions occur first in the cold and low P top of the AUM, which creates a spatially broad fluid plume rising across the mafic and sediment layers (Fig. 2b).This pattern differs in the warm slab of Cascadia, for example, where chlorite destabilizes first in the hot and deeper segment of the AUM, and spatially disconnects it from the serpentine dehydration front (very hydrated top of the AUM, Fig. 2a).The entire Nicaragua slab is dry by 170 km depth, and the entire Cascadia slab is dry by 100 km depth based on the HHS, IHS and LHS model results.By contrast, only the top 500~1000 m of the cold Honshu oceanic mantle (antigorite) may dehydrate at a depth of about 200 km.The fractionated water from AUM then transports to the igneous layer where it temporarily accumulates.Most of the slab water input, locked in chlorite and serpentine (or Phase-A at deeper depth) in the upper mantle or in phengite and lawsonite in the igneous layer (Fig. 2), is predicted to bypass the sub-arc region.
Globally, about 212−218 Mt H 2 O/yr is released to fore-arc depth (20-80 km), and this value is largely independent of the hydration state of the AUM (see Table 2 for LHS, S5 for IHS and S6 for HHS models).By contrast, the water released within the 80 km to 150 km depth shows marked sensitivity to the hydration state of the AUM.Between 561 Mt/yr (η H2O subarc = 48.1 % of trench input, LHS, Table 1) and 1015 Mt/y (η H2O subarc = 41.8 % of trench input, HHS) is released between for-arc to sub-arc depth (from 80 to 150 km).See Table 1 for detailed values.The initial water input in the mantle is assumed to vary linearly along the in-slab direction.The HHS model result is illustrated here.For illustration, the thickness of the sediment is plotted as 1 km so that the phase labels can be seen.

Interslab variability of element release out of slab: carbon
Carbon is present as aragonite, dolomite and magnesite in the sediment layer of our model (Fig. 2).
Figure 3 shows that carbon is released from the rocks in pulses, regardless of the hydration state.For all the slabs, the location of carbon pulses reflects dehydration events (e.g.Fig. 2 and Fig. 3).For example, all the carbon from the Cascadia sedimentary layer is released within a narrow depth range (60 km to 80 km, Fig. 3).
Our sensitivity tests in Fig. 3a to 3e (Cascadia) and Fig. 3k to 3o (Honshu) show that the AUM hydration state does not modify the pulsed nature of carbon release qualitatively.However, the rate and magnitude of this release (and related parameters such as η CO2 ) does prove sensitive to AUM hydration state for some slabs, including Nicaragua (Fig. 3f      A near-linear trend is observed for the global subduction system.The η 2 diagrams for the other slabs can be found in supplementary materials (Fig. S16).

Interslab variability of element release out of slab: non-volatile elements
Most of our qualitative findings above apply to all other rock forming elements.The main observation is that alkali metals, alkali-earth metals, aluminium and silicon are released in discrete pulses controlled by the sequence of dehydration reactions within the slabs (Fig. 6).Similarly, we find that the dependence of mass fluxes to slab hydration state differ between oxides, and between slabs.This means that slabs that are supply limited for CO 2 may be transport limited for SiO 2 or other components.This dependency reflects the element concentrations (Fig. 7).For example, the extent of desilicification and dealuminification of slabs such Costa Rica, Guatemala or Mariana depend heavily on AUM (transport limited), just as the decarbonation did (Fig. 3 and 5).However, the extent of desilicification and dealuminification from the hot slabs of Central Cascadia and Mexico prove sensitive to the AUM hydration state, but this was not the case for their decarbonation (Fig. 5).We also find that individual slab units (sediment, basalts) may fall into a specific transport regime, while the slab as a whole may behave differently.For example, the Nicaragua sedimentary layer, and that of other carbonate-rich subduction zones, are supply-limited with respect to K (i.e.sediment loses all its K by sub-arc depth) shown in Fig. 8n.However, the slab as a whole shows some dependence to its hydration state because K is also distributed within the igneous layer.

Intraslab redistribution of elements
Our model also provides insights into intraslab element redistribution (intraslab metasomatism) which contributes to significantly modify the rock compositions through depth.However, the concentrations of components in a rock are by themselves a poor indicator of metasomatic transport (Aitchison, 1982;Gresens, 1967).This is because concentration changes of a components (Δc i ) is affected by devolatilization reactions.As a result, the molar concentration of a non-volatile component may rise when H 2 O and CO 2 are lost from the system.For example, this phenomenon causes the sharp increase of Δc K2O for the Cascadia and Honshu subduction zones (Fig. 8f and 8v), and the >100% increase of Δc Al2O3 in the lower sedimentary section of Nicaragua (Fig. 8k).Instead, to assess the metasomatic redistribution of non-volatile components throughout the slab, we monitor the relative change in mass of a component scaled by the its initial input, Δm i .Interestingly, we found that most mass changes occur at lithological boundaries, where contrasts in chemical potentials between adjacent lithologies are the largest, in particular between the sediments and the igneous layers.For illustrative purpose,  Most slabs with sediments that are initially low in Ca operate as net Ca sink during subduction: such as Mexico (+120%), Ryukyu (+180%), Kamchatka (70%), Kurile (200%) at ca. 150 km depth at the bottom of the sediment layer (Fig. S11).Slabs that are initially enriched in Ca tend to show net Ca loss.For example, the bottom section of the sediment of Nicaragua, which loses about 50% of its Ca content by 150 km depth (Fig. 8).The Ca loss is almost entirely provoked by infiltration of the H 2 O plume, sourced in the AUM, when it reaches the sediment in a narrow depth interval.All slab sediments lose Na, with Δm Na2O variations ranging from a few percent (e.g.Aegean) to almost 100% (e.g.Ryukyu, Mexico) at subarc depth.Many slabs lose K (carbonated slabs with low K such as Guatemala, Columbia Ecuador, Peru, Costa Rica, Kamchatka), too, but many like North Chile (Δm K2O = +20%), North and South Sumatra, Alaska, Cascadia, Java sediment and top of the Kermadec sediment layer are net sink for (igneous) K, although accumulations never exceed a few relative percent (Fig. S13).
All slab sediments are net source of H 2 O and CO 2 to the mantle wedge.However, locally H 2 O may accumulate, as can be seen in the case of Nicaragua sediments in the 140 km region, where the sediment is nominally dehydrated (Fig 3).Minor phengite may form, transiently, in the sediment layer and this accounts for the transient storage of water.However, this process is difficult to resolve as the amount of K involved are low and close to the fluid-dominated limit.

Metasomatic control on carbonate incongruent dissolution
Figure 9 illustrates the consequence of metasomatic mass redistribution on the release of C from a sedimentary layer.This effect is quantified by computing the difference of chemical potentials (Δ s μ) of component between two models, one where all elements are mobile (μ spec ), one where only molecular CO 2 and H 2 O are mobile (μ nonspec ), hence Δ s μ= μ specμ nonspec .The superscript s indicates that the symbol Δ does not denote a temporal or spatial variation, but a contrast between speciation models.The dissolution of carbonate (in the form of aragonite) in the fluid can be represented by the reaction: where the notations CaCO 3 , CaO and CO 2 refer to thermodynamic components of the system.
We have from Eq. 3: which states that at equilibrium, the sum μ(CaO) + μ(CO 2 ) is fixed if carbonate is saturated.If μ(CaO) declines (i.e. higher compatibility of Ca in the rock and lower solubility in the fluid), then μ(CO 2 ) rises to offset this decline and maintain the equality.This would translate to a greater CO 2 solubility as long as the activity coefficient (nature and diversity of bonding environment in the fluid etc.) in the system remains the same.We observe that Δ s μCaO and Δ s μCO 2 are strictly anticorrelated (Fig. 9), which means aragonite is the main C host and it is always present in the sediments of Nicaragua.Moreover, we observe that the variations of Δ s μCaO track the main pulses of devolatilization event.At the onset of basalt devolatilization, Δ s μ(CaO) increases in the sediment, which reduces the solubility of C in the system compared to the non-speciated model, that is, Δ s μCO 2 declines.Conversely, when the fluid plume sourced in the AUM reaches the sediment at ca. 140 km depth, Δ s μCaO (and Δ s μCaSiO 3 ) drops, thus both Δ s μCO 2 and the solubility of C increases in the fluid.
Of course, when dealing with changes of a fluid-rock system in internal equilibrium, the chemical potentials of all components are tied to each other.Figure 9c illustrates that the metasomatic effect on carbon concentration in the fluid cannot be restricted to an interplay between CaO and CO 2 , but variations in the chemical potential of other components such as MgO and Al 2 O 3 , partly affected by metasomatic mass transfer, participate in the observed variability of Δ s μCO 2 (i.e.C solubility).As could be expected, it also illustrates that the quantitative effect of using a model where all components are partly mobile is most pronounced when fluids are fractionated from the rock, and this occurs at 80 km depth for our Model Nicaragua slab.The energetic deviations between for each component may reach several thousands of joules when elements are strongly soluble (K, Na).Interestingly, the deviation for silica in the Nicaragua sediment is minimal, meaning its metasomatic transport does not affect its energetics in the sediment system.

Return flux of H 2 O and comparison with ocean depth variations
Our predicted devolatilization events reflect the established activity/composition relationships of hydrous phases such as mica, lawsonite, serpentine and chlorite over a wide range of P and T (e.g.Connolly, 2005;Schmidt andPoli, 1998, 1994;Ulmer and Trommsdorff, 1995).The main pulse of slab devolatilization matches the magnetotelluric record of low mantle resistivity in the Central Cascadia subduction zone (Fig. 1).Slab thermo-petrological models may therefore be compatible with a flux-melting origin for the resistivity anomaly at the mantle source of Cascadia arc magmatism, even though water fluxes may be minimal.Schmidt and Poli (1998) Parai and Mukhopadhyay (2012).However, it may be relevant for specific slabs that are particularly cold, such as those of the western Pacific (e.g.Mariana, Honshu) (Cai et al., 2018).

Three regimes of carbonate recycling
Most models of carbonate subduction to date (Gorman et al., 2006;Kerrick and Connolly, 2001a, 2001b, 1998) support the relatively high stability of carbonates in the slabs.Decarbonation efficiencies range from 40% in Gorman et al. (2006) to 18-70% in Johnston et al. (2011), and this is qualitatively supported by field observations (Collins et al., 2015;Cook-Kollars et al., 2014).However, large uncertainties in slab parametrizations and carbon solubilities have led others to suggest that most slab C may be dissolved and released to the mantle wedge (Kelemen and Manning, 2015).Our work helps to clarify this conundrum.We confirm that C solubilities in metamorphic fluids may have been underestimated (Gorman et al., 2006;Kerrick and Connolly, 2001b, 2001a, 1998;Tian et al., 2019), and this discrepancy may reach at least 2 orders of magnitudes at T < 400 °C (see also Galvez et al. (2016)).However, fluid fluxes are vanishingly low at those temperatures, which explains why the discrepancy does not affect much the overall carbon budget.Hence, our results do bring support to the existing notion that C may bypass the sub-arc depth if one assumes it is part of a cohesive slab (Fig. 5a and 6d).Our total output from the slabs (ca. 25 to 31 Mt C/yr) is about half the global trench input (50 to 62% of C input), but is lower than field estimate of passive CO 2 degassing at volcanic arcs of 55 Mt/yr by Fischer et al. (2019).The former discrepancy is well within the estimates of decarbonation efficiencies (15-60 %) used in previous models of global carbon cycling (Hoareau et al., 2015;Volk, 1989).The latter discrepancy is not a surprise, and even reassuring as a fraction of arc carbon degassing may result from the thermal metamorphism and assimilation of carbonate accreted to the overriding crust (Johnston et al., 2011;Lee et al., 2013;Mason et al., 2017).Both imply that the slabmantle wedge-arcs system is imbalanced.
Regardless of our specific model parametrization, our analysis shows that the decarbonation efficiency of a slab is a variable which depends on its composition, hydration state, subduction rate and thermal structure.This means that using a given value through extended periods of time, e.g. when tectonic configurations, sediment composition and/or slab geometries may differ such as in pre-Mesozoic times (Galvez and Pubellier, 2019), would lead to potentially erroneous results.More generally, we have identified "transport" and "supply" limited regimes of element subductions, and this typology depends on which elements are considered.The classification of a subduction zone within one or the other regime for a given element depends on model parametrization (rock compositions, and thermal structures).For example, subduction zones of intermediate thermal structure are transport-limited, that is, sensitive to the hydration state of the AUM (Fig. 5).This regime includes most slabs of the Sumatra, eastern Pacific, Aegean (see Table 1 and Fig. 4).This is where most C is released (Costa Rica and Guatemala are the largest C emitter per km of trench), those slabs are also the most sensitive to the hydration state of the AUM (Fig 5).We therefore need a more accurate knowledge of the hydration state of those specific slabs.Subduction zones that are hot, such as Mexico and Cascadia, are "supply-limited", i.e. not sensitive to the hydration state of the AUM simply because that the slabs are depleted of carbon at shallow depth.In this case, priority should be given to improving knowledge of their carbon content and distribution.Most of the subduction zones of the western Pacific (e.g. Mariana, Izu, Tonga) also do not depend on the AUM hydration state.The absence of AUM dehydration before sub-arc depth prevents C dissolution from the slab.We call this sub-class of transport-limited regime "thermally-frozen".

Metasomatic redistribution of elements and arc magma signatures
We have provided insights into the subsolidus metasomatic redistribution of mass within a subducting slab.This ubiquitous natural process (Angiboust et al., 2014;Piccoli et al., 2016), to our knowledge, had never been addressed quantitatively at the scale of a subduction zone, and should be included in the design of future thermodynamic and thermomechanical models.Our observation that some slab sediments may get quantitatively depleted in some elements (e.g.K, Na), while other slabs' sediments may locally accumulate elements such as alkali, but also Al and Si reflects the contrasting chemical potentials of the oxide components between adjacent lithologies.For example, μAl 2 O 3 is initially lower in the carbonate sediments of Nicaragua than it is in the metabasalt, triggering Al accumulation in the sediment.We demonstrate in particular the importance of silicification and aluminification in driving the decarbonation of carbonate-rich lithologies, typically those of the Eastern Pacific.Field evidence for decarbonation driven by metasomatic mass transfer have long existed in contexts of contact (e.g.Bergell (Bucher, 1999), Adamello (Abart, 1995), and Beinn an Dubhaich (Ferry et al., 2011) aureoles) and regional metamorphism [e.g.Wepawaug schists (Ague, 2003)].We have found that the incongruent dissolution of carbonate in the sediment of Nicaragua (and also, e.g., Tonga, Costa Rica, Mariana) is associated with a synchronous gain in mass of Al 2 O 3 and SiO 2 in the sediment (Fig. 7), and with a modal increase of garnet (grossular).Noting the conspicuous absence of phengite from the paragenesis [as a potential solid source of Al (Galvez and Martinez, 2013;Malvoisin et al., 2012) at this stage points to a scenario where garnet (serving as Ca, Al and Si sink in the process) grows by infiltration of aluminous and siliceous fluids sourced in the dehydrating oceanic lithosphere where Al 2 O 3 (aq) and SiO 2 (aq) denote that Al and Si come, at least in part, from the infiltrating fluid.
Alternative phases serving as Ca, Al and Si traps during slab decarbonation would be zoisite/clinozoisite or calcic-pyroxene (diopside), depending on P, T and composition of the parent rock.
Because garnet is a major host for heavy rare earth elements (HREE) in the slab (Spandler and Pirard, 2013), this observation may imply that sedimentary lithologies that are initially poor in HREE may become important sinks for this suite of elements upon metasomatic growth of garnet.On the contrary, the abrupt destabilization of phengite (loss of K and H 2 O) in metacarbonate sediments would trigger a massive release of large ion lithophile elements (Cs, Rb, Ba and some mica hosted Pb and Sr) to the fluid, and contribute to their premature enrichment in the mantle wedge.More generally, the progressive loss of K from metacarbonate sediments, and enhanced molar K/Na ratio of our fluids with depth, reflects the experimental observations of the high incompatibility of K in high-P metasediments (Schmidt, 2015;Schmidt et al., 2004).
Our results may also have implications for the isotopic signature of the mantle wedge.When an element is partially dissolved from the slab and released to the mantle wedge, which we show is most of the time the case (e.g.Si, Al, Ca, Mg, Table 2, Figure 7), the isotopic signature of subduction fluid departs from that of the bulk slab, and becomes a complex function of fluid/rock ratio, mineral assemblage, temperature and kinetics of dissolution reactions.Assuming fluids tend to incorporate the lighter isotopes, our findings may help understand the cause for the light isotopic signature of Ca (Antonelli et al., 2021) or Mg (Shen et al., 2018) in the metasomatized mantle wedge.

Sensitivity and Thermomechanical limitations of model
The quantitative results of our model depend on our fluid composition and stoichiometry.In  2010), but the orders of magnitude are correct.Increasing divergence with T is caused by the progressive polymerization of alumino-silicate clusters (Manning et al., 2010;Newton and Manning, 2008).The recent experiments on aragonite solubility (Facq et al., 2014) and their integration within the thermodynamic database of Ca-electrolytes, ensures that the Ca solubilities in our system (where the C speciation is controlled by Ca-C complexes) are probably accurate.There is no equivalent data for Mg-carbonates and Mg-electrolytes, which means Mg solubilities and Mg fluxes out of slabs may be underestimated by about 2 orders of magnitude (Galvez et al., 2015).The solubility of Fe-bearing phases, because of their redox sensitivity, may be the least well constrained of all.Manning (2007) have shown that the solubility of corundum (Al 2 O 3 ) in fluid increases by about 2 orders of magnitude in the presence of quartz at 700 °C and 10 kbar.By analogy, we conjecture that Fe 3+ may be readily dissolved and transported as Si or Al-Si-complexes in metasomatic environments.As a divalent cation like Ca 2+ and Mg 2+ , Fe 2+ should readily form complex with dissolved silica and carbonate, for example, but none of those species are currently included in our model.Therefore, Fe concentration in modelled metamorphic fluids may be currently underestimated by about 2 orders of magnitude, and Fe return flux to the mantle overestimated.
As described previously (Connolly and Galvez, 2018;Galvez et al., 2015), our model loses accuracy in predicting fluid and rock compositions when the Gibbs energy of a fluid component represents a non-negligible part of the total Gibbs energy of this component in the fluid-rock system.
In our model setup, this exclusively occurs when a component concentration in the rock is close to zero and a large amount of fluid is involved.This situation arises for K 2 O when the AUM devolatilization plume reaches the K 2 O depleted Nicaragua sediment (Fig. 8).It also occurs for CO 2 when basaltic fluids hit the sedimentary cover of the Mexico slab, which is C-free.In those specific cases, minute amounts of carbonate or phengite may form, transiently, but the quantities involved are so low that they are beyond the reach of our model.This is inconsequential in our interpretations.This essentially suggests the components K 2 O and CO 2 , which are transported in soluble forms, are almost conservatively transferred from the basaltic layer to the top of the slab with minimal interaction with the sediment.Connolly and Galvez (2018) showed that cycles of fluid speciation by backcalculation (as opposed to one backcalculation only) at a given P and T ('lagged-speciation') should improve the resolution during speciation when the fluid-limited regime is reached.However, electrolytes are not included in the energy minimization step, in both variants of the algorithm, and both are therefore primarily applicable to rock-dominated systems.To be fully consistent, although excessively demanding computationally, both solutes and solvents must be part of the Gibbs energy minimization algorithm.The right balance between computational, geological and chemical accuracy and computational efficiency is determined by the nature of the problem at stake.
Our model captures key aspects of the thermodynamics of fluid-rock interactions in subduction zones, some of which had been elusive for decades, but the geometry and reactivity of a natural fluid flow may still differ from our predictions.For example, accounting for pressure gradients in nature (Wang, 2019;Wilson et al., 2014) would require more complex advection-dispersion designs (Ague, 2007).In addition, field constraints suggest channelized fluid flow should be explored as well (Ague and Baxter, 2007;Angiboust et al., 2014Angiboust et al., , 2012;;Herms et al., 2012;John et al., 2012;Plümper et al., 2017).In several instances, e.g.zones of high permeability contrasts (Breeding et al., 2004), lithological interfaces (Angiboust et al., 2014), shear zone (underpressure) (Mancktelow, 2006), (micro) fractures and connected porosity (wave) (Plümper et al., 2017;Tian et al., 2019), it may occur that only a fraction of the fluid reacts with the rock it traverses (as opposed to 100% in our model).We have tested this situation semi-quantitatively in the case of Cascadia, Nicaragua and Honshu (Supplementary Material Fig. S7).The results suggest that reducing the reactive volume by 50% tends to reduce the amount of carbon released from the slabs by ca.20% to 30% before subarc depth.
However, because fluids tend to dominate the rock mass as their reaction volume is diminished, testing this effect typically requires a computational strategy more adapted to fluid-dominated systems.Besides, Penniston-dorland et al., (2015) have recently proposed that slabs may be hotter than predicted by thermomechanical models.Our semi-quantitative tests (Supplementary Material Fig. S4) suggest that increasing the temperature by +50 °C and +100 °C increases the carbon release by ca. 10 % to 30 %.Although this is preliminary, this indicates that increasing temperature and channelization of fluids may have opposite effects.Finally, we model sediments as a compositionally homogeneous layer, which may not reflect the reality of heterogeneous slabs.Yet, this approximation is not likely to significantly alter our budget.These limitations rest on a key, and fundamental assumption: the mechanically cohesive nature of the slab.The geological record and recent theoretical studies suggest this is certainly not the case in general.First, the Jurassic/Cretaceous metasedimentary units (Liguro-Piemontese metasedimentary units) exposed of the Mediterranean region (Apennines, western and central Alps etc), for example, are built on top of Triassic evaporitic successions inherited from the early rifting stages.These weak evaporitic layers [e.g.'Nappe des Gypses' in the western Alps (Barré et al., 2020;Gabalda et al., 2009)] favour detachment and exhumation of overlying rocks (Malavieille and Ritz, 1989), which preclude their subduction.Examples of sediment underplating are ubiquitous, e.g. in the Cascades (Hacker et al., 2011;Matzel et al., 2004), in the Austroalpine domain (Sesia zone) (Manzotti et al., 2014;Vuichard and Ballivre, 1988), the southern Apennines (HP-LT carbonate-evaporite units of Lungro-Verbicaro (Molli et al., 2020)), and in the Cyclades (Ring and Layer, 2003).
These observations are supported by thermomechanical models by (Gerya and Meilick, 2011), and gravitational instability calculations by Behn et al. (2011), who show that subducted sediments detach from the down going slab at temperatures of 500-850 o C to form buoyant diapirs.Sediment may therefore experience much hotter thermal conditions as they rise through the mantle wedge, and their devolatilization may be driven by partial melting of buoyant diapirs (Behn et al., 2011;Poli, 2015), which is not included in our model.Clearly, better constraints on rheological coupling at the subduction interface could make a decisive contribution to understanding the fate of carbonate in subduction zones.

CONCLUSIONS
We have examined the process of slab decarbonation using a petrological-thermodynamical model where all components (volatiles and non-volatiles) are partially soluble in the fluid phase (Galvez et al., 2015).This provided a high-resolution insight into the reaction paths and mutual dependence of element transport and redistribution during slab dehydration.We have identified three regimes of element subduction: (1) transport (water)-limited, (2) supply-limited, and (3) thermally frozen slabs.
Conventional petrological modelling has assumed slabs as consisting of discrete, compositionally homogeneous layers distinguished by sharp compositional transitions.By contrast, our model leads to intriguing patterns where layer boundaries become diffuse once advection of mass (Al, Si, K, Na, Mg, Ca) is considered, and to significant deviation of bulk rock compositions caused by metasomatic mass redistribution.Our results reappraise the need to treat concentration data with care in field studies of metasomatism, as those are a poor indicator of intra-slab mass redistribution.In particular, we confirm the possibility of early slab depotassification and dehydration (Connolly and Galvez, 2018).The effect is pronounced for carbonate rich slabs, but some slab sediments are net accumulator of alkali, e.g.
Alaska, Chile.Similarly, silica and alumina may get trapped in the sediment layers on their way toward the mantle wedge.Calcsilicate formation by Al-Si metasomatism is an active driver, for example, of decarbonation in the Eastern Pacific slabs.This process is conspicuously marked by garnet (grossular) growth associated with fluid infiltration.Growth of metasomatic garnet by Al and Si-rich basaltic and ultramafic fluids suggests that slabs rich in limestones may acquire a high retention capacity for HREE during their subduction.In the current parametrization, our simulations confirm the notion that about half of C subducted should be release to the mantle wedge and ultimately to the atmosphere (eastern pacific).But this flux is lower than arc degassing, implying subduction zones are imbalanced.The decarbonation efficiency is variable across subduction zones, and markedly sensitive to sediment composition and thermal structures.
Although our equilibrium model predicts a high stability for carbonate along typical geotherms, this does not mean much carbon is mechanically returned to the transition zone.The chemical budget of our slabs depends on their mechanical behaviour, i.e., that they move as a single cohesive unit.This common assumption does not fully represent observational constraints from the Tethysian orogen (Barré et al., 2020;Malavieille and Ritz, 1989), as well as petrological-thermomechanical simulations (Gerya and Meilick, 2011), where sedimentary layers (pelagic or continental) have been massively detached from slabs along weak and ductile evaporitic layers, and underplated.

Fig. 1 :
Fig. 1: (a) Thermal structure of Cascadia subduction zone (redrawn from the results in Syracuse et al., 2010).Layers are composed of sediment, basalt, gabbro and peridotite.The Cartesian coordinate system is transformed into in-slab depth coordinate and slab surface depth coordinate for modelling purpose.(b) Schematic depiction of the modelling approach.The rock undergoes dissolution first during dehydration reaction and the volatile and dissolved non-volatile components are fractionated from the bulk rock.The fluid is then transported to the grid above in the in-slab coordinate.(c) Match between magnetotelluric data (McGary et al., 2014) and H 2 O flux based on 2D dissolution-fractionation-transport model for the Juan de Fuca plate below Cascadia.Three fluid pulses are labelled and enlarged.(d) Fluid flux (H 2 O) as a function of the surface depth.The fluxes at the surface of the sediment and Moho are shown with different colours.The subduction length is 850 km for the Juan de

Fig. 2 :
Fig. 2: Water content variations (H 2 O mass subtracted by the initial H 2 O mass divided by the total rock mass) and their associated equilibrium phase assemblages within a hot (Cascadia, Western U.S.), intermediate (Nicaragua, Central America) and cold (Honshu, Japan) subduction end-members.Red colour indicates net gain in mass due to transport of fluid from below (local composition higher than initial composition) and blue colour indicates net loss.The in-slab depth axis (y-axis) is slightly tilted to better visualize the subduction slab layering.
to 3j), Aegean, Costa Rica and Sumatra(Fig 5b).The integrated decarbonation flux out of the slab increases by 20 to 30% from LHS to HHS models.The decarbonation efficiency, η CO2 , is invariably close to 100 % for hot slabs (e.g.Cascadia, Mexico, Fig.4and Table2).Cold slabs (e.g.most slabs in the western Pacific such as Honshu, Ryukyu, Izu-Bonin, Tonga, Java, Fig.4), are also insensitive to hydration states.Those slabs release negligible amounts of H 2 O from the AUM before sub-arc depth (Fig.5c, 5d), regardless of the hydration model, and therefore conserve their carbon, e.g. as magnesite, to beyond sub-arc depth.This contrasted dependency on AUM hydration state is the basis of our classification of slabs into different regimes of carbonate subduction (cf.Fig.5and discussion).

Fig. 3 .Fig. 4
Fig. 3. Pattern of water and carbon release for warm, intermediate and cold slabs as a function of hydration states.From a to d, the H 2 O and CO 2 mass changes, respectively, are scaled by the initial input.The fraction of the cumulated CO 2 loss in the entire slab is shown in e.The same applies to Nicaragua (subplot f to j) and Honshu (subplot k to o).The warm colour indicates net gain (higher in volatile compared to initial composition)and cold colour indicates net loss.For Honshu slab, the first grid undergoes dehydration in the sediment already due to its high initial H 2 O content.For illustration purpose, we compare the mass variation with respect to the state after the fractionation of the first grid.These three slabs have been chosen because they represent three broad regimes of carbon subduction: i.e. 1) supply-limited slab, 2) transport-limited and 3) thermally-frozen, respectively.

Fig. 5 :
Fig. 5: a and c show the cumulated CO 2 and H 2 O output per km of slab from surface to 150 km depth.b and d show the flux difference between HHS and LHS model scaled by the input to visualize the dependence of C extraction due to AUM water.Slabs that are sensitive to the hydration state of AUM are labelled in black (waterlimited regime), slabs that are insensitive because they are fully decarbonated with LHS model are labelled in red (carbon-limited regime), slabs that are insensitive because they do not lose much C even with HHS model are labelled in blue (thermally-frozen regime).The ranking in d (H 2 O) follows that in b (CO 2 ), so that the three regimes are separated into three regions.For example, in Cascadia, although the water flux difference between HHS & LHS models is high as shown in d, the CO 2 flux difference is rather low as shown in b, indicating that carbon is depleted in this slab and extra water will not further increase CO 2 flux.

Fig. 6
Fig. 6. η 2 diagrams showing the decarbonation efficiency ηCO 2 versus dehydration efficiency ηH 2 O for Cascadia (a, b), South Chile (c), and the global subduction system (d) (results from model LHS).Curves, or portion of curves, with highly positive slopes indicate that most of the carbon is removed with minimal H2O release from the slab (i.e.high C solubility), while portions of curves with flat trends express the opposite.Canales et al. (2017) suggested that the H 2 O content in the upper mantle of the Cascadia subduction may be about 5 times lower than our LHS model for Central Cascadia, and 7.5 times lower for North Cascadia (British Columbia) due to limited plate bending and H 2 O penetration within the oceanic mantle.Panel b shows that even with an upper mantle containing 7.5 times lower H 2 O than LHS model in the case of the North Cascadia subduction zone, the system would still lose >90% of its C by subsolidus dissolution ("Dry" UM model).The LHS model is compared with the case where the AUM is 7.5 times drier than LHS model (red) and AUM is 1.5 times wetter than LHS model (blue).The results are similar, suggesting that the AUM hydration state is not important for the C extraction in Cascadia.The arrow in a and b indicates the location where mantle dehydration starts.d.A near-linear trend is observed for the global subduction system.The η 2 diagrams for the other slabs

Fig. 7 .
Fig. 7. Global ranking of the SiO 2 and Al 2 O 3 fluxes for subarc and forearc depth.

Figure 8
Figure 8 focuses on the mass and concentration changes occurring at the bottom of sediment layer for

Fig. 8 .
Fig. 8.The panels a, i and q show the flux of oxides out of the slab surface.The other panels show the concentration change and mass change of the oxides at the bottom of the sediment layer, where the metasomatic effect is strongest.Unlike concentration variation, the mass variation reflects the actual metasomatic effect.An absence of non-volatile element transport would induce no net mass change (flat curve).Concentration changes may always occur due to residual enrichment or dilution when volatiles are lost or added, respectively.Only the results of the LHS model are displayed.Note that the plotted scale is relative variation with respect to initial input of the component, not the absolute change with respect to the whole rock mass.

Fig. 9 .
Fig. 9. a, c.Chemical potential difference for various components between our speciation model (μ spec ) and a non-speciation model (only CO 2 and H 2 O are mobile, other oxides are infinitely compatible with rocks, μ nonspec ) at the bottom of the sediment layer in Nicaragua.The chemical potentials of CaO and CO 2 are anti-correlated, implying stability of aragonite.The chemical potential difference between those two models is a marker of metasomatism.b.Mass flux of dissolved carbon through the bottom of the sediment layer (normalized to CO 2 mass) and normalized fluxes of SiO 2 and Al 2 O 3 .d. Garnet, aragonite and phengite mode percentage.Note the correlation of garnet growth with fluid pulse.

(
Fig 9d) according to:3 Al 2 O 3 (aq) + 9 SiO 2 (aq) + 9 CaCO 3 = 3 Ca 3 Al 2 (SiO 4 ) 3 + 9 CO 2 (Eq.6) general, accurate element solubility predictions ensure accurate prediction of its mobility/compatibility relation in the (open) fluid/rock system.But fluid stoichiometry (i.e.element ratio in the fluid) turns out more critical because it determines the accuracy of mineralogical compositions i.e. (relative ratio of elements within a given phase, regardless of its modal abundance).Our predictions for the solubility of the assemblage albite-paragonite-quartz (Al, Si and Na concentrations in fluid, Supplementary Fig. S6) are slightly underestimated compared to the experimental data from Manning et al. ( , ca. 300 Mt H 2 O/yr.Our low hydration state model (LHS) reflects this consensus with an AUM input of 258 Mt H 2 O/yr, or 0.85 wt% of H 2 O over a 4 km AUM layer using the same linear variation scheme as in IHS and HHS models, for a total trench input of ca.1166 Mt H 2 O/yr (Table Figure6shows that the provenance of H 2 O (sediment/igneous layer and AUM) in mediating carbon dissolution varies widely between subduction zones.For example, AUM dehydration represents less than 20% of the H 2 O released from the warm slab of Northern Cascadia before sub-arc depth (LHS), but it contributes for 60% of the total C released (Fig6a and 6b).Conversely, AUM dehydration (30% slab H 2 O) contributes for less than 10% of the carbon released from the South Chile slab (Fig.6c).The rest of the C dissolution is controlled by water originating in the sediments and igneous layers: 60% of slab H 2 O (sediment + igneous layers) removes 90% of slab carbon for South Chile, while 70% of slab water (sediment + igneous) removes less than 40% of the Cascadian slab C load.Globally, only about 3.7% (ca.7 Mt/yr CO 2 ) of the total trench carbon input is released at forearc depth (Table2).It is mostly released in the hot subduction zones such as Cascadia and Mexico, and this flux is only driven by devolatilization of the igneous layer (Fig.2).The CO 2 extraction efficiency to subarc-depth is constrained between η CO2 subarc = 51.0%(LHS model, 94.8 Mt CO 2 /yr or 25.9 Mt C/yr) and η CO2 subarc = 61.9%(HHS model, 115 Mt CO 2 /yr or 31.4Mt C/yr), corresponding to σ CO2 subarc = 49.0 % (LHS) and σ CO2 subarc = 38.1 (HHS).Overall, AUM dehydration contributes for only ca.25% of global slab C loss before sub-arc (Fig. 6), the majority of the decarbonation before subarc depth in our model is driven by H 2 O originated from the sediment and igneous crust.
Parai and Mukhopadhyay (2012)bility of large fluid flux at forearc depth.We found that about 74 Mt H 2 O/yr based on the sediment composition fromPlank (2013)could be indeed lost before 20 km depth.Only our IHS and LHS models (Table2for LHS and S5 for IHS) yield subsolidus H 2 O releases to sub-arc depth that match observational constraints on global H 2 O degassing at volcanic arcs, i.e. ~650 Mt H 2 O/yr(Fischer, 2008).According toParai and Mukhopadhyay (2012), our return flux of H 2 O to the mantle for LHS (393 Mt H 2 O/yr) and IHS (571 Mt H 2 O/yr) models would imply ca. 100 m and 300 m sea-level decrease over the Phanerozoic, respectively (Table1).This is well within the acceptable range of variations permissible based on observational constraints of sea-level changes [360 m, (Parai and Mukhopadhyay, 2012)].The global return flux of the HHS model (1195 Mt H 2 O/yr) would imply excessive sea-level variations according to