Research ArticleEARTH SCIENCES

Continental igneous rock composition: A major control of past global chemical weathering

+ See all authors and affiliations

Science Advances  08 Mar 2017:
Vol. 3, no. 3, e1602183
DOI: 10.1126/sciadv.1602183

Abstract

The composition of igneous rocks in the continental crust has changed throughout Earth’s history. However, the impact of these compositional variations on chemical weathering, and by extension on seawater and atmosphere evolution, is largely unknown. We use the strontium isotope ratio in seawater [(87Sr/86Sr)seawater] as a proxy for chemical weathering, and we test the sensitivity of (87Sr/86Sr)seawater variations to the strontium isotopic composition (87Sr/86Sr) in igneous rocks generated through time. We demonstrate that the 87Sr/86Sr ratio in igneous rocks is correlated to the epsilon hafnium (εHf) of their hosted zircon grains, and we use the detrital zircon record to reconstruct the evolution of the 87Sr/86Sr ratio in zircon-bearing igneous rocks. The reconstructed 87Sr/86Sr variations in igneous rocks are strongly correlated with the (87Sr/86Sr)seawater variations over the last 1000 million years, suggesting a direct control of the isotopic composition of silicic magmatism on (87Sr/86Sr)seawater variations. The correlation decreases during several time periods, likely reflecting changes in the chemical weathering rate associated with paleogeographic, climatic, or tectonic events. We argue that for most of the last 1000 million years, the (87Sr/86Sr)seawater variations are responding to changes in the isotopic composition of silicic magmatism rather than to changes in the global chemical weathering rate. We conclude that the (87Sr/86Sr)seawater variations are of limited utility to reconstruct changes in the global chemical weathering rate in deep times.

Keywords
  • chemical weathering
  • continental crust
  • strontium isotopes
  • hafnium
  • detrital zircon
  • arcs
  • supercontinent
  • Paleoclimate
  • phanerozoic
  • Neoproterozoic

INTRODUCTION

The chemical weathering of silicates transfers elements from the continental crust to seawater and exerts a direct control on several biogeochemical cycles. For instance, chemical weathering of silicates transfers calcium (Ca) and magnesium (Mg) to seawater and regulates atmospheric carbon dioxide levels and temperature at the surface by controlling the rate of marine carbonate precipitation (1). Despite decades of research, the mechanisms controlling chemical weathering throughout Earth’s history remain highly debated (2). The flux of elements from continental weathering to seawater is thought to be primarily controlled by the rate of chemical weathering. This rate depends mostly on runoff, temperature, and erosion and is thus modulated through time by changes in topography (3, 4), paleogeography (5), climate (57), and biological evolution (8, 9). The type of rock being weathered is also a critical parameter in understanding chemical weathering because minerals have distinct chemical composition and dissolution kinetics (10). However, the chemical composition of rocks subject to weathering is usually not accounted for when reconstructing the long-term changes in global chemical weathering.

The interpretation of the strontium isotope ratio in seawater [(87Sr/86Sr)seawater] is a good illustration of this issue. The (87Sr/86Sr)seawater variations are commonly used to estimate changes in global chemical weathering rates throughout Earth’s history (11). The (87Sr/86Sr)seawater curve displays an overall exponential increase controlled by the progressive differentiation of the Earth’s crust and the associated increase in the rubidium-to-strontium ratio (Rb/Sr) (11, 12). This trend is superimposed by second-order fluctuations at the scale of tens of million years, which are interpreted as reflecting changes in the relative Sr flux from isotopically distinct sources (1115). This interpretation assumes that the 87Sr/86Sr ratios of the radiogenic and unradiogenic Sr sources change too slowly or are not large enough to control the (87Sr/86Sr)seawater variations (11). For instance, the steep rise in (87Sr/86Sr)seawater ratio in the last 40 million years (My) has been explained by several competing “flux”-based hypotheses, including (i) an enhanced radiogenic Sr flux from radiogenic continental surfaces associated with uplift (3, 16) or glacial processes (7, 17, 18) and (ii) a decreased unradiogenic Sr flux from the oceanic crust associated with slower seafloor spreading rate (19) or cooler ocean temperature (20). Rapid changes in the composition of the 87Sr/86Sr ratio of continental surfaces have been invoked to explain some more specific features in the (87Sr/86Sr)seawater curve. For instance, the emplacement of large igneous provinces (LIPs) or the uplift of ophiolites can rapidly decrease the average 87Sr/86Sr ratio from continental weathering (2123), whereas the weathering of Sr-rich radiogenic metalimestones could raise the 87Sr/86Sr of rivers (24, 25). In contrast, little attention has been given to the potential variations in the average 87Sr/86Sr ratio in igneous rocks generated at plate margins. Plate margin magmatism is the dominant contributor to new continental crust (26, 27) and serves as parent rock to a large portion of the siliciclastic sediments of Earth’s surface (28). The 87Sr/86Sr ratios in silicic rocks generated in this tectonic setting are more variable than those in basalts (29, 30), and the impact of these compositional variations on the (87Sr/86Sr)seawater variations has not been explored.

