The Subducting Slab as a Chromatographic Column: Regimes of Sub-Solidus Mass Transport as a Function of Lithospheric Hydration State, With Special Reference to the Fate of Carbonate

, the


Introduction
When oceanic lithospheres sink into the mantle in subduction zones, a fraction of their mass dissolves into mobile fluids that are released at increasingly higher pressure and temperature (P-T) conditions (Frezzotti & Ferrando, 2015;Hermann et al., 2013;Manning, 2004;Peacock, 1990;Schmidt & Poli, 1998;Stolper & Newman, 1994;Tatsumi et al., 1986).The chemical budget of a down-going plate and of the mantle wedge depends on the temporal and spatial sequence of this mass transfer (Kerrick & Connolly, 2001;Manning, 2004;Stolper & Newman, 1994).This applies particularly to the fate of carbonate (Ague & Nicolescu, 2014;Caciagli & Manning, 2003;Collins et al., 2015;Connolly, 2005;Gorman et al., 2006;Kerrick & Connolly, 2001a;Menzel et al., 2019;Sverjensky et al., 2014).But in systems composed of complex fluid and solid solutions, the fate of all Abstract We investigate the chemical budget of subduction zones at sub-solidus conditions using a thermodynamic-numerical simulation in which all major rock components are treated as soluble and potentially mobile in aqueous fluids.This new strategy significantly improves the accuracy of predicted fluid-rock equilibrium compositions in open petrological systems.We show that all slabs release volatiles and nonvolatiles to the mantle wedge, contributing to its refertilization.But some mobile constituents, such as alkali and alumina, may be trapped along layer boundaries or traverse without interaction depending on chemical contrasts between adjacent lithologies.The accumulation of igneous alumina and silica in the limestones of the centraleastern Pacific slabs drives their decarbonation and is marked by metasomatic garnet growth.Those slabs are also predicted to lose much of their alkalis before sub-arc depth.Even when they are produced in the altered mafic and ultramafic layers, fluids reach the slab/mantle wedge interface with distinct compositional signatures that are typical of the sedimentary cover.We distinguished supply and transport limited regimes of element subduction by testing the sensitivity of our mass balance to changes in slab hydration state (HS).Transport limited slabs sensitive to HS include notably a hotspot of carbon release to the mantle wedge (e.g., Costa Rica).Finally, we show that the quantitative budgets do depend on the geometry of fluid flows, and on assuming that slabs are mechanically continuous structures, which is questionable.Taken together, these insights will help better constrain the long-term chemical evolution of the shallow planetary interior, and the thermomechanical behavior of the subduction interface.
Although it is not shown, Equation 1 is catalyzed by water at typical sub-solidus subduction zone conditions (Aranovich & Newton, 1999;Kerrick & Connolly, 2001a).If the silica constituent is assumed to be only a part of the rock, the reaction may either proceed from P-T driven (i.e., prograde) changes in the chemical potentials (μ) of CaCO 3 , CaSiO 3 , and/or SiO 2 (Aranovich & Newton, 1999;Harker & Tuttle, 1956), or via fluid-driven modification of μ(CO 2 ) by metasomatism, which drives Equation 1 to the right or left (Ferry & Dipple, 1991;Rumble, 1982).But real subduction fluids contain dissolved solutes (Frezzotti & Ferrando, 2015;Manning, 1998;Philippot & Selverstone, 1991;Rustioni et al., 2021;Schmidt & Poli, 2003;Schneider & Eggler, 1986;Spandler et al., 2007) and those fluids are therefore neither strictly binary, nor purely molecular.Hence, changes in chemical potentials (μ) of the constituents in Equation 1 may also result from diffusion or advection of SiO 2 , CaO, MgO and/or other major rock components in and out of the system (Abart, 1995;Ague, 2003;Ferry et al., 2011).This can be expressed by rewriting Equation 1 as: (2) where SiO 2 (aq) denotes that silica may also be part of a reactive fluid phase.Accounting quantitatively for this mutual dependence of element mobility has long remained elusive in the numerical modeling of subduction zones processes.This is in large part due to the lack of computational tools to connect the mobility of volatile and nonvolatiles at a more fundamental level.
In this context, Sverjensky et al. (2014) have called attention to the capacity of the HKF equations of state (Helgeson et al., 1981), as implemented within the Eq3/6 software package (Facq et al., 2014(Facq et al., , 2016;;Wolery, 1992), to deal with solvent-electrolyte solutions owing to the knowledge of water dielectric constant over an extended range of P-T conditions (Fernández et al., 1995(Fernández et al., , 1997;;Pan et al., 2013).But extensions of electrolyte chemistry to high-P/T petrological systems must account for two of their most defining characteristics: the existence of mixed volatile fluids characterized by variable water activities, and the prevalence of multicomponent minerals and solid-solutions at those conditions.Galvez et al., (2015) introduced a simple approach to achieve that 10.1029/2021JB023851 3 of 28 goal.It takes advantage of the precision of phase equilibria calculations by Gibbs energy minimization (e.g., Connolly, 1990) to compute the speciation and elemental stoechiometry of mixed-volatile high-P fluids using a limited set of chemical potentials, regardless of mineralogical complexity (Galvez et al., 2015(Galvez et al., , 2016)).Their metric of relative acido-basicity ΔpH reflects the fundamentally hybrid character of high-P/T fluids, simultaneously electrolytic and molecular.But the chief advantage of this approach is that it relaxes the "infinite compatibility assumption" of conventional petrological models.The integration of this new method within the PerpleX software by Connolly and Galvez (2018) provided a flexible numerical plateform for applications related to lithospheric devolatilization/revolatilization (Menzel et al., 2020;Piccoli et al., 2019) or the nature of fluid-rock interactions on Enceladus (Daswani et al., 2021).Connolly and Galvez (2018), in particular, reported that the mobility of potassium in fluids released from the model GLOSS sediment (Plank & Langmuir, 1998) may lead to its early depotassification, and therefore dehydration.But this result may have been amplified by infiltrating the sediment with initially pure H 2 O, which calls for a systematic approach of sub-solidus mass transfers where all lithologies, sediment, mafic, and ultramafic would be partly soluble in the fluid phase.
Building on these premises, we report the first attempt at a quantitative chemical budget of subduction zones using a thermopetrological numerical method in which all rock components have a finite and experimentally constrained solubility in the fluid phase.We pay special attention to the process of slab decarbonation.Because subduction zones differ in composition and thermal structure, we run each individually when relevant data are available.For the sake of simplicity, we set up a model which considers dissolution, fractionation, and vertical transport of aqueous fluids and their dissolved solutes within, and out of, subducting slabs.We test the sensitivity of our chemical budgets quantitatively by adopting a range of upper mantle hydration states.Our model relies on a set of geological, chemical, geodynamical, and computational assumptions that are discussed and tested quantitatively when possible.The concepts and methods discussed here are a step forward in understanding the role of fluid-rock interactions on the composition and mechanical behavior of the Earth's lithosphere, and by extension the chemical dynamics of other rocky planets as well.

Chemical Structure of the Subducting Slab: Model Setup
Our slabs are 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 of the pile (see an example in Figure 1).The composition of the sediment is from Plank (2013), including H 2 O and CO 2 , with a few minor modifications to ensure consistent Ca/C ratio when C abundances are lacking, or particularly uncertain (Figure S1 in Supporting Information S1).The altered basalt (1.5 km thick) and gabbro (5.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 & Honnorez, 1984;Staudigel et al., 1989).To simulate this process, we considered 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 in Supporting Information 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 & Teagle, 2003).The altered upper mantle (AUM) composition is from Hart and Zindler (1986), and its thickness and water content is a variable and discussed below.Kelemen and Manning (2015), that is, 147-242 Mt CO 2 /yr.The Makran subduction zone is not included in our simulations because no thermal model was available (Syracuse et al., 2010).Results of a model with higher C contents in the slab (ca.267 Mt CO 2 /yr from Clift, 2017) are also presented for comparison in the supplementary materials.
Because H 2 O is the main solvent in our model, we have adopted a range of hydration states for our slabs.Our trench sedimentary H 2 O composition is from Plank (2013) (Jarrard, 2003), the sedimentary and igneous flux in our models is 833 Mt H 2 O/yr (Table 2).This is consistent with the general The main uncertainty is associated with the serpentinized AUM on a global scale.In detail, Schmidt and Poli (2003) 1 and 2), that is, 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) 1 and 2).Our global estimates may locally overestimate (Canales et al., 2017) or underestimate (Cai As a first-order approximation, pressure within the subducting slab is lithostatic using an average rock density of 3,300 kg/m 3 , and slab thermal structures are from the numerical results of Syracuse et al. (2010) using input data of plate convergence rate and slab geometry.Slab names, length, and sediment thicknesses are taken from Clift (2017).Most of our article refers to the extraction efficiency of a given oxide i, it is denoted by η i .This variable corresponds to the fraction of the trench input of component i released from the slab within a certain depth interval.The subduction efficiency σ = 1−η, describes the fraction of the total input retained within the slab to a given depth (Galvez & Pubellier, 2019;Volk, 1989).The extraction and subduction efficiency of carbon is a key parameter of long-term climate and biogeochemical models (Caldeira, 1991;Hoareau et al., 2015;Volk, 1989).
For our purpose, the forearc region starts at 20 km, the depth of the first grid of our model.It ends at 80 km, unless stated otherwise in our sensitivity tests.The "sub-arc region" includes the forearc and ends at a depth of 150 km.

Computational Strategy
Simulating the processes of internal thermodynamic equilibration between fluids and solids by Gibbs energy minimization requires that the system is closed, that is, 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 nonvolatile) 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.The main novelty here is that all chemical components have a finite solubility in the fluid phase, with mineral solubility relationships conforming to decades of experimental data (see thermodynamic database and comparison with experimental results in supplementary materials).Therefore, while being computationally   S5 and S6 in Supporting Information S1.We provide budgets at fore-arc (<80 km depth) and sub-arc (<150 km depth), and they are expressed as a function of trench inputs (e.g., %).Slab name and length are based on Clift (2017).Table 2 Global Subduction Input and Extraction Efficiency (Flux/Input in %, ‰, ‱) of Major Elements at Fore-Arc and Sub-Arc Depth efficient, this approach offers greater geological consistency and accuracy in modeling the composition of rocks and fluids in open petrological systems.
We ran a high-resolution variant of the original backcalculation algorithm of Galvez et al. (2015).The core component of all variants of this algorithm, that is, simple- (Galvez et al., 2015), lagged- (Connolly & Galvez, 2018), or micro-backcalculation (this study), relies on finding first the Gibbs energy minimum (equilibrium state at a given P and T) of a system composed of rock and solvent using some appropriate Gibbs energy minimization (GEM) technique, for example the Perplex software here (Connolly, 1990;Galvez et al., 2015).In a second stage, a set of independent chemical potentials derived from this GEM phase is used to estimate (backcalculate) the electrolyte speciation (cf.Galvez et al., 2015 for details).The reason why this basic strategy is so precise in most situations (Connolly & Galvez, 2018) is because the rock and solvents represent the dominant mass (and Gibbs energy) of high-P systems.This is the general rock-dominated criterion which is always realized here.It ensures that the predicted fluid and solid phase assemblage is very close to the Gibbs energy minimum even though the electrolytes are not part of the GEM procedure.A more robust, specific, rock-dominated criterion would require that the mass of every single system's component (e.g., SiO 2 , CaO, K 2 O, and CO 2 etc.) is dominantly partitioned into the solid phase.Violation of this criterion does occur occasionally, but the impact on our conclusions is marginal.Indeed, without getting rid of both rock-dominated constraints, here we significantly enlarge their domain by equilibrating aliquots of the infiltrating fluid and its solutes sequentially.This ensures minimal fluid (solute)-rock ratio at each iteration, which is nearly equivalent to assuming a constantly low porosity regardless of integrated fluid flux.This is geologically realistic at high P.We then take advantage of the fact that cycles of fluid speciation and fractionation conducted in the rock-dominated limit may lead to profound, and so far overlooked, changes in rock compositions.
In detail, we first discretize a column of rock into grids (resolution provided below), and we assign the bulk composition to each grid according to its position within the four distinct horizontal layers (see Supporting Information S1).We then proceed in three steps: rock-solvent equilibration (by GEM), fluid speciation, and fluid fractionation.In the first step, we compute the chemical potential of the rock and solvent components by equilibration of the bottom grid in the column of rock using the Perple_X software (Connolly, 2005).Table S4 in Supporting Information S1 lists the solid solution models used for this GEM step where melts are not considered.In the second step, the chemical potentials obtained by GEM are used to backcalculate the fluid speciation using a non-platform specific code, which we run on Mathlab® (cf.Galvez et al., 2015).We provide details on the protocol and activity models in the Supplementary Information S1.In the third step, the fluid (solvent and solutes) is fractionated from the bulk composition of the local grid, and added to the bulk composition of the grid immediately above before repeating the process.For enhanced thermodynamic consistency and chemical accuracy, we estimate fluid speciations in high-resolution mode (i.e., by microbackcalculation) where aliquots of fluids (as opposed to the bulk fluid) are equilibrated, speciated, and fractionated successively.Discretization of the fluid volume (set here to n = 10 sub-volumes) can be made as dense as desired, at the cost of computation time.This mode converges to the rock/fluid composition asymptotically, unless a rock component is entirely dissolved (i.e., fluid-dominated system), that is, when the amount of a dissolved component is predicted to be higher in the fluid than the leftover in the rock.In this case, only the leftover mass of this specific rock component is allowed to dissolve in the fluid.This is done under mass and charge balance constraints, that is, the dissolved component traverses conservatively the rock grid in which it is temporarily out-of-equilibrium.This minor rule ensures mass conservation during the entire transport and subduction sequence.
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.Although Mottl et al. (2004) supported this model at Mariana fore-arc, other studies (e.g., Angiboust et al., 2014;Galvez et al., 2013;Kawano et al., 2011) show that interface-parallel fluid flux is important in subduction settings, suggesting our model may overestimate vertical fluxes.The dissolution-fractionation-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 AUM layers.Subsequently, the entire column of rock is moved downward according to the thermal model considered.The descending rate of slabs sets the subducted mass flux of rock columns (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

Devolatilization Pattern of Cascadia Versus Magnetotelluric Data
Figures 1c and 1d show 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 (cf. Figure S2 in Supporting Information S1 for details on fluid flux and composition)-and the depth of a low resistivity area for the Cascadia subduction zone (McGary et al., 2014).The water plume linked to upper mantle dehydration matches an extensive zone of low resistivity located between 80 and 110 km depth.A minor fluid plume, sourced 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 (Figure 2 and Figure S2 in Supporting Information S1).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 (Figure 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, Figure 2a).Eclogitization of the sediment and basalt releases a free fluid phase between 40 and 60 km.The entire Nicaragua slab is dehydrated by 170 km depth, and the entire Cascadia slab is dehydrated by 100 km depth.By contrast, only the top 500-1,000 m of the cold Honshu oceanic mantle (antigorite) may dehydrate at a depth of about 200 km.The water released from the AUM is transported to the igneous layer.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 (Figure 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 1, and also Table 2 for LHS, Table S5 in Supporting Information S1 for IHS, and Table S6 in Supporting Information S1 for HHS models).By contrast, the water released within the 80-150 km depth shows marked sensitivity to the hydration state of the AUM (Table 1).Between 561 Mt/yr (η H2O subarc = 48.1% of trench input, LHS, Table 1) and 1,015 Mt/y (η H2O subarc = 41.8% of trench input, HHS) is released in the depth interval.

Interslab Variability of Element Release Out of Slab: Carbon
Carbon is present as aragonite, dolomite, and magnesite in the sediment layer of our model (Figure 2).Figures 3  and 4 show that carbon is released from the rocks in pulses, regardless of the hydration state.For all the slabs, the location of carbon pulses is systematically associated with dehydration events (e.g., Figures 2 and 3).For example, the carbon load from the Cascadia sedimentary layer is released within a narrow depth range (60-80 km, Figures 3a-3d) associated with eclogitization of the sediment and igneous layers and infiltration of H 2 O from the AUM flushing C from reaction sites (Figures 3a-3d).
The AUM hydration state does not modify the pulsed nature of carbon release (Figures 3 and 4).However, key quantitative aspects of this release, for example, its rate and magnitude are sensitive to the hydration state (Figures 5 and 6b).For example, the integrated decarbonation flux out of the Nicaragua slab increases by 20%-30% from LHS to HHS models (Figures 4b and 6b).Because the carbon output is sensitive to the fluid expulsion (transport) from the AUM, we classify this regime as transport (water) limited (Figure 6 and Figure S4 in Supporting Information S1).This is not the case for the hot slabs such as Cascadia and Mexico, where the 10.1029/2021JB023851 9 of 28 decarbonation efficiency (Figures 4 and 5 and Table 2) is invariably close to 100% regardless of the hydration state.Because the carbon output is controlled by the initial carbon supply, we classify those slabs as supply limited (Figures 6 and S4 in Supporting Information S1).Interestingly, cold slabs (e.g., most slabs in the western Red color indicates net gain in mass due to metasomatism and blue color indicates net loss.The in-slab depth axis (y-axis) is slightly tilted to better visualize the subduction slab stratification.The initial water input in the mantle is assumed to vary linearly along the in-slab direction, which explains the drier phase assemblages in the lower part of the ultramafic section (cf.Figures 3a-3c).The thickness of the sediment is plotted as 1 km for Nicaragua and Honshu so that the phase assemblages can be seen.Detailed flux of water and composition of fluid along the slab surface (Figure S2 in Supporting Information S1) and composition of fluid along selected vertical transects of the three systems (Figure S3 in Supporting Information S1) are provided in the supplementary material.Pacific such as Honshu, Ryukyu, Izu-Bonin, Tonga, and Java, Figure 5), are also insensitive to hydration states, but for different reasons.Figures 6c  and 6d show that those slabs release negligible amounts of H 2 O from the AUM before sub-arc depth, regardless of the hydration model, and therefore retain a large fraction of their carbon load to beyond sub-arc depth as magnesite.The feature is therefore primarily caused by the hydrodynamic regime of the slab, that is, the lack of water to mediate the dissolution, and only proximally to the slab's thermal structure.We call this regime thermally frozen.
Overall, AUM dehydration contributes only ca.25% of global slab C loss before sub-arc (Figure 7d), and therefore the majority of the subarc decarbonation is predicted to be driven by H 2 O originated from the sediment and igneous crust.About 40% of the sedimentary and igneous carbon remains after sub-arc depth, while about 80% of the ultramafic C remains (Figure S5 in Supporting Information S1).In other words, our current model design predicts inefficient slab decarbonation, and this inefficiency affects all layers, including the most heavily flushed sediments and igneous layers.Of course, AUM dehydration does promote continued dissolution of the slabs (Figures 3 and 7d) beyond 150 km depths.This flux may flow toward the back-arc region, unless it is channeled back toward the sub-arc laterally.

Interslab Variability of Element Release Out of Slab: Nonvolatile Elements
Our qualitative findings above apply to all other rock forming elements.Most notably, the release of K, Na, Ca, Mg, Al, and Si occurs in discrete pulses that are controlled by the sequence of dehydration reactions within the slabs (Figure 2).The composition of fluids available for mantle wedge metasomatism is generally rich in Na, Si, C, and Ca (ca.0.1 to 1 molal beyond 80 km, Figures 3a and 3i), poor in K until ca.90 km depth, reflecting the stability of white mica.They differ greatly by their Al content (Figure S2 in Supporting Information S1).It is highest at high-T (which varies from slab to slabs) where the Al reaches the top of the slab in the form of soluble alumino-silicate complex, and lowest at lower T where Al is highly insoluble.Unlike C and Si, the hotspots of Al release to sub-arc depth are the trenches of Mexico, South West Japan, or Alaska, broadly reflecting their compositional enrichment in Al.At higher P and T, the metasomatic fluids become more potassic while the slabs, taken as whole, get progressively depleted of this element (Figures 3 and 9, Figure S2 in Supporting Information S1).An interesting feature of the Nicaragua sediment is that it is K free after 90 km depth (Figure 3), and yet the fluid  crossing it and reaching the mantle does contain K inherited from the mafic (not sedimentary) layers (Figures S2 and S3 in Supporting Information S1).This is a peculiar case of partial equilibrium between the fluid and the rock.For all the subduction zones studied, the fluids that reach the top of the slab are much richer in Al, Si, C and much poorer in Mg and Fe than they are at the site of their production (in the mafic and ultramafic sections, Figure S3 in Supporting Information S1).
The sensitivity of mass fluxes to slab hydration state differs between oxides and between slabs (Figures 3, 8 and 9) For example, the hot slabs of Central Cascadia and Mexico were supply limited for C because this element was in sufficiently low concentration to be completely removed from the sediment by sub-arc depth (Figure 6).However, the extent of their desilicification and dealuminification (rich in SiO 2 , Al 2 O 3 ) does prove sensitive to the AUM hydration state, which makes them transport limited for those components.By contrast, the desilicification and dealuminification of slabs such as Costa Rica and Guatemala depends heavily on AUM hydration (transport limited), just as their decarbonation did (Figures 6 and 8).A last salient feature is 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, Figure 9n).However, the release of K from the slab as a whole shows some dependence to hydration state because K is also distributed within the igneous layer.

