Uncovering and quantifying the subduction zone sulfur cycle from the slab perspective

Sulfur belongs among H2O, CO2, and Cl as one of the key volatiles in Earth’s chemical cycles. High oxygen fugacity, sulfur concentration, and δ34S values in volcanic arc rocks have been attributed to significant sulfate addition by slab fluids. However, sulfur speciation, flux, and isotope composition in slab-dehydrated fluids remain unclear. Here, we use high-pressure rocks and enclosed veins to provide direct constraints on subduction zone sulfur recycling for a typical oceanic lithosphere. Textural and thermodynamic evidence indicates the predominance of reduced sulfur species in slab fluids; those derived from metasediments, altered oceanic crust, and serpentinite have δ34S values of approximately −8‰, −1‰, and +8‰, respectively. Mass-balance calculations demonstrate that 6.4% (up to 20% maximum) of total subducted sulfur is released between 30–230 km depth, and the predominant sulfur loss takes place at 70–100 km with a net δ34S composition of −2.5 ± 3‰. We conclude that modest slab-to-wedge sulfur transport occurs, but that slab-derived fluids provide negligible sulfate to oxidize the sub-arc mantle and cannot deliver 34S-enriched sulfur to produce the positive δ34S signature in arc settings. Most sulfur has negative δ34S and is subducted into the deep mantle, which could cause a long-term increase in the δ34S of Earth surface reservoirs.