Erosion and burial remove rocks from the surface through time, which complicates the reconstruction of the compositional evolution of igneous rocks (31). Fortunately, eroded components of rocks, such as zircon, are preserved in the siliciclastic sedimentary record. Zircon grains can survive multiple sedimentary cycles, because of their resistance to physiochemical alteration (32), and contain a wealth of geochemical information on their hosting igneous rock (32). Uranium-lead dating (U-Pb date) provides precise crystallization age of the zircon grains, whereas epsilon hafnium (εHf) records the degree to which a melt incorporates juvenile mantle versus reworked preexisting crust sources (32). Large compilations of integrated U-Pb date and εHf from detrital zircon have been used to reconstruct the evolution of the continental crust throughout Earth’s history (3336). However, generation and preservation biases in the detrital zircon record complicate the interpretation of the secular εHf variations through time (35). Mafic rocks (high-temperature, low-silica dry magmas) have low-zircon fertility in comparison with felsic rocks (low-temperature, high-silica hydrous magmas) (28, 37, 38). Zircon grains generated in intraplate volcanism, late accretionary, and collisional tectonic settings are preferentially preserved in sediments relative to those generated in extensional tectonic setting (3943). Ultimately, the detrital zircon record is likely biased toward silicic igneous rocks generated in intraplate volcanism or convergent margins, although the debate on this question is far from being resolved (28, 33, 35, 37, 38, 41, 42, 44). Acknowledging the biases outlined above, we used the detrital zircon record to reconstruct the initial 87Sr/86Sr ratio of zircon-bearing igneous rocks [(87Sr/86Sr)i-zig] over the last 1000 My. We discuss the strengths and limitations of this record, and we reinterpret the (87Sr/86Sr)seawater variations over the last 1000 My by accounting for the changes in isotopic composition of igneous rocks generated through time.

RESULTS

Here, we present a method to reconstruct the (87Sr/86Sr)i-zig variations for the last 1000 My. We compiled a first database gathering integrated U-Pb date and Hf isotopes on detrital zircon grains from globally distributed siliciclastic sediments (database S1). Detailed metadata are used to filter the data to minimize the issues associated with the use of combined U-Pb and Hf isotopes on detrital zircon grains (see Materials and Methods and the Supplementary Materials). The screened database gathers integrated U-Pb date and Hf isotopes from 24,715 zircon grains from 535 individual Mesozoic and Cenozoic sediments (Fig. 1). A debiasing method based on bootstrap resampling was applied to the compiled data set to correct for geographic and sampling biases in the detrital zircon record (see Materials and Methods and the Supplementary Materials). To account for the uncertainty in εHf and U-Pb date of individual zircon grains, we performed the resampling with introduction of Gaussian noise scaled to the εHf value and U-Pb date of the grains (see Materials and Methods and the Supplementary Materials). Using the resampled data set, we applied a smoothing procedure to obtain an approximately unbiased estimate of the εHf in detrital zircon grains through time (Fig. 1). The sensitivity of the smoothed trend to different weighing schemes, uncertainty levels, and resampling strategies was also tested to verify the robustness of our results (see Materials and Methods and the Supplementary Materials).

Fig. 1 Plot of εHf of detrital zircon grains versus U-Pb date.

Black dots represent the nonresampled εHf data from a screened subset of database S1. Moving median and quartiles at 1-My increments for the resampled database S1 (red line and green lines) and for the nonresampled database S1 (black line). The resampled data set includes both the debiasing resampling and the uncertainty propagation resampling (see Materials and Methods and the Supplementary Materials). The uncertainty propagation resampling uses U-Pb date (±2%) and εHf (±0.62).

We compiled a second database combining whole-rock strontium isotope ratio (87Sr/86Sr) data with integrated U-Pb date and Hf isotopes from magmatic zircon for Phanerozoic igneous rocks (database S2). The epsilon strontium (εSr) of the whole-rock samples and the εHf of the associated zircon grains were calculated. We screened the database to minimize uncertainty in the εHf and εSr calculations (see Materials and Methods and the Supplementary Materials). The screened data set gathers 351 individual igneous rock samples with combined εSr whole-rock and associated average integrated U-Pb date and Hf isotopes from their hosted magmatic zircon grains. The εHf of zircon grains correlates with the εSr of their hosting igneous rock, reflecting the coupling between the Sr and Hf isotope systems in magmatic processes (Fig. 2). However, Hf and Sr can be decoupled, with Hf always being incompatible and Sr becoming compatible once plagioclase crystallization begins (45). This decoupling, as well as the broad range of Sr and Hf content in the parent magma, explains the scatter in the correlation between εHf and εSr. It should be noted that the screening procedure did not significantly change the correlation equation between whole-rock εSr and the average εHf of hosted magmatic zircon grains. We used the relationship between the εHf of zircon grains and the εSr of their hosting igneous rocks to estimate the secular variations in the (87Sr/86Sr)i-zig for the last 1000 My (Fig. 2). Secular changes in the (87Sr/86Sr)i-zig ratio reflect the changing proportion of juvenile and reworked materials generated in orogenies through time and vary with the supercontinent cycle (46).

Fig. 2 Plot of εHf values of inherited zircon grains versus whole-rock εSr values for Phanerozoic igneous rocks.

(A) Nonscreened database S2. (B) Screened database S2. Screening procedure removed rock with crystallization age older than 300 Ma and with 87Rb/86Sr superior to 40. Regression models fitted to the data (red lines) with 95% confidence interval (blue lines). RMSE, root mean square error.

