Research ArticleCLIMATOLOGY

Dole effect as a measurement of the low-latitude hydrological cycle over the past 800 ka

See allHide authors and affiliations

Science Advances  07 Oct 2020:
Vol. 6, no. 41, eaba4823
DOI: 10.1126/sciadv.aba4823


The quest of geological proxies to evaluate low-latitude hydrological changes at a planetary scale remains an ongoing issue. The Dole effect is such a potential proxy owing to its global character. We propose a new approach to recalculate the fluctuation of the Dole effect (∆DE*) over the past 800 thousand years (ka). The ∆DE* calculated this way is dominated by precession cycles alone, with lesser variance in the obliquity bands and almost no variance in the eccentricity bands. Moreover, the ∆DE* is notably correlated with Chinese stalagmite δ18O record over the past 640 ka; simulated terrestrial rainfall changes between 30°N and 30°S over the past 300 ka. Our findings highlight the predominant role of the low-latitude hydroclimate in governing the ∆DE* on orbital time scales, while high-latitude climate impacts are negligible. In turn, we argue that the ∆DE* can be used to indicate low-latitude hydrological changes at a global extent.


A variety of geological proxies has been developed to reconstruct the history of the low-latitude hydrological cycle, which plays a crucial role in the global moisture and heat transport. In the past decade, the strong similarity of δ18O variations between atmospheric O218Oatm) and Chinese stalagmite, either on orbital or on millennial time scales, has raised the debate on the correlation between δ18Oatm and the low-latitude hydrology (15). As a chemical signal of the atmosphere, δ18Oatm is an integrated result of different processes that occurred over the globe, which is able to reflect biogeochemical and thus climatic changes at a planetary scale compared with regional proxies.

Interpreting the climate significance of δ18Oatm, however, turns out not to be straightforward. The δ18Oatm value is determined by the oxygen isotopic fractionation associated with photosynthesis and respiration of biosphere as well as the hydrological cycle, which is today more enriched by 23.5 ± 0.3‰ than oceanic water. This isotopic offset is traditionally referred to as the Dole effect. The respiration of both terrestrial and marine biosphere preferentially consumes 16O over 18O, causing a nearly equivalent isotope effect of ~19‰ [e.g., (6)]. This is the biggest source of the Dole effect. On land, plant transpiration can cause an isotopic enrichment of leaf water relative to soil water, resulting in a globally integrated photosynthetic Dole effect of 4 to 8‰ [e.g., (6, 7)]. In the marine realm, O2 produced by phytoplankton shows an 18O enrichment of less than 6‰ with respect to ambient seawater (8). The hydrological cycle, however, changes δ18Oatm in an opposite direction to the aforementioned two processes. This is because the ocean-to-land moisture transfer can always result in an isotopic depletion of terrestrial rainfall with respect to seawater. Rainwater is subsequently used by the terrestrial biosphere, which thus decreases the Dole effect (2, 6).

The fluctuation of the Dole effect (∆DE) over the past 800 thousand years (ka) has been reconstructed by removing the imprint of seawater δ18O from ice core–derived δ18Oatm records [e.g., (3, 6)]. During the mid-to-late Pleistocene, varying global ice volume is the first-order reason causing mean oceanic water δ18O changes. Accordingly, the reconstructed or the modeled ice-volume isotopic effect (δ18Osw) (9, 10) was subtracted from δ18Oatm to obtain the ∆DE [e.g., (3)]. Using this approach, the reconstructed ∆DE shows strong variance in the precession bands with an amplitude of ~1.2‰, believed related to low-latitude hydrological changes (3, 6). However, the ∆DE record contains additional and important variance at the ~100-ka periodicity, suggesting a substantial influence of ice-sheet volume on the ∆DE, possibly via changing terrestrial productivity at mid-to-high latitudes (3). This, therefore, impedes the use of the ∆DE as a proxy for the low-latitude hydroclimate.

