Carbon Storage in the Forearc Produced by Buoyant Diapirs of Subducted Sediment

Carbonate sediments transported into the mantle at subduction zone settings account for the majority of the carbon flux into the Earth's interior and are thus critical to the deep carbon cycle. Understanding carbon storage volumes in the deep earth requires knowledge of the degree to which carbonate sediments are stored in the arc lithosphere or descend to the deep mantle. Here, we use petrological‐thermomechanical modeling to indicate that solid‐state diapirs dominate the removal of carbon from subducting plates, which may be the principal carbon‐release mechanism for the Cyclades (Greece) and Costa Rican forearcs. We find that forearc diapirs remove up to ∼80% of subducting carbon and develop diagonally upward, resulting in massive carbon storage in the subarc lithosphere. Outgassing from the carbon storage may cause high carbon outputs and explain volcanic gas with high δ13C at some subduction zones, affecting atmospheric CO2 concentration.


Introduction
There is increasing consensus that tectonic processes involving carbon play a crucial role in the modulation of Earth's climate system on geological time scales.Subduction is viewed as the most important physical process linking carbon exchange between the Earth's interior and atmosphere, transporting both organic and inorganic carbon (current carbon flux of ∼68-96 Mt C/yr) into the mantle (Barry et al., 2019;Clift, 2017;Dasgupta & Hirschmann, 2010;Kelemen & Manning, 2015).Carbonate sediments are thought to dominate the subducting carbon flux (current carbon flux of ∼48-57 Mt C/yr), which is more than twice the sum of the input fluxes of the oceanic crust and mantle lithosphere (Clift, 2017;Dutkiewicz et al., 2018;Kelemen & Manning, 2015;Plank & Manning, 2019).The subduction zone carbon flux appears to be greater than the volcanic output plus diffuse outgassing combined (Aiuppa et al., 2019;Kelemen & Manning, 2015).It has been proposed that only a small amount of carbon is recycled into the deep convecting mantle, and thus, shallow Earth reservoirs above subarc depths may act as carbon sinks on geological timeframes (Kelemen & Manning, 2015).However, the subduction and recycling of sedimentary carbon below forearcs remain poorly understood.A substantial amount of carbon is thought to be released during subduction at subarc depths (∼70-150 km) (Plank & Manning, 2019), but the forearc (<70 km) has also been postulated to be a significant reservoir for carbon removal from subducting sediments (Müller et al., 2022;Stewart & Ague, 2020;Zhao et al., 2023).
• Young oceanic plates, slow convergence rates, and thick sediments promote forearc diapirs and cause high decarbonation efficiency • Forearc diapirs dominate removal of carbon from the subducting plate and reduce the proportion of the carbon released at subarc depths Supporting Information: Supporting Information may be found in the online version of this article.
Mechanical removal, metamorphic reactions, melting, and aqueous dissolution are the four carbon-release mechanisms thought to occur at shallow depths (<150 km) in subduction zones.At subarc depths, aqueous dissolution of carbonates in the sediments and crust is thought to carry the majority of the carbon into the source of arc magma in the mantle wedge (Ague & Nicolescu, 2014;W. Chen et al., 2023;Farsang et al., 2021), but this is limited to conditions with high aqueous fluid production and high Ca + content carbonates.During the infiltration of a water-bearing fluid at forearc depths, metamorphic decarbonation reactions have been suggested to be a significant mechanism for releasing carbon from sediments (e.g., carbonated siliciclastic) and altered ocean crust (Gorman et al., 2006;Stewart & Ague, 2020;Tian et al., 2019;Zhong & Galvez, 2022).Diapirs have been proposed to be more efficient than metamorphism and melting during recycling of subducted carbon (Behn et al., 2011;C. Chen et al., 2021;Marschall & Schumacher, 2012).At temperatures of 500-850°C and depths >60 km, buoyant sediments are thought to detach from the downgoing slab and form diapirs (Behn et al., 2011), consistent with numerical modeling results (Currie et al., 2007;Ducea et al., 2022;Gerya & Meilick, 2011;Gonzalez et al., 2016;Klein & Behn, 2021).
Geodynamic studies have suggested that forearc diapirs of carbonate sediments may occur during subduction (Gonzalez et al., 2016), but the fraction of subducted sediments that is detached at subarc depths is debated (Ducea et al., 2022).Here, we present results of high-resolution two-dimensional (2D) petrological-thermomechanical models that explore the geodynamic evolution of carbonate sediment subduction and recycling at depths less than 150 km in subduction zones.The two main carbon-release mechanisms addressed here are diapirism and metamorphism, which cause significant carbon release in the forearc.Because of the low degree of these processes at shallow depths (e.g., W. Chen et al., 2023;Thomson et al., 2016), aqueous dissolution and partial melting are not included in our models.Our results provide insights into the nature of carbonate sediment recycling in subduction zones and have implications for carbon storage in the subarc lithosphere and long-term climate system.
We designed a 2D coupled petrological-thermomechanical numerical model to simulate the sediment carbonrelease processes in subduction zones where the oceanic plate underthrusts the overriding continental plate (Figures 1a and 1b).The model domain is 2,000 km wide × 400 km deep, with 2,001 × 401 nonuniformly spaced nodes.In the study area, which spans 1,200-1,800 km in the horizontal direction and 0-160 km in depth, the spatial resolution is 500 × 500 m.The model uses a 1.7 × 3 km resolution on the outside of the main area of interest (0-1,200 and 1,800-2,000 km in the horizontal direction, and 160-400 km in depth).The subducting oceanic plate consists of hydrated and carbonated sediments of varying thickness, a two-layer crust, and the mantle lithosphere.The effects of the sediment thickness (d ss ), which ranges from 0.5 to 4 km, on the modeled diapirs and the decarbonation efficiency were investigated.The oceanic crust is 7 km thick, with 2 km of hydrothermally altered basalt and 5 km of gabbro.The continental crust of the overriding plate is 20 km thick in the upper layer and 10 km thick in the lower layer.The mantle lithosphere of both the oceanic and continental plates is set as anhydrous and carbonate-free peridotite.The accretionary wedge, which is composed of sediments, and the underlying weak zone decouple the subducting oceanic plate from the overlying continental lithosphere during subduction initiation.The initial subduction dip angle is set to 30°in the models.
The top and side are free slip boundaries, and the bottom boundary is open.A sticky air layer, with a thickness of 10 km and a low viscosity of 10 18 Pa s, is implemented on top to mimic a free surface (Crameri et al., 2012;Gerya, 2019).The convergence rates are applied at x = 1,000 km on the subducting oceanic plate and the thermal age increases from 0 to 20 Myr to drive the model evolution and is imposed at 0-80 km (x-axis) along the left boundary.In our model box, the thermal boundary conditions are fixed at the top (0°C) and bottom (1545°C), and there is no heat flux at the side boundaries.A half-space thermal cooling model determines the temperature field of the oceanic plate, the age of which ranges from 10 to 80 Myr.Different thermal ages of the subducting oceanic plate (T age ) and different convergence rates (V push , which range from 20 to 120 km/Myr) were investigated to test the effects of the slab's thermal structure on the decarbonation efficiency of the model.The continental plate has a linear temperature profile, which increases from the top (0°C) to the base (1350°C) of the lithosphere.The adiabatic thermal gradient is 0.5 K/km in the upper mantle.The erosion and accretion at the model's surface are implemented by solving the diffusion equation at the Eulerian nodes in each time step (Gerya & Yuen, 2003).We use the Perple_X software version 6.9 (Connolly, 2009) with a Gibbs free energy minimization approach and an internally consistent thermodynamic database of Holland and Powell (1998) to conduct thermodynamic models.We determine the stable mineralogical phases and fluid compositions at equilibrium as a function of pressure and temperature on a grid of 313 × 313 nodes with a resolution of 3.8 K and 22 MPa.The global subducting sediment composition (GLOSS) was used for the carbonate sediments, one of the best constrained compositions for materials entering subduction zones.The carbonate sediments were assumed to contain 3.01 wt% CO 2 , and the CO 2 contents of altered basalt and mantle peridotite were not taken into account (Hart & Zindler, 1986;Plank & Langmuir, 1998;Staudigel et al., 1989).The CO 2 content of GLOSS sediment decreases with increasing temperature and pressure due to metamorphic decarbonation (Figure 1c).