We compared the smoothed (87Sr/86Sr)i-zig ratio curve with the variations of the (87Sr/86Sr)seawater ratio through time (Fig. 3). To facilitate the time series analysis, we built upon previous work (12) to normalize the (87Sr/86Sr)seawater variations [N(87Sr/86Sr)seawater] (see Materials and Methods and the Supplementary Materials). To normalize this curve, we removed from the (87Sr/86Sr)seawater curve the signal associated with the radiogenic decay of the crust (Fig. 3A) (47). The N(87Sr/86Sr)seawater curve correlates with the (87Sr/86Sr)i-zig curve over the last 1000 My (Fig. 3). The highest correlation [cross-correlation function (CCF) > 0.65] occurs when the time series are centered on each other (lag = 0 ± 20 My), with the correlation decreasing rapidly for time lags larger than 20 My (Fig. 3C). To test the strength of the correlation at different periodicity, we decomposed both N(87Sr/86Sr)seawater and (87Sr/86Sr)i-zig curves into a slow-varying trend component and a fast-varying component by applying two respective fourth-order Butterworth filters (see Materials and Methods and the Supplementary Materials). The signals show a high correlation (CCF ≈ 0.6) for periodicities above 30 My. The correlation decreases for periodicity shorter than 30 My likely because of the uncertainty in U-Pb date and εHf value, which make the (87Sr/86Sr)i-zig ratio increasingly noisy. We tested the sensitivity of the correlation between these time series to our weighing scheme and normalization procedure and found the correlation between these curves to be robust and independent of our data analysis choices.

Fig. 3 Time series analysis between N(87Sr/86Sr)seawater (blue) and the median (red) and mean (green) (87Sr/86Sr)i-zig, both decomposed to their slow-varying trend and fast-varying residues.

N(87Sr/86Sr)seawater ratio (A) and (87Sr/86Sr)i-zig ratio (B) as in Fig. 1. (C) CCF. The first column corresponds to the original data at a 1-My interpolated sampling interval, the second column corresponds to the signal filtered for 700-My or longer periodicity, and the third column corresponds to the residual filtered between 30- and 700-My periodicity. The range of lags tested was set at ±200 My to include slowly exhumed plutonic rocks. A negative lag corresponds to the N(87Sr/86Sr)seawater ratio lagging the (87Sr/86Sr)i-zig ratio.

DISCUSSION

We interpret the correlation between the N(87Sr/86Sr)seawater and the (87Sr/86Sr)i-zig variations over the last 1000 My as reflecting a direct control of the isotopic composition of silicic igneous rocks on the (87Sr/86Sr)seawater variations. We suggest that when the global isotopic composition of silicic igneous rocks increases, the rapid cycling and subsequent weathering of these rocks lead to an increase in the (87Sr/86Sr)seawater ratio within a relatively short time period (<20 My). The (87Sr/86Sr)i-zig variations probably reflect the changes in the relative proportion of evolved to less evolved magmas generated in orogenies and subduction zones during different stages of the supercontinent cycle (46, 48). Here, we reinterpret the (87Sr/86Sr)seawater variations over the last 1000 My by accounting for this compositional variable. We also propose alternative hypotheses that could explain the correlation between these time series.

The N(87Sr/86Sr)seawater and the (87Sr/86Sr)i-zig ratios correlate strongly over different periodicities (Fig. 3C). Igneous rocks forming during collisional magmatism, accretionary orogenies in advancing phase (for example, Andean orogeny), and mature island arcs are characterized by a high reworking of older crust and have high (87Sr/86Sr)i-zig values (46, 48). In contrast, igneous rocks forming in island arcs, extensional arcs, or accretionary orogenies in retreating phase (rollback) are less evolved and have low (87Sr/86Sr)i-zig values (46, 48). We thus suggest that the overall high N(87Sr/86Sr)seawater ratio between 600 and 300 Ma (million years ago) reflects the continued assembly of the Gondwana-Pannotia and Pangea supercontinents (46). During that period, collisional orogenies are frequent and subduction occurs dominantly in advancing phase, with the opening of several ocean basins forming igneous rocks with high (87Sr/86Sr)i-zig ratios (46). The lower N(87Sr/86Sr)seawater ratios before 700 Ma and between 250 and 100 Ma might reflect the increased proportion of less evolved island and extensional arcs as well as the long-term closing of ocean basins during the assembly of the Gondwana-Pannotia and Amasia supercontinents, respectively. During those periods, subduction zones are more frequently in retreating phase forming immature arcs with lower (87Sr/86Sr)i-zig ratios (Fig. 3A) (46). Interpreting the (87Sr/86Sr)i-zig fluctuations on a shorter time scale is more ambiguous because of the bias in the detrital zircon record, the overlapping of breakup and assembly phases, and the difference in the supercontinent assembly mode between Pangea and Gondwana (that is, introversion versus extroversion) (46). However, we argue that the (87Sr/86Sr)i-zig fluctuations track the global evolution of the isotopic composition of silicic magmatism occurring at plate boundaries through time and that these isotopic fluctuations are rapidly transmitted to seawater through erosion and weathering.

Two prerequisites for this interpretative framework to be valid are as follows: (i) the amplitude of (87Sr/86Sr)i-zig variations and the magnitude of the Sr flux from young zircon-bearing igneous rocks have to be significant enough to affect the marine Sr isotope budget, and (ii) young zircon-bearing igneous rocks have to be cycled (that is, exhumed and weathered) rapidly on the surface (<20 My) for the time series to be synchronized (Fig. 3C). The reconstructed (87Sr/86Sr)i-zig ratio ranges from ~0.704 to ~0.710 for the last 1000 My (Fig. 4C). Given this range, the flux of Sr from zircon-bearing igneous rocks needs to contribute around half of the total Sr flux to control the totality of the N(87Sr/86Sr)seawater variations. Zircon-bearing igneous rocks represent only ~10% of the silicates exposed on Earth’s surface at the present day (49). Young igneous rocks and associated volcaniclastic sediments are usually highly weatherable, exhumed rapidly, and uplifted which favors a very high solute flux from these rocks (5053). At the present day, young igneous rocks (including basalts) contribute more than a third of the total Sr flux to seawater (50, 52, 53). Together, these observations suggest a large contribution of young igneous rocks to the Sr budget in seawater.