In this study, we generated a new ∆DE* (to distinguish from the traditional definition of ∆DE) record by removing the globally stacked sea surface δ18O (δ18Osurf) from δ18Oatm, which shows pronounced precession cycles (~23 and ~19 ka) but almost no ~100-ka periodicity over the past 800 ka. Furthermore, we used the Community Climate System Model version 3 (CCSM3) to perform a 300-ka transient simulation of global climate change. The modeled terrestrial rainfall changes between 30°N and 30°S strongly resemble the ∆DE*. We therefore propose that the low-latitude hydrology is the only important forcing for the ∆DE* on orbital time scales, while the influence of other climate factors, if any, should be minor.

Recalculation of the ∆DE*

We argue that δ18Oatm has a stronger relevance with isotopic values of sea surface compared with those of the entire ocean because (i) the terrestrial biosphere uses waters originally evaporated from sea surface and (ii) nearly 100% photosynthesis and 95% respiration of the marine biosphere occur within ocean mixed layers (7). Therefore, it is more reasonable to use δ18Osurf rather than δ18Osw (or mean oceanic water δ18O) to calculate the ΔDE*. δ18Osurf is different from δ18Osw, because in addition to varying ice volume, hydrological changes at mid-to-low latitudes can also exert substantial influence on δ18Osurf through altering evaporation, precipitation, runoff, and reservoir size of liquid water on land. Moreover, δ18Osurf can be reliably reconstructed by removing temperature effect from planktonic foraminiferal δ18O. Through compiling parallel measurements of planktonic foraminiferal δ18O and sea surface temperature (SST) from marine sediment cores, a global δ18Osurf stack over the past 800 ka has been generated by Shakun et al. (11) (Figs. 1 and 2A). We update the 24 to 0 ka BP (Before Present) interval of this stack by using high-resolution and 14C-constrained reconstructions with Mg/Ca-SST results (Materials and Methods as well as the Supplementary Materials, Fig. 1, and fig. S1).

Fig. 1 Modern global monsoon domain and locations of paleoclimatological reconstructions.

Monsoon domain (green) was defined by Wang and Ding (1). Blue and red lines indicate the mean position of the Intertropical Convergence Zone for August and February, respectively (solid for monsoon trough and dashed for trade wind convergence) (24). Blue (11) and red (this study) dots indicate δ18Osurf reconstructions used for the generation of a global stack. WAIS, West Antarctic Ice Sheet; EPICA, European Project for Ice Coring in Antarctica.

Fig. 2 Comparison of proxy records over the past 800 ka.

(A) A global stack of δ18Osurf (blue, with ±1σ standard error) [this study, (11)] and the modeled δ18Osw (red) (10). (B) Normalized δ18Osurf and δ18Osw are obtained by the z-standard method, which are further used to calculate the difference. (C) A composite δ18Oatm record from Antarctic ice cores, including the WAIS Divide and the Siple Dome (50 to 0 ka BP) (2, 5), the Vostok (100 to 50 ka BP) (12), and the EPICA Dome C (800 to 100 ka BP) (1317). (D) Chinese stalagmite δ18O, compiled by Cheng et al. (18). (E) A compilation of Bornean stalagmite δ18O (1922). The duration of terminations I to V is marked by yellow bars. In (B) to (E), dashed lines indicate the linear regression of each proxy record over the past 430 ka.

Discrepancies between the δ18Osurf stack and the modeled δ18Osw (10) are described in the Supplementary Materials (figs. S1 to S4). Briefly, the δ18Osurf and the δ18Osw show three robust discrepancies over the past 430 ka (Fig. 2). First, with respect to δ18Osw, isotopic values of several interstadials within an interglacial period are more comparable in δ18Osurf [marine isotope stage (MIS) 11c versus 11a, 9e versus 9c and 9a, 5e versus 5c and 5a, Fig. 2A]. Second, at the recent five glacial terminations, the rise of δ18Osurf lags that of δ18Osw by 3 to 4 ka (marked by yellow bars in Fig. 2; also, see figs. S1 to S3). Third, δ18Osurf exhibits a long-term trend toward lighter isotopic values for the past 430 ka, which was also noted by Shakun et al. (11). This secular trend is clearly revealed in the difference between δ18Osurf and δ18Osw, showing a decreasing tendency over this period (Fig. 2B and fig. S4).