Results
The thermal parameter, defined as T age × V push , controls the thermal structure and thus the metamorphic devolatilization processes in subduction zones (Kirby et al., 1996).Similar to previous studies (Johnston et al., 2011;Kirby et al., 1996), the sediment thermal parameter (STP) is defined as to identify our simulated results (See Supporting Information S1).We first describe the evolution of three endmember models, which are recognized based on the computed STP values ranging from 1.38 to 2.57, and additional test models.Then, we focus on the influence of STP on sediment diapirs and decarbonation efficiency.Details on additional models and model uncertainties (including our use of GLOSS sediment and closed-system thermodynamic modeling; e.g., Connolly & Galvez, 2018;Galvez et al., 2015;Gorce et al., 2019;Kerrick & Connolly, 2001) are presented in Supporting Information S1.

Carbon Accumulation Due To Sediment Diapirs
The models with medium STP values of 1.84-2.09indicate sediment diapirs form first, followed by diagonally upward migration.The evolution of a Reference Model with an STP of 1.95 shows that buoyant sediments accumulate at forearc depths and sedimentary diapirs with high CO 2 contents initiate in the solid-state at ∼11 Myr (Figure 2a).As a result of these diapirs together with metamorphic decarbonation, the CO 2 content in subducting plates below ∼70 km depth is significantly reduced.The diapirs in the forearc migrate diagonally upward due to mantle corner flow, and then undergo horizontal motion when approaching the overlying lower crust (Figure 2b).In addition to the carbonates derived from incoming sediments, subarc diapirs with low CO 2 contents rise buoyantly along the surface of the subducting slab and mix with the shallower diapirs at forearc depths, adding additional carbon to the forearc reservoir (Figure 2c).Large amounts of carbon are transported and stored in the forarc lithosphere as a result of these diapirs.
The variations in STP do not change the sediment diapirs at forearc depths, but they do alter their dynamic processes.Compared to the Reference Model, forearc carbon accumulation is dramatically less in the models with high STP (2.09-2.57)such as model MC8S05 (Figure 2d).In model MC8S05, sedimentary diapirs flow vertically upward rather than along the slab surface.Besides the formation of diapirs in the forearc, low STP values (1.38-1.82)facilitate diapir exhumation within the subduction channel (Figure 2e).The subduction channel is increasingly widened and filled by exhumed rocks, and diagonal upwelling of sedimentary diapirs occurs, followed by horizontal migration.
Additional test models show that forearc diapirs occur for a wide range of conditions, which including various density of sediment and serpentine, rheological weakening by fluids, and various rheological strength of both the carbonate sediment and overlying continental lithosphere (see Supporting Information S1; e.g., Arcay et al., 2005;Hirth & Kohlstedt, 2004;Huang et al., 2013;Rudnick & Gao, 2003, 2014;Saffer, 2003;Schmidt & Poli, 1998;Zhang et al., 2020).For example, different carbonate sediment rock types including wet quartzite (Ranalli, 1995), calcite (Schmid et al., 1977), dolomite (Holyoke et al., 2013), and magnesite (Holyoke et al., 2014), each of which has a distinct rheological strength, can influence the dynamic processes of forearc diapirs (Figures S1a and S1b in Supporting Information S1).Calcite is the weakest among them, although experimental data on their rheological strengths have large uncertainties due to grain-size dependence (e.g., Brodie & Rutter, 2000;Holyoke et al., 2013).We find that the calcite and dolomite diapirs in the forearc are comparable to computed results in models with low STP and medium STP, respectively.In contrast, a strong sediment layer like magnesite inhibits the accumulation and migration of forearc diapirs (Figure S1b in Supporting Information S1).However, at the forearc depths, the stable crystalline carbonate is generally dominated by calcite/dolomite (Dasgupta & Hirschmann, 2010).Moreover, due to the prograde metamorphism of the subducting slab, the forearc lithosphere is hydrated and weakened (Figure S2a in Supporting Information S1).Semi-analytical solution suggests that a weak overriding lithosphere could promote the formation of the diapirs (Klein & Behn, 2021;Weinberg & Podladchikov, 1994), consistent with our model results (Figure S2b in Supporting Information S1).Thus, the general dynamics of the formation of diapirs in the forearc region is robust.