Intra-Slab Redistribution of Elements
Our model also provides insights into patterns of element redistribution within the slab (intra-slab metasomatism), which contributes to significantly modify the rock compositions.However, the concentrations of components in a rock are by themselves a poor indicator of metasomatic transport (Aitchison, 1982;Gresens, 1967).This is because the concentration changes of a component (Δc i ) hinges on the behavior of the other rock components as well.For example, the concentration of a nonvolatile component may be residually enriched when H 2 O and CO 2 are massively lost from the system.Here, we use the relative change in mass of a component scaled by its initial input Δm i as a primary metric to assess the metasomatic redistribution of nonvolatile components throughout the slab (Figure 9).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 (Figure 3).Another salient feature is that element concentration (Δc i ) and mass (Δm i ) variations may vary in opposite directions (Figure 9, and exhaustive data set in Figures S6-S12 in Supporting Information S1).
For example, Δc SiO2 (increases) and Δm SiO2 (declines) are anticorrelated in the bottom of the sedimentary layers of Nicaragua, Cascadia and Honshu, which are affected by large volatile loss.Similar anticorrelated trends for Al 2 O 3 (most slabs, Figure S7 in Supporting Information S1) and for K 2 O variations in Honshu (NE Japan), Aegean and North Sumatra (Figure S10 in Supporting Information S1) are caused by net loss of H 2 O and CO 2 from the rock, which counters the net gain in mass of the nonvolatile component.This primarily reflects residual change (Figure 9).Conversely, Δm and Δc may covary perfectly in the case of H 2 O, CO 2 , Al 2 O 3 , and K 2 O.For example, most carbonate-rich sediments (e.g., Nicaragua, Costa Rica, Peru, Columbia, and Guatemala, cf.Supporting Information S1 for exhaustive data set in Figure S10) lose all of their initially low K 2 O (Δm and Δc decline), whereas they accumulate Al and Si in all or parts of their thickness (Δm and Δc rise).This feature is mostly associated with infiltration of fluids produced in the dehydrating AUM, which reach the top of the slab at higher temperatures.In this case, dissolved Al promotes the growth of grossular garnet in the sediments (details in discussion).Other slab sediments, such as those of North Chile (Δm K2O = +20%), Kurile, and Sumatra slightly increase their K mass with depth (Figure S10 in Supporting Information S1), while those richer in Si and Al display net Si and Al mass loss (Figures S6 and S7 in Supporting Information S1).The link between composition and patterns of mass variations by metasomatism is also observed for Ca.Ca-deficient sediments accumulate Ca, and the effect is dramatic for Mexico (+120%), Ryukyu (+180%), and Kurile (200%) (Figure S8 in Supporting Information S1).By contrast, the lower section of limestone-subducting slabs loses about 50% or more of its Ca content by 150 km depth (Figure 9).This effect is not observed at the top of this sedimentary layer.This effect may be caused by insufficient H 2 O infiltration.
We also found that patterns of mass changes may not be monotonous, and the profile of CaO in the Honshu sediment is a striking example.After the carbonate (dolomite and aragonite) is fully dissolved from the sediment by 110 km depth, a process accompanied by net Ca loss (Figure 9t), Ca accumulates in the carbonate-free sediment and is marked by metasomatic garnet growth.Moreover, most slab sediments lose Na due to the Na-deficient character of our ultramafic or igneous fluids, with Δm Na2O variations ranging from a few percent (e.g., Aegean) to almost 100% (e.g., Ryukyu, Mexico) at subarc depth.There are a few exceptions, but in those cases the accumulations never exceed a few relative percent at most (Figure S9 in Supporting Information S1).Mg only accumulates at the base of the gabbro layer after dehydration of the upper mantle (less than a few relative %, Figure S18 in Supporting Information S1 for Nicaragua), and to a lesser extent at the bottom of Mg-depleted sediments.We expect this effect to be more pronounced in case of direct ultramafic/sediment interfaces, which are not simulated here.This intermediate storage tends to lock Mg within the slab, instead of releasing it to the mantle wedge.