The aforementioned features of δ18Osurf are shared by δ18Oatm (1217) and tropical-subtropical stalagmite δ18O (Fig. 2) (1822), including comparable amplitude of interstadials within an interglacial (MIS 11, 9, and 5), positive δ18O excursions at glacial terminations, and the long-term trend over the past 430 ka. Together, the common features among different climate archives can be interpreted as that the distinct isotopic imprint of sea surface was transferred to continental precipitation and leaf water via the hydrological cycle, which was subsequently recorded by stalagmite and δ18Oatm. Therefore, δ18Osurf should be used to adjust precipitation isotopic records (e.g., stalagmite δ18O, plant wax hydrogen isotope) and calculate the ΔDE*.


The new ΔDE* dominated by precession cycles

As illustrated in Fig. 3, compared with the previous estimate of the ΔDE (δ18Oatm − δ18Osw), the new ΔDE* (δ18Oatm − δ18Osurf, Materials and Methods) displays a considerably better match with precession in amplitude over the past 800 ka. An evident mismatch between the previous ΔDE and precession is found between 450 and 360 ka BP (Fig. 3A), which is improved between the ΔDE* and precession but still exists to some extent (Fig. 3C). Note that during this period, δ18Oatm and δ18Osurf contain relatively less variance in the precession bands (Fig. 2, A and C), and astronomically tuned age models of both climate records have relatively large uncertainties (the establishment of chronologies are described in Materials and Methods) (17, 23). Therefore, the ambiguity of precession signals in climate reconstructions together with age model uncertainties may partly explain this mismatch.

Fig. 3 Comparison of two different estimates of the ΔDE.

(A) The previous estimate of ΔDE (red) and precession (gray). (B) Eccentricity cycles. (C) The new estimate of ΔDE* (blue, with ±1σ standard error) and precession (gray). (D) The global stack of benthic foraminiferal δ18O (23). Eccentricity, obliquity, and precession are derived from Berger (41), which are further normalized to calculate the ETP (E + T-P). Cross-spectral analysis results of the ETP with the previous (E) and the new (F) estimate of ΔDE, respectively, over the past 800 ka. In (E) and (F), the spectrum of each record is presented with the 95% confidence level (dashed lines). Below, coherency spectra are indicated by gray-filled curves associated with the 95% Monte Carlo false-alarm level (black dashed lines). The time-series analysis was performed using the REDFIT program (46). The numbers denote primary orbital cycles.

The ΔDE* shows no evident imprint of an individual glacial-interglacial cycle, nor does it respond to an increase in the magnitude of glacial-interglacial cycles that occurred after glacial termination V, referred to as the Mid-Brunhes Event. In contrast to the previous estimate of ΔDE, the spectrum of the ΔDE* shows strong variance in the ~19-, ~23-, and ~41-ka bands but lacks variance in the ~100-ka bands (Fig. 3F). This suggests that the presence of the ~100-ka periodicity in the previous ΔDE (Fig. 3E) (3) is an artificial signal. Together, precession and obliquity to a lesser extent are dominant forcing of the ΔDE* over the past 800 ka, while the influence of 100-ka cycles (ice volume) seems negligible.

The ΔDE* exhibits a remarkable similarity to Chinese stalagmite δ18O (hereafter refers to the one adjusted for δ18Osurf changes, simply by subtracting δ18Osurf from stalagmite δ18O) over the past 640 ka (R = 0.60; P < 0.001), both following the temporal rhythm of the Northern Hemisphere summer insolation (Fig. 4). Their spectra are also well matched, both showing strong variance in the precession bands (Fig. 4D). In addition, the Mid-Brunhes Event is not reflected in Chinese stalagmite δ18O reconstructions, either (Fig. 4A) (18).

Fig. 4 Comparison of the ΔDE*, Chinese stalagmite δ18O, and simulated rainfall changes.