Fig. 4 Interpretation of the marine Sr budget.

Maximum correlation coefficient through time (A) between the N(87Sr/86Sr)seawater ratio (B) and the (87Sr/86Sr)i-zig ratio (green) (C). Curves (B) and (C) are obtained by summing the slow-varying trend and fast-varying residues of the mean value obtained in Fig. 3 (A and B, respectively). The CCF is calculated for a 150-My time window moving in 10-My increments, with lags ranging between 0 and 20 My. (D) Frequency of collisional orogenies (46), major LIPs (23, 314), and major glaciations. Shaded blue bars represent ocean opening stages.

The strength of the correlation between the time series is not uniform through time (Fig. 4A). Periods of high correlation (CCF > 0.7) coincide with supercontinent dispersal and assembly stages (Fig. 4D) (50). This observation suggests that during periods when plate margins are dominated by subduction magmatism, the (87Sr/86Sr)i-zig variations, associated with changes in the relative proportion of evolved versus less evolved arcs, are the dominant factor of N(87Sr/86Sr)seawater variations. One issue with this interpretation is that the detrital zircon record, and hence (87Sr/86Sr)i-zig variations, does not record in equal proportion between evolved and less evolved magmatism because zircons are preferentially formed in silicic magmas (37, 41). Hence, the sensitivity of the (87Sr/86Sr)i-zig variations to magmatism from less evolved arcs remains unknown. At the present day, less evolved arcs contribute significantly to the Sr budget in seawater and are likely to have been an important source of Sr to the ocean in the past (52). Owing to the large volume of magma generated during less evolved arc magmatism (42), zircon grains from this tectonic setting might still be the dominant control of the reconstructed (87Sr/86Sr)i-zig ratio during periods dominated by this style of magmatism. Another possibility to explain the strength of the correlation is that if the rates of mid-ocean ridge magmatism and subduction magmatism are at steady state (54, 55), then the combined Sr flux from these sources should be relatively constant through time. As the 87Sr/86Sr ratio of basalts generated in mid-oceanic ridges and less evolved arcs varies little through time (29), the isotopic composition of the combined subduction and mid-oceanic ridge Sr fluxes should be mostly sensitive to the 87Sr/86Sr variations of silicic magmatism in more evolved arcs. In any cases, the strong correlation between (87Sr/86Sr)i-zig and N(87Sr/86Sr)seawater variations through most of the last 1000 My suggests a strong control of the isotopic composition of silicic arcs on the (87Sr/86Sr)seawater variations.

Periods with a low correlation coefficient (CCF < 0.5) between the time series could be explained by the following: (i) the incomplete geographic coverage, potential biases, and/or uncertainty in the detrital zircon database or (ii) the decoupling between εHf and εSr in magmatic processes (Fig. 2). However, we notice that intervals when the time series are not correlated coincide with specific paleogeographic, tectonic, or climatic events. We suggest that although the relative proportion of evolved versus less evolved arcs is the primary driver of (87Sr/86Sr)seawater variations, the remaining variance could be explained by changes in the relative Sr fluxes from isotopically distinct sources associated with paleogeographic, tectonic, or climatic events. Our reconstructed (87Sr/86Sr)i-zig variations do not account for the radiogenic Sr flux from radiogenic recycled silicates, the unradiogenic Sr flux from mid-oceanic ridges, or the Sr flux from recycled carbonates. Although the relative Sr fluxes and isotopic composition from these sources might be relatively constant over long time periods, a shift in the relative flux from one of these isotopically distinct sources will alter the correlation between (87Sr/86Sr)i-zig and (87Sr/86Sr)seawater variations. The correlation coefficient decreases during the Cryogenian, the Early Ediacaran, and the Neogene periods. During these periods, the (87Sr/86Sr)i-zig values decrease, whereas the (87Sr/86Sr)seawater values increase, indicating a “missing” radiogenic component to the Sr isotope budget in seawater (Fig. 4C). These periods coincide with time when the rate of collisional orogenies was high and/or when the climate was in an icehouse period (Fig. 4). During periods of collisional orogenies or during icehouse periods, increased erosion might increase the radiogenic Sr flux relative to the unradiogenic Sr flux (Fig. 4) (3, 18). The correlation between the time series also decreases during the amalgamation of the Pangea supercontinent. However, during this period, the (87Sr/86Sr)i-zig values increase, whereas the (87Sr/86Sr)seawater values decrease (Fig. 4C). Global runoff has been relatively constant over long time scale, and only a few paleogeographic configurations have led to significant runoff changes in the last 1000 My (5, 56). During Pangea amalgamation, increased continental interior aridity led to a twofold decrease of the global runoff relative to other time periods, which would have largely decreased the radiogenic Sr flux from cratonic areas (5). This progressive decrease in Sr flux from radiogenic continental areas combined with the increase in mid-oceanic ridge magmatism associated with the Tethys opening might explain the lower (87Sr/86Sr)seawater values during the Permian period (57). Changes in the isotopic signature of the Sr flux exported from continents are associated with the uplift of metamorphic rocks, ophiolites, and LIPs (22) can also play a significant role in the Sr isotope budget (56). The rapid weathering of these isotopically distinct rocks in favorable paleogeographic or tectonic configuration likely contributes to some of the (87Sr/86Sr)seawater variations over the last 1000 My (Fig. 4).