Metasomatic Control on Carbonate Dissolution
Metasomatic mass redistribution affects the release of C from the slab (Figure 9).This effect is monitored by computing the system's chemical potential difference (Δ s μ) between two models, one where all components 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 a carbonate (e.g., CaCO 3 ) in a fluid can be described by an equilibrium relating the Ca-carbonate constituent to its thermodynamic components (e.g., CaO and CO 2 ), regardless of their hosts (fluid or solid) and of the silicate environment: (3) Under thermodynamic equilibrium of a closed system, Equation 1gives: which also implies that the two potentials μ(CaO) and μ(CO 2 ) cannot vary independently if a pure CaCO 3 endmember (i.e., aragonite) is stable, regardless of the system's composition.For example, if μ(CaO) declines (i.e., higher compatibility of CaO in the rock and lower solubility in the fluid), μ(CO 2 ) rises to offset this decline and maintain the equality, which practically means greater CO 2 concentration in the fluid if the activity coefficient remains similar.We do observe that Δ s μCaO and Δ s μCO 2 are strictly anticorrelated along the subduction path (Figure 10a), indicating that aragonite is indeed always stable in the sedimentary mineral assemblage of Nicaragua (see also Figure 2).And the marked deviation of Δ s μCO 2 from the non-speciated reference frame (Δ s μCO 2 = 0) is the most direct marker of non-volatile metasomatism-the distinctive feature of our model.As could be expected from a metasomatic effect, the fluctuations of Δ s μCaO track the two main fluid infiltration events into the sediment, at 80 and 140 km (Figures 2b and 10b).These events contribute not only mafic/ultramafic CO 2 and H 2 O (Figure 2b), but also CaO, Al 2 O 3 and SiO 2 into the system (Figure 10b, see also Figures 9i-9p).The energetic deviations between the speciation and non-speciation model for each component may sometimes reach several thousands of joules, especially when elements are particularly soluble (e.g.K, Na) or particularly tieed to the fate of the alkalis (e.g.Al).The silica component is an exception, however, as the deviation of Δ s μSiO 2 for all sediments is generally minimal.This means that the net silicification (or desilicification, N/S Chile, SW Japan) of sediments (Figures 9j and S6 in Supporting Information S1) is insufficient to affect the overall energetics of the SiO 2 component in those systems.