(A) The ΔDE* (blue, with ±1σ standard error) and Chinese stalagmite δ18O (orange, after adjusted for changes in δ18Osurf) (17). (B) Simulated terrestrial rainfall changes between 30°N and 30°S (annual mean, green) and July 21 daily insolation at 20°N (gray) (41). (C) Simulated terrestrial rainfall changes over 0 to 30°N (annual mean, red) and 0 to 30°S (annual mean, blue). Simulation results are from the experiment CCSM3_orb+ghg + ice (Materials and Methods), which are further fitted using the Savitzky-Golay algorithm with 15 points at second order. In (D) and (E), spectra of the 800-ka ΔDE*, the 640-ka Chinese stalagmite δ18O, and the 300-ka simulation output are presented with the 95% confidence level (dashed lines). Below, coherency spectra of their cross-spectral analysis are indicated by gray-filled curves associated with the 95% Monte Carlo false-alarm level (black dashed lines). The time-series analysis was performed using the REDFIT program (46). The numbers denote primary orbital cycles.

Transient simulation of the hydrological cycle

An accelerated simulation was performed using the CCSM3 from 300 to 0 ka BP, forced by orbital-driven insolation, greenhouse gases, and ice sheets (Exp_orb+ghg + ice; Materials and Methods and figs. S5 and S6). To evaluate the effect of different forcing on the hydrological cycle, two additional experiments were run, forced by orbital parameters (Exp_orb) and orbital parameters plus greenhouse gases (Exp_orb+ghg), respectively (Materials and Methods and figs. S5 and S7). Terrestrial rainfall changes in a specific region are indicated by their proportion in the global total.

In the model output of the Exp_orb+ghg + ice (fig. S6), except for those over 60°S to 90°S, terrestrial rainfall changes over other latitudes all contain strong variance in the precession bands. Especially over 0 to 30°N, precession is the predominated cycle for the past 300 ka. Terrestrial rainfall changes over 30°N to 60°N and 60°N to 90°N also include additional but less important variance in the 100-ka bands. Our simulation successfully reproduces a fundamental feature of the global hydrological cycle as revealed in many geological records [e.g., (24)] that terrestrial rainfall changes between Northern and Southern Hemisphere low latitudes are antiphased in the precession bands (Fig. 4C and figs. S6 and S8).

Because terrestrial rainfall changes over 0 to 30°N show stronger variations than over 0 to 30°S (Fig. 4C), the combined rainfall changes between 30°N and 30°S show strong variance only in the precession bands and follow the temporal rhythm of the Northern Hemisphere summer insolation (Fig. 4B). Using cross-spectral analysis, simulated rainfall changes between 30°N and 30°S are strongly coherent with the ΔDE* in the precession bands (Fig. 4E). Negative shifts of the ΔDE* (and Chinese stalagmite δ18O) are always concurrent with large terrestrial rainfall between 30°N and 30°S, and vice versa. Moreover, due to a high fraction of tropical rainfall in the global total, temporal changes of the global terrestrial rainfall are also highly similar to those over 0 to 30°N, showing a domination of precession cycles (fig. S6).

A comparison of model output reveals that all three experiments generate similar results on terrestrial rainfall changes in tropics (fig. S7). This suggests that solar insolation is much more important than greenhouse gases and ice volume in regulating the temporal pattern of tropical hydroclimate. Including the ice-volume forcing in the model can enhance and weaken the variance in the precession bands over 0 to 30°N and 0 to 30°S, respectively (fig. S7).


Causes for the ΔDE* on orbital time scales

Previously, three factors were considered important for causing the Dole effect fluctuation on orbital time scales: changes in the ratio of terrestrial to marine productivity (6, 7), vegetation changes induced respectively by ice volume (3), and the low-latitude hydrological cycle [e.g., (2, 6)]. Estimates of paleoproductivity are still subject to large uncertainties. Nevertheless, it is generally accepted that changes in the global biogenic productivity were significantly affected by glacial-interglacial cycles during the late Quaternary (25). With respect to interglacial periods, a substantial reduction of terrestrial productivity and a slight increase (or unchanged) of marine productivity during glacial time have been found (25, 26). Because no evident imprint of glacial-interglacial cycles (100-ka periodicity) is detected in the ΔDE* as aforementioned, past changes in the land/sea productivity ratio seem not critical for the ΔDE*. This inference is supported by recent studies that suggest a nearly equivalent Dole effect value generated by modern terrestrial and marine biosphere [e.g., (8)]. This is in contrast to the previous estimate, which considered a higher terrestrial Dole effect than the marine value (6, 7).