Decarbonation Efficiency Via Diapirs
We grouped the models according to their STP values to obtain a steady-state decarbonation efficiency for each investigated STP (Table S2 in Supporting Information S1).Our results (Figure 3) reveal that the decarbonation efficiency of forearc diapirs increases linearly with decreasing STP, reaching up to ∼80%, implying that substantial carbon can be stored forearc regions above subduction zones with young oceanic plates, slow convergence rates, and thick layers of sediments.At low STP values (1.38-1.82),these three parameters speed up the migration and accumulation of sedimentary diapirs, resulting in a high decarbonation efficiency.As a result of the removal of carbon-bearing sediments at shallower depths, the decarbonation efficiency of subarc diapirs reaches zero in most of the modeled models.At medium (1.84-2.09) to high STP values (2.09-2.57),as the STP increases, subarc diapirs are enhanced because there are fewer shallow forearc diapirs.However, only a few models produced a high decarbonation efficiency for the subarc diapirs.Therefore, the decarbonation efficiency of the diapirs at subarc depths is positively correlated with the STP, but is generally low.

Evidence for Diapir-Induced Carbon Storage in the Forearc
Subduction zones in which carbonates in subducting sediments dominate the carbon input were chosen to compare modeled and observed decarbonation efficiencies based on carbon input, carbon output, and the STP parameters (T age , V push , and d ss ; Figure 3; Tables S2 and S3 in Supporting Information S1; Barry et al., 2019;Clift, 2017;Fischer et al., 2019;Jarrard, 1986Jarrard, , 2003;;Li et al., 2019;Stewart & Ague, 2020;Syracuse et al., 2010).The forearc is considered to be a carbon reservoir in our geodynamic models.Seismic evidence for diapirs within the mantle wedge has been discovered by recent observations of scattering obstacles located above the Ryukyu slab, ranging from ∼60 to ∼95 km (Lin et al., 2019(Lin et al., , 2021)).Similar to diapirs, accreted underplates have revealed low-seismic velocity anomalies in many Pacific subduction zones, for example, beneath the forearcs of Cascadia, Alaska, and the Aleutian Islands (Scholl, 2019).Forearc carbon reservoir is also supported by geological and geochemical evidence from the subduction zones in the Cyclades (Greece) and Costa Rica regions.Forearc decarbonation in the Cyclades (Greece) accounts for ∼40%-65% of the carbon input of the subducting slab according to geochemical analysis and thermodynamic modeling (Stewart & Ague, 2020).Our results indicate a ∼51.47% decarbonation efficiency via advection of carbon in sedimentary diapirs in model MA8S4, but predict a low metamorphic decarbonation of <0.01%(Table S2 in Supporting Information S1).Thus, rather than the metamorphic reactions proposed by Stewart and Ague (2020), diapirism may be the carbon-release mechanism at forearc depths that explains the massive carbon release in the Cyclades (Greece).Based on helium and carbon isotopes, 0.01%-11.82% of the subducting carbon is released beneath the forearc in Costa Rica, and 91% of the carbon released is stored in the forearc crust via calcite deposition (Barry et al., 2019).In the forearc in Costa Rica, sediment diapirs would be the principal mechanism for removing carbon from the slab, and their corresponding decarbonation efficiency of 3.63% in model MC8S05 (STP of 2.40; Figure 3) is consistent with the observations.