Return Flux of H 2 O and Comparison With Ocean Depth Variations
Our predicted devolatilization events for the Central Cascadia subduction zone match the magnetotelluric record of low mantle resistivity (Figure 1).This supports a flux-melting origin for the resistivity anomaly at the mantle source of Cascadia arc magmatism, even though the water fluxes involved may be limited (Canales et al., 2017).
The limited fore-arc devolatilization that we predict globally, and in the Kamchatka subduction zone in particular, is consistent with findings by Konrad-Schmolke et al. (2016).However, our results also show fore-arc estimates are surprisingly sensitive to model design.For example, shifting the fore-arc/subarc threshold to 90 km depth doubles the predicted forearc H 2 O flux (Figure S13 in Supporting Information S1).This could typically account for interface-parallel upward fluid transfer, in which a fraction of the fluid pulse released from the basaltic layer between 80 and 90 km would be channelized to the fore-arc by migration within the slab, or along the subduction interface (Angiboust, Glodny, et al., 2021;John et al., 2012;Kawano et al., 2011).Only our IHS and LHS models (Table 1) yield sub-solidus H 2 O releases to sub-arc depth that match observational constraints on global H 2 O degassing at volcanic arcs, that is, ∼650 Mt H 2 O/yr (Fischer, 2008).Our corresponding return flux of H 2 O to the mantle (393 Mt H 2 O/yr for LHS, and 571 Mt H 2 O/yr for HHS) imply ca. 100 and 300 m sea-level decrease over the Phanerozoic, respectively (Table 1).This is within the acceptable range of sea-level variations permissible based on observational constraints, that is, 0-360 m (Parai & Mukhopadhyay, 2012).The global return flux of our HHS model (1,195 Mt H 2 O/yr) would imply excessive sea-level variations meaning it is not realistic globally.But it may be relevant regionally when slabs are particularly cold (e.g., Mariana, Honshu; Cai et al., 2018).The impact on C and nonvolatile element fluxes would be sensible only if the AUM would dehydrate to some significant extent before sub-arc depth, which is not observed here (Figure 2c).But the contention by Penniston-Dorland et al. ( 2015) that the thermal models of Syracuse et al. (2010) do not reflect accurately the metamorphic record leaves room for debate.It is important to note, however, that the metamorphic record studied by the former includes rocks detached from the slab, underplated and exhumed, i.e. geodynamic contexts quite different from the subduction zone structures simulated by Syracuse et al. (2010).Both studies may therefore be relevant in different geodynamic contexts.In the present context of a mechanically continuous slab, our work suggests that the low upper-mantle hydration states of Van Keken (2011) or Schmidt and Poli (2003), and even the intermediate hydration state of Hacker (2008), are all compatible with the constraint of limited sealevel variability in the phanerozoic if one considers that sub-arc depths extend to 150 km at least (as opposed to ca. 100 km in Parai & Mukhopadhyay, 2012).

