Oxidation of Archean upper mantle caused by crustal recycling

The redox evolution of Archean upper mantle impacted mantle melting and the nature of chemical equilibrium between mantle, ocean and atmosphere of the early Earth. Yet, the origin of these variations in redox remain controversial. Here we show that a global compilation of ∼3.8-2.5 Ga basalts can be subdivided into group B-1, showing modern mid-ocean ridge basalt-like features ((Nb/La)PM ≥ 0.75), and B-2, which are similar to contemporary island arc-related basalts ((Nb/La)PM < 0.75). Our V-Ti redox proxy indicates a more reducing upper mantle, and the results of both ambient and modified mantle obtained from B-1 and B-2 samples, respectively, exhibit a ∼1.0 log unit increase in their temporal evolution for most cratons. Increases in mantle oxygen fugacity are coincident with the changes in basalt Th/Nb ratios and Nd isotope ratios, indicating that crustal recycling played a crucial role, and this likely occurred either via plate subduction or lithospheric drips.

R edox state of the Archean upper mantle buffered atmospheric composition and influenced the early oceanatmosphere system through the geochemical behavior of volatile elements (e.g., S, C and N) [1][2][3] . Previous investigations of the mantle redox state have focused on the Fe 3+ /ΣFe content in pristine mantle and mantle-derived rocks 1,4,5 , the behavior of redox-sensitive elements (e.g., V) in komatiites 6 , and the V-Sc (Sc is a typical redox-insensitive element) redox proxy in Archean basalts 7,8 . These methods provide an in principle estimate of oxygen fugacity (fO 2 ; relative to the fayalite-magnetite-quartz buffer (ΔFMQ)) of the upper mantle. However, in terms of Archean samples, the Fe 3+ /ΣFe content of minerals and glasses might be readily modified during alteration, metamorphism and syn-eruptive process 6,9,10 . The behavior of V in komatiites shows a significant increase of~1.3 log units in mantle fO 2 betweeñ 3.5 Ga and 1.9 Ga 6 . This change, however, might not reflect the nature of the upper mantle because komatiites constitute less than 10% of exposed Archean volcanic rocks, and their magmatic sources are thought to be complicated with materials from either the core-mantle or the upper-lower mantle boundaries 11,12 . Compared to komatiites, basalts constitute the dominant component of Archean continental crust (~65-75 %) and can provide alternative insights into resolving fO 2 of the mantle 1,8,13 . Some studies indicated that the mantle fO 2 has been at present-day levels since~3.9 Ga 7 , but others inferred a dominantly reducing upper mantle on the basis of the significantly lower V/Sc ratios in Archean mid-ocean ridge basalt (MORB)-like samples relative to their modern counterparts 1,8 . The reason of these contrasting interpretations is attributed to choices of different mantle melting reactions and partition coefficients (Kds) in converting V/Sc ratios to mantle fO 2 1,7,8 . Our recent study found that the Kds of V, Sc and Ti are influenced not only by melting temperaturepressure (P-T) and fO 2 , but also by the compositions of mantle mineral and water contents 14 , which have not been fully considered in previous fO 2 estimations. Thus, suitable Kds and mantle melting reactions are critical for the accurate estimation of fO 2 . In addition, although Archean MORB-like basalts or eclogites are increasingly used in mantle fO 2 calculations, limited investigations of basalts with characteristics indicative of lithospheric recycling (e.g., modern island arc basalt (IAB)) have impeded a comprehensive understanding of the early Earth's mantle-crust interactions [15][16][17][18] .
We have assembled a global whole-rock geochemical database of Archean (~3.8-2.5 Ga) basalts (N = 2304; See more details of data collection and filtration in Methods) from fourteen cratons ( Supplementary Fig. 1). They exhibit a significant compositional diversity, and are divided into Basalt-1 (B-1) and Basalt-2 (B-2) groups with modern MORB-and IABlike chemical features, respectively. We applied an updated V-Ti redox proxy, which may be more sensitive than V-Sc systematics because of the more incompatible character of Ti during mantle melting and the influence of Sc content by garnet residues 14 , to this database. The results from applying this proxy are: (1) The average V/Ti ratios of Archean basalts are relatively lower than those of modern counterparts, and when converting them to mantle fO 2 , the Archean upper mantle was more reducing than it is today; (2) There is a marked difference in redox state between ambient (ΔFMQ −1.31 ± 0.58 (1 SD)) and modified mantle (ΔFMQ −0.88 ± 0.84 at 2 wt.% H 2 O and ΔFMQ −0.39 ± 0.89 at 4 wt.% H 2 O) obtained from B-1 and B-2 magmas, respectively; (3) Both the ambient and modified mantle have commonly undergone significant oxidization with 1.0 log unit increase of most Archean cratons and possible craton groups; and (4) Changes in mantle fO 2 are closely associated with changes in basalt Th/Nb ratios and Nd isotope ratios that are sensitive to crustal recycling.

