Does δ 18 O of O 2 record meridional shifts in tropical rainfall?

. Marine sediments, speleothems, paleo-lake elevations, and ice core methane and δ 18 O of O 2 ( δ 18 O atm ) records provide ample evidence for repeated abrupt meridional shifts in tropical rainfall belts throughout the last glacial cycle. To improve understanding of the impact of abrupt events on the global terrestrial biosphere, we present composite records of δ 18 O atm and inferred changes in fractionation by the global terrestrial biosphere ( (cid:49)ε LAND ) from discrete gas measurements in the WAIS Divide (WD) and Siple Dome (SD) Antarctic ice cores. On the common WD timescale, it is evident that maxima in (cid:49)ε LAND are synchronous with or shortly follow small-amplitude WD CH 4 peaks that occur within Heinrich stadials 1, 2, 4, and 5 – periods of low atmospheric CH 4 concentrations. These local CH 4 maxima have been suggested as markers of abrupt climate responses to Heinrich events. Based on our analysis of the


Introduction
The last glacial cycle shows two modes of abrupt climate variability, both with a center of action in the North Atlantic. The first mode is Dansgaard-Oeschger (D-O) events, most clearly identified in Greenland ice core records (Andersen et al., 2004;Dansgaard et al., 1993). D-O events refer to abrupt Northern Hemisphere warming, on the timescale of decades, followed by more gradual cooling. In the Southern Hemisphere, D-O events exhibit the opposite response, cooling followed by warming, lagging the Northern Hemisphere by several centuries on average (Buizert et al., 2015a). Most D-O events can be linked to abrupt variations in the strength and/or northward heat transport of the Atlantic Meridional Overturning Circulation (Lynch-Stieglitz, 2017). The second mode of abrupt climate change is Heinrich events (HEs), periods of extreme cold in the North Atlantic associated with extensive layers of ice-rafted detritus in ocean. Heinrich events occur within longer Heinrich stadial (HS) periods (Hemming, 2004).
High-resolution records of marine sediments (Peterson et al., 2000), ice core methane (Rhodes et al., 2015), speleothem calcite δ 18 O (Cheng et al., 2016;Kanner et al., 2012;Wang et al., 2017Wang et al., , 2008Wang et al., , 2001, and lake elevations (Goldsmith et al., 2017;Yu et al., 2013) all demonstrate that abrupt meridional shifts in tropical rainfall belts occurred repeatedly during the last glacial cycle in response to both D-O events and HEs. This response is also observed in climate models, which consistently simulate a shift of tropical rainfall towards the warmer hemisphere (Broccoli et al., 2006;Chiang and Bitz, 2005;Cvijanovic et al., 2013) ing the impact of these events on global hydrology, monsoon systems, and terrestrial biosphere remains a key goal of paleoclimate research.
Measurements of long-lived atmospheric gases trapped in ice cores offer an opportunity to better constrain past global processes due to the short atmospheric mixing time and therefore globally integrated signal of these gases. Measurements of concentrations and isotopic ratios of atmospheric gases have provided important clues into the biogeochemical responses to substantial climate changes in the past. Motivated by the recent high-resolution West Antarctic Ice Sheet (WAIS) Divide (WD) ice core methane record (Rhodes et al., 2015), in this study we explore a new record of atmospheric oxygen isotopes from WD and revisit published measurements in the Siple Dome (SD) ice core. A composite WD-SD record of atmospheric oxygen isotopes allows us to investigate past changes in the isotopic fractionation of O 2 on a common timescale (Buizert et al., 2015b) with other recent records.
Past changes in the Dole effect -the amount by which δ 18 O of atmospheric O 2 (δ 18 O atm ) exceeds that of seawater (Dole, 1935;Morita, 1935) -have been shown to be closely linked to terrestrial hydroclimate (Bender et al., 1994;Landais et al., 2010;Reutenauer et al., 2015;Severinghaus et al., 2009). We quantitatively explore a simple relationship between meridional shifts in terrestrial oxygen production and δ 18 O of oxygen production-weighted terrestrial rainfall (δ 18 O precip ) over the modern seasonal cycle using spatially gridded monthly observations. We use these observations as a guide to consider instantaneous changes in O 2 fractionation during HEs and discuss the implications of these changes, and related changes in atmospheric CH 4 , for the meridional distribution of tropical rainfall associated with these events.

The Siple Dome and WAIS Divide ice cores
The 1004 m long SD ice core (81.65 • S, 148.81 • W) was drilled in the late 1990 s at an elevation of 621 m a.s.l. The site has modern-day mean-annual surface temperature and accumulation rates of −24.5 • C and 13.5 cm ice equivalent per year, respectively (Taylor et al., 2004). The youngest 8.2 kyr of the record were dated by visual and electrical counting of annual layers (Taylor et al., 2004), tuned to stratigraphic markers, while the older portion of the ice core was dated by synchronization of CH 4 and δ 18 O atm measurements to the GISP2, Greenland ice core (Brook et al., 2005). The 3404 m long WD ice core (79.48 • S, 112.11 • W) was drilled at an elevation of 1766 m between 2007 and 2011. Present mean-annual air temperature and accumulation rate at the WD ice core site are −30 • C and 22 cm ice equivalent per year (Banta et al., 2008;Morse et al., 2002). The youngest 31.2 kyr (2850 m) was dated by annual-layer counting (Sigl et al., 2016), while the oldest portion of the core (up to ∼ 68 ka) was dated by CH 4 synchronization to the Greenland NGRIP ice core, scaled linearly to match the timing of abrupt events recorded in U/Th-dated speleothem records (Buizert et al., 2015b).

δ 18 O atm and the Dole effect: a brief overview
The enrichment of δ 18 O atm relative to δ 18 O of mean seawater (δ 18 O sw ) is known as the Dole effect (Dole, 1935;Morita, 1935) and has a modern-day value of ∼ 23.88 ‰ . Most of the Dole effect is attributable to the isotopic discrimination of marine and terrestrial respiration, both of which preferentially utilize 16 O, thus enriching atmospheric O 2 in 18 O (Bender et al., 1994;Luz and Barkan, 2011). Photosynthesis does not fractionate the 18 O/ 16 O ratio of chloroplast water (e.g., Helman et al., 2005), and therefore the δ 18 O of O 2 produced by photosynthesis is equal to the δ 18 O of the plant source water (ultimately derived from precipitation) plus fractionation due to evapotranspiration. Past changes in the Dole effect estimated from measurements of atmospheric oxygen isotopes in ice core bubbles and reconstructions of δ 18 O sw from marine sediment cores have been shown to be small over the past 130 kyr, during which the SD of the Dole effect was only 0.24 ‰ (Bender et al., 1994).
Variability in the Dole effect throughout the past glacialinterglacial cycle has been hypothesized to be related to terrestrial hydrology based on coherence with atmospheric CH 4 measured in ice core bubbles (Bender et al., 1994), comparison with Chinese cave records (Severinghaus et al., 2009), and climate modeling experiments (Reutenauer et al., 2015). To investigate instantaneous changes in the isotopic fractionation of atmospheric O 2 , Severinghaus et al. (2009) used a one-box model deconvolution to account for the smoothing effect of the ∼ 1000-year residence time of atmospheric O 2 . The parameter ε LAND , associated with this box model, approximately represents changes from the present in globally integrated (mostly) terrestrial O 2 isotopic fractionation: where ε L refers to the change from present in the fractionation of terrestrially produced oxygen relative to seawater, ε RL and ε RO respectively refer to changes from present in terrestrial and marine effective respiratory fractionation, and f L is the fraction of oxygenesis occurring on land. Because the three terms in Eq. (1) are not independently discernible, Severinghaus et al. (2009) derive a formula for ε LAND in terms of knowable quantities: δ 18 O sw , δ 18 O atm , and time derivative of δ 18 O atm (dδ 18 O atm dt −1 ). The calculation of ε LAND in terms of these variables is described in Sect. 2.1. ε LAND , as estimated from SD δ 18 O atm measurements, is strongly correlated with Dongge cave δ 18 O (Wang et al., 2005(Wang et al., , 2001) over the past 12 kyr (Severinghaus et al., 2009). Before this period, uncertainty in the original SD gas timescale prevented reliable high-precision comparison with other paleoclimate records. Improving the timescale for SD δ 18 O atm measurements is therefore a key goal of this work.
Modeling results from a recent freshwater hosing experiment suggest that changes in δ 18 O atm over a HS are dominated by changes in δ 18 O precip over the terrestrial biosphere (Reutenauer et al., 2015). In terms of Eq. (1), this finding suggests that ε LAND may be driven by ε L . Below we build on the conclusion of this recent modeling study by exploring a composite WD-SD ε LAND record on an improved timescale and drawing insights from the modern seasonal cycles of terrestrial oxygenesis and δ 18 O precip . In order to evaluate a smooth derivative for the calculation of ε LAND , discrete WD δ 18 O atm values were fit to a Fourier series ranging from 1/697 to 3.2 cycles kyr −1 (spacing: 1/697 cycles kyr −1 ) by a Bayesian, weighted, linear least-squares technique, which a priori assumes a red spectrum and factors this assumption into the cost function (Severinghaus et al., 2009; Appendix A). The time derivatives of SD (on its original timescale) and WD δ 18 O atm fitted curves were used to determine depth-gas age tie points between the two cores.
After synchronization of SD gas ages to the WD2014 timescale (Sect. 2.2), a composite fitted curve and time derivative of δ 18 O atm were calculated using all discrete SD and WD δ 18 O atm younger than 50.12 ka BP. ε LAND was calculated independently in three ways, from fitted δ 18 O atm curves and derivatives in (1) SD only, (2) WD only, and (3) the composite SD-WD record, each time using the same smoothed, interpolated δ 18 O sw record (Waelbroeck et al., 2002): where f L is the fraction of photosynthesis occurring on land (assumed to be 0.65) and τ is the residence time of atmospheric O 2 (1 kyr).

Synchronization of Siple Dome gas ages to WAIS Divide
Siple Dome gas ages between 50 ka and 1950 CE were synchronized with the WAIS Divide WD2014 timescale (Buizert et al., 2015b) using the smooth annual layer thickness (ALT) method of Fudge et al. (2014). The ALT method, as applied in this study, is a standard regularized linear least-squares method (Aster et al., 2013) that estimates an annual layer thickness history constrained by depth-gas age tie points and smoothness of the ALT time series. Integration of the resulting time series of annual layer thicknesses yields a gas agedepth chronology. Formally, the ALT method seeks an optimal annual layer thickness history, m, that minimizes the cost function C as defined: where G is a matrix that maps annual layer thicknesses m to depths d corresponding to tie points of a given gas age. L is a second-derivative operator and α is a trade-off parameter between smoothness and tie-point agreement (the notation || || 2 2 indicates the squared L2 norm of the argument). Depth constraints are weighted based on individual estimates of tiepoint uncertainty, σ , as discussed in detail in Appendix B.
In total, 56 tie points were used to constrain the ALT inverse problem: 36 from SD and WD records of dδ 18 O atm dt −1 , and 20 from SD and WD CH 4 records. Because δ 18 O atm is smoothly varying on millennial timescales due to the ∼ 1 kyr atmospheric residence time of O 2 , WD and SD dδ 18 O atm dt −1 time series (calculated from fitted curves of δ 18 O atm at 10-year resolution, Appendix A) were instead used for synchronization. Tie points based on dδ 18 O atm dt −1 were determined by calculating the midpoint times of abrupt transitions common to both records. CH 4 -based tie points were determined by value-matching unambiguous abrupt transitions identified in the two records, within a prescribed error range of ±15 ppb. More complete explanations of the tie-point selection process, an estimate of its uncertainty, the application of the ALT method, and sensitivity tests are provided in Appendix B. Figure 1 shows the fitted SD δ 18 O atm curve on its original (Brook et al., 2005) and WD2014synchronized timescales alongside the WD δ 18 O atm discrete measurements and fitted curve. Figure 2 shows the synchronization results as a comparison between the original and new SD gas timescales, CH 4 and dδ 18 O atm dt −1 tie points, and the annual (gas) layer thickness profile of the new and original timescales.

Modern seasonal cycle analysis
To explore seasonal changes in the precipitation supplied to terrestrial oxygen production (Sect. 3.2), modern mean monthly gridded datasets of gross primary productivity (GPP; Jung et al. 2011) and δ 18 O precip (Terzer et al., 2013) were analyzed. The δ 18 O precip isoscape product used    (Brook et al., 2005) and WD2014-synchronized timescale (this study), and gas age differences between the original and WD2014-synchronized Siple Dome gas timescales (b) constrained by tie points selected from abrupt transitions in the CH 4 and dδ 18 O atm dt −1 records at SD and WD. for our analysis (Sect. 3.2) was created from regionalized cluster-based water isotope prediction (RCWIP) constrained by 57 000 Global Network of Isotopes in Precipitation (GNIP) measurements from 1960 to 2009 (IAEA, 2017; IAEA/WMO, 2017; Terzer et al., 2013). The gridded monthly mean GPP product used in our study was determined from upscaled eddy covariance observations (FLUXNET) from 1982 to 2011 (Jung et al., 2011).
To estimate monthly terrestrial GPP-weighted δ 18 O precip , monthly-mean isoscapes (originally 1/6 • × 1/6 • ) were regridded using MATLAB's geoloc2grid function to 0.5 • × 0.5 • resolution to match the resolution of the GPP dataset. GPP data were converted from carbon (or, equivalently, oxygen after accounting for a photosynthetic quotient) fluxes to total emissions by multiplying by grid cell area. In Eqs. (4) and (5) below, GPP refers to an emission, rather than a flux, and is therefore independent of latitude-varying grid cell area. Only land area between 60 • S and 84 • N was included due to data availability, with several gaps over low-GPP regions like the Sahara and Greenland. Monthly global terrestrial GPP-weighted mean δ 18 O precip was calculated as follows: where δ 18 O precip (t) is terrestrial GPP-weighted mean δ 18 O precip during month t, and ϕ and λ refer to center latitudes and longitudes of terrestrial grid cells with spacing dλ = dϕ = 0.5 • . The centroid latitude of terrestrial oxygenesis (ϕ TOE ) is defined as the latitude at which f GPP , the meridionally cumulative fraction total monthly GPP, equals 0.5. Integrating northward from 60 • S, we define ϕ TOE such that In practice, we estimate ϕ TOE at finer resolution than the 0.5 • spacing by first calculating zonal integrals around each band of grid cells and then linearly interpolating such that dϕ in Eq. (5) equals 0.006.

∆ε LAND variations over past 50 kyr
The WD-SD composite ε LAND record is shown in Fig. 3 alongside a composite East Asian speleothem record of calcite δ 18 O (Cheng et al., 2016) and discrete and continuous CH 4 records from WD (Rhodes et al., 2015;Sowers, 2012). The speleothem δ 18 O record has been corrected for seawater δ 18 O changes (Waelbroeck et al., 2002) and is shown both unsmoothed and smoothed (200-year boxcar filter). There are several notable features of the composite ε LAND time series and its comparison to these other records.  (Rhodes et al., 2015;Sowers, 2012). Shaded regions indicate periods of climate impact due to Heinrich events 1, 2, 4, and 5 identified by Rhodes et al. (2015), coincident with the highest observed glacial values of ε LAND .
First, ε LAND is positively and significantly correlated with the seawater-corrected East Asian speleothem δ 18 O record (r = 0.70, p < 0.0001, both records interpolated linearly to 200-year resolution between 50 and 0 ka BP). The regression coefficient, β, over the entire 50 kyr time window is 0.215 ‰ atm ‰ −1 pdb . β varies between 0.387 ‰ atm ‰ −1 pdb during the Holocene (0-12 ka) to 0.168 ‰ atm ‰ −1 pdb during the last glacial period (12-50 ka). Although the reason for the glacial-interglacial change in β is unclear and a worthy topic of future study, the high correlation (r = 0.62, p < 0.0001) beyond the ∼ 0-12 kyr window originally considered by Severinghaus et al. (2009) confirms that the ε LAND and the Asian cave records share common variability throughout the last glacial period in addition to the Holocene. We note that the misalignment of ε LAND and Chinese cave δ 18 O excursions during HS2 may indicate a dating error in the WAIS Divide timescale. The remainder of our analysis is focused on comparison with the WAIS Divide CH 4 record which shares a common timescale (WD2014) with the composite ε LAND record. Absolute dating errors in WD2014 therefore do not hamper comparison of these two records.
Second, over the portion of the last glacial period covered by this record (∼ 50-12 ka BP), ε LAND only consistently exceeds 0 ‰ (for more than a century) during the Younger Dryas and HS 1, 2, 4, and 5. ε LAND exhibits pronounced abrupt increases (downward in Fig. 3) during HEs. Following HE 1, 2, 4, and 5 -all of predominant Hudson Strait origin (Hemming, 2004) -ε LAND abruptly increases by ∼ 0.1-0.3 ‰ over 200 to 300 years and remains generally high for ∼ 1-2000 years. The timings of ε LAND maxima during each of HS 1, 2, 4, and 5 all occur within the proposed time windows of abrupt, large-scale climate responses to each HE (Rhodes et al., 2015;Fig. 4). Similarly, these proposed windows of HE climate correspond to periods of sustained above-average ε LAND values. Uncertainty in the measurement of δ 18 O atm and the computation of its time derivative preclude meaningful analysis of centennial-scale changes in ε LAND of below ∼ 0.05 ‰.
Increases in ε LAND over HSs stand in contrast to repeated decreases during D-O warmings. Notably, atmospheric CH 4 increases during both D-O warmings and HE 1, 2, 4, and 5 (Northern Hemisphere cooling events). Thus, the relationship between ε LAND and atmospheric CH 4 reverses in sign between Northern Hemisphere warming and cooling. In Sect. 4 we discuss the implications of this reversal for inwww.clim-past.net/13/1323/2017/ Figure 4. Atmospheric CH 4 measured at WAIS Divide and composite WD-SD ε LAND throughout and following Heinrich events 1, 2, 4, and 5. ε LAND is highlighted green during periods of climate impact from Heinrich events 1, 2, 4, and 5 proposed by Rhodes et al. (2015), beginning with an abrupt CH 4 rise (blue dashed line with shaded ±1σ dating uncertainty) and ending with another CH 4 rise of a presumed northern source (red dashed line with shaded ±1σ dating uncertainty).
ference of the position of the thermal equator and tropical rainfall belts.

Seasonal cycle of terrestrial photosynthetic source water δ 18 O
We analyzed monthly-mean datasets of δ 18 O precip and GPP on land (Sect. 2.3) to test the sensitivity of GPP-weighted mean δ 18 O precip to the median latitude of oxygenesis over the seasonal cycle, which serves as a rough analogue for past shifts in the thermal equator. Figure 5 shows December-January-February (DJF) and June-July-August (JJA) mean spatial distributions of δ 18 O precip and GPP. Throughout all seasons, δ 18 O precip generally trends isotopically lighter northward between the southern mid-latitudes and northern high latitudes. The meridional trend in GPP, however, exhibits a strong seasonal dependence in its sign. As the northern terrestrial biosphere grows from boreal spring to summer, a northern maximum develops in the mid-to high latitudes (Fig. 6, Movie S2), shifting the centroid of terrestrial oxygen production northward. During boreal fall and winter, Northern Hemisphere productivity decreases, shifting terrestrial oxygenesis southward toward the relatively seasonally invariable southern tropical maximum near 5 • S. For the purpose of quantifying meridional shifts in terrestrial GPP, we define the terrestrial oxygenesis equator (TOE) as the latitude of the terrestrial GPP centroid (Fig. S1,  (Terzer et al., 2013). Data gaps (e.g., over oceans, Greenland, Sahara) are left unfilled in this plot. Figure 6. Zonally integrated terrestrial gross primary productivity as a function of latitude for monthly-mean data averaged over years from 1982 to 2011 (Jung et al., 2011). The dashed red line shows monthly-mean TOE latitude.
Sect. 2.3). In simpler terms, half of all terrestrial oxygen production occurs to the north and half to the south of the seasonally varying TOE. The TOE shifts from ∼ 5 • S in boreal winter to ∼ 35 • N in boreal summer (Fig. 7, Movie S1).
To determine the direct influence of δ 18 O precip on the δ 18 O of O 2 produced by the terrestrial biosphere, we calculated monthly GPP-weighted mean δ 18 O precip . Throughout the entire seasonal cycle, GPP-weighted mean δ 18 O precip varies by ∼ 2.6 ‰, with the isotopically lightest and heaviest GPP-weighted precipitation occurring during boreal sum-mer and winter, respectively. GPP-weighted mean δ 18 O precip was found to be strongly and significantly correlated with TOE latitude over the seasonal cycle (r = 0.95, p < 0.0001), with a mean sensitivity of −0.066 ‰ degree −1 TOE latitude. Applying a last glacial ice mask (Argus et al., 2014) such that GPP was set to zero of glaciated land areas at 26 ka, we find the mean sensitivity of GPP-weighted mean δ 18 O precip remains similar (−0.072 ‰/ • latitude).

Discussion
Considering the observed seasonal-cycle relationship between the TOE and global terrestrial GPP-weighted mean δ 18 O precip , here we propose a mechanism by which meridional shifts in the TOE are unambiguously recorded in the sign of ε LAND changes. Although the seasonal cycle is an imperfect analog for abrupt climate change during the last glacial period, the analysis presented in Sect. 3.2 highlights the importance of the meridional distribution of oxygen production for the GPP-weighted mean δ 18 O of photosynthetic source water. Specifically, a northward shift of the TOE (e.g., during a D-O warming) lowers global GPPweighted mean δ 18 O precip , while a southward shift (e.g., during a HS) does the opposite. We suggest that the sign of this relationship is robust over long timescales although the magnitude likely varies. A simple experiment in which the modern terrestrial GPP climatology (with a 26 kyr ice mask) is perturbed by reducing Northern Hemisphere GPP and increasing Southern Hemisphere GPP each uniformly by 25 % suggests that several-degree-latitude shifts in the TOE may produce changes in terrestrial GPP-weighted mean δ 18 O precip that are comparable in magnitude (order 0.1 ‰) to ε LAND changes over HSs (Fig. S4). To further explore this notion, we consider our original definition of ε LAND (Eq. 1). In this definition, ε L encompasses all changes in global leaf water isotopic composition, which may be divided into separate source water and evapotranspiration components: δ 18 O S is the change from present in oxygen productionweighted source-water δ 18 O (having removed the influence of changes in seawater δ 18 O) and ε ET is the change from present in oxygen production-weighted evapotranspiration enrichment of leaf water δ 18 O. Many studies of Dole effect variability over the last glacial cycle (e.g., Bender et al., 1994;Severinghaus et al., 2009;Landais et al., 2010) have identified low-latitude precipitation as a key driver, based on coherence with orbital precession and atmospheric methane as well as high correlation with Chinese stalagmite δ 18 O. A recent quantitative, earth-system-model-based analysis of Dole effect variations over a simulated HS corroborates the hypothesized influence of low-latitude precipitation and finds that the increase in δ 18 O atm over a HS can be explained almost entirely by changes in terrestrial oxygen production-weighted δ 18 O precip (Reutenauer et al., 2015). Reutenauer et al. (2015) find that oxygen production-weighted relative humidity stays effectively constant between HS and background glacial conditions. Because relative humidity governs leaf water fractionation due to evapotranspiration (Craig and Gordon, 1965), this suggests negligible change in ε ET and implies (in terms of Eq. 6) that ε L ≈ δ 18 OS The latitude dependence of global GPP-weighted δ 18 O precip over the modern seasonal cycle implies that ε LAND should be sensitive to the position of the TOE over a HE. At present, however, it is impossible to determine the exact TOE sensitivity of ε LAND because it is modulated by changes in the underlying spatial distribution of δ 18 O precip . For instance, Chinese speleothem calcite records (Cheng et al., 2016;Wang et al., 2008Wang et al., , 2001 as well as isotope-enabled climate model simulations (Battisti et al., 2014;Liu et al., 2014) demonstrate that δ 18 O precip in southeastern China changes with EASM (East Asian summer monsoon) strength, such that local δ 18 O precip decreases during periods of an intense monsoon and increases with a weaker monsoon. Since this region is important for global oxygen production, it follows that a decrease in EASM intensity after a HE would act to increase ε LAND . However, concurrent, smaller-amplitude decreases in southern tropical rainfall δ 18 O precip during a HS (Kanner et al., 2012), particularly in the productive Amazon lowlands (Wang et al., 2017), may somewhat offset this increase in ε LAND Reutenauer et al. (2015) find that decreased southern tropical δ 18 O precip does not fully compensate for increased Northern Hemisphere δ 18 O precip , and conclude that spatial changes in δ 18 O precip over the terrestrial biosphere impart a net increase on δ 18 O atm over a simulated HS.
Clearly, to understand δ 18 OS one must consider both (1) changes in regional δ 18 O precip ( δ 18 O precip ), and (2) changes in the distribution of oxygen production (the TOE, to first order). We suggest that δ 18 OS is driven by the superposition of these two components, both of which likely operate in the same direction. For a northward migration of the thermal equator and tropical rainfall belts, rainfall over productive northern regions like Southeast Asia becomes more isotopically depleted (Battisti et al., 2014;Cheng et al., 2016;Liu et al., 2014;Wang et al., 2008Wang et al., , 2001. Assuming Northern Hemisphere δ 18 O precip outweighs Southern Hemisphere δ 18 O precip , a northward migration of the tropical rainfall belts imparts a negative signature on δ 18 O s . Similarly, a northward shift of the thermal equator and therefore a warmer and wetter Northern Hemisphere (and cooler, drier Southern Hemisphere) shifts the TOE northward as well. By increasing the relative contribution of the lowδ 18 O precip Northern tropics and mid-latitudes to global terrestrial O 2 , global δ 18 OS will also decrease. Each of these effects on (1) and (2) -and especially their combinationwould lead to a decrease in δ 18 OS and thus a decrease in ε LAND for an abrupt Northern Hemisphere warming (e.g., www.clim-past.net/13/1323/2017/ a D-O warming). In the opposite scenario (southward shift of thermal equator) during a HE, for example, the superposition of (1) and (2) would act to increase δ 18 OS because of a strong increase in Northern Hemisphere δ 18 O precip and a southward shift of the TOE. Therefore, we suggest that the sign of a ε LAND change faithfully reflects the direction of meridional shifts in the thermal equator and tropical rainfall belts.
While changes in respiratory fractionation, ε RL , are poorly known, two leading suggestions for the response of ε RL to large-scale hydroclimate change each reinforce the notion that the sign of ε LAND change is a robust indicator of meridional shifts in tropical rainfall. First, Reutenauer et al. (2015) find that changes in terrestrial respiration over a modeled HS are uncertain but small, one order of magnitude less important than GPP-weighted δ 18 O precip in explaining the modeled change in δ 18 O atm . If true, this is consistent with the notion that ε LAND is largely controlled by ε L , which mostly reflects δ 18 OS. Second, if instead global ε RL does appreciably vary with shifts in tropical rainfall, a candidate mechanism is the "monsoon rectifier effect", which proposes that changes in ε RL may oppose and partially counteract changes in δ 18 OS at high latitudes but amplify changes in δ 18 OS in monsoon regions (Luz and Barkan, 2011;Severinghaus et al., 2009). This notion hinges on the close coupling between local photosynthesis and respiration, along with the temperature dependence of soil respiration fractionation, which is stronger (i.e., increased discrimination of 18 O) in cold high-latitude soils (Angert et al., 2012). If, for example, northern high-latitude GPP were to increase (e.g., for a northward shift of the thermal equator during a D-O warming), so too would northern highlatitude soil respiration. In terms of Eq. (1), then, the lowering of ε L (due to contribution of the isotopically light highlatitude source water) would be somewhat offset by a more negative ε RL (stronger fractionation), leading to attenuation of any change in ε LAND . In monsoon regions, however, a stronger monsoon accompanying a northward shift of the thermal equator would result in both lower δ 18 O precip and weaker (less negative) soil respiration fractionation due to slow diffusion of oxygen in wet soils (Angert et al., 2001). In these regions, a lowering of ε L would be accompanied by less negative ε RL , both of which act to lower ε LAND . Therefore, in either case -whether ε RL is negligible or amplifies monsoon δ 18 O precip signals -one would still expect a net decrease (increase) in ε LAND for a northward (southward) shift of the thermal equator and tropical rain belts.
The TOE-latitude dependence of past changes in ε LAND makes the composite SD-WD ε LAND record a useful complement to other records of past tropical hydroclimate change. A recent analysis of the high-resolution WD CH 4 record suggests that abrupt peaks in CH 4 during Hudson Strait HS (HS 1, 2, 4, and 5) are imprints of HE (Rhodes et al., 2015). Because past changes in tropical wetland CH 4 emissions are thought to be a bimodal function of the latitude of the thermal equator (Rhodes et al., 2015), an abrupt Northern Hemisphere cooling and southward shift of the thermal equator and land extension of the Intertropical Convergence Zone (ITCZ) would intensify southern tropical precipitation and thus enhance CH 4 emissions. However, because Northern Hemisphere warming (e.g., D-O events) also increases atmospheric CH 4 emissions by stimulating northern wetland emissions, CH 4 measurements alone cannot definitively reveal the direction of meridional shifts in tropical rainfall.
ε LAND , which increases during and after CH 4 peaks representing HE (Fig. 4), implies a southward shift of the TOE and therefore provides strong evidence supporting the interpretation that these CH 4 peaks have a Southern Hemisphere source.
Future higher-resolution measurements of δ 18 O atm from WD alone may yield further insight into large-scale hydroclimate changes. A WD-only record would remove any uncertainty in ε LAND induced by the WD-SD synchronization process and allow for direct comparison with WD CH 4 . Centennial-scale WD CH 4 variability may indicate important tropical dynamics (Rhodes et al., 2017) that could be elucidated by a companion high-resolution ε LAND record. For instance, the centennial-scale oscillations in ε LAND following HE 1, which appear independently in both the SD and WD records, possibly hint that ε LAND may have spectral power in the centennial band.
Finally, we caution against over-interpretation of small changes in ε LAND as clear indications of shifts in TOE latitude. Although global GPP-weighted mean δ 18 O precip is a first-order control on ε LAND , small, under-constrained changes over time in respiratory fractionation, evapotranspiration fractionation, the marine Dole effect, the terrestrial fraction of global oxygen production, and the slope of the TOE-δ 18 O precip relationship preclude meaningful analysis of 0.01 ‰-order ε LAND variability.

Conclusions
We analyzed a composite record of δ 18 O atm from the Siple Dome and WAIS Divide ice cores to compare synchronous instantaneous changes in the terrestrial fractionation of δ 18 O atm to records of Chinese speleothem δ 18 O and WAIS Divide atmospheric methane. Based on an analysis of the modern seasonal cycles of terrestrial oxygen production and production-weighted δ 18 O precip , we propose a simple relationship between spatial shifts in terrestrial oxygenesis and δ 18 O atm . Specifically, we identify a strong negative correlation between the centroid latitude of terrestrial oxygenesis and GPP-weighted mean δ 18 O precip . This relationship suggests that positive δ 18 O atm (or ε LAND ) anomalies should accompany southward shifts of terrestrial oxygen production. Because productivity has strong hydroclimate controls, we propose that this mechanism may explain much of the shared variability between past changes in the Dole effect and proxies for low-latitude hydrological changes. Finally, positive excursions in ε LAND during HS 1, 2, 4, and 5 culminate in local ε LAND maxima which shortly follow abrupt increases in atmospheric CH 4 (of several tens of ppb) that occurred within these stadial periods (times of otherwise low CH 4 concentrations). These changes in ε LAND imply southward shifts of terrestrial rainfall and therefore strongly support that interpretation that the near-contemporaneous spikes in methane -proposed to be evidence of a rapid atmospheric teleconnection during HE 1, 2, 4, and 5 -had a Southern Hemisphere source.  Following Severinghaus et al. (2009), WD δ 18 O measurements were corrected for enrichment due to gravitational settling and fractionation during gas loss. Gravitational settling leads to an increase in heavy to light gases, scaling nearly linearly with their mass difference (Craig et al., 1988). Gravitational corrections were made to each measurement (δ 18 O, δO 2 /N 2 , δ 15 N, δAr/N 2 ) as follows: where δ meas is a raw measured gas ratio, δ gravcorr is its gravitationally corrected value, and m is the mass difference between heavy and light gases in amu (4 for δO 2 /N 2 , 12 for δAr/N 2 , and 2 for δ 18 O).
To determine δ 18 O atm , δ gravcorr was corrected empirically for fractionation due to gas loss. Slow loss of gases from the lattice to microcracks causes mass-independent fractionation which affects δO 2 /N 2 and δAr/N 2 , due to their different diameters, but does not affect isotopes (Bender et al., 1995). Mass-dependent fractionation, which affects elemental and isotopic ratios, occurs due to diffusive loss of gases through the microcracks (Kobashi et al., 2008). Following Severinghaus et al. (2009) we find a plane of best fit (Fig. S3) by regressing gravitationally corrected δ 18 O replicate pair differences ( δ 18 O gravcorr ) against gravitationally corrected δO 2 /N 2 and δAr/N 2 pair differences ( δO 2 /N 2,gravcorr and δAr/N 2,gravcorr ): The regression coefficients governing this plane of best fit (R 2 = 0.31), A and B, are nearly equal and opposite in sign (A = −0.00454; B = 0.00447) similar to the finding of Severinghaus et al. (2009) in SD samples, implying that fractionation of δ 18 O gravcorr is indeed due to mass-dependent fractionation process (proportional to δAr/N 2,gravcorr − δO 2 /N 2,gravcorr ), as the mass-independent fractionation process equally fractionates δAr/N 2 and δO 2 /N 2 . Finally, δ 18 O atm is calculated by removing this empirically determined mass-dependent fractionation and fixing the mean δ 18 O atm value over the last 1000 years to 0 ‰ by choosing C (below) to be 1.357: In order to reduce the influence of noise in calculating time derivatives of δ 18 O atm from SD, WD and the composite record, each record was fit to a Fourier series by weighted linear least squares. Following Severinghaus et al. (2009), Fourier amplitudes were assumed to be inversely proportional to their squared frequency (i.e., red spectrum). A Fourier series with frequencies ranging from 1/697 to 3.2 cycles kyr −1 with spacing of 1/697 cycles kyr −1 was used for curve fitting. δ 18 O atm data were weighted according to measurement uncertainties (pooled SDs of 0.0095 and 0.0085 ‰ for SD and WD, respectively). SD outliers with anomalous δ 15 N were downweighted as in Severinghaus et al. (2009). After determining Fourier amplitudes, continuous fitted curves and analytical derivatives were evaluated at 10-year resolution.

Appendix B: Details of synchronization method
The gas timescale synchronization method implemented in this study (ALT; Fudge et al., 2014) optimizes error-weighted tie-point agreement and smoothness of the annual layer thickness profile. Note that because this study only concerns gas ages, annual layer thicknesses are defined as vertical distances between two ice layers containing gas bubbles differing in age by 1 year. Here we present (1) the formally defined linear least squares problem and its uncertainty estimate, (2) methods for tie-point selection and uncertainty analysis, (3) an evaluation of the WD2014-synchronized SD gas age timescale, and (4) the results of a sensitivity test using an alternate initial SD timescale for tie-point selection.
The generalized ALT cost function (Eq. 3) is given here in the full form used in our analysis (bold-and-capitalized, boldand-lower-case and unbolded letters denote matrices, vectors and scalars, respectively): where (as in Eq. 3) G maps annual layer thicknesses, m, to a depth-gas age relationship evaluated at tie-point depths, d.
L is a second derivative operator, α is a trade-off parameter between smoothness and tie-point agreement, and W is a weighting matrix with inverse tie-point depth uncertainties (1/σ 2 ) along its diagonal. The time series of annual layer thicknesses that minimizes Eq. (B1),m, found at the zero crossing of ∂C ∂m : Equation B2 follows the common notation of Lundin et al. (2012) and Aster et al. (2013), such that G # is defined as We note that an alternate formulation of G # is possible by using generalized singular value decomposition (Aster et al., 2013;Lundin, 2012). We solved form via singular value decomposition and found that the maximum absolute disagreement of the resulting annual layer thickness profile with the one calculated using B2 and B3 was negligible (on the order of 10 −8 m).
The 1σ uncertainty ofm (for a fixed value of α) is initially estimated by taking the square root of the diagonal values on the covariance matrix, V m : However, uncertainty in α also contributes to uncertainty in m We account for this additional source of uncertainty by determining two candidate optimal values of α and evaluating resulting timescales using each value. These two candidate values of α minimize root mean square deviations between WD fitted-curve and SD discrete (1) δ 18 O atm values, or (2) CH 4 values, with SD gas ages shifted in each case to a new estimate of the WD2014-synchronized timescale. We evaluatem separately for each value of α and take the absolute difference between these two time series ofm as an estimate of the 1σ uncertainty in ourm due to uncertainty in α.
The final estimate of 1σ uncertainty inm is then given by the quadrature sum of the covariance-matrix (fixed α) estimate (Eq. B3) and the estimated uncertainty due to α described above. Tie points were estimated by matching abrupt transitions in dδ 18 O atm dt −1 and CH 4 in the WD and SD records. To find dδ 18 O atm /dt-based tie points, continuous time derivatives of the fitted δ 18 O atm curves were first evaluated as described in Appendix A. All local minima and maxima in dδ 18 O atm dt −1 were then identified and 36 abrupt transitions identifiable in each record were selected (each one between an identified local minimum and maximum). The midpoint time of a transition was found by determining when dδ 18 O atm dt −1 was halfway between the bounding extrema. The times at which dδ 18 O atm dt −1 was 25 and 75 % between bounding extrema were taken as ±1σ uncertainty estimates. Figure B1 shows abrupt dδ 18 O atm dt −1 transition midpoints, bounding extrema, and 25-75 % markers identified in the SD (original timescale) and WD records. For each dδ 18 O atm dt −1 midpoint, corresponding SD depths and WD gas ages were used as tie points for the synchronization method.
A similar approach was employed to determine CH 4 -based tie points using CH 4 records at SD and WD. Since WD CH 4 records (Rhodes et al., 2015;Sowers, 2012) are far higher resolution than SD CH 4 records (Ahn et al., 2012Brook et al., 2005;Ahn and Yang, unpublished), discrete SD CH 4 values were selected during abrupt transitions for which the magnitude of SD and WD CH 4 values were in good agreement. Twenty of these discrete midtransition values were identified and used as tie points. To find a WD gas age matching the SD depth of a discrete CH 4 measurement, the gas age of the nearest value in the 2-year resolution WD record was determined. The 1σ uncertainty range for CH 4 tie points was estimated from the range of gas ages corresponding to all nearby WD CH 4 values that agreed with the selected SD CH 4 value to within ±15 ppb. Figure B2 illustrates this CH 4 tie-point value matching and uncertainty process.    (Buizert et al., 2015b). Root mean squared deviations between WD fitted curve and SD discrete measurements are 0.029 ‰ on the original timescale and 0.022 ‰ on the new timescale. Figure B4. CH 4 in SD (Ahn et al., 2012Brook et al., 2005;Ahn and Yang, unpublished) and WD (Rhodes et al., 2015;Sowers, 2012): discrete SD measurements on original gas timescale (blue markers; Brook et al., 2005) and on new gas timescale (green markers, this study) are superimposed on WD interpolated curve (blue line) on WD2014 gas timescale (Buizert et al., 2015b). Root mean squared deviations between WD fitted curve and SD discrete measurements are 33.7 ppb on the original timescale and 18.3 ppb on the new timescale.
We also tested the sensitivity of the resulting timescale to the chosen CH 4 matching value range (±15 ppb) by varying it from ±5 to ±25 ppb to construct alternate timescales. Fig-ure S5 shows deviations of these alternate timescales from the final timescale used in this study. The results of this test suggest that the choice of ±15 ppb for the final timescale is reasonable, because alternate timescales constructed with ranges of ±10 and ±20 ppb agree with the ±15 ppb-based timescale within ±1σ uncertainty at all times, with a maximum absolute deviation of ∼ 58 years.
Together, dδ 18 O atm dt −1 and CH 4 tie points were used to constrain equation B2 via the weighting matrix W to produce a best estimate of SD annual layer thicknesses at 10year resolution from 50.12 to 0 ka BP (before 1950). Integrating this annual layer thickness time series yields a depthgas age relationship. Figure 2 shows the annual layer thickness estimate, its uncertainty, and a comparison between our new ALT-based and the original Brook et al. (2005) SD timescales. To assess the validity of the new SD gas timescale estimate, we calculated root mean squared deviations (RMSDs) between WD (fitted curve or interpolated) and SD (discrete) δ 18 O atm and CH 4 on the original and new timescales. Although a non-ideal evaluation metric for relatively flat portions of these records, RMSDs penalize timescale errors during transitions. The substantial RMSD reductions between the original and new timescale add confidence to the synchronization method. CH 4 RMSDs decrease from 33.7 to 18.3 ppb on the original and new timescales, respectively. δ 18 O atm RMSDs decrease from 0.029 to 0.022 ‰ on the original and new timescales, respectively. Figures S5 (in the Supplement) and B3 show fitted and interpolated WD δ 18 O atm and CH 4 curves, respectively, alongside discrete SD δ 18 O atm and CH 4 measurements on the original and new timescales.
To investigate the sensitivity of our SD depth-gas age to the initial SD timescale, we repeated the SD tie-point selection procedure and least squares minimization using an alternate SD timescale. This sensitivity test was motivated by concern that error in the original timescale would lead to error in the magnitude of the time derivative of δ 18 O atm . The tie-point method implemented in this study (aligning the midpoints of dδ 18 O atm dt −1 transitions rather than value-matching) partially alleviates this concern. Nonetheless, for completeness, we carried out a sensitivity test using an alternate initial SD timescale for curve fitting, differentiation, and tie-point determination to constrain synchronization method (Eq. 3). Because the Brook et al. (2005) timescale was constructed by tying SD gas measurements to the Meese/Sowers GISP2 (Greenland) ice core timescale (Meese, 1999), we added the difference in gas ages between the Meese/Sowers and more recent GICC05 timescales (Rasmussen et al., 2014) to the Brook et al. (2005) SD timescale. This resulting alternate timescale differs from the original Brook et al. (2005) timescale by up to ∼ 2 kyr in the interval from ∼ 50 ka to present. This alternate timescale is notable older before ∼ 45 ka, as shown clearly by comparison of SD dδ 18 O atm dt −1 records on the original and alternate GISP2-tied timescales (Fig. S6). The close similarities in the magnitude of SD dδ 18 O atm dt −1 on the original and alternate timescales strongly indicates that the small dating errors are not significant enough to substantially affect the magnitude of the δ 18 O atm derivative. Each of the 36 transitions associated with dδ 18 O atm dt −1 tie points was identifiable on this alternate timescale (Fig. S7). After determining the depths corresponding to these midpoints on the alternate timescale (in order to determine tie-point array d for equation B2), midpoint depths were compared between the original and alternate timescales and found to agree within the 25-75 % transition value derived 2σ uncertainty range for 35 of 36 tie points. The same CH 4 tie-point method was implemented and the resulting set of 46 depth-gas age tie points was used to estimate an alternate annual layer thickness profile via the ALT method (Eq. 3). This alternate estimate was found to agree the original estimate within 2σ uncertainty throughout the entire time range considered (Fig. S8). We therefore conclude that, given these relatively small dating errors of below ∼ 2 kyr, the δ 18 O atm derivative-based tie-point determination method is justified.