An alternative and/or complementary interpretation to explain the correlation between the (87Sr/86Sr)i-zig and (87Sr/86Sr)seawater ratios is that the generation and preservation biases of the detrital zircon record correlate to changes in Sr flux from radiogenic and unradiogenic sources of Sr. The (87Sr/86Sr)seawater variations are thought to be highly dependent on the rate of collisional orogenies because of the potential large increase in the Sr flux from uplifted recycled silicates with radiogenic isotopic signatures (3, 16, 25). During periods when both collisional orogenies and subduction magmatism occur, the global detrital zircon record will display a high global zircon count associated with generation and/or preservation biases (42, 43). The detrital zircon record will be biased toward magmatism occurring during collisional orogenies and/or arc magmatism preserved in this convergent tectonic setting (43). The zircon grains generated in this tectonic setting might have a more radiogenic isotopic signature that will bias the (87Sr/86Sr)i-zig variations toward higher values. These high (87Sr/86Sr)i-zig values could coincide with an increase in the weathering rate of uplifted radiogenic recycled silicates (Fig. 4). Conversely, during periods with no collisional orogenies, the global detrital zircon record will display a low global zircon count associated with generation and/or preservation biases (42, 43). The reconstructed (87Sr/86Sr)i-zig ratio might be lower during these periods, reflecting the higher contribution of zircon grains with low (87Sr/86Sr)i-zig ratio (46). These low (87Sr/86Sr)i-zig values could coincide with increased Sr flux from unradiogenic basalts at island arcs or mid-oceanic ridges as observed for the Mesozoic (Fig. 4). However, this framework requires that magmatism associated or preserved during collisional orogeny has a distinct isotopic signature relative to other magmatism preserved in the detrital zircon record (43, 46, 48, 58). This might be true for some supercontinents and some collisional orogenies, but it also depends on the mode of supercontinent assembly and the type of arcs preserved during collision (43, 46, 58).

In summary, we present a new method to reconstruct the evolution of the strontium isotopic composition of silicic igneous rocks through time. We argued that the (87Sr/86Sr)i-zig variations are the dominant control of (87Sr/86Sr)seawater evolution over the last 1000 My. Instead of interpreting (87Sr/86Sr)seawater variations as reflecting changes in the chemical weathering rate of radiogenic continental surfaces, we suggest that they primarily track the relative proportion of evolved versus less evolved magmas in the continental crust. These changes in the global isotopic composition of silicic magmas are likely controlled by the different types of orogeny and the variations in subduction phases occurring with distinct supercontinent cycle stages. The remaining (87Sr/86Sr)seawater variations are probably explained by changes in the relative Sr flux from isotopically distinct sources associated with specific paleogeographic configurations, mountain-building events, emplacement of LIPs, or climate variations. We conclude that the (87Sr/86Sr)seawater variations are of limited utility to reconstruct the long-term chemical weathering rate.

MATERIALS AND METHODS

Database compilation and screening procedures

Database S1. We assembled a geochemical database with integrated U-Pb date and Hf isotopes from a total of 183 individual studies (34, 59241), including some from preexisting compilations (database S1) (35, 36, 46, 90, 242). We added a large number of Mesozoic and Cenozoic samples from regions that were underrepresented in previous databases, such as the North American Cordillera, the South American Cordillera, Antarctica, the Middle East, Turkey, and Australia. Each sample was paired with a geospatial location. Additional data from rivers in Europe and western Africa would be required to fill geographic gaps. When location was not given, we used Google Earth and metadata from the study to estimate the approximate latitude and longitude of the samples. Other metadata compiled include the lithostratigraphic description of the sediment, age of deposition, instrument used for analysis, and range of Th/U ratio and cathodoluminescence description of the samples. For each zircon grain present in the databases, we also calculated εHf using Eq. (1). The result is a data set including up to 37 defined variables for each of the compiled samples from varied geographic locations. These detailed metadata were used to filter the database and to minimize the issues associated with the use of combined U-Pb and Hf isotope on detrital zircon grains (36, 243). To minimize crystallization age and analytical uncertainty on εHf calculation, we only selected the following: (i) zircon grains analyzed using in situ analysis, (ii) zircon grains with high U-Pb date concordance (<10%), and (iii) zircon grains with low analytical uncertainty and isobaric interferences 176Yb/177Hf < 0.2 or 176Lu/177Hf < 0.005. To minimize the possible impact of mixed sampling of complex zircon grains, we screened zircon grains that were reportedly affected by metamorphism using Th/U ratio (Th/U < 0.05) and/or visual occurrence of metamorphic rims from cathodoluminescence. The resulting database contains 54,406 zircon grains with integrated U-Pb date and Hf isotopes from 1263 individual siliciclastic sediments.

Database S2. We compiled a second database combining whole-rock 87Sr/86Sr ratio data and average integrated U-Pb date and Hf isotopes from hosted magmatic zircon grains from 60 individual studies (database S2) (244304). For each whole rock, we calculated εSr using Eq. (1). We compiled igneous rocks data only from the Phanerozoic to minimize the uncertainty in the εSr calculations. The data set contains 441 analyses of 87Sr/86Sr data from whole igneous rocks coupled with the average U-Pb date and Hf isotopes from their hosted magmatic zircon grains. We used this database to relate the εHf of zircon grains with the εSr of their hosting igneous rock (Fig. 2). We first attempted to minimize uncertainty in εSr calculations by filtering data with the highest analytical or calculation uncertainty. At equal age, whole rock with very high Rb/Sr ratio will have a higher uncertainty in their calculated εSr. Similarly, at equal Rb/Sr ratio, older rocks will have a higher uncertainty in the εSr calculation. We tested several thresholds of Rb/Sr and age uncertainty and found that keeping whole-rock data younger than 300 Ma and with 87Rb/86Sr inferior to 40 led to the best improvement in correlation coefficient and root mean square error. This screening effectively removes 90 εSr from the data set or ~20% of the data.