Regimes of Carbonate Subduction
Another implication of our work concerns the fate of carbonate.Most models of carbonate subduction to date (Gorman et al., 2006;Kerrick and Connolly, 1998, 2001a, 2001b) support their relatively high thermodynamic stability in subducting slabs.Predicted 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;Lefeuvre et al., 2020).But large uncertainties in slab parameterizations and carbon solubilities have led others to suggest that most slab C may be dissolved, released, and stored temporarily in the fore-arc and mantle wedge (e.g., Barry et al., 2019;Kelemen & Manning, 2015).Our work helps to clarify this conundrum.We confirm that C solubilities in metamorphic fluids may have been locally underestimated in previous models of slab decarbonation (Gorce et al., 2019;Gorman et al., 2006;Kerrick and Connolly, 1998, 2001a, 2001b;Tian et al., 2019).By comparing the C concentration and Ca/C ratio of fluids predicted by speciated (hybrid electrolytic/molecular) and non-speciated (molecular) models, Galvez et al. (2016) have shown that this discrepancy in C concentration may reach at least 2 orders of magnitudes at T < 400-500°C in pelitic or basaltic lithologies because the dissolution of carbonate may be mostly congruent at those temperatures (Caciagli & Manning, 2003;Fein & Walther, 1987).However, fluid fluxes are also vanishingly low at those temperatures (Figures 1d and 2) such that the effect of enhanced CO 2 solubility is countered and turns out inconsequential in the predicted slab carbon budget, by 150 km depth at least.In other words, if one assumes carbon is part of a cohesive slab, and fluxes are dominantly vertical, then our results do support the existing notion that about half of subducted C (ca. 25 to 31 Mt C/yr, 50%-62% of C input) may be released before sub-arc depth (Caldeira, 1991;Gorman et al., 2006;Hoareau et al., 2015;Volk, 1989), and the other half subducted beyond (Figures 6a and 7d).
But this does not entirely solve the above conundrum as we noted that our fore-arc H 2 O and mass fluxes may be underestimated by a factor of 4 at least if the hydrodynamic regime is only slightly modified, ie.if lateral transport of fluid over 10 km is considered.Globally, this is a shift from 7 to 32 Mt CO 2 /yr at fore-arc, or 17% of the input (Figure S13 in Supporting Information S1).This sensitivity may be quite dramatic locally depending on slab composition and temperature profile.For example, our predicted fore-arc C flux in Costa Rica is negligible, 0.0015 Mt CO 2 /yr (Figure 6a and Table 1), which matches the lower-end isotope-based estimate of Barry et al. (2019), who report a flux of 0.005-0.57Mt CO 2 /yr in this region.But if one allows a fraction of the fluids released at slightly greater depths (e.g., between 80 and 90 km) to be channelized from the sub-arc to the forearc by migration along the subduction interface (e.g., Angiboust et al., 2014;Kawano et al., 2011), or along the impermeable carbonate sediment and its mafic (or ultramafic [e.g., Galvez et al., 2013]) basement, this estimate rises to 0.4 Mt CO 2 /yr.This happens to match the upper estimate of Barry et al. (2019) for Costa-Rica.Although this match may well be fortuitous, the sensitivity we report shows that it will be difficult to draw a definitive mass-balance estimate for the forearc, ie.one of the most dynamic and compositionally complex environment of the subduction zone.Therefore, our study supports the contention by Kelemen and Manning (2015) and Barry et al. (2019) that carbon release and storage toward the forearc may be greater than previously thought, as well as the notion of a high thermodynamic stability of carbonate beyond forearc and to subarc depth at least (Collins et al., 2015;Cook-Kollars et al., 2014;Gorman et al., 2006;Lefeuvre et al., 2020).The complex hydrodynamic regime of fluids in the slab may be the primary cause for this conundrum, more than the existence of possible inaccuracies in thermodynamic models of subduction zones.
Plank and Manning (2019) and others (Galvez & Pubellier, 2019) have stressed the importance of compositional heterogeneities between subduction zones.Our results reinforce this notion, as the decarbonation efficiencies of the western Pacific slab is much lower than that in the Eastern Pacific.This difference reflects fundamental compositional and thermal controls.Because subduction efficiencies are typically used as a parameter in longterm climate models such as those of Hoareau et al. (2015); Volk (1989), our work stresses that it should be treated as a variable, in space and certainly in time as well.Caution is therefore needed in applying numeric values of subduction efficiencies to situations where tectonic configurations, sediment composition, and/or slab geometries may have differed from today (e.g., before the Mesozoic).
Our global output from the slabs is lower than field estimate of passive CO 2 degassing at volcanic arcs of 55 Mt/ yr by Fischer et al. (2019), but the latter may include a fraction released by thermal metamorphism and assimilation of carbonate accreted (with both pelagic and shelf origins) to the overriding crust (Johnston et al., 2011;Lee et al., 2013;Mason et al., 2017).More importantly, our sensitivity test for hydration state revealed transport and supply limited regimes of element subductions, as well as thermally frozen slabs (Figure S4 in Supporting Information S1) which informs future research targets.For example, subduction zones of intermediate thermal structure are transport-limited, that is, C fluxes are sensitive to the hydration state of the AUM (Figure 5).This regime includes most slabs of the Sumatra, eastern Pacific, and Aegean (see Table 1 and Figure 5).This is where most C is released (Costa Rica and Guatemala are the largest C emitter per km of trench; Figure 6).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 (C)-limited that is, C fluxes are 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.In that case, the absence of AUM dehydration before sub-arc depth prevents C dissolution from those slabs which are coined frozen.