S ulfur is one of the most common volatiles on Earth. It plays key roles in, for example, the redox evolution of the sub-arc mantle 1,2 , the formation of ore deposits 2 , and the composition of the atmosphere through volcanic SO 2 degassing 3 . Subduction zones are the primary locations for the global sulfur cycle, transporting sulfur to the deep mantle via the descending slab or returning it to the surface by arc magmatism 2,4,5 . Compared to fresh MORB, the relatively high sulfur concentrations ([S], up to 3000 µg g −1 ) and positive δ 34 S values (+5 to +11‰) of volcanic rocks and melt inclusions in some arcs (e.g., Western Pacific) [4][5][6][7] , and the presence of sulfate in mantle xenoliths 8 , have been attributed to the addition of slab-derived sulfate to arc magmas by fluids 8,9 . Alternatively, some deep arc cumulates (e.g., Eastern Pacific) with mantle-like δ 34 S values suggest more limited slabderived sulfate contributions to arc lavas and that the positive δ 34 S signature of the lavas results from crustal assimilation 10 .
The role of slab fluids in delivering sulfur species to the mantle wedge is central to this debate. Experimental results suggest that slab-derived aqueous fluids are an effective agent for transporting sulfur from the slab to the mantle wedge 11,12 . In addition, some studies predict that sulfates are likely the dominant sulfur species in slab-derived fluids 9,13 . On the other hand, sulfate is relatively rare in high-pressure (HP) rocks [13][14][15][16][17][18] , and experimental studies have proposed that reduced sulfur species are dominant in slab fluids 11 . Furthermore, in situ measurements of the δ 34 S compositions of sulfides from HP eclogites and serpentinites reveal significant isotopic heterogeneity and complicated sulfur behavior during slab metamorphism and metasomatism [13][14][15][16] .
Clearly, large gaps in our knowledge of the speciation, flux, and isotopic composition of sulfur in slab fluids remain. Understanding these is of utmost importance for addressing slab-arc sulfur recycling, and has global geochemical significance for deciphering the redox state of the mantle 2 and constraining the formation of arcrelated ore deposits 12 . Direct examination of devolatilization pathways in exhumed HP rocks is essential to provide independent new perspectives critical to resolving this debate, as it provides the necessary field-based evidence for the sulfur redox state and δ 34 S signature of fluids released from subducted slabs.
Sulfur is transported into the subduction zone by sediments, variably altered oceanic crust (AOC), and hydrated slab mantle (serpentinites) [19][20][21][22] . Sulfides are commonly observed in exhumed fragments of the oceanic lithosphere such as eclogites, blueschists, HP-metapelites, and serpentinites, as well as related HP veins [13][14][15][16][17][18]23 . Such vein systems represent fossilized pathways for channelized flow of dehydration-related slab fluids and, thus, directly record fluid geochemical signatures 17,24,25 . Consequently, the study of HP vein-rock systems provides important information regarding sulfur behavior during slab dehydration and fluid transfer 17 . Isotopic constraints on S-bearing HP rocks and veins linked to the sequence of slab dehydration allow quantification of sulfur release during subduction of oceanic lithosphere.
Here, we report bulk-rock and in situ sulfur isotope compositions for sulfide-bearing HP rocks and veins from the late Paleozoic southwestern Tianshan (ultra-)high-pressure/low-temperature ((U)HP/LT) metamorphic belt (China). The sulfides in these HP rocks and veins 17 provide an exceptional window into the fate of subducted sulfur. Analytical data and thermodynamic calculations point to low sulfur concentrations in slab fluids, which have negative δ 34 S values and are dominately composed of reduced sulfur species. Hence, we determine modest slab-to-arc sulfur transport, and find neither significant slab sulfate flux to the mantle wedge nor a direct link between slab-derived sulfur and the positive δ 34 S signature of arc settings.
In addition, three representative sulfide-bearing dehydrationrelated veins in blueschists or eclogites were investigated to reconstruct the sulfur behavior in subduction fluids derived from different sources (Fig. 1). Vein_1 (JTS) consists of a well-studied wallrock-selvage-vein system 25,29 formed by fluid-rock interaction during prograde metamorphism (Fig. 1a) dominated) and a blueschist-eclogite transition zone due to reaction with an external fluid (Fig. 1a). The wallrock-selvage-vein system equilibrated at peak metamorphic conditions of ∼510°C and 2.1 GPa 29 . Sr and Ca isotope compositions trace the fluid source to seawater-altered lithospheric slab-mantle and/or oceanic crust 25 . Vein_2 (L1422) is a 2-cm-wide dolomite-quartz-epidote-dominated vein crosscutting a massive host blueschist (Fig. 1b). Similar to Vein_1, the blueschist-eclogite transition zone and eclogite selvage formed due to interaction with Ca-rich fluid along the conduit (Fig. 1b). The similar structure, mineral assemblages, and compositions of Vein_1 and Vein_2 indicate that they formed at similar P-T conditions. Vein_3 (L1013) is a 1-3 cm wide dolomite-quartz-epidote-dominated vein cutting massive eclogite (Fig. 1c). Occurrences of high-pressure minerals such as omphacite and rutile in Vein_3 along with the observation that no reaction halo occurs between the vein and host eclogite (Fig. 1c) indicate that the vein also formed at eclogite-facies conditions. Considering most eclogite samples in the Tianshan HP metamorphic belt were exhumed from ∼80 km depth 26 , all three HP veins are thought to represent the fluid activity that took place at 70-90 km depths in the subduction zone. Sulfides are found in all vein samples. In sample JTS, pyrite is the dominant sulfide inclusion in garnet, but pyrrhotite dominates in the matrix (Fig. 2a, b)  during garnet growth due to changing sulfur fugacity-oxygen fugacity (fS 2 -fO 2 ) conditions associated with fluid metasomatism. Sulfide in the vein is mostly pyrite but also includes minor pyrrhotite (Fig. 2c). Along the traverse of samples L1422, sulfide (mainly pyrite) abundances increase toward the vein. The Co-Ni element distribution maps (Method 1) reflect distinct differences between host blueschist pyrite ( Fig. 2d-f) and selvage-vein pyrites ( Fig. 2g-i). Both selvage and vein pyrites display multiple growth generations ( Fig. 2g-i). In sample L1013, Co-Ni element distribution maps and contents show core-rim textures in both eclogite (Fig. 2m) and vein pyrite ( Fig. 2n-p).
In most samples, fine fractures were occasionally observed surrounding sulfide grains, reflecting rigidity contrasts between sulfide and matrix minerals. These fractures are usually filled with albite + magnetite + calcite ± chalcopyrite ± barite due to latestage fluid infiltration, accompanying variable retrogression of neighboring omphacite and glaucophane. Sulfate and magnetite were observed only in these late-stage, retrograde fractures.  Table 1) and the metabasites have [S] WR = 841-3978 µg g −1 and a range in δ 34 S WR of −7.2 to +3.6‰, averaging −2.7‰ (n = 5). In contrast, all measured serpentinites ([S] WR = 124-422 µg g −1 ) have positive δ 34 S WR values (+3.6 to +12‰), signifying high-temperature waterrock interaction during oceanic serpentinization 30 . These results are similar to unmetamorphosed oceanic lithosphere [19][20][21]30 and suggest that the main slab sulfur reservoirs have distinct δ 34 S WR compositions, which are generally consistent with previous studies of exhumed slab rocks 14,15 (Fig. 3 Sulfur geochemistry in HP vein systems. In the host blueschist to blueschist-eclogite transition zone of Vein_1 the [S] WR varies between 698 µg g −1 and 841 µg g −1 , but toward the vein and veinlike eclogite selvage, the [S] WR increases up to 2183 µg g −1 (Fig. 4a). In contrast, the δ 34 S WR values decrease gradually from +0.43‰ to −0.98‰ toward the vein (Fig. 4a). The in situ δ 34 S values of sulfides in garnet show little variation; however, matrix sulfides display a decreasing trend from host rock toward the vein ( Fig. 4b; Method 1). The local bulk isotopic compositions of sulfides (δ 34 S sulfide ), calculated using mean in situ δ 34 S values of individual sulfides (Fig. 4b) and their mineral volume ratios ( Supplementary Fig. 1), have a narrow range of +0.34-+0.72‰ (mean +0.60‰) in the garnet along the traverse (Fig. 4b). This reflects shielding of the sulfide inclusions by garnet during fluid-rock interaction. In contrast, the δ 34 S sulfide in the matrix decreases gradually from about 0.00‰ to −1.35‰ toward the vein (Fig. 4b). The δ 34 S WR and δ 34 S sulfide display similar decreasing trends along the traverse. Vein sulfides have uniform δ 34 S sulfide values of about −1.0‰ (Fig. 4b).
The vein pyrite of Vein_3 has a thick Co-poor but Ni-rich core (δ 34 S = −8%) and a thin Co-rich but Ni-poor rim (δ 34 S = −5‰) (Figs. 2n-p, 5f). In contrast, pyrite in the host eclogite contains a Co-Ni-rich core with MORB-like δ 34 S values (−1.3 to −0.5‰), but its rim is analogous to the vein pyrite core in Co-Ni-As contents, δ 34 S values (about −8%), and mineral inclusions (omphacite and rutile) (Figs. 2m, 5e). This indicates that the vein-forming fluid also altered the immediate eclogite and caused pyrite regrowth surrounding the cores.
Sulfur concentrations in aqueous fluids from DEW modeling. The sulfur concentration in fluids ([S] fluid ) is the most important factor determining the slab sulfur output, as aqueous fluids are thought to be the major agent for slab-mantle sulfur transfer 11 . We use the DEW (Deep Earth Water) model 31,32 to calculate subduction zone [S] fluid (Method 2), as this allows a quantitative prediction of speciation and solubilities in fluids at upper mantle conditions. Because sulfur solubility and speciation is redoxdependent, an estimate of the fO 2 is required prior to calculation. The fO 2 of subducted AOC is FMQ + 1 (ref. 33 ) (one log unit above Fayalite-Magnetite-Quartz buffer) at the trench and decreases gradually with increasing depth (below FMQ at depths corresponding to eclogite-facies conditions) 17 , as generally reducing fluids (<FMQ) are generated 17,34 . In contrast, the redox  Errors (2σ) of δ 34 S WR and S WR calculated from measurement reproducibility are less than the symbol size. Sulfur isotope values of HP Tianshan (TS) rocks are from this study, Franciscan (FS) eclogites and blueschists from ref. 10 , high-pressure serpentinites of the Voltri Massif (VM, Italy) from refs. 47,48 , high-pressure serpentinites of the Cerro del Almirez (CA, Spain) from ref. 21 , and the mantle value (−0.91 ± 0.50‰) from ref. 45 . The shaded area indicates δ 34 S range of metabasites, which is comparable to the mantle value. Source data are provided in Supplementary Data. state of slab serpentinite is more complicated (either above or below FMQ) and is suggested to vary due to different degrees of pre-subduction serpentinization 35 , producing both highly oxidizing or reducing fluids 35 . Dehydration of incompletely serpentinized rocks (usually those beneath oceanic crust) in which awaruite is present produces reducing H 2 -bearing fluids, whereas deserpentinization of completely serpentinized rocks (usually those once directly exposed to seawater) in which awaruite is absent produces oxidizing fluids in the subduction zone 35 . The former is applicable in our case, as the majority of slab mantle occurs beneath oceanic crust and is not fully serpentinized. For details regarding slab fO 2 estimates see the Supplementary Note 1.
Hydrothermal sulfur isotope fractionation. Effects of sulfur isotope fractionation during hydrothermal processes are largely influenced by pressure, temperature, fO 2 and the pH of the fluids 40,41 . We thermodynamically calculated fO 2 -pH diagrams 42 (Fig. 7a-d) (Method 3) for different P-T conditions to reveal δ 34 S fractionation in subduction zones. Our DEW calculations suggest that slab fluids are generally alkaline and the pH value ranges from neutral (pH n ) to pH n + 2, consistent with previous work 43 . According to fO 2 (Supplementary Table 2) and pH estimates of subduction zone fluids, all the fO 2 -pH conditions plot in fields dominated by the species H 2 S or HS − (yellow area, Fig. 7a-e), in agreement with DEW results (Supplementary Fig. 2). The fO 2 -pH diagram indicates limited sulfur isotope fractionation  In particular, at the vein-forming P-T conditions of this study, fO 2 (<FMQ) 17 and pH range (pH n to pH n + 2) suggest sulfur isotope fractionations <1.3‰ (Fig. 7e). In addition, precipitation styles of sulfide in hydrothermal settings (closed-or open-system) may also influence the sulfur isotope fractionation 41 . 3) are similar to their protoliths, the young marine sedimentary rocks (Phanerozoic) that mostly have δ 34 S values of −24 to −8‰ due to the presence of biogenically produced sulfide 19,44 . It is suggested that sulfide in sediments retains its δ 34 S characteristics during subduction metamorphism 13 and that metasediments may act as a negative δ 34 S reservior in subducting slabs. The vein pyrite core with negative δ 34 S (−8‰) in sample L1013 (Fig. 2n) thus likely represents the δ 34 S signature of fluids derived from the abundant subducted metasediments in the Tianshan HP belt 16,26 . Pristine oceanic crust is typically within the range of the average mantle δ 34 S of −0.91 ± 0.50‰ 45 . However, presubduction seafloor processes lead to considerable δ 34 S heterogeneities in the upper crust (Fig. 3, Supplementary Table 1 15 . However, these sulfides are associated with blueschist/greenschist retrogression 15 , which may record the oxidizing fluids at shallow depth that usually cause retrogression of exhumed eclogites/blueschists 17 . The positive δ 34 S may originally come from seawater hydrothermal alteration, or represent the fluids derived from serpentinite dehydration (see below). Therefore, although pre-subduction seafloor processes cause negative or positive δ 34 S shifts in metabasites, bulk-rock geochemistry shows that many mafic eclogites and blueschists still retain their mantle-like sulfur isotope signature throughout HP metamorphism (Fig. 3). Furthermore, the lower crust (dike and gabbro) retains its mantle-like δ 34 S values as well 46 .  (-1)   (Fig. 3) and in situ sulfide δ 34 S compositions from Corsican serpentinites 14 . We suggest that variably serpentinized slab mantle beneath oceanic crust is characterized by positive δ 34 S compositions as a result of sulfide addition via sulfate reduction at high-temperatures 30

Discussion
Thermodynamic modelling shows that at subduction zone P-T-fO 2 -pH conditions sulfur in fluids is dominated by the reduced H 2 S and HS − species, whereas sulfate species (e.g. SO 4 2− , HSO 3 − ) are rare ( Supplementary Fig. 2). This is consistent with our petrological evidence for the occurrence of sulfide, but not sulfate, in the veins (Figs. 1-2), the very low sulfate concentrations in rocks and veins ( Fig. 4a; Supplementary Table 1), and previous experimental results 11 . If slab fluids are dominated by sulfate as some recent studies propose 13 , several predictions follow. First, the oxidizing fluid will produce a redox gradient in the immediate wallrock, but this is not recorded in the selvages we examined. Second, reduction from S 6+ to S 2− will cause oxidation within the immediate rock and vein (in particular during vein crystallization) in the form of hematite or magnetite 17 , which, however, are absent from the veins or selvages. Third, complete sulfate-sulfide transformation will produce very high δ 34 S values in the product phases, which is also not observed in the veins or   selvages. The sulfate introduced during pre-subduction hydrothermal seafloor alteration 20,21,30 may have been lost or converted to sulfide at early stages of subduction, for example, at fore-arc depths 10,12 . Thus, we conclude that the dehydration-related slab fluids likely transport reduced sulfur species such as aqueous H 2 S and HS − at sub-arc depths. Mass-balance calculations were used to estimate the sulfur influx (F S ) and outflux (f S ) of subduction zones, as well as their δ 34 S values. The sulfur input estimate was computed (Method 4) based on the average [S] and δ 34 S compositions of our best current understanding of oceanic lithosphere stratigraphy (Supplementary Fig. 3, Note 2) in combination with the global length of subduction zones and their average convergence rate, sequence thickness, and density. The resulting subduction zone sulfur influx is estimated to be 4.65 × 10 13 gram per year (g yr −1 ) with a bulk negative δ 34 S value of −3.60‰. Gabbro (49%) and sediment (23%) are two important sulfur reservoirs (Fig. 8a), whereas serpentinite is insignificant due to its low [S]. This influx is slightly lower than, but generally within the same order of magnitude of, previous estimates 2, 19,21,22 (Fig. 8b).
The sulfur output can be calculated from the product of [S] fluid and the fluid fluxes from the dehydrating slab. Previous thermodynamic constraints 37-39 on [S] fluid predict a rather low proportion of H 2 S (<0.01 mol.% in equilibrium with H 2 O + pyrite + pyrrhotite) 38,39 . This low [S] fluid , in turn, predicts negligible slab sulfur output (f S /F S < 1%), which is not compatible with the high sulfur contents observed in arc settings. Based on our novel [S] fluid results from the DEW calculations (Fig. 6a) and depthdependent fluid fluxes from the dehydrating slab 49 (Fig. 6b), a total sulfur outflux (Method 5) of 2.91 × 10 12 g yr −1 (6.3% of F S ) is estimated via fluids derived from the subducting slab between 30-230 km (Fig. 8a). The sheeted dikes (31%) and the upper volcanic layer (26%) contribute most of the sulfur release, but the sediment (15%) and the lower volcanic layer (15%) are also important (Fig. 8a). This sulfur outflux from the slab is about 1/5 to 1/3 of the sulfur output estimates from arcs 2,50,51 (Fig. 8b).
Importantly, our sulfur output estimate shows a major sulfur release of 2.46 × 10 12 g yr −1 (5.3% of F S ) to the mantle wedge at depths of 70-100 km (Fig. 9) due to both elevated [S] fluid (Fig. 6a) and released H 2 O flux 49 (Fig. 6b). The volcanic layers, dikes, gabbro, and sediments contribute to the major sulfur release at this depth interval (Fig. 8a). This sulfur release window coincides with pyrite-to-pyrrhotite breakdown 12 (releasing H 2 S) and the major slab fluid release (∼32% H 2 O/H 2 O total ) 49 , which subsequently acts as the trigger for partial melting of the mantle wedge and ultimately arc magmatism 25,52,53 . The release of slab sulfur is minor at other depths (Fig. 9). The DEW model may underestimate [S] fluid due to unknown sulfur species. In addition, considering the uncertainties on the thermodynamic data (±0.3 to ±0.5 of logK) 31 , the slab sulfur loss estimate may extend to a maximum of ∼20%. This is comparable to the estimate obtained from natural rocks and experiments (30%, Supplementary Note 3), and previous estimates (8-18%, Fig. 8b) from arcs 2,50,51 .
Knowing [S] fluid (Fig. 6a) and water fluxes 49 (Fig. 6b) liberated from slab sequences as well as the δ 34 S fluids of sediments (−8‰), oceanic crust (−1‰) and serpentinites (+8‰), mass-balance calculations (Method 5) predict a δ 34 S value of slab fluids of −1.8 ± 3‰ at depths of 70-100 km. Both the gradually decreasing δ 34 S WR and δ 34 S sulfide values (sample JTS, Fig. 4a) and the trend of increasing δ 34 S WR values within the studied profile (sample L1422, Fig. 5a) provide clear evidence for isotope exchange during fluid-rock interaction, despite rapid fluid transport rates 24,25 . The sulfur isotope composition of channelized fluids generated in the lower parts of a slab could be slightly altered within the upper parts during fluid transport. Assuming a 10% sulfur isotopic contamination rate by sediment and 20% by oceanic crust for underlying sulfur sources, the net δ 34 S value is calibrated to −2.5 ± 3‰ at depths of 70-100 km (Fig. 9). If δ 34 S fluids values from different reservoirs are kept constant at different depths, δ 34 S of slab fluids at shallower depths remains negative but shifts to positive at greater depths, and a net δ 34 S of −2.1‰ is obtained for slab fluids between the whole 30-230 km (Fig. 9). Regardless of flux uncertainties, we emphasize that the negative δ 34 S value of slab fluids is robust and insensitive to a range of different model scenarios (Method 6), including thick serpentinized slab-mantle scenarios 54 , which may supply more fluids with positive δ 34 S values.
This study provides the first comprehensive and quantitative view of the flux and isotope compositions of sulfur-bearing slab fluids, which likely mirror the slab sulfur contributions to arcs 25,55 . We show that dehydration-related fluids transfer modest amounts of sulfur (6.4% of total subducted sulfur, up to 20% maximum)  from the slab to the mantle wedge. This maintains elevated sulfur contents of the mantle source for arc magmas (250-500 µg g −1 ) 6,7 in comparison to MORB (80-300 µg g −1 ) 56 . Additional significant release by, for example, slab melting is unlikely, as [S] in melts is much lower than in aqueous fluids (D S fluid/melt usually >200) 11 . This slab-arc sulfur cycle is operated by fluid-mediated H 2 S and/or HS − transport with negative δ 34 S composition, which has no direct links to the high oxygen fugacity and heavy δ 34 S signature observed in arc volcanic rocks. Our work also sheds further light on the nature of arc magmas. The reason for the higher fO 2 of arc magmas 57,58 relative to MORBs is still debatable. Intraoceanic or rare continental arcs, like those of the Western Pacific, may record flux melting; mantle peridotites with elevated fO 2 values in these settings have been thought to be influenced by slab-derived oxidizing agents 59,60 . In contrast, continental arcs, like those of the Eastern Pacific arcs where crust thickness may modulate the melting degree 61 , may represent a complicated melting mode involving decompression and mantle peridotite that is not necessarily oxidized 59,62 . Direct comparisons of Western and Eastern Pacific arcs may be challenging due to their different melting modes and arc maturity 59,60 . The P-T evolution of the Tianshan eclogites (representing cold/old subducted oceanic slab) corresponds more closely to the thermal structure of subduction zones beneath the Western Pacific arcs 59 . However, our finding of negligible sulfate in the slab fluids indicates that slab SO 4 2− was unlikely to be the main oxidizing agent during South Tianshan Ocean subduction. In such environments, the high fO 2 of sub-arc mantle may instead result from addition of slab H 2 O and CO 2 (refs 63,64 ), instead of oxidized sulfur species. Processes including incorporation of H 2 into orthopyroxene 63 and the formation of diamond 64 and CH4 (ref. 65 ) in the mantle wedge may produce oxidized melts that elevate the fO 2 of Western Pacific arc magmas .
The calculated negative δ 34 S (−2.5 ± 3‰) released from the subducted slab (Fig. 9) contrasts with the positive δ 34 S values found in the inclusions of Western Pacific arc rocks 4,5,8 . In general, the mantle wedge should have mantle-like δ 34 S values of ∼0‰. Therefore, the positive δ 34 S signature in arc-related rocks requires additional sulfur sources or processes for 34 S enrichment. Volcanic degassing effects on melt δ 34 S are highly dependent on redox state 66,67 . But even under oxidizing conditions (>FMQ + 2), increases in melt δ 34 S caused by degassing are modest (∼1.5‰) 66 . Therefore, the negativeto-positive shift in the δ 34 S composition of melts must happen as part of the partial melting processes, such that significant sulfur isotope fractionation accompanies the melt oxidization. For example, 32 S may be scavenged into surrounding mantle to form sulfides while H 2 is incorporated into orthopyroxene 63 , producing 34 S-rich sulfate in oxidizing melt and finally isotopically heavier arc magmas. Thus, ARTICLE further studies will be necessary to assess the processes that may lead to the positive δ 34 S compositions in arc magmas.
Our comparison of subduction input with output fluxes indicates that most of the sulfur (>80%) with negative δ 34 S values (<−3.7‰) is retained in the descending slab and recycled to the deep mantle (Fig. 9). This may have resulted in a progressive 34 Senrichment of Earth's surface sulfur reservoirs 19 , and can explain the negative δ 34 S values of alkaline magmas related to ocean island basalts (OIBs) since the Phanerozoic 68 .

Methods
Analytical methods. Bulk-rock sulfur contents and isotope compositions were measured at the Geological Institute at the Freie Universität Berlin. Extraction of the bulk-rock sulfur was performed by extracting the acid volatile sulfide (AVS), chromium reducible sulfide (CRS), and the sulfate fraction 69 . Sulfur isotope measurements of AVS, CRS, and sulfate fractions were done on a Thermo Fisher Scientific MAT 253 mass spectrometer combined with a Eurovector elemental analyzer. The [S] WR of individual samples were calculated by summing sulfur amounts of measured AVS, CRS, and sulfate. The δ 34 S WR was calculated by measured δ 34 S values of AVS, CRS, and sulfate in combination with their amounts. In situ sulfur isotopes of sulfides on epoxy discs were analyzed via Secondary Ionization Mass Spectrometry (SIMS) using a Cameca IMS 1280 instrument located at the Swedish Museum of Natural History, Stockholm, Sweden (NORDSIM facility) 70 for sample JTS and at the Institute of Geology and Geophysics, Chinese Academy of Sciences (IGGCAS, Beijing, China) 71 for other samples. Measurements were conducted over a rastered 10 × 10 μm area using a 133 Cs + primary beam with 20 kV incident energy (10 kV primary, −10 kV secondary) and a primary beam current of ∼1.0 nA. All δ 34 S results are reported with respect to the V-CDT standard 72 . Detailed descriptions of δ 34 S measurement parameters and standard references are given in the Supplementary Note 4. Elemental Co and Ni X-ray maps of pyrite were made in wavelength-dispersive spectrometer mode by electron microprobe (CAMECA SXFive FE) at the IGGCAS. An acceleration voltage of 20 kV, beam current of 100 nA, 3-5 μm pixel size, and dwell time of 50 ms were used. In situ trace-element analyses by laser ablation inductively-coupled plasma mass spectrometry of sulfides were made on thin sections in the GeoZentrum Nordbayern of the University Erlangen-Nürnberg, Erlangen, Germany.
Reactions and equations of sulfur species include H 2 S (aq) , HS − , HSO 4 − , and SO 4 2− , and the relative isotopic fractionation (Δ i = 34 S i -34 S H2S ) for sulfur species were calculated at higher temperatures according to the equation 40 : where a, b, and c are empirically-determined constants, and T is temperature in Kelvin. The equilibrium constants for reactions and activity coefficients of aqueous species were recalculated for higher P-T conditions based on the DEW model 31,32 . The abundance of sulfur species and contours of sulfur isotope fractionation (compared to initial sulfur isotope of hydrothermal fluid at δ 34 S ∑S = 0‰) as functions of fO 2 and pH were calculated using Eqs. (17)(18)(19)(20)(21)(22)(23)(24)(25) listed in ref. 42 .
The changes in isotope fractionation during pyrite crystallization from slab fluids (H 2 S dominated or SO 4 2− -dominated) in closed system and open-system (Rayleigh fractionation) processes were modeled at the vein-formation temperature 550°C. In a hydrothermal system, the isotope composition of an instantaneously separated solid phase i from a fluid is: where δ 34 S 0 is the initial fluid isotope composition (set as 0‰) and F is the fraction of sulfur remaining in the fluid.
Estimate of global sulfur input into subduction zones. The input sulfur flux F S and its isotope composition into subduction zones are: where L is the global length of subduction zones, R is the convergence rate, t is the thickness of the sequence layers in the slab, ρ is the density of the sequence layers, C S is the sulfur concentration [S], and δ 34 S t is the sulfur isotope composition of the sequence layers. The total effective length of subduction zones is~38,500 km, which covers more than 90% of global trench length 49 . The convergence rate of 6.2 cm yr −1 used here is taken from an average rate of 17 active oceanic subduction zones 36 . Based on the oceanic lithosphere stratigraphy (Penrose style) and its average [S] and δ 34 S composition from the best current understanding (Supplementary Fig. 3), the calculated global sulfur input via subducting slabs is estimated to be 46.5 × 10 12 g yr −1 . The bulk slab sulfur isotope composition of this sulfur input is estimated at −3.60‰. Using the same method, the calculated global water flux (1.06 × 10 15 g yr −1 ) of subducted slabs is very close to previous estimates (1.0 × 10 15 g yr −1 ) 49 .
Sulfur output and net δ 34 S released by slab fluids. The output sulfur flux released from the subducted slab via fluids (f S ) is: f S ¼ C SÀfluid Á f fluid and the net sulfur isotope composition of fluids released from the subducted slab (δ 34 S net ) is: where C S-fluid refers to [S] fluid and f fluid to the fluid flux released from the subducting slab. Based on water flux (0.32 × 10 15 gyr −1 ) and [S] fluid from the DEW model, the calculated sulfur output at 70-100 km is 2.46 × 10 12 g yr −1 (5.3% of total input F S ) with a δ 34 S value of −1.84 ± 3 ‰. The net δ 34 S value of slab fluids released at 70-100 km depths is further adjusted to −2.54 ± 3 ‰ (Fig. 9) considering fluid-rock isotopic exchange.
Uncertainties on output sulfur δ 34 S estimates. The estimates of sulfur fluxes released from the slab have significant uncertainties. However, our study provides a robust isotopic signature for the slab fluids. The δ 34 S estimate remains at slightly negative values in all of the following scenarios: The uncertainty of δ 34 S net is mostly dependent on the δ 34 S value of fluids released by the AOC at 70-100 km, which provides the major fluid flux and has a relatively high sulfur concentration (0.74 wt.%). Although we consider a large δ 34 S range of fluid AOC (−6 to +4‰), the errors on the δ 34 S of slab fluid released at 70-100 km are all less than ±2.5‰ (2σ). Changes in other parameters and assumptions cause variations of less than ±1‰ in δ 34 S. Hence, we estimate ±3‰ as a reasonable uncertainty.
Our study is based on the best current knowledge of slab structure and water budget 49 . However, new research based on ocean-bottom seismic data reports that mantle hydration may extend up to 24 km beneath the Moho 54 , which indicates that the subducting plate may contain much more water than previously thought 49 . If we adopt this assumption of a thicker serpentinized upper mantle 54 and recalculate the water and sulfur fluxes (i.e., enlarged the serpentinite-dehydrated water amounts in the subduction zone accordingly), the sulfur input increases to 7.6 × 10 13 g yr −1 and sulfur output increases to 3.93 × 10 12 g yr −1 but the sulfur productivity (5.2%) of the subducting slab shows little variation. More importantly, the net δ 34 S displays almost no change at 70-100 km (−2.4‰). This consolidates our prediction of slab-released sulfur regarding the δ 34 S signatures of arc settings.
Subducted sediment types and their redox state may have a potential effect on our results, even though there currently is no firm consensus about how much sediment is subducted. Metasediment in the Franciscan complex contains red ferruginous chert, but its proportion is subordinate compared to greywackes 77 . Moreover, at the major sulfur release window (70-100 km), the sediment contribution to the total sulfur loss is small (less than 20%). The channelized fluids 24,25 in subducted slabs (instead of pervasive fluids) prevent intensive sediment contamination of deeply-derived fluids. Here we conclude that the variable sediment protoliths and redox state will not influence our results significantly.

Data availability
The source data underlying Figs

Code availability
The computer code including the Deep Earth Water model (2019)