Epsilon value calculation. The formula used to calculate the epsilon value is given below for the Hf isotope system, but the formula is also valid for the strontium isotope systemEmbedded Image(1)where (176Hf/177Hf)sample and (176Hf/177Hf)CHUR are the Hf isotope ratio of the zircon and chondritic uniform reservoir (CHUR) at crystallization age, respectively. These values were calculated using the radiogenic equation, the present-day measured 176Hf/177Hf and 176Lu/177Hf, and a radiogenic decay constant of 1.867 × 10−11 years (305). (176Hf/177Hf)CHUR and (176Lu/177Hf)CHUR at the present day are 0.282772 ± 29 and 0.0332 ± 2, respectively (306). (87Sr/86Sr)CHUR and (87Rb/86Sr)CHUR at the present day are 0.7045 and 0.0824 (307, 308).

Debiasing the detrital zircon record

Debiasing method. In accordance with the Horvitz-Thompson estimation theory, we used weighted resampling to correct for biases in the detrital zircon record and to obtain the best estimate of the εHf ratio of zircon-bearing igneous rocks. The Horvitz-Thompson estimation theory dictates that an unbiased estimate of εHf at any time can be obtained by weighting each observation with the inverse of its inclusion probability in the sample (309). We resampled with replacement (bootstrap) from our database with resample probabilities equal to the Horvitz-Thompson weights and performed a smoothing procedure on the resampled data set to obtain an approximately unbiased estimate of the εHf trend through time. Note that our resampling procedure is similar to that of Keller and Schoene (31), although we made the connection between weighting and unbiasedness explicit via the Horvitz-Thompson estimation theory.

We note three key sources of preferential sampling (bias) in our database. First, sediments present in the database integrate different proportions of the Earth’s surface. Some sediments in our database drain very large catchment (for example, the Mississippi River), whereas some others represent much smaller catchments. The sediments draining large catchments tend to have a broader U-Pb date distribution than those draining smaller catchments. The zircon grains collected in small catchments are more likely to be biased toward a more local magmatic signature than those in the larger catchments. We did not identify a proper way to correct for this bias, but to minimize it, we used only sediments with Cenozoic and Mesozoic depositional age, which represent 24,715 zircon grains from 535 individual Mesozoic and Cenozoic sediments. The idea is that older zircon grains in those young sediments likely represent a more integrated picture of Earth’s surface because they had time to mix during multiple sedimentary cycles.

Second, a disproportionate number of samples were observed from similar geospatial locations (geographic bias) (fig. S1). If uncorrected, Asian samples would dominate the estimation of trend εHf. We assumed that the present-day sediment geographic distribution was a sufficient estimate for geographic debiasing because of the relatively similar paleogeography throughout this period. We weighted the geographic locations of each sediment sample usingEmbedded Image(2)where dist(i,j) represents the haversine distance between the location of sediment i and the location of sediment j, and var denotes the sample variance of the distances.

Third, to correct the bias introduced by a larger number of zircon grains sampled from the same sediment (sampling bias) (fig. S2), we introduced the weighting byEmbedded Image(3)where Ni is the number of zircon grains sampled for the sediment i. To avoid overweighing poorly sampled sediments, we screened out sediments with less than 10 zircon grains sampled. The final database used during the debiasing procedure contained 24,377 individual zircon grains from 470 individual sediments deposited during the Cenozoic or Mesozoic.

To combine these sampling mechanisms into a resample probability that places equal importance on each mechanism, it is necessary to rescale each component. We chose to do this by matching the first and second sample moments via a location- and scale-preserving correction function, z(·), which ensures that the collection of weights has a mean of 100 and an SD of 10 for each of the two mechanisms. Thus, our overall resample probabilities are given byEmbedded Image(4)

We chose a resample of size 106 to balance computational constraints and ensure that actual resampled proportions were approximately equal to the desired resample proportion.

Response of the εHf trend to resampling procedure. The spatial distribution of geographic weights shows higher weights for isolated samples and lower weights for clustered samples, as we intended our weighted resampling procedure to induce (fig. S3). The prior and posterior density distributions were compared to the present-day distribution of continental surfaces (fig. S4). The results suggest substantial improvement in geographic sampling distribution. However, as observed with the relatively bimodal distribution of the geographic weights (fig. S3), the geographic declustering procedure could benefit from applying more advanced three-dimensional declustering methods. Prior and posterior resampling distributions of zircon U-Pb date as well as secular εHf curve do not change markedly after resampling (fig. S5), although undersampled and isolated samples have a much higher weight (figs. S6 and S7). The U-Pb date distribution of the resampled database resembles the U-Pb date distributions of previous databases, with a series of U-Pb peaks centered on supercontinent periods (fig. S5). However, we used only sediments deposited during the Mesozoic and Cenozoic, which explains the progressive increase in zircon frequency throughout Earth’s history that is not observed in other databases compiling sediments of all depositional ages. Previous work also noted that the secular trends in εHf differed between compilations for the last 250 Ma of Earth’s history (46). We suggest that these differences are associated with significant geographic biases in these databases. The database compiled in this work markedly extends the geographic coverage of Mesozoic and Cenozoic zircon grains. Our debiasing approach also limits the influence of sampling and geographic biases in investigating the resulting εHf trend.

Propagating U-Pb date and εHf uncertainties in the smoothed εHf trend