Metasomatic Redistribution of Elements and Arc Magma Signatures
Knowledge of the absolute mass change of an element typically matters in interpreting the chemical budget and isotopic signature of a rock (Angiboust et al., 2014;Bebout, 2007;Bebout & Barton, 1993;Galvez et al., 2013;Miller et al., 2009;Penniston-Dorland et al., 2012;Piccoli et al., 2016;Vitale Brovarone & Beyssac, 2014), whereas concentration changes typically matter for deriving the geophysical properties of a rock, including its rheology, mechanical behavior, and seismic properties (e.g., Connolly, 2005).Isolating both may be cumbersome, as it is well-known that compositional data may not reliably reflect metasomatic mass transfers (Ague, 2011;Aitchison, 1982;Gresens, 1967).Our results provide a first numerical demonstration that this phenomenon is at play in subducting slabs.We find that patterns of mass and concentration changes may be anticorrelated or correlated, reflecting the antagonistic role of metasomatism versus residual (passive) enrichment caused by volatile loss.By monitoring mass and concentration changes simultaneously, we confirm the possibility of early slab depotassification and dehydration as reported by Connolly and Galvez (2018), but this seems to only apply to certain slab compositions.The effect is most dramatic for carbonate rich and H 2 O deficient slabs (e.g., Nicaragua, Costa Rica), while some clay-rich slab sediments are net accumulator of alkali, for example, Alaska and Chile.
More importantly, all chemical elements are redistributed within the slabs before being released from them, and this redistribution reflects local gradients in chemical potential.As a result, most compositional and mass changes occur at the sediment/basalt and basalt/ultramafic interfaces (Figure 3).We note that our 4-layer slab structure is most representative of fast-spreading ''pacific-type'' subduction settings.Atlantic-type environments that lack a continuous mafic layer are relevant as well (Agard, 2021) and they would include highly contrasted sediment/ serpentinite interfaces (e.g., Galvez et al., 2013;Malvoisin et al., 2012;Scambelluri et al., 2022) which are not represented here.Interestingly, the first fluids produced in the altered ultramafic and mafic layers may be rich in both K and Mg, ie."shoshonitic".This is due to the highly incompatible nature of K, a minor component in those lithologies (Figure S3 in Supporting Information S1).However, regardless of where the fluids are produced in the rock column (e.g., basalt or ultramafic layers), those have generally acquired a clear sedimentary character by the time they reach the top of the slab.They have become enriched in Si, Al, alkali, carbon, and orders of 10.1029/2021JB023851 20 of 28 magnitude poorer in Mg than ultramafic fluids (Figure S3 in Supporting Information S1).In other words, we have demonstrated that slabs operate as "chromatographic" columns, and this effect should be considered in discussing the role that slab fluids may play in modifying the composition of the mantle wedge through time.
The fluids reaching the slab/mantle wedge interface are enriched in C, Si and Na (cf. Figures S2 and S3 in Supporting Information S1 for details).However, only hot fluids deliver Al in appreciable quantity to the slab surface, that is, at different depths depending on slabs thermal structure.This not only supports a metasomatic origin for jadeites (Angiboust, Muñoz-Montecinos, et al., 2021;Harlow & Sorensen, 2005), but it suggests that the alumina content (e.g., their aluminum saturation index) of subduction zone fluids may be used as chemo/geothermometer.The experimental evidence for the high incompatibility of K in high-P/T metasediments (Schmidt, 2015;Schmidt et al., 2004) suggests that the mobility of this element may be elevated at depth.The increasing K/Na ratio of our fluids with depth and temperature supports this notion (Figure S2 in Supporting Information S1).But this incompatibility also implies that dry segments of the slab maintaining K, Na in their compositions may be metastable with respect to a fluid phase.For example, we attribute the flat profile of sodium in the Cascadia sediment after 100 km (Figure 9e) to this metastability effect, caused by the lack of water influx.This may, in turn, cause extreme fluid peralkalinity if a fluid traverses the rock later in the subduction process.
Metasomatic decarbonation reactions have been described in the field (Abart, 1995;Ague, 2003;Bucher, 1999;Ferry et al., 2011), but never tested quantitatively in computational models of slabs compositions.We found that the dissolution of carbonate in the sediment of Nicaragua (and e.g., Tonga, Costa Rica, and Mariana) is associated with a synchronous gain in masses of Al 2 O 3 (Figure 10) and increased concentration of SiO 2 (Figure S6 in Supporting Information S1), reflected in a modal increase of garnet (grossular).This suggests garnet grows by infiltration of aluminous and siliceous fluids sourced in the dehydrating oceanic lithosphere (Figure 10d) according to the simplified reaction: where Al 2 O 3 (aq) and SiO 2 (s/aq) denotes the presence of reactive Al and Si in the infiltrating fluid.The concentration of Al in those high-P fluids is greatly enhanced by forming (alkali)-Al/Si clusters (Anderson & Burnham, 1983;Manning, 2007;Manning et al., 2010), here represented by the species AlO 2 [SiO 2 ] − , cf.Table S4 in Supporting Information S1).This illustrates that the fate of carbonate, which is here directly linked to the potentials μ(CaO) and μ(CO 2 ) (e.g., Equation 4), in fact hinges on multicomponent interactions in the bulk system.This is because CaO associates with SiO 2 and a range of other components such as MgO, FeO and Al 2 O 3 in minerals such as epidote, lawsonite, clinopyroxene, and garnet, as well as in carbonate solid solutions such as dolomite or ankerite (Lefeuvre et al., 2020) depending on P, T, and bulk composition.And since all those rock components are affected, to some degree, by metasomatic mass transfer (Figures 3 and 9) they all indirectly contribute to regulate μ(CO 2 ) (i.e., C solubility), and therefore the fate of carbon in those lithologies.We have not yet examined whether Mg and/or Fe cation substitution in carbonate plays an important role in the fluid-driven remobilization of CaO, as reported by Lefeuvre et al., 2020 in the western Alps.
An important implication of metasomatic growth of garnet in the subducting limestones of the central-eastern Pacific slabs, along with their progressive depletion in alkali and volatile (Figure 10), is that it would make those rocks much denser, and less buoyant, than previously thought (Behn et al., 2011;Ducea et al., 2022;Gerya & Meilick, 2011).In addition, assuming fluids tend to incorporate the lighter isotopes, our findings may also help understand the light isotopic signature of Ca (Antonelli et al., 2021) or Mg (Shen et al., 2018) in the metasomatized mantle wedge.Indeed, 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, and Mg, Table 2, Figure 8), 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.