The absence of 100-ka glacial-interglacial cycles in the ΔDE* also excludes an important role of the ice volume. In high latitudes, the photosynthetic Dole effect of boreal biomes is partially canceled out by the respiratory effect in soils (27). Therefore, pronounced changes in the amount of boreal biomes due to the expansion and the retreat of ice sheets are inferred to exert limited influence on the global ΔDE* (8).

The dominance of precession cycles in the new ΔDE*, however, highlights the role of the low-latitude hydrology and vegetation as previously thought [e.g., (2, 6)]. Tropical and subtropical vegetation comprises the major part of the global productivity on land (4). On one hand, the respiratory Dole effect in low latitudes is considered lower than the global mean value due to slow diffusion of O2 in wet soils [e.g., (28)]. A strengthened hydrological cycle and subsequently an expansion of wetlands can therefore decrease the global respiratory Dole effect. On the other hand, a strengthened hydrological cycle is associated with the enhanced convection intensity over oceanic moisture source areas, the increased isotopic fractionation during the ocean-to-land moisture transfer, and the expansion of monsoon realms [e.g., (29)]. All these processes can result in the isotopic depletion of meteoric waters on land. The substantial impact of precipitation isotope composition on δ18Oatm has been clearly revealed by the strong similarity between Chinese stalagmite δ18O and the ΔDE* (Fig. 4). Furthermore, elevated humidity levels over vegetated areas can decrease the evaporative enrichment of leaf water δ18O [e.g., (6)]. Together, the respiratory and the photosynthetic Dole effect of the low latitude both decrease during periods with a strengthened hydrological cycle.

Apparently, our simulation supports the correlation between the tropical hydroclimate and the ΔDE* on orbital time scales. Given its high proportion in the global terrestrial rainfall and its spectral distribution, the tropical hydrological cycle can most likely explain the ΔDE* (fig. S7). In addition to stronger rainfall changes over 0 to 30°N than over 0 to 30°S, generally stronger continentality and higher altitudes in the Northern Hemisphere can decrease meteoric water δ18O values to a greater extent due to the moisture isotopic fractionation along transport paths (2). Therefore, the northern low-latitude Dole effect signal should overwhelm the southern counterpart, and the global ΔDE thus follows the temporal pattern of the Northern Hemisphere summer insolation (Fig. 4).

Obliquity cycles of the ΔDE*

The obliquity cycle is recognizable in the ΔDE* but nearly undetectable in Chinese stalagmite δ18O (Fig. 4D), suggesting that meteoric water δ18O is a very important but not the only cause for the ΔDE* as aforementioned. Although terrestrial rainfall changes over northern and southern high latitudes include weak variance in the obliquity bands as shown in our model output (figs. S6 and S9), they unlikely can account for the obliquity cycle of the ΔDE* because of their small contribution to the global terrestrial rainfall (figs. S6). We thus exclude high latitude-sourced causes for the obliquity cycle in the ΔDE*.

In contrast to our transient simulation, some equilibrium simulations have revealed that large obliquity can enhance summer monsoon intensity in both hemispheres, particularly in the Northern Hemisphere [e.g., (30)]. Although it is model dependent, the effect of obliquity has been clearly observed in some regional monsoon reconstructions [e.g., (30)]. Therefore, the obliquity cycle of the ΔDE* can also be interpreted in the context of the low-latitude hydrology. High-obliquity values should correspond to relatively low Dole effect values, as confirmed by our data analysis (fig. S9). To summarize, the hydroclimate over low-latitude vegetated areas should be the only dominant cause for the orbital-scale ΔDE* over the past 800 ka.

The ΔDE* as a proxy for the low-latitude hydrology

Among a few different but interrelated forms of the low-latitude hydrological cycle (e.g., the El Niño–Southern Oscillations and the Intertropical Convergence Zone), monsoons are recognized as the most mutable component and represent the dominant mode of annual low-latitude variation (Fig. 1) (1). Moreover, a compilation of geological records has confirmed the coherent variability of regional monsoons across orbital to millennial time scales (Fig. 1) (24). This has further led to an evolving concept, the geological evolution of “global monsoon” [e.g., (24)]. Because monsoons are viewed as a planetary-scale rather than a regional phenomenon, it would be significant if a common proxy could be found to describe the global monsoon or low-latitude hydrological changes. Because of its global character and strong correlation with simulated tropical rainfall changes (Fig. 4), the ΔDE* is suggested as such a potential proxy on orbital time scales.