Gaussian noise resampling. On the basis of the resampled data set, the εHf trend was calculated using a rolling median or rolling mean at 1 My, chosen to balance smoothness and informativity (Fig. 1). However, although we tried to minimize uncertainty through our screening procedure, both the U-Pb date and εHf value of individual zircon grains contain significant uncertainty. To test the sensitivity of the smoothed εHf trend to these uncertainties, we resampled the debiased data set by introducing Gaussian noise around individual U-Pb and εHf values using a range of uncertainties for each variable. We calculated a median age discordance of 2% (Q1, 0.7%; Q3, 4.1%) using the zircon grains for which the absolute age discordance was compiled (21% of the data). Similarly, we calculated a median uncertainty for εHf values of ±0.62 (Q1, ±0.45; Q3, ±1).

Response to Gaussian noise resampling. We found that the smoothed εHf trend is mostly sensitive to uncertainty in U-Pb date (fig. S8). The smoothed εHf trend shows very little change when testing a range of uncertainty for the εHf value (fig. S9). Conversely, the overall amplitude of variations decreases progressively as larger U-Pb date discordances are considered (fig. S8). For U-Pb date discordance superior to 10%, the smoothed εHf trend loses most of its structure, suggesting that lower discordance threshold should be considered when trying to reconstruct the long-term evolution of εHf through time (fig. S8). We choose the εHf trend accounting for the median uncertainty in U-Pb date (±2%) and εHf (±0.62) as our reference smoothed curve (Fig. 1).

Normalized (87Sr/86Sr)seawater curves

We used existing compilations of 87Sr/86Sr analysis in carbonates to generate the (87Sr/86Sr)seawater ratio curve over the last 850 My (13, 15). The quality and uncertainty associated with these data vary (310). The age model of Phanerozoic carbonates is usually more precise (<1 Ma) than that of Proterozoic carbonates (>5 Ma). The sampling density is much higher in the Cenozoic and Mesozoic and decreases progressively with time because of the scarcity of exposed carbonates with ages older than the Mesozoic (fig. S10). The sampling density remains higher than 1 My for the entire Phanerozoic, whereas several data gaps exist throughout the Ediacaran and Cryogenian. This low sampling density complicates the comparison of the Proterozoic trend with the (87Sr/86Sr)i-zig ratio. We applied a linear resampling to fill data gaps in the Proterozoic and calculated the median (87Sr/86Sr)seawater ratio at a 1-My time step.

To facilitate the time series analysis between (87Sr/86Sr)seawater and (87Sr/86Sr)i-zig curves, we detrended the (87Sr/86Sr)seawater variations [N(87Sr/86Sr)seawater]. (87Sr/86Sr)seawater ratio displays an exponential increase throughout Earth’s history primarily driven by the radiogenic decay of 87Rb into 86Sr in the silicate crust (12). On the basis of the radiogenic equation, the rate of radiogenic decay through time is controlled by the Rb/Sr ratio of the silicate crust. Rb and Sr are incompatible elements and are progressively concentrated into crustal silicates throughout Earth’s history, leaving the mantle depleted in these elements. Rb, being less compatible than Sr, further concentrates into crustal silicates, which leads to an increase in the Rb/Sr ratio of the crust through time. If the rate of crustal growth was constant throughout Earth’s history, the Rb/Sr ratio of crustal silicates would exponentially increase as a mirror image of the Rb/Sr ratio of the depleted mantle, but with a steeper gradient due to the small volume of the crust in comparison with the mantle (12). However, the rate of crustal growth was not constant through time and altered this idealistic model (35). On the basis of geochemical reconstructions, the time-integrated 87Rb/86Sr ratio could have remained relatively constant (47) or increased slightly (31). Consequently, instead of using an idealized growth curve with a continuous enrichment of the crust (47), we used a constant time-integrated 87Rb/86Sr ratio for silicates over the last 1000 My as our reference curve. However, we also tested the sensitivity of our conclusions to a potential twofold increase or decrease of the time-integrated 87Rb/86Sr ratio. We calculated the evolution of 87Sr/86Sr ratio in seawater associated with radiogenic decay over the last 1000 Ma asEmbedded Image(5)where [87Sr/86Sr(t)]seawater-decay is the 87Sr/86Sr ratio of seawater associated with radiogenic decay through time, [87Sr/86Sr(t1000 My)]seawater is the 87Sr/86Sr of seawater at t = 1000 My, (87Rb/86Sr)cs is the 87Rb/86Sr of the crust silicates for the last 1000 My, and λ is the decay constant of the parent isotope (1.42 × 10−11 year−1).

N(87Sr/86Sr)seawater is calculated asEmbedded Image(6)

We calculated the present-day value of the time-integrated 87Rb/86Sr ratio of the continental crust asEmbedded Image(7)where (87Sr/86Sr)riverine ratio is the 87Sr/86Sr ratio of the Sr exported from rivers at the present day (0.71144) (311, 312); (87Sr/86Sr)crust and (87Sr/86Sr)carb are the 87Sr/86Sr ratio of the continental crust silicates and carbonates, respectively, with (87Sr/86Sr)carb equal to 0.708 at the present day (52, 311, 313); and fcs and fcarb are the relative fluxes of Sr from the continental crust silicates and carbonates at the present day (0.37 and 0.63, respectively) (52, 311, 313).Embedded Image(8)where λ is the decay constant of the parent isotope (1.42 × 10−11 year−1) and (87Sr/86Sr)i is the initial 87Sr/86Sr of the bulk silicate Earth 4550 Ma (0.69897) (307).

We applied Eq. (7) and found a 87Sr/86Sr ratio of continental crust silicates of 0.7178 at the present day, with a corresponding 8Rb/86Sr ratio of silicates of 0.283. This value is used in Eqs. (6) and (7) to calculate the N(87Sr/86Sr)seawater ratio.

Time-series analysis