Sensitivity and Thermomechanical Limitations of Model
Our model captures key aspects of the thermodynamics of fluid-rock interactions in subduction zones, but a vertical transport of fluid does not adequately reflect the natural environment in multiple ways.Future works may explore the effect of pressure gradients (Wang, 2019;Wilson et al., 2014), more complex advection-dispersion designs (Ague, 2007), and channelized fluid flow which may apply to zones of high permeability contrasts (Breeding et al., 2004), lithological interfaces (Angiboust et al., 2014), shear zone (Mancktelow, 2006), (micro) fractures, and connected porosity (wave; Plümper et al., 2017;Tian et al., 2019).While lateral fluid flow may significantly modify our fore-arc budgets, they may also create situations where only a fraction of the fluid reacts with the rock it traverses (as opposed to 100% in our model) while moving along a preferential flow path.We have tested this semi-quantitatively in the case of Cascadia, Nicaragua, and Honshu (Figure S14 in Supporting Information S1).The results suggest that reducing the reactive volume by 50% tends to reduce the amount of carbon released from the slabs by ca.20%-30% before subarc depth.As mentioned previously, Penniston-Dorland et al. ( 2015) argued that slabs may be hotter than predicted by thermomechanical models (by > 300°C).
Our semi-quantitative tests (Table S7 and Figures S14, S21 and S22 in Supporting Information S1) show that increasing the temperature by +50 and +100°C increases the carbon release by ca.10%-30%.In other words, temperature effects may counter those of a channelization of fluid flow.
Computational aspects are relevant as well.By definition, none of the existing variants of the backcalculation algorithm (Galvez et al., 2015;Connolly & Galvez, 2018) include electrolytes in the search for the Gibbs energy minimum of the rock/solvent system.Although it is highly precise in most rock-dominated situations (i.e., where rock/electrolyte ratio is very high), a backcalculation strategy necessarily loses precision when specific components reach vanishingly low concentrations in the rock while it is present in the low volume of infiltrating fluid (e.g., specific fluid-dominated condition only occurs here in the case of K 2 O, H 2 O and CO 2 for Nicaragua and Mexico sediments).Practically, this means our algorithm may sometimes fail in refining low amounts of phengite and carbonate that may form transiently in these situations.However, this is inconsequential here because (a) the mass involved are always negligible, and, (b) more importantly, mass balance is maintained at all times.That is, the absence of suitable solid hosts for a specific component results here in its conservative, and out-of-equilibrium, transfer from the basaltic layer to the top of the slab in soluble form (e.g., potassic fluid in K-free sediment, Figures S2 and S3 in Supporting Information S1), with minimal interaction in the sediment.The implication is that the composition of a natural rock may not faithfully record the composition of an infiltrating fluid if some fluid components do not have a thermodynamically stable host in the rock.Thermodynamically more rigorous alternatives to our approach exist and would include all solids, solutes, and solvents in finding the Gibbs energy minimum of the system.The softwares HCh/HighPGibbs (Shvarov, 2008;Zhong et al., 2020) or GEM-Selektor (Kulik et al., 2013) implement such comprehensive GEM approach, but their algorithmic advantage (in part designed to address low-P systems were fluid-rock ratio are typically high) comes at the cost of a more limited solid solution package, which currently restricts their applicability to high-P systems were multicomponent minerals and solid-solutions predominate.The steep computational cost of such strategies is deemed unjustified for the present purpose.Indeed, we believe the main obstacle to accurate prediction of high-P fluid compositions is not algorithmic anymore, but experimental.
The accuracy of our fluid composition and stoichiometry depends on the quality of our thermochemical databases and of the experiments used to build them (Methods, Table S4 in Supporting Information S1).Our predictions for the solubility of the assemblage albite-paragonite-quartz (Al, Si, and Na concentrations in fluid, Figure S16 in Supporting Information S1) are slightly underestimated compared to the experimental data from Manning et al. (2010), but the orders of magnitude are correct over the P-T range considered.Increasing difference with T is caused by the progressive polymerization of alkali-alumino-silicate clusters (Anderson & Burnham, 1983;Newton & Manning, 2008;Pokrovski, 1998).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 C speciation is controlled by Ca-C complexes (cf. Figure S2 in Supporting Information S1)-are also probably relatively accurate.However, there is no equivalent data for Mg-carbonates and for the thermodynamics of Mg-bearing electrolytes.This means that unlike Ca, Mg concentration and Mg fluxes out of slabs may be underestimated by about 2 orders of magnitude or more (Galvez et al., 2015), and that our predicted solid solutions may bias toward their Mg-endmembers.The solubility of Fe-bearing phases, with their well-known redox sensitivity, may be the least well constrained of all.In this respect, it is worth noting that the solubility of corundum (Al 2 O 3 ) increases by about 2 orders of magnitude in the presence of quartz, at 700°C and 10 kbar (Manning (2007)).By analogy with the Al 3+ cation, we conjecture that Fe 3+ may be readily dissolved and transported as Si or Al-Si-complexes in high P-T metasomatic environments.The experiments of Pokrovski et al. (2003) support this hypothesis.Similarly, 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.Fe concentration in modeled metamorphic fluids may therefore be currently underestimated by about 2 orders of magnitude at least, and our predicted solid solutions may bias toward Fe-rich compositions.This hints to a set of relevant experimental data that would enhance the predictive accuracy of thermopetrological models, and allow exploring more complex redox systems (e.g., Evans & Frost, 2020;Galvez et al., 2013;Galvez & Jaccard, 2021).More generally, we note that accurate element solubility predictions ensure accurate prediction of its mobility/compatibility relation in the (open) fluid/rock system.But accurate fluid stoichiometry (i.e., element ratio in the fluid) is often more important because it determines the accuracy of mineralogical compositions (relative ratio of elements within a given phase, regardless of the modal abundance of the phase).
As for the accuracy of global elemental budgets, the above limitations appear relatively minor in front of the most critical assumption of our simulation, which is mechanical.We have made the common assumption that slabs are mechanically continuous, but the geological record and theoretical studies suggest this is certainly not the case in general (Malavieille & Ritz, 1989;Platt, 1986) as there are many non-advective (fluid-mediated) ways to remove material from the slab and export it to the wedge.Possibilities include the buoyant rise of subducted slabs of crustal material relative to higher density mantle rocks (Ernst, 1975;Gerya & Meilick, 2011;Molnar & Gray, 1979;Platt, 1986), and underplating of material at the base of accretionary wedge or of the continental mantle (Angiboust et al., 2022;Chakraborty et al., 1995;Cowan & Silling, 1978;Platt, 1986).For example, the metasedimentary units exposed in the Mediterranean region (Apennines, western and central Alps etc.) are built on top of Triassic evaporitic successions.These weak evaporitic layers in the western Alps (Barré et al., 2020;Gabalda et al., 2009), and elsewhere (Molli et al., 2020;Ring & Layer, 2003), favor detachment and exhumation of the sedimentary cover, and preclude its subduction (Malavieille & Ritz, 1989;Menant et al., 2020;Platt, 1986).
Evidence for active underplating in the 25-45 km depth range have also been reported in today's Cascadia, Northern Chile, Alaska, and other subduction zones of the circumpacific area, and reviewed recently (Angiboust et al., 2022).Within the context of this study, this means that an unquantified fraction of mass transport to the mantle wedge may occur mechanically/gravitationally, possibly promoted by hybridization at the slab/mantle interface (Bebout & Barton, 1993;Sorensen et al., 2010), and not by advective transport.

Conclusions
Conventional petrological modeling has assumed slabs as consisting of discrete, compositionally homogeneous layers distinguished by sharp compositional transitions and exchanging H 2 O and CO 2 only.By contrast, our model leads to intriguing patterns where layer boundaries become diffuse once advection of mass is considered.This is made possible by considering both volatiles and nonvolatiles as mobile components in the fluid phase.We found that the mobility of volatile and non-volatile elements is interdependent due to the formation of multicomponent aqueous species, and to the importance of multicomponent mineral and solid-solutions at elevated pressure and temperature.
Our slab operates as a chromatographic column in which elements are redistributed according to their respective affinity, or not, with the local rocks.As a result, rock compositions deviate significantly by metasomatism.Al, Si metasomatism, for example, drives the decarbonation of subducting limestones of the Central-eastern Pacific, and the synchronous growth of metasomatic garnet and loss of alkali and volatiles from these systems should form denser and less buoyant assemblages than previously thought.Conversely, some slabs accumulate alkali, and maintain them across dry segments of their trajectory, suggesting they may be metastable with respect to a fluid phase until deep in the subduction.
The fluids produced in the altered oceanic lithosphere gradually loose their ultramafic/mafic character as they move up in the rock column, and they reach the top of the slab with distinct sedimentary signatures.This means slab sediments are overall efficient at buffering fluid compositions before they reach and metasomatize the mantle wedge, at least until ca.150 km depth.Of course, the present modelling framework may be applied to other types of lithological interfaces of local or global geological relevance (e.g., sediment-ultramafic interfaces).However, the constraints that our study brings on the compositional evolution of slabs is key because the composition of subducted rocks influences their rheology and mechanical properties.
Our forearc carbon flux ranges from 3% to 17% of the trench input, in large part due to its pronounced sensitivity to the choice of forearc-subarc threshold depth.This suggests that the existing discrepancy between model predictions and field observations of carbon release at fore-arc depths may be caused by the complex hydrodynamic regime of the shallow slabs rather than by inaccurate thermodynamics.The subarc flux represents 51%-62% of the input, or 26-32 Mt C/y, which supports the existing notion of a high carbonate stability in subduction zones today.But we found that local heterogeneities between subduction zones are the rule rather than the exception.Moreover, this result does not mean much carbon is mechanically subducted beyond subarc depth, let alone the transition zone.This is because the chemical budget of our slabs depends on their structure and mechanical behavior, that is, that they move as a single mechanically continuous unit.This common assumption of thermodynamic-petrological simulations does not fully represent the observational, analogical, or thermomechanical evidences suggesting that sedimentary layers (pelagic or continental) may be detached from slabs along weak and ductile evaporitic layers, and underplated.
From a broader perspective, this work poses the question as to what part of the slab contribution to the mantle wedge is caused by fluid-driven mass advection, as simulated here, or by gravitational transport, for example by mantle hybridization, buoyancy forces and/or underplating.These issues must await future theoretical and experimental efforts.