The low-latitude hydrological cycle as a whole, as revealed by the ΔDE*, appears to be governed by insolation forcing only. This differs from regional hydroclimate reconstructions, which show substantial response to a set of other factors besides solar insolation, including changing ice volumes, greenhouse gas concentrations, SST, sea-land distribution, and orography (31). Therefore, a global proxy is in favor of going beyond regional character and complexity and is advantageous to identify common and most fundamental mechanisms for the operation of the low-latitude hydroclimate at a large spatial extent.

The idea of questing for a global monsoon proxy is similar to using benthic foraminiferal δ18O as an indicator for the global glacial-interglacial cycles. Although regional ice sheets can have different timing of advance and retreat, benthic foraminiferal δ18O has been used as a globally integrated proxy to depict the cyclicality, duration, and transition of Quaternary glacial-interglacial climate. Likewise, the ΔDE* may become a benchmark to measure the dynamics of the low-latitude hydrological cycle during the Quaternary.


The δ18Osurf stack over the past 24 ka

Paired measurements of planktonic foraminiferal δ18O (δ18Oc) and Mg/Ca-SST with 14C-constrained chronology from 36 cores were compiled to generate a global stack of δ18Osurf (Fig. 1, fig. S1, and table S1). We excluded published alkenone-SST reconstructions. Although, in many regions, Mg/Ca-SST and alkenone-SST records are largely consistent on orbital time scales, they show discrepancies on millennial time scales mainly due to their different seasonal preferences (32). δ18Oc and Mg/Ca-SST are both sourced from calcareous shells, which can therefore provide more robust reconstructions of δ18Osurf than using alkenone-SST. Through applying the Marine13 dataset to recalibrate 14C dates (33), the age model of nearly all records was readjusted. One record was left on its original 14C and ice-core tuned age model (MD97-2120). The time resolution of most reconstructions is <500 years (n = 32), while four other records range between 520 and 580 years. In each core, δ18Osurf [‰, SMOW (Standard Mean Ocean Water)] was calculated by (SST-16.5)/4.8 + δ18Oc [‰, PDB (Pee Dee Belemnite)] + 0.27 (34). An evenly spaced δ18Osurf with a resolution of 500 years was further produced by averaging data in 1000-year bins. For stacking, each δ18Osurf record was shifted to a mean of zero and combined as unweighted global averages. Last, because only the amplitude of δ18Osurf changes is concerned, both δ18Osurf stacks produced in this study and by Shakun et al. (11) were shifted with a late Holocene value of zero (Fig. 2). Considered propagating uncertainties of δ18Oc and Mg/Ca measurements, the Mg/Ca-SST calibration, and the seawater δ18O-paleotemperature equation, the standard error (±1σ) of a single δ18Osurf estimate is approximately ±0.30‰ (35). The standard error (±1σ) of the δ18Osurf stack is between ±0.06 and ± 0.08‰, depending on the number of available data points.

Chronologies of climate reconstructions

The chronology of Chinese stalagmite records was established by radiometric dating, with a typical error of <1 and 1 to 3.5 ka (1σ), respectively, for the interval 400 to 0 ka BP and 640 to 400 ka BP (18). Age models of marine and ice-core records are basically derived from the astronomical tuning, with an uncertainty less than 4 ka for most periods of the past 800 ka (23, 36, 37). The global stack of δ18Osurf (800 to 23.5 ka BP) (11) and the modeled δ18Osw (800 to 0 ka BP) (10) are both plotted against the LR04 age model (23). The δ18Oatm record is plotted against the AICC2012 age model (36, 37), with that between 640 and 100 ka BP being retuned to the Chinese stalagmite chronology according to the similarity between δ18Oatm and stalagmite δ18O (17). Applying the AICC2012 age model to the entire δ18Oatm record, however, would not change the spectral analysis results of the ΔDE and the ΔDE* and thus the conclusions of this study. Because the δ18Oatm record contains strong variance in the precession bands, and the recognition of precession cycles in δ18Oatm is additionally constrained by a few discrete absolute dating results (36, 37), the spectral analysis of δ18Oatm and the ΔDE* does not seem to be substantially influenced by the tuning method.