Both (87Sr/86Sr)seawater and (87Sr/86Sr)i-zig curves have been decomposed into a slow-varying trend component and a fast-varying component by applying two respective fourth-order Butterworth filters. A Butterworth filter is a widely used technique to decompose a time series into components with a given periodicity range. The slow-varying component was isolated from the original record by using a cutoff corner frequency at 1/700 My−1 with a low-pass Butterworth filter. The filter results in a slow-varying component that has a time scale of 700 My or longer. The fast-varying component was calculated using a band-pass Butterworth filter with a cutoff frequency between 1/700 and 1/30 My−1. These cutoff frequencies were selected so that the half periodicity covers the range of subduction and orogeny durations during Earth’s history (46). Cross-correlation coefficient analysis was applied to investigate the potential linear dependence between the pair of slow-varying components and also that between the pair of fast-varying components. In addition to correlation coefficient, which shows the linear dependence between two time series, the cross-correlation can further reveal possible lead-lag relationship, which depends on the sign of lag that corresponds to the maximum coefficient.

Sensitivity analysis

In the first step, we verified whether the smoothed (87Sr/86Sr)i-zig trend holds for different debiasing and resampling procedures. This is an important question considering the scatter in εHf values (fig. S5). We have already shown that the smoothed εHf trend [hence, the (87Sr/86Sr)i-zig ratio] is sensitive to the uncertainty in U-Pb date (fig. S8) but not sensitive to the uncertainty in εHf value. We further tested the sensitivity of the time series analysis to our screening and debiasing procedures. In the first simulation, we tested the impact of our debiasing procedure choices. We changed the sampling weight calculation from 1/n to 1/n2. The resulting smoothed (87Sr/86Sr)i-zig ratio and the time series analysis show almost no difference with our reference case (Fig. 3 and fig. S11). In the second simulation, we tested the sensitivity of our screening procedure to the time series analysis. Instead of using a data set with only sediments with Cenozoic and Mesozoic depositional ages, we used a data set with sediment of all depositional ages. We applied the same screening and debiasing procedure as in our reference case. In this scenario, the resulting smoothed (87Sr/86Sr)i-zig ratio shows some significant differences with our reference data set (Fig. 3 and fig. S12). The correlation between the time series is slightly lower than in our reference case, but the conclusions remain unchanged: The CCF between the two time series remains significant and centered on zero (fig. S12). We suggest that the slight decrease in correlation observed in this scenario reflects the integration of zircon data from small rivers that oversample local magmatism and bias the smoothed εHf values.

In the second step, we tested the sensitivity of the time series analysis when using a range of calibration procedures for the N(87Sr/86Sr)seawater ratio. In the first scenario, we assumed that the time-integrated 87Rb/86Sr ratio increased, following an ideal growth model increasing by a factor of 2 throughout the last 1000 Ma (fig. S13) (310). In the second scenario, we assumed that the time-integrated 87Rb/86Sr ratio decreased by a factor of 2 throughout the last 1000 Ma (fig. S14) (47). The modifications in the parameterization of the time-integrated 87Rb/86Sr ratio calibration lead to some significant changes in the slow-varying trend of (87Sr/86Sr)i-zig ratio but changes little the fast-varying (87Sr/86Sr)i-zig trend. In both scenarios, the correlation between (87Sr/86Sr)i-zig ratio and N(87Sr/86Sr)seawater ratio remains significant and centered on zero (figs. S12 and S13).

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/3/3/e1602183/DC1

fig. S1. Map of sediment locations.

fig. S2. Sampling bias: Histogram of zircon grain count per sediment from the screened data set (24,377 individual zircon grains from 470 individual sediments).

fig. S3. Geographic bias and geographic weights: Map of sediment sample locations for the screened data set with associated geographic resampling weights.

fig. S4. Geographic debiasing validation.

fig. S5. Resampled εHf using combined weights (red) superimposed on nonresampled εHf (black) from detrital zircons from Cenozoic and Mesozoic sediment data compiled in database S1.

fig. S6. Map of sediment sample locations with associated final resampling weights combining geographic and sampling weights.

fig. S7. Resampling weights visualization.

fig. S8. Sensitivity of the smoothed εHf trend to the uncertainty in U-Pb date.

fig. S9. Sensitivity of the smoothed εHf trend to the uncertainty in εHf values.

fig. S10. Evolution of the (87Sr/86Sr)seawater ratio in the last 1000 My from the analysis of 87Sr/86Sr ratio of carbonates.

fig. S11. Sensitivity of time series analysis to debiasing procedures.

fig. S12. Sensitivity of time series analysis to data screening.

fig. S13. Sensitivity of time series analysis to N(87Sr/86Sr)seawater curve parameterization by assuming a linear increase in Rb/Sr ratio of the crust in the last 1000 My.

fig. S14. Sensitivity of time series analysis to N(87Sr/86Sr)seawater curve parameterization by assuming a linear decrease in Rb/Sr ratio of the crust in the last 1000 My.

database S1. Excel data table with worldwide U-Pb and Hf isotopes in detrital zircons compiled from the literature.

database S2. Excel data table with 87Sr/86Sr and 143Nd/144Nd ratio data from whole igneous rock and U-Pb and Hf isotope data from their hosted zircons compiled from the literature.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank G. J. Bowen, C.-T. Lee, and R. Mills for helpful discussions. We thank B. Schoene, C. Spencer, and F. Macdonald for their constructive review comments. Funding: C.P.B. and X.-M.L. acknowledge funding from X.-M.L.’s University of North Carolina at Chapel Hill startup fund. Author contributions: C.P.B. designed the project and compiled the data sets. C.P.B., A.W., and X.Y. performed the statistical analysis and model development. All authors contributed to the interpretation of the results and writing of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Related Content

More Like This

Navigate This Article