Figure 1 .
Figure 1.(a) Thermal structure of the Cascadia subduction zone (redrawn from 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 modeling purpose.(b) Schematic depiction of the modeling approach.The rock undergoes dissolution first during dehydration reaction and the volatile and dissolved nonvolatile 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 highlighted.(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 colors.The subduction length is 850 km for the Juan de 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.
values are based on low upper mantle hydration state model (LHS) (ca.1,166 Mt/yr H 2 O/yr global water input with ca.250 Mt/yr from mantle peridotite).The full mass balance results from IHS and HHS models are provided in Tables computation time while offering sufficient resolution to monitor the mass transport process.It limits overstepping of the dehydration reactions, and ensures minimal fluid/rock ratio.Temperature at each grid point is interpolated based on the thermal model fromSyracuse et al. (2010).Hotter alternatives(Penniston-Dorland et al., 2015) are tested(Figures S21 and S22  in Supporting Information S1) and discussed.Pressure is related to depth assuming lithostatic pressure gradient.The P-T field is therefore specific for each slab.

Figure 2 .
Figure 2. Relative water mass variations per 100 g of rock (H 2 O mass subtracted by the initial H 2 O mass divided by the total rock mass × 100).Each phase field present the associated equilibrium phase assemblages within a hot (North Cascadia, Western U.S.), intermediate (Nicaragua, Central America), and cold (Honshu, Japan) subduction end-members (HHS model).Red color indicates net gain in mass due to metasomatism and blue color indicates net loss.The in-slab depth axis (y-axis) is slightly tilted to better visualize the subduction slab stratification.The initial water input in the mantle is assumed to vary linearly along the in-slab direction, which explains the drier phase assemblages in the lower part of the ultramafic section (cf.Figures3a-3c).The thickness of the sediment is plotted as 1 km for Nicaragua and Honshu so that the phase assemblages can be seen.Detailed flux of water and composition of fluid along the slab surface (FigureS2in Supporting Information S1) and composition of fluid along selected vertical transects of the three systems (FigureS3in Supporting Information S1) are provided in the supplementary material.

Figure 3 .
Figure 3. Pattern of water, carbon, Al 2 O 3, and K 2 O redistribution and release for warm, intermediate, and cold slabs as a function of hydration states.Relative H 2 O, CO 2 , Al 2 O 3, and K 2 O mass changes (Δm i ) of components are given in % of the total rock mass at first grid.The percentage values displayed as inset within the panels express this variation in % of the input mass of the component at first grid.The warm color indicates absolute mass gain and cold color indicates absolute mass loss.For Honshu slab, the first grid undergoes dehydration in the sediment because of its excessively 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.

Figure 4 .
Figure 4. Cumulated fraction of CO 2 lost from top of the slab for the three subduction zones, highlighting contrasted sensitivity to altered upper mantle hydration state for the carbon limited (Cascadia), water limited (Nicaragua), and frozen (Honshu) regimes.

Figure 5 .
Figure 5.Total top-of-the-slab CO 2 flux (left panel) and decarbonation efficiency (right panel, ηCO 2 ) to sub-arc depth 150 km for all the studied subduction zones, ranked geographically.The left panel shows the integrated flux of CO 2 in Mt/yr from surface to sub-arc depth (150 km).It shows that the canonical value of ca.ηCO 2 = 0.5(Volk, 1989) is only relevant in average, and today(Galvez & Pubellier, 2019).Note that the E-Pacific slabs have the highest extraction efficiency.The decarbonation efficiency of slabs generally increases with the geothermal gradient (FigureS4in Supporting Information S1).

Figure 6 .
Figure 6.Panels (a and c) show the cumulated CO 2 and H 2 O output per km of trench from surface to 150 km depth.Panels (b and d) show the flux difference between high hydration state (HHS) and low hydration state (LHS) model.Slabs that are sensitive to the hydration state of AUM are labeled in black (water-limited regime), slabs that are insensitive because they are fully decarbonated with LHS model are labeled in red (carbon-limited regime), slabs that are insensitive because they do not lose much C even with HHS model are labeled in blue (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 and 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 increase the CO 2 flux further.Note that the forearc flux are particularly sensitive to the choice, semi-arbitrary, of the fore-arc/sub-arc threshold.Here it is placed at 80 km (see text for discussion).

Figure 7 .
Figure 7. η-η diagrams showing the decarbonation efficiency ηCO 2 versus dehydration efficiency ηH 2 O for (a and b) Cascadia, (c) South Chile, and (d) the global subduction system (results from model low hydration state [LHS]).Curves, or portion of curves, with highly positive slopes indicate that most of the carbon is removed with minimal H 2 O release from the slab (i.e., high C solubility), while portions of curves with flat trends express the opposite.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 η-η diagrams for the other slabs can be found in supplementary materials (Figure S19 in Supporting Information S1).Panel b (red) is the result from runs where H 2 O content in the upper mantle of North Cascadia (British Columbia) is 7.5 times lower than even our driest LHS model, due to limited plate bending and H 2 O penetration within the oceanic mantle (cf.Canales et al., 2017 for details), or 1.5 wetter than our LHS model (blue).Regardless of the hydration state, the similarity shows that the slab is supply limited.

Figure 8 .
Figure 8. Global ranking of the SiO 2 , Al 2 O 3 , Na 2 O, and K 2 O fluxes for subarc and forearc depth with the same setup as in Figure 5.Total fluxes and total flux per km of trench are given.

Figure 9 .
Figure 9.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.Relative concentration (Δc i ) and mass variations (Δm i ) of component i are given in % of the input mass of the component at first grid.Unlike concentration variation, the mass variation reflects the actual metasomatic effect.Models assuming infinite compatibility of nonvolatile major elements in the hydrous fluid phase would typically show no net mass change (i.e., flat blue curve at 0% mass change for nonvolatile components).Concentration changes may always occur due to residual enrichment or dilution when volatiles are lost or added, respectively.Only the results of the low hydration state model are displayed.

Figure 10 .
Figure 10.(a and 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, and 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.
The global input of CO 2 associated with the basalt + gabbro (igneous) layer is ca.85.8 Mt CO 2 / yr (23 Mt C/yr), and it is around 7.5 Mt CO 2 /yr (2.1 Mt C/yr) for the altered upper mantle.The total CO 2 trench input is ca.186 Mt CO 2 /yr with 92.6 Mt CO 2 /yr from the sedimentary layer, which reflects the global estimate of Keken et al., 2011)timate AUM input of ca.860 Mt H 2 O/yr, while the upper-estimate of Rüpke et al. (2004) is ca.1,200MtH 2 O/yr (see Table 3 in VanKeken et al., 2011).We opted for a conservative value, and our high hydration state (HHS) model assumes an AUM input flux of ca.
Parai and Mukhopadhyay (2012)ent compositions ofPlank (2013).bFromJarrard(2003).cFromSchmidtand Poli (2003).dFromHacker(2008).eCorresponds to the low end fromRüpke et al. (2004);Van Keken et al. (2011)and low end ofSchmidt and Poli (2003).fThisrefers to the H 2 O lost at the first grid during phase equilibrium computation.gThisindicatesthereturn fluxes associated with 3 different sea-level decline scenario in the past 542 Myr based onParai and Mukhopadhyay (2012), and 0 m indicate no sea-level change.The sea-level change is calculated based on our computed return flux and the return flux vs. sea level change relationship fromParai and Mukhopadhyay (2012).Global Water Budget Derived From Our Models and Implications for Sealevel Variations in the Phanerozoic a