Results
Geochemical characteristics of Archean basalt. Based on the established tectonic settings, modern basalts are commonly divisible into non-subduction and subduction-related types. The former records the ambient mantle conditions via adiabatic decompression melting, and the latter indicates fluid flux melting due to crustal/lithospheric recycling 19 . Our statistical analysis of modern basalts suggests that (Nb/La) PM ≥ 0.75 (PM, primitive mantle normalized 20 ) is a powerful criteria identifying primitive melts that were generated in the non-subduction settings (e.g., midocean ridge, oceanic island and oceanic plateau), whereas the crustal recycling-related basaltic melts generally have (Nb/La) PM < 0.75 ( Supplementary Fig. 2). In previous studies, the differentiated degree of rare earth elements (REEs; e.g., La/Sm and La/Yb ratio) was taken as the main criteria to distinguish Archean basalts with crustal recycling-related origins 21 . However, on modern Earth, some initial arc basalts and most of the backarc basin basalts (BABBs) also show low tholeiite-like La/Yb ratios and similar REE patterns to those of modern tholeiitic MORBs (Data from the GEOROC and Schmidt's database) 19 . Therefore, on this basis, we subdivide the Archean basalts into B-1 and B-2 groups with (Nb/La) PM ≥ and < 0.75, respectively.  22 , clearly distinct from B-1 samples ( Supplementary Fig. 3). Some B-2 samples have LREE depleted to flat patterns, similar to the modern basalts in the initial arc and back-arc basin settings 19 , but most of them display LREE enriched patterns with average (La/Yb) N of 2.19 (Supplementary Table 1).
Redox state of Archean upper mantle. The principle of the V-Ti redox proxy is that V is a multivalent element (2+, 3+, 4+, and 5+), and its Kds vary with fO 2 , melting P-T conditions, water contents and mineral compositions, in contrast the Kds of Ti are independent of fO 2 7,8,14 . In addition, both V and Ti are generally immobile during dehydration and melting of oceanic crust (including sediments), and are not strongly enriched in continental crust 7,14 . Together with their highly incompatible features among most mantle minerals, the V-Ti systematics can 'see through' early magmatic differentiation. Based on a thorough consideration of peridotite melting reactions and water contents, the V-Ti redox proxy has successfully estimated the fO 2 of both modern ambient and modified mantle 14 , providing a viable solution to the mantle fO 2 debate in the early Earth.
Recent studies highlight that degassing or interaction with polyvalent gas species (e.g., S) of basalts can potentially change the Fe 3+ /ΣFe contents, leading to an incorrect estimation of redox state inferred for their magmas 9,10 . However, there is no evidence that these changes would affect the V/Ti ratio, which is an indicator to the mantle source of less-or undifferentiated magmas 1 . Therefore, in this study, we only use basalts with MgO ≥ 8.5 wt.% to avoid possible effects from the late overprint of fractional crystallization and contamination from continental material 23 .
Sources of uncertainties that arise from trace element analyses, mantle compositional heterogeneity, P-T determinations and systematic biases for the Kds of V in mantle minerals and the fO 2 are incorporated into the application of V-Ti redox proxy and conversion of V/Ti ratio into the mantle fO 2 . These uncertainties are incorporated as the following: (1) Analytical uncertainties of V and Ti are set to 5 %, assuming that all of samples collected in this study were analyzed by the inductively-coupled plasma mass spectrometry (ICP-MS) technique 6,14 ; (2) Uncertainties of mantle heterogeneity were taken as 7 % for V and 12 % for Ti 14,31 ; (3) Uncertainties of P-T determinations are based on an absolute value of 0.20 GPa and 3 %, respectively 23 , assuming that there were no additional uncertainties in the water content assessment; (4) Estimate of uncertainties and systematic biases for the functions between Kds of V in mantle minerals and the fO 2 were incorporated into further propagations 14 . In this study, uncertainties (1) and (2) were firstly propagated into the calculation of the compositions of V and V/Ti in the melts and their Kds in bulk rock compositions, which were then propagated to the Kds in mantle minerals and fO 2 on the basis of a full account of uncertainties (3) and (4). On this basis, the average propagated 1 SD uncertainties of the mantle fO 2 revealed by B-1 and -2 samples are 0.40 and 0.32 log units, respectively (see the black error bar in Figs. 1-3).
During the evolution of the Superior Province between~3.0 and 2.7 Ga, the fO 2 values of its modified mantle increase from a low threshold (ΔFMQ −1.31 ± 0.66) to a high threshold (ΔFMQ −0.47 ± 0.54) at 2 wt.% H 2 O (Fig. 2a). Similarly, the ambient mantle beneath both the Superior Province and the Yilgarn Carton have been oxidized between~2.9 Ga and 2.7 Ga from ΔFMQ −1.79 ± 0.48 to ΔFMQ −1.16 ± 0.52 and ΔFMQ −1.76 ± 0.42 to ΔFMQ −1.29 ± 0.41, respectively (Fig. 1a). Although the amount of data is relatively limited, the fO 2 values for the Tanzania and Zimbabwe cratons fall within this~3.0-2.7 Ga mantle oxidation trend (Figs. 1 and 2), suggesting that along with the Superior and Yilgarn cratons, these four cratons may also constitute a craton group and probably share a similar Archean evolutionary history.
Therefore, the upper mantle (both ambient and modified mantle) beneath most Archean cratons, involving three temporal groupings has undergone significant oxidation during their evolutionary histories. Notably, the starting values of mantle oxidation beneath each craton or craton group are roughly the same (ambient mantle at approximate ΔFMQ −1. Fe 3+ with a corresponding release of O 2 at lower mantle 2,6 . It is believed that all these mechanisms are mutually responsible for Archean mantle oxidation. However, if mechanism (2) is dominant, the starting mantle fO 2 values of most cratons should gradually rise over time, distinct to our discoveries that initial fO 2 remain unchanged (Figs. 1-3). Furthermore, Archean mantle convection is commonly considered to be sluggish 40,41 , suggesting that mechanism (3) was not a dominant factor. Therefore, it is crucial to explore the relationship between crustal recycling and mantle oxidation.
Basalt Th/Nb ratios are generally taken to reflect the interaction between mantle-derived magmas and crustal materials, namely, crustal recycling, due to higher elemental abundances of Th against Nb in the continental crust [42][43][44] . Notably, the Th/Nb ratios of B-2 samples from each craton show a covariant elevated trend with the increase of mantle fO 2 . For example, the North Atlantic Craton increases from 0.28 ± 0.14 (1 SD) to 0.37 ± 0.12 between~3.8 and 3.0 Ga, the Superior Province from 0.16 ± 0.02 to 0.30 ± 0.22 between~3.0 and 2.7 Ga and the Dharwar Craton from 0.23 ± 0.19 to 0.37 ± 0.09 between~2.7 and 2.6 Ga (Supplementary Table 1; Fig. 2b-d). These increasing values suggest a significant increase in crustal recycling. Similarly, the εNd (t) values of B-2 samples from the Superior, Yilgarn, Tanzania and Zimbabwe cratons exhibit a significant mantle enrichment process between 3.0 Ga and 2.7 Ga 45 (Supplementary Fig. 6a; Data from the GEOROC database). Given that the Kds of Th are significantly lower than those of Nb within most mantle minerals, there should be a sharp decrease in Th/Nb ratio of mantle residues and hence B-1 samples during the melting and depletion of the ambient mantle ( Supplementary Fig. 6b). However, the B-1 samples from most cratons show a nearly constant Th/Nb trend under the background of gradual oxidation of ambient mantle, for example, the~3.5-3.2 Ga Pilbara Craton at~0.13, thẽ 2.9-2.7 Ga Yilgarn Craton at~0.14 and the~2.7-2.5 Ga North China Craton at~0.10 (Supplementary Table 1; Fig. 1b-d), indicating that recycled crustal materials can also influence the chemical compositions of ambient mantle perhaps via increasing mantle convection with craton maturation 46 . These processes would offset the decrease of basalt Th/Nb ratio during the depletion of the ambient mantle.
Two main mechanisms have been proposed to convey crustal materials into the hotter mantle of the early Earth: a form of plate subduction (or hot subduction) 17,47-54 and lithospheric drips 55,56 . Although our data alone cannot differentiate between these mechanisms, it does establish that Archean crustal recycling was widespread and ongoing within each craton over hundreds of millions of years 34,37 .

Methods
Data collection and filtration. The whole-rock geochemical database of thẽ 3.8-2.5 Ga basalts is a combination of the GEOROC database (http://georoc. mpch-mainz.gwdg.de/georoc/) and a subset of the North China Craton volcanic database assembled by our research group. In this study, the principles of data filtration are listed as follows: (1) We define basalts as volcanic samples with SiO 2 = 45-54 wt.% (including some basaltic andesites, because the primitive melts generated by melting of mantle peridotite may be more siliceous, especially in modern arc-related settings 19 ) and MgO < 18 wt.%, distinct to those of komatiites with MgO ≥ 18 wt.% 11,12 ; (2) They should have reliable formation ages traced by isotopic dating methods (mainly via whole-rock Sm-Nd, Re-Os and magmatic zircon U-Pb isotopic systematics, which are commonly used to mafic rock dating, and believed to be relatively accurate and not easily susceptible to late thermal events) or limited by geological relationships, for example, the interlayered relationship with felsic volcanic rocks or the intrusive relationship with granitoids and mafic dykes, providing the lower limit to their respective formation age); (3) The least-altered Archean basalts are reported to have low loss on ignition (LOI < 6 wt.%) and negligible Ce anomalies (0.9 < Ce/Ce* < 1.1, calculated as Ce N / Sqrt(La N × Pr N )) 57 ; and (4) Basalt samples that have been reported in the literature with significant crustal contamination are excluded. After these processes, this database comprises a total of 2304 basalt samples that cover fourteen Archean cratons (Supplementary Fig. 1; Supplementary Table 1).
Near fractional melting modeling. In the process of modelling, the content of trace element i in the instantaneous melt is calculated by C nþ1 i;m ¼ C n i;r D nþ1 i;bulk þ Fð1 À D nþ1 i;bulk Þ , where C nþ1 i;m is the element content in the melt at step n + 1, C n i;r is the element content in the residue after melt extraction at step n, D nþ1 i;bulk is the bulk partition coefficient at step n + 1, and F is the melt fraction of each step 14 14 . The initial V and Ti contents are set to be 79 ppm and 798 ppm, respectively, according to the depleted mantle (DM) compositions 31 . The initial mineral compositions apply orthopyroxene (Opx) Wo# = 4.38 (Wo# = X Wo / (X Wo + X En + X Fs ), in which X Wo , X En and X Fs are fractions of wollastonite, enstatite and ferrosilite, respectively), clinopyroxene (Cpx) Al T = 0.17 (Al in the tetrahedron-coordinated site), Opx Al T = 0.15 and spinel (Sp) Cr# = 10.65 (Cr# = Cr 2 O 3 / (Cr 2 O 3 + Al 2 O 3 ) on a molar basis), respectively 30 .
Parameter selection of V-Ti redox proxy. In model A, the initial mineral assemblages of peridotite constitute by 57% olivine (Ol), 28% Opx, 13% Cpx and 2% Sp 30 . The melting reactions of anhydrous depleted peridotite at 1.0 GPa were used in this model 27 . The melt productivity of 0.23°C for peridotite and 0.14°C for harzburgite were used to calculate the melting degrees and temperatures 14 . The recent empirical equations were used to calculate the Kds of V in Ol, Opx, Cpx and Sp, and Ti in Opx and Cpx at each melt fraction 14 . The initial mineral assemblages in model B consists of 52.3% Ol, 17.4% Opx, 26.2% Cpx and 4.1% Grt 28 , and the melting reactions were performed by anhydrous peridotite at~3.0 GPa 28 . In models C and D, the initial mineral assemblages are 53% Ol, 37% Opx, and 9.3% amphibole (Amph), and 52% Ol, 40% Opx, 0.6% Cpx, 7.7% Amph and 0.1% Grt, respectively 29 . The melting reactions were performed on a depleted peridotite, which is metasomatized by a MORB-derived hydrous silicate melt, at low pressure (~1-2 GPa) for modal C and high pressure (~3 GPa) for modal D 29 . The Kds of V in Ol, Opx and Cpx, and Ti in Opx and Cpx at each melt fraction were calculated by the above-mentioned empirical equations with variables of mineral compositions and melting P-T conditions 14 . In addition, the Kds of V in Grt and Amph, and Ti in Ol, Grt and Amph are assumed to be constant and used the values from http://earthref.org/GERM/. All of the calculations are compiled into a macro document and listed in Supplementary Table 2.

Data availability
The whole-rock geochemical, Nd-isotope and locational data for basalt samples are provided in Supplementary Data 1, and can be also downloaded from the GEOROC database (http://georoc.mpch-mainz.gwdg.de/georoc/). Data for mineral partition coefficients are provided in Supplementary Data 2, and can be obtained from the GERM database (http://earthref.org/GERM/). The authors declare that all data supporting the findings of this study are available online.