Recalculation of the ΔDE

We followed the approach of Extier et al. (17). Series of δ18Osw, δ18Osurf, and δ18Oatm were first interpolated to 100-year resolution and then fitted using the Savitzky-Golay algorithm with 31 points at second order. Fitted curves of δ18Osw and δ18Osurf were extracted from that of δ18Oatm, respectively, to calculate the ΔDE and the ΔDE*. Given the standard error of the δ18Osurf stack (0.05 to 0.12‰, 1σ) [this study, (11)] and the δ18Oatm measurement (±0.04‰, 1σ) (13), the standard error of the ΔDE* was estimated to be around 0.07 to 0.15‰ (1σ).

Model simulations

The CCSM3 includes dynamically coupled atmosphere, ocean, land, and sea ice components (38, 39). The atmosphere model has a ~3.75° horizontal resolution (T31, varying from 340 km at the equator to 40 km around Greenland) and 26 hybrid coordinate levels in the vertical. The land model has the same T31 resolution. The ocean model has a 3° horizontal resolution with 25 vertical levels. The resolution of the sea ice model is identical to that of the ocean model, including a subgrid-scale ice thickness distribution.

To distinguish the effect of each forcing on the global climate, three experiments were designed, forced by orbital parameters (Exp_orb), orbital parameters combined with greenhouse gases (Exp_orb+ghg), and orbital parameters combined with greenhouse gases and ice volume (Exp_orb+ghg + ice), respectively (figs. S5 and S7). In the first two experiments, orbital parameters and greenhouse gases were advanced by 100 years at the end of each model year. In the third experiment, the continental ice-sheet volume and the position of coastlines were altered at steps of an equivalent 40-m sea level rise or fall to reconfigure and restart the model. The upper ocean and the atmosphere are considered to reach quasi-equilibrium in these accelerated simulations (40). Annual-mean rainfall changes over lands of six latitude bins, 90°N to 60°N, 60°N to 30°N, 30°N to 0°, 0° to 30°S, 30°S to 60°S, and 60°S to 90°S, were calculated and shown in this study. Terrestrial rainfall changes in a specific region were represented by their proportion in the global total (Fig. 4 and figs. S6 and S7).

As shown in fig. S5, the orbital and greenhouse gases forcing were reconstructed from Berger (41) and the Antarctic ice core records (42, 43), respectively. Continental ice-volume and sea-level changes were calculated by scaling the ICE-5G ice distributions for the past 21 ka (44) to the globally stacked benthic foraminiferal δ18O record (23). A part of simulation results has been published by Lu et al. (45).

We note that the Exp_orb+ghg + ice output contains some suspicious signals, including pulsed changes during the recent three glacial terminations (fig. S7, A and B), and relatively weak terrestrial rainfall during interglacial peaks (MIS 7e, 5e, and 1) over 0 to 30°N (fig. S7A). The former is induced by a rapid jump in the ice-volume forcing (fig. S5), and the atmosphere and surface ocean in the model need time to reach quasi-equilibrium. The latter seems to contradict with climate reconstructions, which indicate that interglacial peaks were always characterized by very humid conditions over the Northern Hemisphere tropical lands [e.g., (36)]. In this study, we mainly focus on astronomical cycles of simulation output, which is therefore not subject to these unreal signals.


Supplementary material for this article is available at

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.


Acknowledgments: Funding: This research is funded by the National Science Foundation of China (nos. 41776054, 41525020, 91128208, 41606045, and 41976047) and the National Key Research and Development Program (2018YFE0202401). Author contributions: E.H. conceived the study, compiled all data shown in this study, and wrote the manuscript. M.Y. and Y.W. analyzed the CCSM3 simulation results. S.L. generated the global stack of sea surface oxygen isotope record. All authors discussed the results and commented on 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.

Stay Connected to Science Advances

Navigate This Article