Implications for Massive Carbon in the Subarc Lithosphere
We propose that carbon stored in the subarc lithosphere can be derived from subducted sediments via transport in forearc diapirs.Our results indicate that forearc sedimentary diapirs initiate in the solid-state at the slab surface (Behn et al., 2011;Gerya & Yuen, 2003;Marschall & Schumacher, 2012).In addition to their movement along a subduction channel, they migrate diagonally upward, as induced by mantle corner flow (Freda et al., 2008), and flow laterally into the subarc mantle lithosphere (Figure 4).This dynamic transport is one form of relamination (Hacker et al., 2011(Hacker et al., , 2015;;Kelemen & Behn, 2016) and is compatible with the conceptual model that proposes a diagonal rather than vertical direction for diapiric rocks (C.Chen et al., 2021;Marschall & Schumacher, 2012).Furthermore, the efficiency of carbon extraction from the subducting plates via forearc diapirs was calculated to reach ∼80% in this study.This result is consistent with thermodynamic calculations, that indicating the carbon outflux from subducting sediments above a depth of 125 km can account for 83% (up to 49 Mt C/yr) of the sedimentary carbon input (up to 59 Mt C/yr) in some arcs (Müller et al., 2022).In contrast, the overlying crust beneath arcs (Figure 4) could be one of shallow carbon reservoirs (e.g., Lee et al., 2013).However, the mechanism for the emplacement of sedimentary carbonates from the surface into the deep crust remains unresolved.
Carbonate sediments in the subarc lithosphere can be remobilized and may dominate volcanic arc emissions (Lee et al., 2013;Mason et al., 2017;Shen et al., 2023;Wang et al., 2023).Geochemical and geophysical evidences for the transport of carbonate sediment diapirs to the surface have been suggested in the subduction systems including northern China (Liu et al., 2015), southern China (Liu et al., 2021), northern California (Tewksbury-Christle et al., 2021), Mexico (Gómez-Tuena et al., 2018), and central Alps (Condit et al., 2022).Anomalously 13 Cenriched volcanic gas compositions of the Kusatsu-Shirane and Izu-Oshima arcs in SW Japan; Guadeloupe and Dominica arcs in the Antilles; and the Java arc suggest a sedimentary source for the CO 2 (Mason et al., 2017).For example, the decarbonation efficiency beneath the SW Japan arc is up to 52.05% ± 30.92% (Table S3 in Supporting Information S1; e.g., Clift, 2017;Stewart & Ague, 2020), which is approximately equal to the decarbonation efficiency of ∼46.39% due to forearc diapirs in our models (model MST2; Figure 3).And the Antilles and Java subduction zones have STP parameters similar to our models MA8S2 and MA8C8, and high decarbonation efficiencies of 35.15% ± 20.80% and 9.33% ± 5.43% (Table S3 in Supporting Information S1; e.g., Clift, 2017;Stewart & Ague, 2020), respectively.Thus, these observations are consistent with our suggestion that the abnormally high magmatic carbon recycling efficiency in these arcs is caused by remobilization of sedimentary carbonates stored in the subarc lithosphere (Wang et al., 2023).In addition, the fate of forearc diapirs has crucial implications for the carbon reservoirs in the deep convecting mantle in addition to the shallow reservoirs and carbon outputs (Figure 4).The presence of massive sedimentary carbon sinks at forearc depths would reduce the recycling of carbon into the convecting mantle (Kelemen & Manning, 2015;Sieber et al., 2020;Stewart & Ague, 2020).We suggest that the sedimentary carbon entering the deep mantle could be <20% of the carbon input for small STP subduction zones, since the size of shallow carbon reservoirs may be underestimated in our study, which did not include the carbonate dissolution mechanism.

Conclusions
Systematic petrological-thermomechanical models are conducted to investigate the geodynamic evolution of sediment subduction and recycling at depths less than 150 km in subduction zones.We find that forearc diapirs dominate the carbon release from subducting plates and that metamorphic decarbonation is a relatively minor process during long-term subduction of carbonate-bearing sediments.Young oceanic plates, slow convergence rates, and a thick layer of sediments (corresponding to low STP values) promote forearc diapirs, thereby causing high decarbonation efficiency.It is important to notice that diapirs are transported diagonally upward, resulting in massive carbon storage in the overlying subarc lithosphere, reducing the proportion of the carbon released at subarc depths and recycled into the convecting mantle reservoir.We suggest that high carbon outputs and volcanic gas with high δ 13 C values may be caused by remobilization of sedimentary carbonates stored in the subarc lithosphere.

Figure 1 .
Figure 1.Numerical model configuration.(a) Compositional setup of the full model domain.The thick black line and arrow indicate the location and direction of the externally imposed velocity (V push ), respectively.(b) Enlargement of the area in the dashed box in panel (a) shows the geometry of the subducting oceanic plate and the weak zone prescribed for subduction initiation.The black lines show the initial temperature fields (°C).(c) Pressure-temperature diagrams of the CO 2 content for average global subducting sediment composition sediments derived using Perple_X and Gibbs free energy minimization(Connolly, 2009).

Figure 2 .
Figure 2. Model evolution results showing the dynamics of diapirs and the migration of subducting sediments (a-c) in the Reference Model with sediment thermal parameter (STP) of 1.95, (d) in model MC8S05 with STP of 2.40, and (e) model MC2S4 with STP of 1.38.The first and third columns are high-resolution composition and temperature fields of the studied domain.The second column shows the CO 2 content retained in the sediments, which is enlargement of the area in the dashed box in the first column.The blue arrows indicate the migration direction of sedimentary diapirs.The white dashed lines constrain the forearc depths (15-70 km) and subarc depths (70-150 km).The left boundary of the subarc lithosphere is limited by the yellow dashed lines.The closed yellow curves represent diapiric sediments/CO 2 migration to the subarc lithosphere.

Figure 3 .
Figure 3.Effect of sediment thermal parameter (STP) on decarbonation efficiency at forearc and subarc depths for the modeled and observed data.The figure includes two observations in the forearc, and eight studies reporting on a large proportion of sedimentary carbon in arc volcanism.The uncertainties correspond to the range of the decarbonation efficiency obtained from the observations.The blue and yellow lines are linear fits for our modeled decarbonation efficiency of forearc diapirs and subarc diapirs, respectively.The numbers (in order of increasing STP) correspond to the following subduction zones: 1 = Cascadia, 2 = Cyclades (Greece), 3 = SW Japan, 4 = Alaska, 5 = Sumatra, 6 = Antilles, 7 = Aleutians, 8 = Guatemala-EI Salvador, 9 = Java, and 10 = Costa Rica.The lowercase letters indicate the additional test models: a = model MST2, b = model MA8S2, c = model MA8C8, and d = model MC8S05.

Figure 4 .
Figure 4. Schematic diagram showing the transport of forearc diapirs and the possible carbon reservoirs in the forearc mantle-lithosphere.The white and black arrows indicate the subduction direction of the oceanic slab and the migration direction of the carbonate sediments, respectively.The boldness of the black arrows denotes the proportion of carbonate rocks.① Shows the transport of forearc diapirs to the subarc lithosphere and exhumation of the diapiric sediments within the subduction channel in the left corner; the subarc mantle-lithosphere, and crust beneath arc represented by ② are two possible carbon reservoirs.