# Quantifying molecular oxygen isotope variations during a Heinrich stadial

C. Reutenauer, A. Landais, T Blunier, C Bréant, M Kageyama, M.-N Woillez, C Risi, V Mariotti, P Braconnot

► To cite this version:

C. Reutenauer, A. Landais, T Blunier, C Bréant, M Kageyama, et al., Quantifying molecular oxygen isotope variations during a Heinrich stadial. Climate of the Past, 2015, 11 (11), pp.1527-1551. 10.5194/cp-11-1527-2015. hal-01259465

HAL Id: hal-01259465

https://hal.sorbonne-universite.fr/hal-01259465

Submitted on 20 Jan 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Quantifying molecular oxygen isotope variations during a Heinrich stadial

C. Reutenauer<fnref n="1" />, A. Landais<fnref n="2" />, T. Blunier<fnref n="1" />, C. Bréant<fnref n="2" />, M. Kageyama<fnref n="2" />, M.-N. Woillez<fnref n="2" />,ª, C. Risi<fnref n="3" />, V. Mariotti<fnref n="2" />, and P. Braconnot<fnref n="2" />

*Correspondence to:* C. Reutenauer (creuten@nbi.ku.dk)

Received: 22 May 2015 – Published in Clim. Past Discuss.: 23 June 2015

Revised: 1 November 2015 – Accepted: 5 November 2015 – Published: 19 November 2015

**Abstract.** δ¹⁸O of atmospheric oxygen (δ¹⁸O_atm) undergoes millennial-scale variations during the last glacial period, and systematically increases during Heinrich stadials (HSs). Changes in δ¹⁸O_atm combine variations in biospheric and water cycle processes. The identification of the main driver of the millennial variability in δ¹⁸O_atm is thus not straightforward. Here, we quantify the response of δ¹⁸O_atm to such millennial events using a freshwater hosing simulation performed under glacial boundary conditions. Our global approach takes into account the latest estimates of isotope fractionation factor for respiratory and photosynthetic processes and make use of atmospheric water isotope and vegetation changes. Our modeling approach allows to reproduce the main observed features of a HS in terms of climatic conditions, vegetation distribution and δ¹⁸O of precipitation. We use it to decipher the relative importance of the different processes behind the observed changes in δ¹⁸O_atm. The results highlight the dominant role of hydrology on δ¹⁸O_atm and confirm that δ¹⁸O_atm can be seen as a global integrator of hydrological changes over vegetated areas.

## 1 Introduction

Oxygen is one of the most abundant species in atmospheric air. As oxygen is produced by photosynthesis and consumed by respiration, a record of oxygen concentration in the past should help us to constrain these two major biospheric fluxes on Earth and potentially provide information on their link with the carbon cycle.

Changes in the O₂/N₂ ratio can be measured in air trapped in ice cores back to 800 kyr (Bender, 2002; Kawamura et al., 2007; Landais et al., 2012; Lemieux-Dudon et al., 2015). Unfortunately the O₂/N₂ ratio in ice cores does not provide a direct information on the true atmospheric variations because it is affected by permeation through the ice lattice during bubble formation at pore close-off, roughly 100 m below the ice sheet surface, and by gas loss during ice core storage. These effects have less impact on the isotopic composition of oxygen. These isotopic compositions have thus been explored as possible constraints on biospheric productivity (Luz et al., 1999).

When dealing with isotopes, it is standard to use the isotope ratio, *R*, defined as the fraction of the abundance of the rare isotope over the dominant one in a substance. Since changes in isotope ratios through natural processes are very small, they are expressed in relation to a standard (recent air and Vienna Standard Mean Ocean Water (VSMOW) being

<footnote n="1">¹Centre for Ice and Climate, Niels Bohr Institute, University of Copenhagen, Copenhagen, Denmark</footnote>

<footnote n="2">²Institut Pierre Simon Laplace (IPSL), Laboratoire des Sciences de Climat et de l'Environnement (LSCE), UMR8212 (CEA-CNRS-UVSQ), CE Saclay, Orme des Merisiers, Bat. 701, 91191 Gif-sur-Yvette Cedex, France</footnote>

<footnote n="3">³Laboratoire de Météorologie Dynamique UMR8539, IPSL/CNRS/UPMC, 4, place Jussieu, 75252 Paris Cedex 05, France</footnote>

<footnote>ªnow at: IFPEN, 14 avenue de Bois Préau, 92852 Rueil-Malmaison Cedex, France</footnote>
used for O₂ and H₂O, respectively) using the δ notation,

$$ \delta^{18}\text{O} = \frac{^{18}\text{R}_{\text{sample}}}{^{18}\text{R}_{\text{standard}}} - 1. \qquad (1) $$

δ¹⁸O and δ¹⁷O of atmospheric oxygen have been measured for the period of the past 800 kyr with a mean resolution of about 1500 years (e.g., Landais et al., 2010, and references therein; Blunier et al., 2012). As shown by Bender et al. (1994), δ¹⁸O of atmospheric oxygen, noted δ¹⁸O_atm hereafter, cannot easily be related to biospheric productivity through photosynthesis and respiration fluxes. δ¹⁸O_atm variations actually reflect for a large part the isotopic composition of the meteoric water. The latter is transmitted to the plant through its roots and stems to the leaves, where photosynthesis produces oxygen with an isotopic composition close to the isotopic composition in leaf water. Respiration modifies the isotopic composition of atmospheric oxygen in a complex way. While the processes consuming oxygen enrich atmospheric oxygen through a preferential consumption of the lightest molecules, individual biological pathways are associated with a wide range of oxygen fractionations (Helman et al., 2005).

Based on fractionation factors available at the time, Bender et al. (1994) established that the relative proportion of oceanic vs. terrestrial biospheric productivities together with the difference in isotope fractionation over land and ocean were driving the δ¹⁸O_atm budget. Several studies have built on this idea and interpret the δ¹⁸O_atm variations mainly as variations in the oceanic vs. terrestrial biospheric productivities (Hoffmann et al., 2004; Ciais et al., 2012). However, recent measurements have revealed that, overall, fractionation associated with oceanic productivity is very similar to its terrestrial counterpart (Hendricks et al., 2004; Eisenstadt et al., 2010), questioning the interpretation of δ¹⁸O_atm as an indicator of the relative proportion of oceanic vs. terrestrial biosphere productivity (Luz and Barkan, 2011b).

Despite the complex interpretation of δ¹⁸O_atm, several robust features have already been observed that highlight the potential of these measurements. At the orbital scale, δ¹⁸O_atm is showing clear variations at a 23 kyr periodicity (Dreyfus et al., 2007). This strong link with precession is probably related to the variations in the hydrological cycle at low latitudes (Bender et al., 1994). Indeed, variations related to the monsoon regime strongly imprint the isotopic composition of meteoric water as observed in speleothem records (e.g., Wang et al., 2008). They are easily transmitted to the isotopic composition of atmospheric oxygen because the major part of the biospheric productivity, and hence photosynthesis, is occurring in the tropics and subtropics.

At the millennial scale, it has recently been shown that δ¹⁸O_atm is responding to the abrupt climate changes of the last glacial period (Landais et al., 2007a; Severinghaus et al., 2009). Millennial-scale climate variability is perhaps best known from the Greenland ice cores, where it is manifested

in the stable water isotopes of ice. During the last glacial period, these cores show 25 Dansgaard-Oeschger (DO) events (NGRIP members, 2004). A DO event typically exhibits a sawtooth pattern: (i) a cold phase (Greenland stadial, noted GS hereafter) lasting from centuries to millennia, followed by a warm phase (Greenland interstadial, GI) starting with (ii) a rapid transition (a few decades) with an amplitude of up to 16 ± 2.5°C (Landais et al., 2004; Huber et al., 2006; Kindler et al., 2014), and ending with (iii) a gradual cooling before an abrupt decrease towards cold, stadial values.

During the last decade, mechanisms of glacial abrupt events have been investigated using coupled ocean-atmosphere models of varying complexity (e.g., Kageyama et al., 2010, 2013; Stouffer et al., 2006). Recent hypotheses often invoke internal variability (Kleppin et al., 2015; Dokken et al., 2013), involving sea ice-atmosphere interactions (e.g., Li et al., 2005, 2010), through ice-albedo feedback and the impact of sea ice cover on regional temperatures by preventing heat exchange between the ocean and atmosphere. There remains robust evidence from multiple lines of paleoceanographic information and modeling that millennial-scale variability is linked to changes in the Atlantic meridional overturning circulation (AMOC) intensity (e.g., Mc Manus et al., 1998), potentially initiated by large freshwater input in the North Atlantic (e.g., Broecker et al., 1990). The presence of ice rafted debris (IRD; Ruddiman, 1977; Heinrich, 1988) in marine sediments from the North Atlantic region during the largest GS document episodes of massive iceberg discharge in the North Atlantic (Heinrich events) mainly from the Laurentide (H2, H4, H5) and Fennoscandian (H3, H6) ice sheets (Grousset et al., 1993; Guillevic et al., 2014 and references therein). Even though IRD is present in each GS (Elliot et al., 2002), not all GSs contain a Heinrich event. Heinrich stadials (noted HSs hereafter) are GSs associated with a Heinrich event (Barker et al., 2009; Sanchez Goni and Harrison, 2010).

Several aspects of the observed patterns during DO events can be captured through the response of the Earth system to imposed freshwater perturbations in the North Atlantic (Liu et al., 2009; Otto-Bliesner and Brady, 2010; Kageyama et al., 2010; Roche et al., 2010), mimicking Heinrich events. Depending on the background state of the climate (glacial or interglacial, orbital context) and the AMOC, as well as on the magnitude of the freshwater forcing, these models produce a complete shutdown of the AMOC (HS-like state) or a reduction of the strength of the AMOC (GS-like state; e.g., Menviel et al., 2014). The injection of freshwater produces in all models a significant cooling of the North Atlantic region. The amplitude of the associated temperature change is probably affected by the simulated change in sea ice extent and feedbacks between sea ice and temperature that vary in the different models (Kageyama et al., 2013). These housing experiments also produce an interhemispheric see-saw temperature pattern, associated with a southward shift of the Intertropical Convergence Zone (ITCZ) (e.g., Dahl et al.,
2005; Broccoli et al., 2006; Krebs and Timmermann, 2007; Swingedouw et al., 2009; Cvijanovic and Chiang, 2013). Abrupt climate variation associated with the Greenland signal is found down to low latitudes in numerous terrestrial and marine archives (e.g., Clement and Peterson, 2008). Its climatic impact is recorded in large parts of the North Atlantic region, both in marine cores (e.g., Bond et al., 1993; Broecker, 2000) and in speleothems (Fleitmann et al., 2009). Concomitant methane excursions and variations in the isotopic composition of the calcite of speleothems in eastern Asia (e.g., Wang et al., 2001; Cheng et al., 2012) strongly support the fact that these DO events are associated with major reorganization of the tropical water cycle and hence monsoon intensity through a shift in the ITCZ and its terrestrial equivalent, the tropical rain belt (Chappellaz et al., 2013; Wang et al., 2008; Pausata et al., 2011).

For this period of millennial-scale variability, high-resolution measurements of δ¹⁸O_atm have been obtained in Greenland and Antarctic ice cores (e.g., Guillevic et al., 2014; Landais et al., 2007a, 2010; Severinghaus et al., 2009). In Fig. 1 we present a synthesis of δ¹⁸O_atm evolution from the Siple Dome ice core over HSs displayed on Greenland Ice Core Chronology 2005 (GICC05) timescale, using definitions of Rasmussen et al. (2013) of the onset of GS. The δ¹⁸O_atm records show a systematic increase in a few thousand years following the onset of a HS (Fig. 2) by around 0.13‰, from +0.08‰ (HS1) to +0.18‰ (HS5). The difference in the slope inflection at the onset of HS4 and HS5, more pronounced that for HS1, HS2 and HS3, may be due to the long-term trend observed in δ¹⁸O_atm. Indeed, from 35 to 15 kyr, δ¹⁸O_atm exhibits a constant increase, consistent with the build-up of polar ice sheet, and hence enrichment of ocean water in ¹⁸O, but remains relatively stable over MIS3 (Fig. 1).

Because of its global character, δ¹⁸O_atm should provide added value compared to the different local records of hydrological cycle variations in different continental and marine archives. However, until now, no quantitative, robust interpretation of past variations in δ¹⁸O_atm has been established, which limits the use of δ¹⁸O_atm as a quantitative indicator for past biospheric production or variations in the hydrological cycle. The aim of this modeling study is thus to provide a quantitative interpretation for the systematic increase in δ¹⁸O_atm over HSs. To reach this objective, we propose a global approach incorporating outputs from a general circulation model implemented with water isotopes and focus on the millennial variability in the last glacial period. We follow a modeling approach already developed by Hoffmann et al. (2004). We combine climatic parameters (temperature and humidity), isotopic composition of meteoric water, vegetation distribution and productivity simulated by different models with monthly mean temporal resolution.

In the following section, we describe the general method used to simulate a global δ¹⁸O_atm signal. Section 3.1 is dedicated to model validation and Section 3.2 proposes to quan-

tify the different contributions (hydrology, vegetation, climatic conditions) to the δ¹⁸O_atm signal over a HS equivalent.

## 2 Method

According to Landais et al. (2007a, 2010) and Severinghaus et al. (2009), the millennial variations in δ¹⁸O_atm during the last glacial period are driven by shifts in the tropical rain belt modifying the relative humidity distribution and the isotopic composition of meteoric water consumed by terrestrial biosphere. The isotopic content of atmospheric oxygen is controlled by numerous processes, so we must consider (i) the worldwide meteoric water isotopic composition, from which ground water is derived; (ii) the worldwide temperature and humidity, from which evaporative enrichment of leaf water δ¹⁸O is calculated; (iii) the worldwide vegetation cover and gross primary productivity, defining the photosynthetically and respiratory active areas that contribute to δ¹⁸O_atm; and (iv) respiratory processes.

### 2.1 Oxygen isotopes mass balance model

Oxygen is exchanged with the terrestrial and marine biospheres as well as with the stratosphere. Assuming a steady state, δ¹⁸O_atm can thus be expressed as follows:

$$ \delta^{18}\text{O}_{\text{atm}} = \frac{(F_{\text{terr}} \cdot \delta^{18}\text{O}_{\text{terr}} + F_{\text{mar}} \cdot \delta^{18}\text{O}_{\text{mar}})}{(F_{\text{terr}} + F_{\text{mar}})} - 18\epsilon_{\text{strat}}, \quad (2) $$

where ¹⁸ε_strat represents the stratospheric isotope fractionation caused by photochemical reaction in the stratosphere involving O₂, O₃ and CO₂. F_terr and F.mar denote O₂ fluxes of gross terrestrial and oceanic productivity, respectively. δ¹⁸O_terr and δ¹⁸O.mar are the isotopic composition arising from the terrestrial and oceanic realms, respectively. ¹⁸ε_strat is a small term, 0.4‰ compared to ~23.8‰ for δ¹⁸O_atm with reference to V-SMOW (Luz et al., 2014) and is not assumed to change significantly over a HS because CO₂ level remains relatively stable. We assume a constant CO₂ level between the LGM and HS in our study. Ahn and Brook's (2014) study shows that variations over HSs are small (increase of less than 20 ppm). Effect of isotopic exchange between CO₂ and O₂ in the stratosphere on δ¹⁸O_atm is expected to be proportional to CO₂ mixing ratio. Following the calculation of Bender et al. (1994), which estimates a δ¹⁸O_atm depletion of 0.4‰ for a CO₂ concentration of 353 ppm, we can estimate that a 20 ppm increase between the LGM and HS can modify δ¹⁸O_atm by -0.023‰. The sign of this change is actually opposite to the sign of the observed δ¹⁸O_atm signal. We focus mainly on the millennial-scale variations in the terrestrial contribution to δ¹⁸O_atm signal, i.e., $F_{\text{terr}} \cdot \delta^{18}\text{O}_{\text{terr}}/(F_{\text{terr}} + F_{\text{mar}})$.

We do not consider the marine influence, in this first approach, for the following reasons. First, our aim is to test the hypothesis of Landais et al. (2007b) and Severinghaus et al.
<figure><img src="image_3.png" /></figure>

**Figure 1.** Greenland stadials and Heinrich stadials during the last glacial period. Red numbers indicate GSs and black labels HSs. (a) Black line: NGRIP δ¹⁸O, ‰, on the GICC05 timescale back to 60 ka b2k. (b) Red line: Siple Dome atmospheric δ¹⁸O (Severinghaus et al., 2009) on GICC05 timescale. Red dots denote intervals in the core where deep air convection or cracked firn layer may have thermally fractionated these samples, as shown by anomalous δ¹⁵N in Figs. S1 and S8 of Severinghaus et al. (2009). The transfer of Siple Dome atmospheric δ¹⁸O on the GICC05 chronology is achieved by using the Siple Dome gas age scale compatible with the GICC05 chronology (Guillevic, 2013), based on match points between Siple Dome methane variations (Brook et al., 2005; Ahn et al., 2012; depth point) and NEEM methane variations (Chappellaz et al., 2013; GICC05 gas age point). A linear interpolation is then performed between match points to calculate the Siple Dome gas age. Green line: NEEM atmospheric δ¹⁸O, ‰ (±0.03 ‰; Guillevic et al., 2014). Colored areas: GSs. Grey: GS with no major Heinrich event. Orange: HS1, HS2, HS3, HS4, HS5 and end of HS6 – GSs with a major Heinrich event. Black error bars indicate HS onset and end uncertainty (2σ), based on the Rasmussen et al. (2013) maximum counting error (MCE). Top (right to left): black and white horizontal bar indicate Marine Isotope Stage 1 to 4.

<figure><img src="image_4.png" /></figure>

**Figure 2.** Evolution of Siple Dome atmospheric oxygen δ¹⁸O (Severinghaus et al., 2009) during Heinrich stadials on the GICC05 timescale.

(2009) that δ¹⁸O_atm millennial-scale variations are largely driven by changes in the low-latitude hydrological cycle through changes in the δ¹⁸O of precipitation. Second, Hendricks et al. (2004) and Luz and Barkan (2011a) have shown that the difference between δ¹⁸O_terr and δ¹⁸O-mar is not significant. Finally, the spatial and temporal variations in water δ¹⁸O and respiration pathways in the ocean are expected to be relatively small compared to the variations on land, which renders their integration for δ¹⁸O_atm modeling less crucial, as illustrated in the study of Hoffmann et al. (2004).

## 2.2 Calculation of δ¹⁸O_terr

The major source of atmospheric oxygen from the terrestrial biosphere is the oxygen produced during photosynthesis. The fractionation associated with photosynthesis is small (Guy et al., 1993; Eisenstadt et al., 2010). The oxygen produced by this process thus has almost the same isotopic composition as the leaf water. Consumption of oxygen is also associated with biosphere productivity through different pathways (dark respiration, photorespiration, Mehler reaction). δ¹⁸O_terr thus results from isotope fractionation associated with photosynthesis.
thesis and oxygen uptake. Following Bender et al. (1994), Blunier et al. (2002), Hoffmann et al. (2004), and Landais et al. (2007b), we assume a steady state, where photosynthesis equals respiration. δ¹⁸Oₜerl calculates as

$$ \delta^{18}\text{O}_{\text{terr}} = \frac{(\delta^{18}\text{O}_{\text{lw}} + 1)}{^{18}\alpha_{\text{resp}}} - 1, \quad (3) $$

where δ¹⁸Oₜlw is the global production-weighted average isotopic composition of leaf water and ¹⁸αₜresp is the global apparent respiratory isotope fractionation factor associated with global oxygen uptake, i.e., oxygen consumption weighted average of fractionation factors associated with specific respiratory pathways.

## 2.3 Photosynthetic oxygen

As classically done, we estimate the value of δ¹⁸O of leaf water, hereafter δ¹⁸Oₜlw, based on the Craig and Gordon (1965) equation (C&G) of evaporation applied to leaf transpiration (Dongmann, 1974; Flanagan et al., 1991b). It is thus calculated in the following way:

$$ \delta^{18}\text{O}_{\text{lw}} = h \cdot (\delta^{18}\text{O}_{\text{vap}} + ^{18}\epsilon_{\text{eq}}) + (1 - h) \cdot (\delta^{18}\text{O}_{\text{gw}} + ^{18}\epsilon_{\text{eq}} + ^{18}\epsilon_{\text{kin}}), \quad (4) $$

where h is the relative humidity at the site of photosynthesis, ¹⁸εₜeq is the temperature-dependent liquid-vapor equilibrium isotope effect (Majoube, 1971), ¹⁸εₜkin is the kinetic isotope effect occurring when humidity is below saturation, δ¹⁸Oₜgw is the isotopic composition of soil water and δ¹⁸Oₜvap is the water vapor δ¹⁸O near the surface. ¹⁸εₜkin is deduced from the ratio of the diffusion coefficient associated with H₂¹⁶O (D) and H₂¹⁸O (D*). Several values for the ratio D/D* can be found in the literature (Merlivat, 1978; Cappa et al., 2003; Luz et al., 2009), varying from 1.028 to 1.032. For leaf water evaporation, many studies have reported lower enrichment in δ¹⁸Oₜlw than that predicted by Eq. (4) with ¹⁸εₜkin = D/D* - 1 as classically assumed (e.g., Allison et al., 1985; Bariac et al., 1989; Walker et al., 1989; Walker and Brunel, 1990; Yakir et al., 1990; Flanagan et al., 1991b, a, 1993, 1994). Farquhar et al. (1989) suggested that ¹⁸εₜkin depends on the importance of either stomatal or boundary layer resistances. In moist conditions, stomata resistance is low and boundary layer resistance high, leading ¹⁸εₜkin to values as low as 19‰ when using the Merlivat (1978) value for D/D*. In this study, we have imposed a mean value for ¹⁸εₜkin of 20‰ because higher values led to too high a global value for δ¹⁸Oₜatm.

The calculation of δ¹⁸Oₜlw using Eq. (4) requires spatial and temporal variations in temperature and relative humidity as well as the variations in the isotopic composition of water vapor and meteoric water, from which δ¹⁸Oₜgw will be deduced (Sect. 2.3.2). These variables are obtained from outputs of modeling experiments.

### 2.3.1 Simulated climatic variations over an abrupt cooling

Temperature and relative humidity variations over a HS are inferred from simulations with the atmosphere-ocean general circulation model (AOGCM) IPSL_CM4 (Marti et al., 2010) with a horizontal resolution for the atmosphere grid of 3.75° × 2.5° (latitude × longitude) and with a 19 vertical layer atmosphere. To model the δ¹⁸Oₜatm variations over a HS, we have used a glacial simulation perturbed by a freshwater hosing experiment. We will compare in the following the outputs of two simulations: one for the Last Glacial Maximum (LGMCtrl) and one for the Heinrich stadial (HS_exp).

The LGMCtrl boundary conditions are as follows (see Kageyama et al., 2009, for a detailed presentation of the climate setup): orbital parameters for 21 ky BP, CO₂, CH₄ and N₂O levels set to 185 ppm, 350 and 200 ppb, respectively (Monnin et al., 2001; Dällenbach et al., 2000; Flückiger et al., 1999), ICE-5G ice sheet reconstruction and land-sea mask (Peltier, 2004).

The first experiment is an equilibrated glacial run (LGMCtrl) used as a reference run (see LGMb in Kageyama et al., 2009). The second experiment (HS_exp) is a water hosing experiment, where an additional freshwater flux of 0.1 Sv (1 Sv = 10⁶ m³s⁻¹) is imposed instantaneously in the Atlantic north of 40°N and the Arctic (see LGMc in Kageyama et al., 2009) from year 150 for 400 years. The input of freshwater in HS_exp, mimicking a Heinrich event, leads to an AMOC collapse in 250 years (see Fig. 1 in Kageyama et al., 2009). We selected this HS_exp experiment since the most efficient way to simulate the climate state during an Heinrich event with a model is to add freshwater in the high latitudes of the Atlantic ocean, even though results from recent studies (Marcott et al., 2011; Guillevic et al., 2014; Rhodes et al., 2015; Alvarez-Solas et al., 2013) suggest that this does not satisfactorily explain the observed sequences of events (freshwater discharges from ice sheets might not be the initial trigger of Heinrich events) and especially the decoupling between Greenland and low latitudes. We therefore center our study on the mean state of two contrasted periods rather than investigate the dynamics of the transition from the LGM to HS. It also explains why we only focus on HSs and not DO events.

The climate response to the AMOC collapse in the HS_exp is of global extent and qualitatively agrees with paleoclimate reconstructions for the North Atlantic cooling, southwards Atlantic ITCZ migration and weakening of Indian and African monsoons (Kageyama et al., 2009, 2013). However, the model does not simulate an Antarctic warming or weakened East Asian monsoon (Kageyama et al., 2009). In the following we used the monthly averaged spatial fields of temperature (Tm) and humidity calculated on the first layer of the atmosphere grid. From these data we followed the Lloyd and Farquhar (1994) approach to link leaf temperature during photosynthesis, Tp, to Tm, through
<figure><img src="image_5.png" /></figure>

**Figure 3.** Amount-weighted annual-mean δ¹⁸O of precipitation for (a) LGMCtrl and (b) HS_exp experiments obtained with LMDZ-iso. Note that the anomaly can be seen in Fig. 4a. (c) HS_exp–LGMCtrl annual-mean anomaly of groundwater δ¹⁸O. Groundwater δ¹⁸O represents the isotopic value of the substrate water for photosynthesis (see text for details). (d) HS_exp–LGMCtrl annual-mean anomaly of rainfall amount.

$T_p = 1.05 \cdot (T_m + 2.5)$. An additional 2.5°C is added to $T_m$ to account for the daytime increase in air temperature at the time of photosynthesis. The 5% increase allows for net canopy to air heat fluxes (Farquhar and Lloyd, 1993). The relative humidity from the first layer is not modified (Farquhar et al., 2007).

### 2.3.2 Modeling of δ¹⁸O of meteoric water and groundwater

For our estimate, we also need the distribution of the oxygen isotopic composition of meteoric water. We extract it from the isotopic version of the atmospheric general circulation model developed at the Laboratoire de Météorologie Dynamique (LMDZ4; Risi et al., 2010). LMDZ is the atmospheric component of the IPSL-CM4 model used above. The physical package is described in detail by Hourdin et al. (2006). It includes in particular the Emanuel convective parameterization (Emanuel, 1991; Grandpeix et al., 2004) coupled to the Bony and Emanuel (2001) cloud scheme. Each grid cell is divided into four subsurfaces: ocean, land, ice sheet and sea ice (Risi et al., 2010). The monthly sea surface temperature and sea ice fields obtained from the two aforementioned experiments at equilibrium (LGMCtrl and HS_exp) have been used as surface boundary conditions for the isotopic simulations. Monthly mean outputs of the IPSL-CM4 are imposed to the LMDZ4 model, so there is no coupling between ocean and atmosphere, and nonlinear sub-

monthly scale processes are thus not taken into account here. This choice should not alter our results. Indeed, LeGrande and Schmidt (2008) analyze changes in water isotopes following the 8.2 ka event's meltwater pulse (meltwater δ¹⁸O: -30‰) in a fully coupled AOGCM (Goddard Institute for Space Studies ModelE-R) and show that the effect of the freshwater impulse on δ¹⁸O of precipitation, noted δ¹⁸Oₚ hereafter, can be neglected because the signal is very short-lived, only a few decades, before the climatic component dominates. Furthermore, the main changes are constrained to the northern North Atlantic and its surrounding regions (Fig. 6 of LeGrande and Schmidt, 2008). Those regions only have a limited contribution to the δ¹⁸Oₐtmos signal, as most of the terrestrial photosynthesis occurs in the tropics. Figure 7 of LeGrande and Schmidt (2008) shows that two decades after a meltwater pulse, the ensemble mean (five simulations) anomaly of δ¹⁸Oₚ calculates to -0.01‰, confirming the small impact of δ¹⁸O depleted meltwater.

Figure 3 shows the mean annual δ¹⁸Oₚ, simulated for the LGMCtrl and the HS_exp experiments. The δ¹⁸Oₚ distribution for the LGMCtrl experiment has already been confronted to observations of water isotopes in vapor and precipitation (Risi et al., 2010) and seasonal patterns are well captured. For validation purposes, δ¹⁸Oₚ changes from the LGM to HS are compared with changes in (i) calcite δ¹⁸O in speleothems, noted δ¹⁸Oₑ hereafter, and in (ii) ice δ¹⁸O
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) Model-data comparison of δ¹⁸O precipitation anomaly during HS compared to the LGM. Data represent speleothem's calcite δ¹⁸O from various locations (see Table 3 for details). (b) Comparison of reconstructed HS precipitation anomaly Δδ¹⁸O from selected proxies shown in (a) and simulated Δδ¹⁸Oₚ anomaly (R = 0.89, n = 13). Note that the correlation is done with point XIII corrected, as we assume a bias in the model. Refer to Table 3 for details on reconstructed precipitation. Points falling on the line depict the same anomaly in the reconstruction and the simulation. Note that Timta and Dongge (Wang et al., 2005) Cave δ¹⁸Oₚ values were estimated from the Younger Dryas excursion, sometimes called H0 and characterized by a large freshwater input in the North Atlantic (Pausata et al., 2011).</figcaption></figure>

in Greenland ice cores over Heinrich events (Sect. 3.1.2; see Fig. 4a for δ¹⁸Oₚ anomaly).

Finally Eq. (4) also requires an estimate of δ¹⁸O₉ₗ. Groundwater pumped through the plant's roots represents a mixture of stored water and incoming precipitation water (McGuire et al., 2002). Indeed, during spring/summer, when maximal productivity occurs, the groundwater is composed of significant amounts of fall/winter precipitation. In order to take the mixing into account we use amount-weighted annual-mean δ¹⁸O of precipitation. This approach has been shown to be realistic in a field experiment in Kenya (Wang et al., 2012) and has been implemented in a global coupled

<table><caption>**Table 1.** Plant functional types (PFT) in ORCHIDEE, abbreviation used and mega-biome assignment in this study.</caption><thead><tr><th>PFT</th><th>Abbreviation</th><th>Mega-biome</th></tr></thead><tbody><tr><td>Bare soil</td><td>Bare soil</td><td>Bare soil</td></tr><tr><td>Tropical broadleaf evergreen trees</td><td>TrBE</td><td>Tropical trees</td></tr><tr><td>Tropical broadleaf raingreen trees</td><td>TrBR</td><td>Tropical trees</td></tr><tr><td>Temperate needleleaf evergreen trees</td><td>TempNE</td><td>Temperate trees</td></tr><tr><td>Temperate broadleaf evergreen trees</td><td>TempBE</td><td>Temperate trees</td></tr><tr><td>Temperate broadleaf summergreen trees</td><td>TempBS</td><td>Temperate trees</td></tr><tr><td>Boreal needleleaf evergreen trees</td><td>BoNE</td><td>Boreal trees</td></tr><tr><td>Boreal broadleaf summergreen trees</td><td>BoBS</td><td>Boreal trees</td></tr><tr><td>Boreal needleleaf summergreen trees</td><td>BoNS</td><td>Boreal trees</td></tr><tr><td>C3 grass</td><td>C3 grass</td><td>C3 grass</td></tr><tr><td>C4 grass</td><td>C4 grass</td><td>C4 grass</td></tr></tbody></table>

climate model of intermediate complexity (Caley et al., 2014). When implementing this, we neglect the fractionation effects that can significantly affect the soil water isotopic composition, especially in dry regions (Kanner et al., 2014).

## 2.4 Oxygen uptake in respiratory processes

### 2.4.1 Global oxygen production

The geographical distribution of respiratory O₂ fluxes (noted GPP_O₂ hereafter from the hypothesis of equilibrium between oxygen production and consumption) is computed from the vegetation cover and gross primary productivity (GPP) distribution provided by the ORCHIDEE model, which simulates the vegetation distribution of 10 natural plant functional types (PFTs) and bare soil (Krinner et al., 2005). ORCHIDEE is the land surface component of the IPSL-CM4 model. It is used here with the same spatial resolution as the aforementioned models (latitude × longitude 3.75° × 2.5°) and is run in offline mode, forced by the high-frequency outputs from IPSL-CM4 for the two experiments LGMCtrl and HS_exp. The prescribed forcing has a 6h time step temporal resolution, and thus takes into account daily variability and diurnal cycle simulated by the IPSL model (Woillez, 2012). Each grid cell simulates the vegetation cover by splitting the cell into fractions of 10 PFTs (see Table 1 for names) and provides GPP fluxes (g C m⁻² yr⁻¹) for each PFT. Validation of the vegetation cover can be found in Woillez et al. (2011) for LGMCtrl and Woillez et al. (2013) for HS_exp over western Europe. Here we extend this work and compare the LGM and HS simulation with worldwide vegetation reconstructions selected from marine (coastal) and terrestrial paleoarchives with high temporal resolution (Sect. 3.1.3).

In this study, following Hoffmann et al. (2004), we calculate the terrestrial biosphere's O₂ fluxes in three steps. First, the outputs from ORCHIDEE provide the GPP expressed in g C m⁻² yr⁻¹ for each of the 10 PFTs on each model grid point. Second, simulated carbon molar fluxes for each PFT are converted to oxygen molar fluxes, based on the biochemical model of photosynthesis from Farquhar et al. (1980). The
<table><caption>**Table 2.** Parameters involved in the calculation of δ¹⁸O_atm. Uncertainties are given for most of the parameters except for those derived from the ORCHIDEE model.</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Definition</th><th>LGMCtrl</th><th>HS_exp</th></tr></thead><tbody><tr><td>[CO<sub>2</sub>]</td><td>ppm</td><td>carbon dioxide mixing ratio in the troposphere</td><td>190</td><td>190</td></tr><tr><td>t</td><td>°C</td><td>temperature at the site of photosynthesis</td><td>21.08</td><td>21.41</td></tr><tr><td>h</td><td>%</td><td>relative humidity at the site of photosynthesis</td><td>66.09</td><td>66.12</td></tr><tr><td>GPP_C</td><td>Pmol C yr<sup>-1</sup></td><td>gross photosynthetic molar carbon flux from the terrestrial biosphere</td><td>6.758</td><td>6.450</td></tr><tr><td>GPP_O<sub>2</sub></td><td>Pmol O<sub>2</sub> yr<sup>-1</sup></td><td>gross photosynthetic molar oxygen flux from the terrestrial biosphere</td><td>11.768</td><td>11.410</td></tr><tr><td>f<sub>C4</sub></td><td>%</td><td>C4 fraction (in terms of GPP_C)</td><td>36.92</td><td>35.59</td></tr><tr><td>f<sub>photo</sub></td><td>%</td><td>fraction of photorespiration</td><td>28.55</td><td>29.52</td></tr><tr><td>f<sub>soildark</sub></td><td>%</td><td>soil fraction of dark respiration</td><td>62<sup>a</sup></td><td>62<sup>a</sup></td></tr><tr><td>f<sub>darksoil</sub></td><td>%</td><td>fraction of soil respiration</td><td>38.10</td><td>37.50</td></tr><tr><td>f<sub>darkleaves</sub></td><td>%</td><td>fraction of leaf respiration</td><td>23.35</td><td>22.98</td></tr><tr><td>f<sub>Mehler</sub></td><td>%</td><td>fraction of Mehler respiration</td><td>10<sup>b</sup></td><td>10<sup>b</sup></td></tr><tr><td>¹⁸ε<sub>darksoil</sub></td><td>%</td><td>global isotopic fractionation associated with dark soil respiration</td><td>16.242 ± 0.5<sup>c</sup></td><td>16.056 ± 0.5<sup>c</sup></td></tr><tr><td>¹⁸ε<sub>darkleaves</sub></td><td>%</td><td>mitochondrial (AOX + COX) isotopic fractionation in leaves</td><td>19 ± 1<sup>c</sup></td><td>19 ± 1<sup>c</sup></td></tr><tr><td>¹⁸ε<sub>Mehler</sub></td><td>%</td><td>global Mehler respiration isotopic fractionation</td><td>10.8 ± 0.2<sup>d</sup></td><td>10.8 ± 0.2<sup>d</sup></td></tr><tr><td>¹⁸ε<sub>photo</sub></td><td>%</td><td>global photorespiration isotopic fractionation</td><td>21.4 ± 1<sup>d</sup></td><td>21.4 ± 1<sup>d</sup></td></tr><tr><td>¹⁸ε<sub>resp</sub></td><td>%</td><td>global terrestrial respiration isotopic fractionation</td><td>17.83</td><td>17.80</td></tr><tr><td>δ¹⁸O<sub>p amount</sub></td><td>%</td><td>global precipitation water isotope delta</td><td>-6.689</td><td>-6.781</td></tr><tr><td>δ¹⁸O<sub>p</sub></td><td>%</td><td>global photosynthesis precipitation water isotope delta</td><td>-5.530</td><td>-5.289</td></tr><tr><td>δ¹⁸O<sub>vap amount</sub></td><td>%</td><td>global water vapor isotope delta</td><td>-12.648</td><td>-12.653</td></tr><tr><td>δ¹⁸O<sub>vap</sub></td><td>%</td><td>global photosynthesis water vapor isotope delta</td><td>-12.483</td><td>-12.295</td></tr><tr><td>δ¹⁸O<sub>leafwater</sub></td><td>%</td><td>global leaf water isotope delta</td><td>5.164 ± 1<sup>e</sup></td><td>5.301 ± 1<sup>e</sup></td></tr><tr><td>δ¹⁸O<sub>terr</sub></td><td>%</td><td>global terrestrial tropospheric isotope delta</td><td>23.407 ± 1</td><td>23.516 ± 1</td></tr><tr><td>δ¹⁸O<sub>mar</sub></td><td>%</td><td>global marine tropospheric isotope delta</td><td>25.3 ± 2<sup>f</sup></td><td>25.3 ± 2<sup>f</sup></td></tr><tr><td>δ¹⁸O<sub>atm</sub></td><td>%</td><td>global tropospheric isotope delta</td><td>23.88 ± 2</td><td>23.95 ± 2</td></tr></tbody></table>

### 3.1.1 Simulated humidity validation

According to climate reconstructions, during HS, wet periods in northeastern Brazil are synchronous with periods of weak East Asian summer monsoons (Wang et al., 2001) and with cold periods in Greenland (Grootes and Stuiver, 1997) and Europe (Genty et al., 2003). Reorganization in tropical rainfall patterns leads to wetter conditions in southwestern North America (Asmerom et al., 2010) and southern South America (Kanner et al., 2012), and to dryer conditions in the Australian–Indonesian monsoon region (Mohtadi et al., 2011), wide parts of Asia (Wang et al., 2008), northern South America (Peterson and Haug, 2006), Mediterranean region (Fleitmann et al., 2009), and equatorial western Africa (e.g., Weldeab, 2012).

During HS, the model simulates similar humidity patterns as reconstructed ones, with dryer conditions over Europe, the Mediterranean region, northern and equatorial Africa, southern and eastern Asia, Middle East, India, southern Australia and parts of Indonesia. In South America, a region of particular interest, where major simulated changes in vegetation and oxygen production occur, the model captures well the observed contrast with increased moisture in northeast-

ern Brazil, and drying in northern South America and Central America.

These rapid comparisons show that there is a good general agreement between modeled changes in humidity over an HS and climatic reconstructions over the different regions.

### 3.1.2 Simulated amount-weighted δ¹⁸O<sub>p</sub> validation

#### Tropics

Comparisons of modeled hosing-driven amount-weighted Δδ¹⁸O<sub>p</sub> anomalies with reconstructed Δδ¹⁸O<sub>p</sub> of speleothem's calcite during HSs are presented in Table 3 and Fig. 4. Thirteen Heinrich Δδ¹⁸O<sub>p</sub> proxy reconstructions arise from Lewis et al. (2010) (n = 11) and Pausata et al. (2011) (n = 4). They are located in the eastern Mediterranean, and in the regions of the South American monsoon, East Asian monsoon, Indian summer monsoon, North American monsoon and the Australian–Indonesian monsoon. Those regions represent the most productive ones and therefore carry a substantial part of the δ¹⁸O<sub>terr</sub> signal. Although reconstructed Δδ¹⁸O<sub>p</sub> anomalies from Lewis et al. (2010) and Pausata et al. (2011) studies were estimated differently, they are consistent and common reconstructed Δδ¹⁸O<sub>p</sub> estimates (for Hulu and Songjia Cave) are similar in both

<footnote><sup>a</sup>(Schlesinger and Andrews, 2000). Note that this estimation is for present day, and here we assume it was similar during the last glacial period. <sup>b</sup>(Badger et al., 2000). <sup>c</sup>(Landais et al., 2007a). <sup>d</sup>(Helman et al., 2005). <sup>e</sup>(Gillon and Yakir, 2001). <sup>f</sup>Note that the increase of 1‰ compared to the Luz et al. (2014) value accounts for the 1‰ enrichment of the glacial ocean (Waelbroeck et al., 2002).</footnote>
studies (Table 3, this study; method section of Pausata et al., 2011).

The dominant hydrological controls on reconstructed δ¹⁸Oₚ are site-specific and are described by Lewis et al. (2010). Figure 4 demonstrates the ability of the AOGCM LMDZ-iso to reproduce the observed Δδ¹⁸Oₚ spatial pattern for most of the sites, particularly in regions strongly affected by ITCZ (and its land extension) variations and hence by changes in the water cycle, regardless of the processes at play. There are two regions where the model does not properly reproduce the observed signal over HS. The first one occurs in the Indian summer monsoon domain (Timta Cave). Model and observation would reconcile two grid cells south of Timta Cave, as it is located just at the transition between a positive and negative simulated δ¹⁸Oₚ anomaly. This disagreement can be due to a model bias. The IPSL model indeed does not simulate the monsoonal signal at the right place, with an Indian monsoon located too far south even for modern climate (Marti et al., 2010). In the LGM simulation also, the IPSL model predicts that it mostly takes place over the ocean (−0.5 to −2 mm day⁻¹ is only simulated over the ocean; see Fig. 9 (lower panel) of Kageyama et al., 2009), while there is evidence for a monsoonal signal over land. In northern India, i.e., the Timta Cave site, the model does not simulate any significant rainfall change between the two periods. A more intense weakening of the Indian monsoon over land in the HS run, and hence less rainfall, would have helped in reconciling the model and data at Timta Cave, since δ¹⁸Oₚ would have been enriched through the amount effect. Pausata et al. (2011) recently suggested that change in rainfall amount associated with Indian monsoon rather than in southeastern Asia explains changes observed in calcite δ¹⁸O in Chinese stalagmites (in southeastern Asia). As in Pausata et al. (2011), a freshwater impulse was applied to the control simulation with LGM background climate. Rainfall amount drops in eastern Asia and northwestern India, mostly over the ocean, but increases in southeastern India, as shown in Fig. 3d. Values of −0.17 and −0.13 mm day⁻¹ are simulated at Hulu and Songjia Cave during HS, respectively. The enrichment in δ¹⁸Oₚ observed in Chinese caves is reproduced by the model, but the latter fails to capture the enrichment in Timta Cave. Overall, δ¹⁸Oₚ is enriched over the whole of India (with an abrupt change south of Timta Cave) and southern Asia. The possible role of the Indian monsoon in the oxygen isotopic enrichment of Chinese stalagmites is limited in our simulation, probably because the monsoonal signal is located too far south in the IPSL model. The increase in δ¹⁸O over southeastern Asia is consistent with local amount effect.

Another mismatch occurs in the North American monsoon domain (Cave of the Bells), where the observed Δδ¹⁸Oₚ (−0.8‰) and the modeled Δδ¹⁸Oₚ (0.9‰) are of opposite sign. The elevation of the site (1700 m a.s.l.) might explain the disagreement between model and data, in a region where the coarse model resolution does not allow for the role of orography to be properly represented. At Timta Cave and

Cave of the Bells, our model fails to capture the calcite δ¹⁸O anomaly recorded in speleothems. These two sites are located at high altitude and do not correspond to the regions where most of the oxygen is produced.

As shown in Table 3, the modeled increase in δ¹⁸Oₚ quantitatively agrees with data δ¹⁸Oc increase during HS in most of the compared sites (Fig. 4). In conclusion, the key features of HS precipitation inferred from speleothem's δ¹⁸Oc, i.e., a low-latitude interhemispheric see-saw pattern (Cheng et al., 2012), are generally well captured by the LMDZ model.

### High latitudes

In Greenland, HS can hardly be distinguished from the GS or from the mean LGM state. The only clear δ¹⁸Oice signal is observed from GI to GS (or HS) with an approximately 4‰ decrease in central Greenland sites (GRIP, GISP2, NGRIP). The depletion simulated in Greenland, with a 1.6‰ decrease at the GRIP site, by the model for a HS compared to a glacial background state does not compare well with available data. However, it is difficult to compare the δ¹⁸O change simulated by a freshwater input (the most efficient way to model a Heinrich event, as mentioned in Sect. 2.3.1) and the δ¹⁸O depletion between a GI and a GS. Indeed, there is more and more evidence that the δ¹⁸Oₚ depletion at the end of a GI is not due to the same freshwater discharge than the one associated with a Heinrich event. It can well be due to a threshold in the extent in sea ice or an atmospheric heat transport. Therefore our choice of modeling approach may potentially explain some of the discrepancies observed in the low latitudes, but our approach is the best we can realize today.

### 3.1.3 Validation of simulated vegetation

In order to compare model and data easily, simulated PFTs are gathered into five mega-biomes (boreal, temperate and tropical trees, C3 and C4 grasses) as well as bare soil. We distinguish between C3 and C4 plants as their partitioning has a strong impact on photorespiration fraction. The simulated dominant vegetation fraction is shown for LGMCtrl (Fig. 5) and HS_exp (Fig. 6), together with pollen-based reconstructed mega-biomes. Given its domination, we display bare soil fraction only if it covers more than 80% of a grid cell.

### Global oxygen production

Present-day carbon and oxygen productions amount to 10.5 Pmol C yr⁻¹ and 17.95 Pmol O₂ yr⁻¹ (taking into account photorespiration) in the ORCHIDEE model, respectively. This is in line with other estimates, e.g., Angert et al. (2003) or Welp et al. (2011), estimating 8 to 13 Pmol C yr⁻¹ and 12.5 to 14.2 Pmol O₂ yr⁻¹, respectively. For the LGM and HS, land carbon production estimates from the ORCHIDEE model are rather low, 6.8 and 6.5 Pmol C yr⁻¹,
<table><caption>**Table 3.** Comparison of isotopic proxy records (speleothem's calcite δ¹⁸O) with annual average modeled amount-weighted δ¹⁸O<sub>p</sub>. Note that anomalies from Pausata et al. (2011) are calculated from H1 and YD, while anomalies from Lewis et al. (2010) arise from all identifiable δ¹⁸O<sub>c</sub> excursions.</caption><thead><tr><th>Core</th><th>ID</th><th>Region</th><th>Latitude</th><th>Longitude</th><th>Data Δδ<sup>18</sup>O</th><th>Model Δδ<sup>18</sup>O</th><th>Reference</th></tr></thead><tbody><tr><td>Hulu Cave</td><td>XI</td><td>China</td><td>32.5</td><td>119.2</td><td>1.4</td><td>1.1</td><td>Pausata et al. (2011); Lewis et al. (2010)</td></tr><tr><td>Songjia Cave</td><td>XII</td><td>China</td><td>32.3</td><td>107.2</td><td>1.4</td><td>1.1</td><td>Pausata et al. (2011); Lewis et al. (2010)</td></tr><tr><td>Dongge Cave</td><td>IX</td><td>China</td><td>25.3</td><td>108.8</td><td>1</td><td>0.8</td><td>Pausata et al. (2011)</td></tr><tr><td>Timta Cave</td><td>XIII</td><td>India</td><td>29.8</td><td>80.0</td><td>3</td><td>-3.7</td><td>Pausata et al. (2011)</td></tr><tr><td>Sanbao Cave</td><td>X</td><td>China</td><td>31.7</td><td>110.5</td><td>1.2</td><td>1.1</td><td>Lewis et al. (2010)</td></tr><tr><td>Borneo</td><td>VII</td><td>Indonesia</td><td>4.0</td><td>114.0</td><td>0.8</td><td>0.6</td><td>Lewis et al. (2010)</td></tr><tr><td>Moomi Cave</td><td>VIII</td><td>Yemen</td><td>12.5</td><td>54.3</td><td>0.9</td><td>0.6</td><td>Lewis et al. (2010)</td></tr><tr><td>Soreq Cave</td><td>VI</td><td>Israel</td><td>31.5</td><td>35.0</td><td>0.5</td><td>0.5</td><td>Lewis et al. (2010)</td></tr><tr><td>Rio Grande do Norte</td><td>II</td><td>northeastern Brazil</td><td>-5.7</td><td>-37.7</td><td>-1.6</td><td>-0.1</td><td>Lewis et al. (2010)</td></tr><tr><td>Santana Cave</td><td>V</td><td>southern Brazil</td><td>-24.5</td><td>-48.7</td><td>-0.8</td><td>-0.9</td><td>Lewis et al. (2010)</td></tr><tr><td>Botuvera Cave</td><td>IV</td><td>southern Brazil</td><td>-27.2</td><td>-49.2</td><td>-1.1</td><td>-1.1</td><td>Lewis et al. (2010)</td></tr><tr><td>Cave of the Bells</td><td>III</td><td>North America</td><td>31.7</td><td>-110.8</td><td>-0.8</td><td>0.8</td><td>Lewis et al. (2010)</td></tr><tr><td>Poleva Cave</td><td>I</td><td>Europe</td><td>44.7</td><td>21.8</td><td>-2</td><td>-1.0</td><td>Lewis et al. (2010)</td></tr></tbody></table>

which translates into 11.8 and 11.4 Pmol O₂ yr⁻¹ for the LGM and HS, respectively. This is up to a factor of 2 lower than model-based LGM estimates from Joos et al. (2004), Hoffmann et al. (2004) or Bender et al. (1994), ranging from 23 to 16.7 Pmol O₂ yr⁻¹.

The ORCHIDEE model is known to underestimate LGM productivity at both low latitudes (too low productivity in tropical forests, especially Amazonia) and high latitudes (too low productivity in the absence of permafrost modeling).

Photorespiration fraction (see Sect. 2.4.2) may also be invoked to explain part of the model-data discrepancy for the LGM and HS. Underestimation of photorespiration may arise from uncertainties related to the time of photosynthesis. In the real world, plants must reduce their CO₂ uptake under water stress, as stomata close to preclude water loss. This leads to a higher proportion of photorespiration, not necessarily considered during experiments performed under ideal hydric conditions, whose results are used in the classical Farquhar parameterization (Farquhar et al., 1980).

The classical scaling factor between carbon uptake and oxygen production (Keeling, 1988) of 1.07 used in our study may also have been underestimated. Indeed, plants can produce oxygen without involving carbon uptake during times of stress, which is not considered in experiments run under ideal conditions. The oxygen production calculated here from the ORCHIDEE model seems to be biased toward too low values for the LGM. The same bias is true for HS. Still, it should be noted that the δ¹⁸Oₜerrett calculations of our study and hence the final results do not depend on the absolute value of oxygen production at the LGM and HS.

## LGM\_ctrl vegetation

The main features of the glacial vegetation are correctly reproduced by the ORCHIDEE model, as briefly presented by Woillez et al. (2011): reduced fractions of tropical forest, particularly in Amazonia, and high grass fractions in Siberia,

Alaska, and western North America. Main biases are an overestimation of the tree fractions over western Europe, eastern Eurasia and eastern North America, as well as an overestimation in bare soil fractions over India, southern Africa, Siberia and South America (Woillez et al., 2011), leading to an underestimation of the global carbon production, as mentioned in the previous section.

The model simulates temperate trees in southeastern Asia (Vietnam, southern China, Cambodia), tropical trees and grasses over the western pacific warm pool (Malaysia, Thailand, Indonesia) and over southern Africa, in agreement with BIOME6000 reconstructions (Prentice et al., 2000). The model underestimates temperate trees in Asia and overestimates bare soil in South Africa.

A more detailed comparison shows that the important (bare) tree fraction over southwestern Europe differs from palynological reconstructions depicting an important grass fraction, but this bias mainly comes from the overestimation of Boreal broadleaf summergreen trees, which is a common feature in the version of ORCHIDEE used here, also found in present-day vegetation simulations. The bias might also be the sign that the LGM climate simulated by IPSL-CM4 over western Europe is too warm and wet (Woillez et al., 2011). In a few regions, ORCHIDEE correctly simulates the presence of forest, but the dominant type of tree disagrees with pollen reconstructions: tropical trees over Papua New Guinea and western Indonesia, while reconstructions reveal the presence of temperate trees over these regions. Simulated forests over southern Australia (a thin coastal band in the southeast) are composed of temperate and tropical trees, while reconstructions rather indicate the presence of a few tropical trees. The model simulates mixed vegetation composed of grasses, bare and temperate trees in eastern North America, consistent with pollen data, but the spatial distribution is incorrect.

It is important to keep in mind that model-data comparison of vegetation can only remain qualitative given the coarse resolution of the vegetation model, related to the model res-
<table><caption>**Table 4.** Comparison of mega-biomes during Heinrich stadials between pollen reconstructions (references are included in the table) and simulated vegetation (compiled from HS_exp using the ORCHIDEE vegetation model). Note that simulated C3 and C4 grasses are merged into one mega-biome because pollen-based biome reconstructions do not allow us to distinguish between the two PFTs.</caption><thead><tr><th rowspan="2">Core</th><th rowspan="2">ID</th><th rowspan="2">Region</th><th rowspan="2">Latitude</th><th rowspan="2">Longitude</th><th rowspan="2">Resolution<sup>a</sup></th><th rowspan="2">Which HS?</th><th colspan="3">Mega-biome distribution</th><th rowspan="2">Agreement</th><th rowspan="2">Reference</th></tr><tr><th>(yr/sample)</th><th>pollen data</th><th>biome(s) designation</th><th>model results dominant, subdominant biome</th></tr></thead><tbody><tr><td>Kashiru Bog</td><td>1</td><td>equatorial Africa</td><td>-3.47</td><td>29.57</td><td>410</td><td>HS1</td><td>grassland and dry shrubland, savannah and xerophytic scrubland</td><td>grasses</td><td>tropical forest, grasses</td><td>fair</td><td>Hessler et al. (2010); Handiani et al. (2012)</td></tr><tr><td>Lake Tanganyika</td><td>2</td><td>equatorial Africa</td><td>-8.5</td><td>30.85</td><td>610</td><td>HS1</td><td>warm temperate mixed forest, savannah and xerophytic scrubland</td><td>temperate forest, grasses</td><td>grasses, tropical forest</td><td>fair</td><td>Hessler et al. (2010); Handiani et al. (2012)</td></tr><tr><td>Lake Masoko</td><td>3</td><td>equatorial Africa</td><td>-9.33</td><td>33.75</td><td>550</td><td>HS1</td><td>warm temperate mixed forest, savannah and xerophytic scrubland</td><td>temperate forest, grasses</td><td>grasses, tropical forest</td><td>fair</td><td>Hessler et al. (2010); Handiani et al. (2012)</td></tr><tr><td>Lake Malawi</td><td>4</td><td>equatorial Africa</td><td>-11.29</td><td>34.44</td><td>200</td><td>HS1</td><td>savannah and xerophytic scrubland, tropical forest</td><td>grasses, tropical forest</td><td>grasses</td><td>good</td><td>Hessler et al. (2010); Handiani et al. (2012)</td></tr><tr><td>Barombi Mbe</td><td>5</td><td>equatorial Africa</td><td>4.51</td><td>9.4</td><td>590</td><td>HS1</td><td>savannah and xerophytic scrubland, tropical forest</td><td>grasses, tropical forest</td><td>grasses, tropical forest</td><td>good</td><td>Hessler et al. (2010); Handiani et al. (2012)</td></tr><tr><td>KS 84-063</td><td>6</td><td>equatorial Africa</td><td>4.4</td><td>-4.18</td><td>450</td><td>HS1</td><td>tropical forest, warm temperate mixed forest</td><td>tropical forest, grasses</td><td>tropical forest, grasses</td><td>good</td><td>Hessler et al. (2010); Handiani et al. (2012)</td></tr><tr><td>ODP 1078-C</td><td>7</td><td>equatorial Africa</td><td>-11.92</td><td>13.4</td><td>140</td><td>HS1</td><td>warm temperate mixed forest, temperate montane forest</td><td>temperate forest, boreal forest</td><td>bare soil, tropical forest</td><td>bad<sup>b</sup> (soil > 90 %)</td><td>Hessler et al. (2010); Handiani et al. (2012)</td></tr><tr><td>GEOB 1023 – Cunene River Mouth</td><td>8</td><td>equatorial Africa</td><td>-17.15</td><td>11.02</td><td>185</td><td>HS1</td><td>savannah and xerophytic scrubland, grassland and dry shrubland</td><td>grasses</td><td>none</td><td>none</td><td>Hessler et al. (2010); Handiani et al. (2012)</td></tr><tr><td>Lake Caco</td><td>9</td><td>South America</td><td>-2.97</td><td>-43.42</td><td>80</td><td>HS1</td><td>warm temperate mixed forest, tropical forest</td><td>temperate forest, tropical forest</td><td>bare soil, tropical forest</td><td>moderate</td><td>Hessler et al. (2010); Handiani et al. (2012)</td></tr><tr><td>Colonia</td><td>10</td><td>South America</td><td>-23.87</td><td>-46.71</td><td>710</td><td>HS1</td><td>savannah and xerophytic scrubland, grassland and dry shrubland</td><td>grasses</td><td>grasses, temperate forest</td><td>good</td><td>Hessler et al. (2010); Handiani et al. (2012)</td></tr><tr><td>La Laguna, Bogota</td><td>11</td><td>South America</td><td>4.92</td><td>-74.03</td><td>670</td><td>HS1</td><td>savannah and xerophytic scrubland, grassland and dry shrubland</td><td>grasses</td><td>bare soil</td><td>bad (soil = 100 %)</td><td>Hessler et al. (2010); Handiani et al. (2012)</td></tr><tr><td>Fuquene</td><td>12</td><td>South America</td><td>4.92</td><td>-74.03</td><td>520</td><td>HS1</td><td>savannah and xerophytic scrubland, temperate montane forests</td><td>grasses, temperate forest</td><td>bare soil</td><td>bad (soil = 100 %)</td><td>Hessler et al. (2010); Handiani et al. (2012)</td></tr><tr><td>GEOB 3104</td><td>13</td><td>South America</td><td>-3.67</td><td>-37.72</td><td>670</td><td>HS1</td><td>temperate montane forest, warm temperate mixed forest</td><td>temperate forest, grasses</td><td>tropical forest, bare soil</td><td>bad<sup>b</sup></td><td>Hessler et al. (2010); Handiani et al. (2012)</td></tr><tr><td>GEOB 3910-2</td><td>14</td><td>South America</td><td>-4.15</td><td>-36.21</td><td>125</td><td>HS1</td><td>savannah and xerophytic scrubland, warm temperate mixed forests</td><td>temperate forest</td><td>tropical forest, bare soil</td><td>bad<sup>b</sup></td><td>Hessler et al. (2010); Handiani et al. (2012)</td></tr><tr><td>MD03-2622</td><td>15</td><td>South America</td><td>10.71</td><td>-65.17</td><td>420</td><td>HS3, HS4, HS5</td><td>montane forest, semi-deciduous forest, savannah (except HS4)</td><td>temperate forest, grasses</td><td>none</td><td>none</td><td>Hessler et al. (2010); Handiani et al. (2012)</td></tr><tr><td>17 962</td><td>16</td><td>Australasia</td><td>7.18</td><td>112.08</td><td>370</td><td>HS4</td><td>tropical forest</td><td>tropical forest</td><td>tropical forest, grasses</td><td>good</td><td>Harrison and Sanchez-Goni (2010)</td></tr><tr><td>18 300</td><td>17</td><td>Australasia</td><td>4.35</td><td>108.65</td><td>526</td><td>HS4</td><td>tropical forest</td><td>tropical forest</td><td>tropical forest, grasses</td><td>good</td><td>Harrison and Sanchez-Goni (2010)</td></tr><tr><td>18 323</td><td>18</td><td>Australasia</td><td>2.78</td><td>107.88</td><td>420</td><td>HS4</td><td>tropical forest</td><td>tropical forest</td><td>tropical forest, grasses</td><td>good</td><td>Harrison and Sanchez-Goni (2010)</td></tr><tr><td>Lake Wangoom</td><td>19</td><td>Australasia</td><td>-38.35</td><td>142.6</td><td>362</td><td>HS4</td><td>herbaceous and shrublands</td><td>grasses</td><td>temperate forest, tropical forest</td><td>bad</td><td>Harrison and Sanchez-Goni (2010)</td></tr><tr><td>Tyrendarra Swamp</td><td>20</td><td>Australasia</td><td>-38.2</td><td>141.76</td><td>337</td><td>HS4</td><td>herbaceous and shrublands</td><td>grasses</td><td>temperate forest, tropical forest</td><td>bad</td><td>Harrison and Sanchez-Goni (2010)</td></tr><tr><td>Lake Surprise</td><td>21</td><td>Australasia</td><td>-38.06</td><td>141.92</td><td>345</td><td>HS4</td><td>herbaceous and shrublands</td><td>grasses</td><td>temperate forest, tropical forest</td><td>bad</td><td>Harrison and Sanchez-Goni (2010)</td></tr><tr><td>Kohuora</td><td>22</td><td>Australasia</td><td>-36.57</td><td>174.52</td><td>375</td><td>HS4</td><td>herbaceous and shrublands</td><td>grasses</td><td>temperate forest, tropical forest</td><td>bad</td><td>Harrison and Sanchez-Goni (2010)</td></tr><tr><td>Native Companion Lagoon</td><td>23</td><td>Australasia</td><td>-27.68</td><td>153.41</td><td>655</td><td>HS4</td><td>tropical forest and open forest, woodland</td><td>tropical forest, temperate forest</td><td>bare soil, temperate forest</td><td>moderate</td><td>Harrison and Sanchez-Goni (2010)</td></tr><tr><td>Ioannina 284</td><td>24</td><td>Europe</td><td>39.75</td><td>20.85</td><td>325</td><td>HS4</td><td>grassland and dry shrubland</td><td>grasses</td><td>boreal forest, grasses</td><td>fair</td><td>Fletscher et al. (2010)</td></tr><tr><td>Megali Limni</td><td>25</td><td>Europe</td><td>39.1</td><td>26.32</td><td>150</td><td>HS4</td><td>grassland and dry shrubland with 40 % xerophytic steppe elements</td><td>grasses</td><td>boreal forest</td><td>good</td><td>Fletscher et al. (2010)</td></tr><tr><td>Lago Grande di Monticchio</td><td>26</td><td>Europe</td><td>40.93</td><td>15.62</td><td>210</td><td>HS4</td><td>grassland and dry shrubland with 40 % xerophytic steppe elements</td><td>grasses</td><td>boreal forest, grasses</td><td>fair</td><td>Fletscher et al. (2010)</td></tr><tr><td>MD04-2845</td><td>27</td><td>Europe</td><td>45.35</td><td>-5.22</td><td>540</td><td>HS3</td><td>grassland and dry shrubland</td><td>grasses</td><td>none</td><td>none</td><td>Fletscher et al. (2010)</td></tr><tr><td>MD99-2331</td><td>28</td><td>Europe</td><td>41.15</td><td>-9.68</td><td>390</td><td>HS4</td><td>grassland and dry shrubland</td><td>grasses</td><td>none</td><td>none</td><td>Fletscher et al. (2010)</td></tr><tr><td>MD95-2039</td><td>29</td><td>Europe</td><td>40.58</td><td>-10.35</td><td>300</td><td>HS4</td><td>grassland and dry shrubland</td><td>grasses</td><td>none</td><td>none</td><td>Fletscher et al. (2010)</td></tr><tr><td>MD95-2042</td><td>30</td><td>Europe</td><td>37.8</td><td>-10.17</td><td>360</td><td>HS4</td><td>grassland and dry shrubland with 40 % xerophytic steppe elements</td><td>grasses</td><td>none</td><td>none</td><td>Fletscher et al. (2010)</td></tr></tbody></table>
<table><caption>**Table 4. Continued.**</caption><thead><tr><th rowspan="2">Core</th><th rowspan="2">ID</th><th rowspan="2">Region</th><th rowspan="2">Latitude</th><th rowspan="2">Longitude</th><th rowspan="2">Resolution<sup>a</sup></th><th rowspan="2">Which HS?</th><th colspan="3">Mega-biome distribution</th><th rowspan="2">Agreement</th><th rowspan="2">Reference</th></tr><tr><th>(yr/sample)</th><th>pollen data</th><th>biome(s) assignment</th><th>model results dominant, subdominant biome</th></tr></thead><tbody><tr><td>ODP site 976</td><td>31</td><td>Europe</td><td>36.2</td><td>-4.3</td><td>240</td><td>HS4</td><td>grassland and dry shrubland with 40% xerophytic steppe elements</td><td>grasses</td><td>bare soil, boreal forest</td><td>bad (soil &gt; 90%)</td><td>Fletscher et al. (2010)</td></tr><tr><td>MD95-2043</td><td>32</td><td>Europe</td><td>36.13</td><td>-2.62</td><td>260</td><td>HS4</td><td>grassland and dry shrubland with 40% xerophytic steppe elements</td><td>grasses</td><td>bare soil, boreal forest</td><td>bad (soil &gt; 90%)</td><td>Fletscher et al. (2010)</td></tr><tr><td>Khoe</td><td>33</td><td>Japan</td><td>51.34</td><td>142.14</td><td>750</td><td>HS4</td><td>cold deciduous and evergreen conifer forest</td><td>boreal forest</td><td>boreal forest, grasses</td><td>good</td><td>Takahara et al. (2010)</td></tr><tr><td>Kenbuchi</td><td>34</td><td>Japan</td><td>44.05</td><td>142.38</td><td>250</td><td>HS1, HS2</td><td>cold deciduous forest</td><td>boreal forest</td><td>boreal forest, grasses</td><td>good</td><td>Takahara et al. (2010)</td></tr><tr><td>MD01-2421</td><td>35</td><td>Japan</td><td>36.02</td><td>141.77</td><td>150</td><td>HS4</td><td>cold evergreen conifer forest</td><td>boreal forest</td><td>boreal forest, temperate forest</td><td>good</td><td>Takahara et al. (2010)</td></tr><tr><td>Lake Nojiri</td><td>36</td><td>Japan</td><td>36.83</td><td>138.22</td><td>100</td><td>HS4</td><td>increase in cold evergreen conifer forest within cool conifer forest</td><td>temperate forest, boreal forest</td><td>none</td><td>none</td><td>Takahara et al. (2010)</td></tr><tr><td>Lake Biwa</td><td>37</td><td>Japan</td><td>35.22</td><td>136</td><td>300</td><td>HS4</td><td>increase in cool conifer forest within temperate conifer forest</td><td>temperate forest</td><td>boreal forest, temperate forest</td><td>fair</td><td>Takahara et al. (2010)</td></tr><tr><td>Kamiyoshi Basin</td><td>38</td><td>Japan</td><td>35.1</td><td>135.59</td><td>800</td><td>HS4</td><td>increase in cool conifer within temperate conifer forest, and deciduous broadleaf forest</td><td>temperate forest</td><td>boreal forest, temperate forest</td><td>fair</td><td>Takahara et al. (2010)</td></tr><tr><td>Toushe Basin</td><td>39</td><td>Japan</td><td>23.82</td><td>120.88</td><td>300</td><td>HS4</td><td>temperate deciduous or warm temperate evergreen forest</td><td>temperate forest</td><td>boreal forest, temperate forest</td><td>good</td><td>Takahara et al. (2010)</td></tr><tr><td>Fargher Lake</td><td>40</td><td>North America</td><td>45.88</td><td>-122.58</td><td>270</td><td>HS4</td><td>boreal forest</td><td>boreal forest</td><td>none</td><td>none</td><td>Jimenez-Moreno et al. (2010)</td></tr><tr><td>Carp Lake</td><td>41</td><td>North America</td><td>45.91</td><td>-120.88</td><td>630</td><td>HS4</td><td>open temperate and pine forest</td><td>temperate forest</td><td>none</td><td>none</td><td>Jimenez-Moreno et al. (2010)</td></tr><tr><td>Little Lake</td><td>42</td><td>North America</td><td>44.16</td><td>-123.58</td><td>260</td><td>HS4</td><td>boreal-temperate forest</td><td>boreal forest, temperate forest</td><td>none</td><td>none</td><td>Jimenez-Moreno et al. (2010)</td></tr><tr><td>W8709A-13PC</td><td>43</td><td>North America</td><td>42.25</td><td>-127.66</td><td>430</td><td>HS4</td><td>boreal forest with decrease in heterophylla</td><td>boreal forest</td><td>none</td><td>none</td><td>Jimenez-Moreno et al. (2010)</td></tr><tr><td>EW-9504-17PC</td><td>44</td><td>North America</td><td>42.23</td><td>-125.81</td><td>460</td><td>HS1, HS2, HS3</td><td>warm temperate</td><td>temperate forest, tropical forest</td><td>none</td><td>none</td><td>Jimenez-Moreno et al. (2010)</td></tr><tr><td>ODP 893A</td><td>45</td><td>North America</td><td>34.28</td><td>-120.03</td><td>220</td><td>HS4</td><td>open temperate forest</td><td>temperate forest</td><td>temperate forest, bare soil grasses</td><td>good</td><td>Jimenez-Moreno et al. (2010)</td></tr><tr><td>Bear Lake</td><td>46</td><td>North America</td><td>41.95</td><td>-111.3</td><td>680</td><td>HS4</td><td>xerophytic shrubland</td><td>temperate forest</td><td>grasses</td><td>good</td><td>Jimenez-Moreno et al. (2010)</td></tr><tr><td>Camel Lake</td><td>47</td><td>North America</td><td>30.26</td><td>-85.01</td><td>300</td><td>HS4</td><td>temperate forest with increase in southeastern pine forest</td><td>temperate forest</td><td>boreal forest, temperate forest</td><td>fair</td><td>Jimenez-Moreno et al. (2010)</td></tr><tr><td>Lake Tulane</td><td>48</td><td>North America</td><td>27.58</td><td>-81.5</td><td>480</td><td>HS4</td><td>southeastern pine forest, florida scrub</td><td>grasses, temperate forest</td><td>bare soil, temperate forest</td><td>moderate</td><td>Jimenez-Moreno et al. (2010)</td></tr></tbody></table>

olution of the climatic forcing fields. Furthermore, pollen records represent the surrounding vegetation distribution at different altitudes, while the ORCHIDEE model does not account for elevation changes within a grid cell (Woillez et al., 2013).

## HS_exp vegetation

To validate the simulated HS vegetation, we compare the millennial-scale changes in selected high-resolution (< 800 years; mean resolution is 400 years) pollen records of 48 sites described for HS1-GS2 (n = 16) in South America and southern Africa (Hessler et al., 2010; Handiani et al., 2012), and for HS4-GS9 (n = 31) in Europe (Fletcher et al., 2010), North America (Jimenez-Moreno et al., 2010), Japan (Takahara et al., 2010) and Australasia (Harrison and Goni, 2010). Figure 7 displays the location of paleorecords discussed in this study. Table 4 summarizes the model-data comparison at a grid cell level and provides additional information revealed by palynological reconstructions. The sampling resolution for the analyzed period (MIS2 for HS1, MIS3 for HS4), the other biomes represented for a given site, and the potential

occurrences of similar reconstructed vegetation changes over other HSs are presented.

The model-data comparison has been performed as follows: the two dominant reconstructed biomes are compared with the two dominant simulated biomes over the grid cell covering the site where the proxy originates. Among the 48 sites with pollen reconstructions, 12 were discarded because of absence of vegetation on the considered grid cell. Among the remaining 36 terrestrial and coastal sites, 11 disagree (30%) and 25 (70%) display moderate to good agreement (Table 4). Good agreement (n = 13) is obtained when reconstructed and simulated dominant biomes are alike, fair agreement (n = 9) when a subdominant biome agrees with a dominant one, and moderate agreement (n = 3) when subdominant biomes only are similar.

Simulated vegetation in regions associated with high oxygen productivity agrees well with pollen reconstruction. This is the case for South America, where a strong increase in tropical forest at the expense of bare soil is simulated in eastern Brazil, and in the West Pacific Warm Pool region, where tropical forest represents the dominant biome. The slight southward shift of the southern border of desert areas

<footnote><sup>a</sup> Sampling resolutions of the MIS where vegetation changes occur. Mean sampling resolution is 393 years. <sup>b</sup> Similar to Handiani et al. (2012) model results.</footnote>
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

**Figure 5.** (a) Model-data comparison of mega-biome distribution for LGMCtrl based on dominant PFT type simulated by ORCHIDEE. For each grid cell, the fraction of bare soil, tropical forests, temperate forests, boreal forests, and C3 and C4 grasses is considered. The type covering the greatest cell fraction is the dominant type. Note that dominant bare soil fraction denotes more than 80%. Circles denote LGM mega-biomes inferred from pollen and plant macrofossil records compiled by the BIOME6000 project. Refer to Table 4 to see how PFTs simulated by ORCHIDEE have been assigned to the mega-biomes mapped in this figure. (b) Detail of the averaged vegetation composition in grid cells occupied by a dominant mega-biome for LGMCtrl.

in equatorial Africa is also well captured by ORCHIDEE. In the Indian summer monsoon region, the simulated weakened monsoon (Kageyama et al., 2009) leads to the appearance of desert areas south and east of India, consistent with a dryer climate revealed by a core from the Indus region (Deplazes et al., 2014).

From this analysis it appears that sites showing a disagreement between model and data are coastal sites and/or have a very high bare soil fraction. Coastal sites do not necessarily only represent the vegetation in the coastal region, but they offer numerous records and allow high-resolution analysis thanks to their high sedimentation rate, so it is crucial to include them for millennial-scale analysis. Five (50%) of the sites showing no agreement present a very high simulated bare soil fractions (> 90%). Woillez et al. (2013) already

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

**Figure 6.** (a) Model-data comparison of mega-biome distribution for HS_exp based on dominant PFT type simulated by ORCHIDEE. For each grid cell, the fraction of bare soil, tropical forests, temperate forests, boreal forests, and C3 and C4 grasses is considered. The type covering the greatest cell fraction is the dominant type. Note that dominant bare soil fraction denotes more than 80%. Circles denote HS_exp mega-biomes inferred from pollen and plant macrofossil records compilation. Refer to Table 4 to see how PFTs simulated by ORCHIDEE and reconstructed vegetation have been assigned to the mega-biomes mapped in this figure. (b) Detail of the averaged vegetation composition in grid cells occupied by a dominant mega-biome for HS_exp.

pointed out the overestimation of the bare soil fractions by the ORCHIDEE model. We argue that this bias might partly explain the observed discrepancy between model and data. Furthermore, the other sites showing a disagreement, over Europe and Australia, underestimate grass fraction. For Europe, the bias is already present in the LGM simulation and is probably partly due to the systematic overestimation of forest by ORCHIDEE in this region (Woillez et al., 2011). Given the scarcity of data offering a time resolution high enough to catch millennial-scale vegetation variability, further testing of the simulated vegetation remains challenging. In conclusion, HS_exp vegetation agrees reasonably well with available pollen-based vegetation reconstruction.

Finally, based on the reasonable agreement of the simulated changes in vegetation, humidity and precipitation with
<figure><img src="image_11.png" /></figure>

**Figure 7.** Location of marine and terrestrial paleoarchives sites included in the model-data comparison. Blue diamonds denote speleothem's calcite δ¹⁸O and red circles denote pollen records. Arabic numbers and Roman numerals displayed on the map identify the location of the paleoarchives listed in Tables 3 (hydrology) and 4 (vegetation) for site names, references and further details.

observations depicted in this section, we can rely on the va-
lidity of the model to simulate δ¹⁸Oₜerr over a HS.

3.2 Global increase in δ¹⁸Oₜerr during a HS

The model calculates δ¹⁸Oₜerr for LGMCtrl and HS_exp as 23.41 and 23.52‰, respectively (Table 2). This average δ¹⁸Oₜerr value is coherent with the δ¹⁸Oₐtmosphere value of 23.8‰ with respect to V-SMOW and the finding that terrestrial and marine contribution to δ¹⁸Oₐtmosphere are similar (Luz et al., 2014). Moreover, the global increase in δ¹⁸Oₜerr of 0.11‰ (Fig. 8) can quantitatively explain most of the 0.1‰ δ¹⁸Oₐtmosphere increase over HS (Severinghaus et al., 2009; Guillevic et al., 2014; Fig. 1). In the following, we use the different model outputs to decipher the main influences on δ¹⁸Oₜerr and hence on δ¹⁸Oₐtmosphere.

By construction (Eq. 3), δ¹⁸Oₜerr is linearly dependent on both δ¹⁸Oᵢₘ and ¹⁸εₜresp; we discuss these two effects below. Figure 8 (upper panel) details the different contributions to δ¹⁸Oₜerr change over a HS and demonstrates the dominant role of δ¹⁸Oᵢₘ (Sect. 3.2.1) compared to ¹⁸εₜresp (Sect. 3.2.2). Indeed, the 0.11‰ increase in δ¹⁸Oₜerr exclusively stems from δ¹⁸Oᵢₘ increase (+0.14‰ over HS_exp), while respiratory fractionation leads to a negative anomaly (−0.03‰) over HS_exp. We explore in more details below the origin of the relative changes in ¹⁸εₜresp and δ¹⁸Oᵢₘ as calculated by our modeling approach. In particular, we look at the different regional contributions to the global ¹⁸εₜresp and δ¹⁸Oᵢₘ signals since low-latitude regions are associated with the largest GPP_O₂ (Fig. 9a for the whole latitudinal range and Fig. 10a for a closeup of the tropics) and hence have the strongest influence on the global δ¹⁸Oₜerr signal (Fig. 9c).

3.2.1 Leaf water

We find global δ¹⁸Oᵢₘ values of 5.16 and 5.30‰ for LGMCtrl and HS_exp, respectively. The 0.14‰ difference is similar to the δ¹⁸O increase observed in ice cores during HS (Fig. 2). The increase in δ¹⁸Oᵢₘ is clearly visible in the low-latitude regions. It corresponds to an increase in δ¹⁸Oₚ and a decrease in GPP_O₂ weighted relative humidity, both effects leading to a global δ¹⁸Oᵢₘ increase (Fig. 9b).

Still, when looking at the whole latitudinal range, the GPP_O₂ weighted relative humidity is not significantly different in HS and in LGM state. This is due to the decrease in relative humidity during HS in the extratropical regions (Fig. 9b). The net effect of relative humidity on δ¹⁸Oᵢₘ is thus zero. As a consequence, the main driver of δ¹⁸Oᵢₘ (and hence δ¹⁸Oₜerr) increase is the increase in GPP_O₂-weighted δ¹⁸Oₚ (Fig. 9b) by 0.18‰ (Fig. 8). This increase is linked to the southward shift of the tropical belt occurring during HS, as suggested from the speleothem data compilation (see Sect. 3.1.2).

Figure 10 clearly shows how rainfall amount and δ¹⁸Oₚ
are anticorrelated as expected on most of the intertropical
band. During a HS, δ¹⁸Oₚ is enriched in the Northern Hemi-
sphere down to 14°S. A particular pattern occurs between
the Equator and 14°S, where oxygen production is most
enhanced at HS, as precipitation is more abundant but also
heavier in δ¹⁸Oₚ.

3.2.2 Respiration

Respiratory processes lead to a 0.03‰ decrease in δ¹⁸Oₜerr
in HS_exp compared to δ¹⁸Oₜerr in LGMCtrl (Fig. 8a). This
variation is too small to challenge δ¹⁸Oₜerr enrichment caused
by hydrological processes, but the sign of its anomaly raises
questions. Here we explain the stability of ¹⁸εₜresp on millen-
nial timescales by a compensatory effect taking place be-
<figure><img src="image_12.png" /></figure>

**Figure 8.** Evolution of the main simulated factors controlling atmospheric δ¹⁸O. Note that all variables are oxygen-production-weighted, i.e., integrated over vegetated areas, if not stated with “non w.”. (a) Left panel: HS_exp – LGMCtrl anomalies of temperature and relative humidity. Right panel: (left to right) HS_exp – LGMCtrl anomalies of seawater δ¹⁸O, amount-weighted precipitation δ¹⁸O, precipitation δ¹⁸O (δ¹⁸O<sub>p</sub>), water vapor δ¹⁸O, leaf water δ¹⁸O (δ¹⁸O<sub>lw</sub>), respiratory isotope fractionation (–¹⁸ε<sub>resp</sub>) and terrestrial contribution to atmospheric δ¹⁸O (δ¹⁸O<sub>terr</sub>). (b) LGMCtrl values of same relevant factors as in (a) in δ¹⁸O<sub>terr</sub> budget. Note that plotted respiratory isotope fractionation anomaly is inverted as respiration is an oxygen uptake process. Note also that the –¹⁸ε<sub>resp</sub> anomaly totals 0.02‰ when soil aeration influence on soil respiratory isotope fractionation is considered. δ¹⁸O<sub>lw</sub> is controlled by δ¹⁸O<sub>p</sub>, temperature and relative humidity as described by Eq. (4). Combined with ¹⁸ε<sub>resp</sub> as described in Eq. (3), one obtains δ¹⁸O<sub>terr</sub>.

tween the main respiratory pathways. We then carry out a
simple sensitivity experiment to get a better understanding
of the causes of ¹⁸αdark soil negative anomaly during HS.

¹⁸ε resp is classically separated into four contributions as
given in Eq. (6): soil respiration, leaf respiration, photorespi-
ration and Mehler reaction.

Soil respiration with associated fractionation factor
¹⁸α<sub>dark_soil</sub> represents 63 % of dark respiration (Schlesinger
and Andrews, 2000). It represents 39 % of global terrestrial
respiration in our simulation, in agreement with estimates
from Raich and Potter (1995) of 30 to 40 %. ¹⁸α<sub>dark_soil</sub>
is temperature-dependent (Angert et al., 2003), with higher
fractionation associated with colder soils, causing a high
zonal contrast. This temperature effect leads to a change in
¹⁸α<sub>dark_soil</sub> by 0.19 ‰ at HS, increasing ¹⁸ε<sub>resp</sub> by 0.12 ‰ and
in turn depleting δ¹⁸O<sub>terr</sub>.

The photorespiration fraction, strongly discriminative
against ¹⁸O (¹⁸α<sub>photo</sub> = –21.4‰), represents 28.56 % of the
global terrestrial oxygen uptake in LGMCtrl. The propor-
tion of photorespiration is mainly linked to the change in the
C3 vs. C4 plant proportions and in temperature during pho-
tosynthesis (Von Caemmerer, 2000; Hoffmann et al., 2004).
During HS_exp, photorespiration fraction increases by 1 %
(29.52 %), driven by a slight decrease in C4 grass (1.3 %)
and a slight photosynthesis temperature increase (0.3 °C). As
a result, change in photorespiration fraction and soil respira-
tion lead ¹⁸ε<sub>resp</sub> to decrease by 0.03 ‰ in HS_exp, thus caus-
ing δ¹⁸O<sub>terr</sub> enrichment.

In summary, in our model, a weaker ¹⁸α<sub>dark_soil</sub> during HS cancels out the effect of increased photorespiration on ¹⁸ε<sub>resp</sub>. Note that this compensation between the two main respiratory processes explains not only the temporal stability of ¹⁸ε<sub>resp</sub> but also its zonal stability. Indeed, in high latitudes, cold temperatures lead to a weak photorespiration but strong soil isotope fractionation. In low latitudes, despite a weak soil isotope fractionation, high temperatures and variable C4 fraction lead to a highly variable photorespiration.

The 0.19‰ increase in ¹⁸α<sub>dark_soil</sub> during HS needs to be considered carefully. Indeed, following the southward migration of the tropical rain belt during HS, tropical soils generally dry out during HSs, as inferred from the lower atmospheric nitrous oxide concentration and its isotopic composition during HS1 (Schilt et al., 2014). Well-aerated soils are associated with a strong discrimination (Angert et al., 2003). Well-aerated tropical soils of the Northern Hemisphere should thus lead to a stronger ¹⁸α<sub>dark_soil</sub>, owing to the greater area of the northern vs. southern tropics. We perform a sensitivity study to assess the effect of soil moisture content on ¹⁸α<sub>dark_soil</sub> by allowing the latter to vary in the intertropical band according to the amount of precipitation simulated by the IPSL model, used as a proxy for soil aeration. Concretely, ¹⁸α<sub>dark_soil</sub> is allowed to vary from its initial value (e.g., –10.2‰ for tropical waterlogged soils) up to a factor of 2 (–20.2‰ for well-aerated tropical soils). Note that not only tropical but also temperate soils are simulated in the 30° S–30° N latitudinal band, as the type of soil is related to the vegetation cover in our model (Sect. 2.4.3).
<figure><img src="image_13.png" /></figure>

**Figure 9.** (a) Zonal annual mean of gross primary productivity expressed in terms of carbon (GPP_C) and oxygen (GPP_O₂) annual molar fluxes for LGMCtrl and HS_exp, respectively. (b) Zonal annual-mean anomalies of δ¹⁸O_lw, δ¹⁸O_p and relative humidity (note its inverted x axis). (c) Zonal annual-mean anomalies for -¹⁸ε resp, δ¹⁸O_lw and δ¹⁸O_terr. Note that all variables of panels (b) and (c) are oxygen-production- (GPP_O₂) weighted.

¹⁸αdark_soil totals -22.36‰ in the modified run, 6.12‰ stronger than in the LGM control run (-16.24‰). δ¹⁸O_terr, increases subsequently by 2.91‰. The picture is similar for HS_exp run, where ¹⁸αdark_soil strengthens by 6.17‰ from -16.06 to -22.23‰, leading δ¹⁸O_terr to increase by 2.96‰ at HS.

While the absolute values of ¹⁸αdark_soil, ¹⁸εresp and δ¹⁸O_terr are significantly modified by accounting for tropical soil aeration, this does not modify the sign of ¹⁸αdark_soil anomaly between the LGM and HS. However, its magnitude is slightly reduced by 0.05‰, which causes the ¹⁸εresp anomaly to vary from -0.03 to +0.02‰. As a result, δ¹⁸O_terr is enriched by 0.16‰ during HS. This 0.05‰ increase in δ¹⁸O_terr anomaly may give an estimate of the magnitude of the uncertainty associated with ¹⁸αdark_soil when considering soil wetness. This sensitivity test does not fundamentally affect the conclusion of the present study, as the ¹⁸εresp anomaly, although becoming positive, remains very small (+0.02‰) but underlines the limitations of our approach. Why does a ¹⁸αdark_soil negative anomaly persist when soil aeration is considered? First, in our sensitivity test, we use the amount of precipitation as an index for soil aeration. This approach may be too simple and could be improved by quantitatively relating the soil aera-

<figure><img src="image_14.png" /></figure>

**Figure 10.** Intertropical profiles of (a) zonal annual mean of gross primary productivity expressed in terms of oxygen (GPP_O₂) annual molar fluxes for LGMCtrl and HS_exp. (b) Zonal annual mean of oxygen-production-weighted δ¹⁸O_p for LGMCtrl and HS_exp. (c) Zonal annual mean of rainfall amount for LGMCtrl and HS_exp.

tion to the model's rainfall amount and land surface slope by using existing parameterizations employed in models that predict nitrous oxide production (e.g., LPX-Bern, which is a state-of-the-art bottom-up dynamic global vegetation and land surface process model; Stocker et al., 2013). Second, soil respiration is closely linked to the vegetation cover in our model, as types of soil (boreal, temperate, tropical) and their associated fractionation factor are related to PFTs rather than land area (Sect. 2.4.3). Third, aboveground oxygen productivity controls the rate of soil respiration. As Southern Hemisphere tropics dominate over Northern Hemisphere tropics in terms of GPP_O₂ for both LGM and HS runs, ¹⁸αdark_soil anomaly remains negative even when considering soil aeration. The distribution of vegetation simulated by ORCHIDEE, favoring Southern Hemisphere tropics, combined with the GPP_O₂ weighting of soil respiration, explains why the drying-out of the Northern Hemisphere soils, albeit covering a greater land area, does not lead to a stronger ¹⁸αdark_soil during HS in our model.

## 4 Discussion

Our results suggest a strong control of tropical hydrology on δ¹⁸O_terr through changes in δ¹⁸O_p. It suggests that δ¹⁸O_atm is related to tropical hydrology and may be a good tracer for global monsoon signal. The aim of the following discussion is to evaluate these results by (i) providing some insights on δ¹⁸O.mar estimate and (ii) testing the robustness of our con-
clusion on the driver of δ¹⁸Oₜerf changes through three sen-
sitivity experiments separating the different parameters (hy-
drology, climate and vegetation).

4.1 Estimate of δ¹⁸Omar over a Heinrich stadial

δ¹⁸Oₜar has been recently estimated as 24.3 ± 2.0‰ for present day (Luz et al., 2014). In order to estimate δ¹⁸Oₜar for the LGM and HS, we assumed that fractionation during oxygen uptake by marine biosphere remained constant between the LGM/HS and present day and used a mean ocean δ¹⁸O enriched by 1‰ at the LGM (Waelbroeck et al., 2002). This results in a value of δ¹⁸Oₜar of 25.3 ± 2.0‰ for the LGM and HS. It is important to note that a rise in sea level during a HS would lead to a depleted mean ocean δ¹⁸O, as polar ice sheets accumulate ¹⁶O, and can therefore not explain the increased δ¹⁸Oₜatm observed during stadials.

Because of the spatial limitation of paleorecords to provide
a global picture of marine primary productivity, we have esti-
mated the marine productivity for the LGM and HS using the
Pelagic Interaction Scheme for Carbon and Ecosystem Stud-
ies (PISCES) model. The PISCES model is a biogeochemi-
cal model of the global ocean including a simple representa-
tion of marine ecosystem and forced offline by the AOGCM
IPSL-CM4 (Aumont and Bopp, 2006; Mariotti et al., 2012).
The model PISCES has already been compared under glacial
conditions with observations (Mariotti et al., 2012; Tagli-
abue et al., 2009; Bopp et al., 2003) and reproduces roughly
the paleoproductivity reconstruction of Kohfeld et al. (2005).
Using the same forcings as for our simulations, Mariotti
et al. (2012) simulate a global decrease in oceanic primary
productivity of 16% during a Heinrich event, in agreement
with independent modeling studies (Schmittner, 2005; Men-
viel et al., 2008; Schmittner and Galbraith, 2008) and more
important than the one in terrestrial GPP (3.5%). Because
δ¹⁸Oₜar is larger than δ¹⁸Oₜerr, this decrease in marine pro-
ductivity would lead to a decrease in global δ¹⁸Oₜatm during
the HS, opposite to the observation. We simulated a change in
the marine production to assess its impact on δ¹⁸Oₜatm signal.
With a change of 10% in marine export, δ¹⁸Oₜatm varies by
0.05‰. However this result needs to be treated with caution
for two main reasons. (i) The fraction of land versus ocean
production is strongly affecting δ¹⁸Oₜatm if the isotope frac-
tionation factors associated with the terrestrial and marine
production are not similar. Pioneer studies on the Dole effect
often invoked the marine-to-terrestrial production ratio to ex-
plain the observed variations. Though it is not the case in
our model, recent studies suggest their magnitude to be very
close, with δ¹⁸Oₜar 1.8‰ higher than δ¹⁸Oₜerr. The change
in δ¹⁸Oₜatm is thus likely to be smaller in the real world. (ii)
Recent studies (eg. Mariotti et al., 2012) rather suggest a de-
crease in marine export after a Heinrich event. This would
lead to a decrease in δ¹⁸Oₜatm signal that is in opposition to
the observations as mentioned in Section 4.1. We conclude

that marine productivity is not the driver for δ¹⁸Oₜatm increase
during HS.

4.2 Disentangling the influences of climate, hydrology
and vegetation on δ¹⁸Oₜatm: sensitivity experiments.

In order to assess the robustness of our conclusion stating that
the low-latitude hydrological cycle is the driver of δ¹⁸Oₜatm
changes, we have run three different experiments:

- In HSclim we test the impact of the climatic conditions. This simulation is similar to LGM-ctrl, except that the temperature and relative humidity from HS_exp are prescribed as boundary conditions.

- In HShydro we test the impact of the hydrological cycle. This simulation is similar to LGM-ctrl, except that δ¹⁸Oₚ and δ¹⁸O_vap from HS_exp are prescribed as boundary conditions.

- In HSveget we test the impact of the vegetation. This simulation is similar to LGM-ctrl, except that the vegetation production and distribution from HS_exp are prescribed as boundary conditions.

Figure 11 summarizes the results of the sensitivity analy-
sis in terms of temperature, humidity, and oxygen isotopic
composition of δ¹⁸Oₚ, δ¹⁸O_Iw and δ¹⁸Oₜerf. The sensitivity
tests show that the implementation of one parameter from
HS_exp (HSclim, HShydro and HSveget) leads to a simu-
lated δ¹⁸Oₜerf anomaly similar to or higher than in the full
HS_exp.

In HSclim, δ¹⁸Oₜerf enrichment is mostly caused by the
0.3% decrease in relative humidity over the LGM veg-
etated areas after the AMOC collapse, since δ¹⁸Oₚ and
δ¹⁸O_Iw are not modified by definition. This global decrease
in GPP_O₂-weighted relative humidity is not visible in the
global HS_exp (Fig. 8) and hence does not explain the
δ¹⁸Oₜerf increase in HS_exp. In HSveget, the southward shift
of HS_exp vegetation leads to a global GPP_O₂ weighted
relative humidity decrease by 0.5%. As in HSclim, this leads
to a high δ¹⁸Oₜerf positive anomaly. The effect of relative
humidity on δ¹⁸Oₜerf is minimized in HS_exp, as the south-
ward shift in vegetation counterbalances the change in cli-
matic conditions. This compensation explains why the final
GPP_O₂ weighted relative humidity does not vary.

In HShydro, the δ¹⁸Oₚ increase explains the whole
δ¹⁸Oₜerf increase. This increase in GPP_O₂-weighted δ¹⁸Oₚ
is similar in HShydro and HS_exp, which confirms that
the δ¹⁸Oₜerf simulated by the model in the full experiment
HS_exp is arising from a change intrinsic to the hydrological
cycle, only slightly affected by vegetation distribution.

In addition, note that a weakening of δ¹⁸O_Iw only oc-
curs with HS vegetation (HS_exp, HSveget), regardless of
the climatic conditions. Moreover, HSclim depicts a stronger
δ¹⁸O_Iw, caused by a higher temperature increase over HS,
<figure><img src="image_15.png" /></figure>

**Figure 11.** Summary of sensitivity study experiments. HSfull uses the same settings as HS_exp, while HSclim, HShydro and HSveget experiments are identical to LGMCtrl, except for meteoric water isotopic composition, climatic conditions, or vegetation production and distribution, respectively, originating from HS_exp. For each of the experiments, an annual-mean anomaly (experiment – LGMCtrl) of oxygen-production-weighted temperature, relative humidity, δ¹⁸O~p~, water vapor δ¹⁸O, δ¹⁸O~lw~, ¹⁸ε~resp~ and δ¹⁸O~terr~ is presented.

leading to an enhanced photorespiration fraction among C3 plants. However, the temperature effect on δ¹⁸O~lw~ remains minor.

The global impact of ¹⁸ε resp is negligible in all experiments, with variations 1 order of magnitude lower than δ¹⁸O p.

Finally, these tests confirm the strong control of hydrological processes on δ¹⁸O terr, and highlight the role of the vegetation distribution in defining δ¹⁸O p and climatic conditions recorded by δ¹⁸O terr.

**5 Conclusions and perspective**

Our study first aimed at quantitatively testing the driving of δ¹⁸O_atm by tropical hydrology as suggested by the strong correlation between local records of δ¹⁸O_c and global record of δ¹⁸O_atm on the millennial scale. For this we used a HS-type simulation under LGM background conditions with an oxygen isotope mass balance model using spatial and temporal fields of (i) temperature and relative humidity from the AOGCM IPSL-CM4, (ii) PFT distribution and GPP provided by the dynamic global vegetation model ORCHIDEE, (iii) oxygen isotope composition of water vapor and precipitation from the AOGCM LMDZ-iso, and (iv) the latest isotope fractionation factor measurements involved in respiratory and photosynthetic processes.

Validation of AOGCM outputs feeding the oxygen isotope mass balance model was performed through a model-data comparison of the main drivers of δ¹⁸O_terr: (i) simulated

δ¹⁸O~p~ was compared to speleothem's calcite δ¹⁸O anomalies, and in most sites showed excellent agreement despite the complexity of the δ¹⁸O~c~ signal; (ii) simulated HS humidity was compared to the reconstructed ones, broadly agreeing with paleodata; and (iii) simulated vegetation was compared with palynological reconstructions for LGM and HS, and was qualitatively consistent.

The model simulates a terrestrial enrichment of δ¹⁸O_terr of 0.11‰, which mostly arises from the δ¹⁸O_p signal. On a global scale, respiration fractionation only plays a minor role in the anomaly observed during HS, and slightly decreases δ¹⁸O_terr in our simulation, driven by a weaker isotope fractionation of soil respiration during HS that masks the effect of increased photorespiration. Accounting for the effect of soil aeration on ¹⁸α-dark_soil modifies the ¹⁸ε resp anomaly by 0.05‰, a change too small to challenge δ¹⁸O_p main control on δ¹⁸O_terr. However, the simplicity of our approach to model soil moisture content and the GPP_O₂ weighting of ¹⁸α-dark_soil favoring Southern Hemisphere tropics are likely to explain this result. In order to precisely evaluate the uncertainty associated with ¹⁸ε resp, an improved parameterization of soil aeration will be needed. This limitation of the model must thus be kept in mind when considering the conclusion of this study.

The strong control of the low-latitude hydrological cycle on simulated δ¹⁸O_atm on millennial timescales suggests that δ¹⁸O_atm records, at first order, changes in monsoonal activity on millennial timescales, in agreement with CH₄ mixing ratio variations. Indeed, rapid CH₄ variations during the last glacial period are generally attributed to changes in the low-latitude water cycle (Baumgartner et al., 2014; Brook et al., 2000; Chappellaz et al., 1993) driven by latitudinal shifts of the ITCZ and the monsoon systems (e.g., Chiang, 2009). Such a signal can also be used for exploring the low-latitude hydrological cycle characteristic of Heinrich events. Indeed, the recent study of Rhodes et al. (2015) suggests that observed CH₄ spikes in WAIS Divide ice core during the cold phases of HSs represent the hydrological signature of Heinrich events, through activation of Southern Hemisphere wetlands. The Guillevic et al. (2014) multi-ice-core proxies approach over GS9-HS4 also suggests a decoupling between changes in Greenland temperatures and low-latitude hydrology identified in both CH₄ and δ¹⁸O_atm, and demonstrates the need for high-resolution data with common precise chronology to explore submillennial variations.

δ¹⁸O_atm is a valuable tool to assess the validity of Earth system model simulations, as it integrates a combination of hydrological, climatic and biological processes. Furthermore, δ¹⁸O_atm is a global signal, which mostly arises from the tropics and integrates all vegetated areas. Therefore, the ability of a model to catch δ¹⁸O_atm millennial-scale variations implies a correct spatial representation of an ensemble of processes. The comparison of Earth system model outputs with global proxies such as δ¹⁸O_atm, involving the main components of the climatic system, is crucial for gaining
confidence in their ability to represent the real world. Our approach is mainly restricted to terrestrial contribution, but future modeling exercises should also include the oceanic δ¹⁸O_atm signal.

**Acknowledgements.** Thanks go to the Marie Curie Initial Training Network INTRAMIF (FP7), which has funded C. Reutenauer's PhD at CIC; to the CIC, for their full support; and to the ANR EL PASO (no. 2010 BLANC 608 01). The climate model simulations were run on CEA supercomputers, and we thank GENCI and TECC for computer time and storage. We are also very grateful for the detailed and constructive comments of Jeff Severinghaus and an anonymous reviewer, which greatly helped in improving the manuscript.

Edited by: E. Wolff

## References

Ahn, J. and Brook, E. J.: Siple Dome ice reveals two modes of millennial CO₂ change during the last ice age, Nat. Commun., 5, 3723, doi:10.1038/ncomms4723, 2014.

Ahn, J., Brook, E. J., Schmittner, A., and Kreutz, K.: Abrupt change in atmospheric CO₂ during the last ice age, Geophys. Res. Lett., 39, 18, L18711, doi:10.1029/2012GL053018, 2012.

Allison, G., Gat, J., and Leaney, F.: The relationship between deuterium and oxygen-18 delta values in leaf water, Chem. Geol., 58, 145–156, doi:10.1016/0168-9622(85)90035-1, 1985.

Alvarez-Solas, J., Robinson, A., Montoya, M., and Ritz, C: Iceberg discharges of the last glacial period driven by oceanic circulation changes, P. Natl. Acad. Sci. USA, 110, 16350–16354, doi:10.1073/pnas.1306622110, 2013.

Angert, A., Barkan, E., Barnett, B., Brugnoli, E., Davidson, E. A., Fessenden, J., Manepong, S., Panapitukkul, N., Rander son, J. T., Savage, K., Yakir, D., and Luz, B.: Contribution of soil respiration in tropical, temperate, and boreal forests to the ¹⁸O enrichment of atmospheric O₂, Global Biogeochem. Cy., 17, 1089, doi:10.1029/2003GB002056, 2003.

Asmerom, Y., Polyak, V. J., and Burns, S. J.: Variable winter moisture in the southwestern United States linked to rapid glacial climate shifts, Nat. Geosci., 3, 114–117, doi:10.1038/NGEO754, 2010.

Aumont, O. and Bopp, L.: Globalizing results from ocean in situ iron fertilization studies, Global Biogeochem. Cy., 20, GB2017, doi:10.1029/2005GB002591, 2006.

Badger, M. R., von Caemmerer, S., Ruuska, S., and Nakano, H.: Electron flow to oxygen in higher plants and algae: rates and control of direct photoreduction (Mehler reaction) and rubisco oxygenase, Philos. T. Roy. Soc. B, 355, 1433–1446, doi:10.1098/rstb.2000.0704, 2000.

Bariac, T., Rambal, S., Jusserand, C., and Berger, A.: Evaluating water fluxes of field-grown alfalfa from diurnal observations of natural isotope concentrations, energy budget and physiological parameters, Agr. Forest Meteorol., 48, 263–283, doi:10.1016/0168-1923(89)90073-7, 1989.

Baumgartner, M., Kindler, P., Eicher, O., Floch, G., Schilt, A., Schwandt, J., Spahni, R., Capron, E., Chappellaz, J., Leuen-

berger, M., Fischer, H., and Stocker, T. F.: NGRIP CH₄ concentration from 120 to 10 kyr before present and its relation to a δ¹⁵N temperature reconstruction from the same ice core, Clim. Past, 10, 903–920, doi:10.5194/cp-10-903-2014, 2014.

Barker, S., Diz, P., Vautravers, M. J., Pike, J., Knorr, G., Hall, I. R., and Broecker, W. S.: Interhemispheric Atlantic seesaw response during the last deglaciation, Nature, 457, 1097–1102, doi:10.1038/nature07770, 2009.

Bender, M. L.: Orbital tuning chronology for the Vostok climate record supported by trapped gas composition, Earth Planet. Sc. Lett., 204, 275–289, doi:10.1016/S0012-821X(02)00980-9, 2002.

Bender, M. L., Sowers, T., and Labeyrie, L.: The Dole Effect and its variations during the last 130,000 years as measured in the Vostok Ice Core, Global Biogeochem. Cy., 8, 363–376, doi:10.1029/94GB00724, 1994.

Blunier, T., Barnett, B., Bender, M. L., and Hendricks, M. B.: Biological oxygen productivity during the last 60,000 years from triple oxygen isotope measurements, Global Biogeochem. Cy., 16, 3–1–3–13, doi:10.1029/2001GB001460, 2002.

Blunier, T., Bender, M. L., Barnett, B., and von Fischer, J. C.: Planetary fertility during the past 400 ka based on the triple isotope composition of O₂ in trapped gases from the Vostok ice core, Clim. Past, 8, 1509–1526, doi:10.5194/cp-8-1509-2012, 2012.

Bond, G., Broecker, W., Johnsen, S., McManus, J., Labeyrie, L., Jouzel, J., and Bonani, G.: Correlations between climate records from North Atlantic sediments and Greenland ice, Nature, 365, 143–147, 1993.

Bony, S. and Emanuel, K. A.: A Parameterization of the Cloudiness Associated with Cumulus Convection; Evaluation Using TOGA COARE Data., J. Atm. Sci., 58, 3158–3183, doi:10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2, 2001.

Bopp, L., Kohfeld, K. E., Le Quere, C., and Aumont, O.: Dust impact on marine biota and atmospheric CO₂ during glacial periods, Paleooceanography, 18, 1046, doi:10.1029/2002PA000810, 2003.

Broccoli, A. J., Dahl, K. A., and Stouffer, R. J.: Response of the ITCZ to Northern Hemisphere cooling, Geophys. Res. Lett., 33, L01702, doi:10.1029/2005GL024546, 2006.

Broecker, W. S.: Abrupt climate change: causal constraints provided by the paleoclimate record, Earth-Sci. Rev., 51, 137–154, 2000.

Broecker, W. S., Bond, G., Klas, M., Bonani, G., and Wolfli, W.: A salt oscillator in the glacial Atlantic? 1. The concept, Paleooceanography, 5, 469–477, doi:10.1029/PA005i004p00469, 1990.

Brook, E. J., Harder, S., Severinghaus, J., Steig, E. J., and Sucher, C. M.: On the origin and timing of rapid changes in atmospheric methane during the Last Glacial Period, Global Biogeochem. Cy., 14, 559–572, doi:10.1029/1999GB001182, 2000.

Brook, E. J., White, J. W., Schilla, A. S., Bender, M. L., Barnett, B., Severinghaus, J. P., Taylor, K. C., Alley, R. B., and Steig, E. J.: Timing of millennial-scale climate change at Siple Dome, West Antarctica, during the last glacial period, Quaternary Sci. Rev., 24, 1333–1343, doi:10.1016/j.quascirev.2005.02.002, 2005.

Caley, T., Roche, D. M., Waelbroeck, C., and Michel, E.: Oxygen stable isotopes during the Last Glacial Maximum climate: perspectives from data–model (iLOVECLIM) comparison, Clim. Past, 10, 1939–1955, doi:10.5194/cp-10-1939-2014, 2014.
Cappa, C. D., Hendricks, M. B., DePaolo, D. J., and Cohen, R. C.: Isotopic fractionation of water during evaporation, J. Geophys. Res.-Atmos., 108, D16, doi:10.1029/2003JD003597, 2003.

Chappellaz, J., Blunier, T., Raynaud, D., Barnola, J., Schwander, J., and Stauffer, B.: Synchronous changes in atmospheric CH₄ and Greenland Climate between 40-kyr and 8-kyr BP, Nature, 366, 443–445, doi:10.1038/366443a0, 1993.

Chappellaz, J., Stowasser, C., Blunier, T., Baslev-Clausen, D., Brook, E. J., Dallmayr, R., Faïn, X., Lee, J. E., Mitchell, L. E., Pascual, O., Romanini, D., Rosen, J., and Schüpbach, S.: High-resolution glacial and deglacial record of atmospheric methane by continuous-flow and laser spectrometer analysis along the NEEM ice core, Clim. Past, 9, 2579–2593, doi:10.5194/cp-9-2579-2013, 2013.

Cheng, H., Sinha, A., Wang, X., Cruz, F., and Edwards, R.: The global paleomonsoon as seen through speleothem records from Asia and the Americas, Clim. Dynam., 39, 1045–1062, doi:10.1007/s00382-012-1363-7, 2012.

Chiang, J. C.: The tropics in paleoclimate, Annu. Rev. Earth Pl. Sc., 37, 263–297, doi:10.1146/annurev.earth.031208.100217, 2009.

Ciais, P., Tagliabue, A., Cuntz, M., Bopp, L., Scholze, M., Hoffmann, G., Lourantou, A., Harrison, S., Prentice, I., Kelley, D., Koven, C. and Piao, S. L.: Large inert carbon pool in the terrestrial biosphere during the Last Glacial Maximum, Nat. Geosci., 5, 74–79, doi:10.1038/ngeo1324, 2012.

Clement, A. C. and Peterson, L. C.: Mechanisms of abrupt climate change of the last glacial period, Rev. Geophys., 46, RG4002, doi:10.1029/2006RG000204, 2008.

Craig, H. and Gordon, L.: Deuterium and oxygen 18 variations in the ocean and the marine atmosphere, in: Stable Isotopes in Oceanographic Studies and Paleotemperatures, edited by: Tongiorgi, E., 26–30 July 1965, Spoleto, Italy, 9–130, 1965.

Cvijanovic, I. and Chiang, J.: Global energy budget changes to high latitude North Atlantic cooling and the tropical ITCZ response, Clim. Dynam., 40, 1435–1452, doi:10.1007/s00382-012-1482-1, 2013.

Dahl, K., Broccoli, A., and Stouffer, R.: Assessing the role of North Atlantic freshwater forcing in millennial scale climate variability: a tropical Atlantic perspective, Clim. Dynam., 24, 325–346, doi:10.1007/s00382-004-0499-5, 2005.

Dällenbach, A., Blunier, T., Flückiger, J., Stauffer, B., Chappellaz, J., and Raynaud, D.: Changes in the atmospheric CH₄ gradient between Greenland and Antarctica during the last glacial and the transition to the Holocene, Geophys. Res. Lett., 27, 1005–1008, doi:10.1029/1999GL010873, 2000.

Deplazes, G., Lückge, A., Stuut, J.-B. W., Pätzold, J., Kuhlmann, H., Husson, D., Fant, M., and Haug, G. H.: Weakening and strengthening of the Indian monsoon during Heinrich events and Dansgaard-Oeschger oscillations, Paleooceanography, 29, 99–114, doi:10.1002/2013PA002509, 2014.

Dokken, T. M., Nisancioglu, K. H., Li, C., Battisti, D. S., and Kissel, C.: Dansgaard-Oeschger cycles: interactions between ocean and sea ice intrinsic to the Nordic seas, Paleooceanography, 28, 491–502, doi:10.1002/palo.20042, 2013.

Dongmann, G.: The contribution of land photosynthesis to the stationary enrichment of ¹⁸O in the atmosphere, Radiat. Environ. Bioph., 11, 219–225, doi:10.1007/BF01323191, 1974.

Dreyfus, G. B., Parrenin, F., Lemieux-Dudon, B., Durand, G., Masson-Delmotte, V., Jouzel, J., Barnola, J.-M., Panno, L., Spahni, R., Tisserand, A., Siegenthaler, U., and Leuenberger, M.: Anomalous flow below 2700 m in the EPICA Dome C ice core detected using δ¹⁸O of atmospheric oxygen measurements, Clim. Past, 3, 341–353, doi:10.5194/cp-3-341-2007, 2007.

Eisenstadt, D., Barkan, E., Luz, B., and Kaplan, A.: Enrichment of oxygen heavy isotopes during photosynthesis in phytoplankton, Photosynth. Res., 103, 97–103, doi:10.1007/s11120-009-9518-z, 2010.

Elliot, M., Labeyrie, L., and Duplessy, J.-C.: Changes in North Atlantic deep-water formation associated with the Dansgaard-Oeschger temperature oscillations (60–10 ka), Quaternary Sci. Rev., 21, 1153–1165, doi:10.1016/S0277-3791(01)00137-8, 2002.

Emanuel, K. A.: A scheme for representing cumulus convection in large-scale models, J. Atm. Sci., 48, 2313–2329, doi:10.1175/1520-0469(1991)048<2313:ASFRCC>2.0.CO;2, 1991.

Farquhar, G. D. and Lloyd, J.: 5 – Carbon and oxygen isotope effects in the exchange of carbon dioxide between terrestrial plants and the atmosphere, in: Stable Isotopes and Plant Carbon-Water Relations, edited by: Farquhar, G. D., Ehleringer, J. R., Hall, A. E., Academic Press, San Diego, 47–70, doi:10.1016/B978-0-08-091801-3.50011-8, 1993.

Farquhar, G. D., von Caemmerer, S., and Berry, J.: A biochemical model of photosynthetic CO₂ assimilation in leaves of C3 species, Planta, 149, 78–90, 1980.

Farquhar, G. D., Hubick, K., Condon, A., and Richards, R.: Carbon isotope fractionation and plant water-use efficiency, in: Stable Isotopes in Ecological Research, Springer, 21–40, 1989.

Farquhar, G. D., Cernusak, L. A., and Barnes, B.: Heavy water fractionation during transpiration, Plant Physiol., 143, 11–18, 2007.

Flanagan, L. B., Bain, J. F., and Ehleringer, J. R.: Stable oxygen and hydrogen isotope composition of leaf water in C3 and C4 plant species under field conditions, Oecologia, 88, 394–400, 1991a.

Flanagan, L. B., Comstock, J. P., and Ehleringer, J. R.: Comparison of modeled and observed environmental influences on the stable oxygen and hydrogen isotope composition of leaf water in Phaseolus vulgaris L., Plant Physiol., 96, 588–596, 1991b.

Flanagan, L. B., Marshall, J. D., and Ehleringer, J. R.: Photosynthetic gas exchange and the stable isotope composition of leaf water: comparison of a xylem-tapping mistletoe and its host, Plant Cell Environ., 16, 623–631, doi:10.1111/j.1365-3040.1993.tb00480.x, 1993.

Flanagan, L. B., Phillips, S. L., Ehleringer, J. R., Lloyd, J., and Farquhar, G. D.: Effect of changes in leaf water oxygen isotopic composition on discrimination against C¹⁸O¹⁶O during photosynthetic gas exchange, Funct. Plant Biol., 21, 221–234, 1994.

Fleitmann, D., Cheng, H., Badertscher, S., Edwards, R. L., Mudelsee, M., Gktürk, O. M., Fankhauser, A., Pickering, R., Raible, C. C., Matter, A., Kramers, J., and Tüysüz, O.: Timing and climatic impact of Greenland interstadials recorded in st mee from northern Turkey, Geophys. Res. Lett., 36, L19707, doi:10.1029/2009GL040050, 2009.

Fletcher, W. J., Goni, M. F. S., Allen, J. R., Cheddadi, R., Combourieu-Nebout, N., Huntley, B., Lawson, I., Londeix, L., Magri, D., Margari, V., Müller, U. C., Naughton, F., Novenko, E., Roucoux, K., and Tzedakis, P.: Millennial-scale variability dur-
ing the last glacial in vegetation records from Europe, Quaternary Sci. Rev., 29, 2839–2864, doi:10.1016/j.quascirev.2009.11.015, 2010.

Flückiger, J., Dällenbach, A., Blunier, T., Stauffer, B., Stocker, T. F., Raynaud, D., and Barnola, J.-M.: Variations in atmospheric N₂O concentration during abrupt climatic changes, Science, 285, 227–230, doi:10.1126/science.285.5425.227, 1999.

Genty, D., Blamart, D., Ouahdi, R., Gilmour, M., Baker, A., Jouzel, J., and Van-Exter, S.: Precise dating of Dansgaard-Oeschger climate oscillations in western Europe from stalagmite data, Nature, 421, 833–837, 2003.

Gillon, J. and Yakir, D.: Influence of carbonic anhydrase activity in terrestrial vegetation on the ¹⁸O content of atmospheric CO₂, Science, 291, 2584–2587, doi:10.1126/science.1056374, 2001.

Grandpeix, J.-Y., Phillips, V., and Tailleux, R.: Improved mixing representation in Emanuel's convection scheme, Q. J. Roy. Meteorol. Soc., 130, 3207–3222, doi:10.1256/qj.03.144, 2004.

Grootes, P. M. and Stuiver, M.: Oxygen 18/16 variability in Greenland snow and ice with 10³ to 10⁵ year time resolution, J. Geophys. Res.-Oceans, 102, 26455–26470, doi:10.1029/97JC00880, 1997.

Grousset, F., Labeyrie, L., Sinko, J., Cremer, M., Bond, G., Duprat, J., Cortijo, E., and Huon, S.: Patterns of ice-rafted detritus in the glacial North Atlantic (40–55 N), Paleooceanography, 8, 175–192, 1993.

Guillevic, M.: Characterisation of rapid climate changes through isotope analyses of ice and entrapped air in the NEEM ice core, PhD thesis, Versailles-St Quentin en Yvelines, 180 pp., 2013.

Guillevic, M., Bazin, L., Landais, A., Stowasser, C., Masson-Delmotte, V., Blunier, T., Eynaud, F., Falourd, S., Michel, E., Minster, B., Popp, T., Prié, F., and Vinther, B. M.: Evidence for a three-phase sequence during Heinrich Stadial 4 using a multi-proxy approach based on Greenland ice core records, Clim. Past, 10, 2115–2133, doi:10.5194/cp-10-2115-2014, 2014.

Guy, R. D., Fogel, M. L., and Berry, J. A.: Photosynthetic Fractionation of the Stable Isotopes of Oxygen and Carbon, Plant Physiol., 101, 37–47, doi:10.1104/pp.101.1.37, 1993.

Handiani, D., Paul, A., and Dupont, L.: Tropical climate and vegetation changes during Heinrich Event 1: a model-data comparison, Clim. Past, 8, 37–57, doi:10.5194/cp-8-37-2012, 2012.

Harrison, S. and Goni, M. S.: Global patterns of vegetation response to millennial-scale variability and rapid climate change during the last glacial period, Quaternary Sci. Rev., 29, 2957–2980, doi:10.1016/j.quascirev.2010.07.016, 2010.

Heinrich, H.: Origin and consequences of cyclic ice rafting in the northeast Atlantic Ocean during the past 130,000 years, Quaternary Res., 29, 142–152, 1988.

Helman, Y., Barkan, E., Eisenstadt, D., Luz, B., and Kapan, A.: Fractionation of the three stable oxygen isotopes by oxygen-producing and oxygen-consuming reactions in photosynthetic organisms, Plant Physiol., 138, 2292–2298, doi:10.1104/pp.105.063768, 2005.

Hendricks, M. B., Bender, M. L., and Barnett, B. A.: Net and gross O₂ production in the southern ocean from measurements of biological O₂ saturation and its triple isotope composition, Deep-Sea Res. Pt. I, 51, 1541–1561, doi:10.1016/j.dsr.2004.06.006, 2004.

Hessler, I., Dupont, L., Bonnefille, R., Behling, H., Gonzalez, C., Helmens, K. F., Hooghiemstra, H., Lebamba, J., Ledru, M.-P.,

Lezine, A.-M., Maley, J., Marret, F., and Vincens, A.: Millennial-scale changes in vegetation records from tropical Africa and South America during the last glacial, Quaternary Sci. Rev., 29, 2882–2899, doi:10.1016/j.quascirev.2009.11.029, 2010.

Hoffmann, G., Cuntz, M., Weber, C., Ciais, P., Friedlingstein, P., Heimann, M., Jouzel, J., Kaduk, J., Maier-Reimer, E., Seibt, U., and Six, K.: A model of the Earth's Dole effect, Global Biogeochem. Cy., 18, GB1008, doi:10.1029/2003GB002059, 2004.

Hourdin, F., Musat, I., Bony, S., Braconnot, P., Codron, F., Dufresne, J.-L., Fairhead, L., Filiberti, M.-A., Friedlingstein, P., Grandpeix, J.-Y., Krinner, G., LeVan, P., Li, Z.-X., and Lott, F.: The LMDZ4 general circulation model: climate performance and sensitivity to parametrized physics with emphasis on tropical convection, Clim. Dynam., 27, 787–813, doi:10.1007/s00382-006-0158-0, 2006.

Huber, C., Leuenberger, M., Spahni, R., Flückiger, J., Schwander, J., Stocker, T., Johnsen, S., Landais, A., and Jouzel, J.: Isotope calibrated Greenland temperature record over marine isotope stage 3 and its relation to CH₄, Earth Planet. Sc. Lett., 243, 504–519, 2006.

Jimenez-Moreno, G., Anderson, R. S., Desprat, S., Grigg, L. D., Grimm, E. C., Heusser, L. E., Jacobs, B. F., Lopez-Martinez, C., Whitlock, C. L., and Willard, D. A.: Millennial-scale variability during the last glacial in vegetation records from North America, Quaternary Sci. Rev., 29, 2865–2881, doi:10.1016/j.quascirev.2009.12.013, 2010.

Joos, F., Gerber, S., Prentice, I. C., Otto-Bliesner, B. L., and Valdes, P. J.: Transient simulations of Holocene atmospheric carbon dioxide and terrestrial carbon since the Last Glacial Maximum, Global Biogeochem. Cy., 18, GB2002, doi:10.1029/2003GB002156, 2004.

Kageyama, M., Mignot, J., Swingedouw, D., Marzin, C., Alkama, R., and Marti, O.: Glacial climate sensitivity to different states of the Atlantic meridional overturning circulation: results from the IPSL model, Clim. Past, 5, 551–570, doi:10.5194/cp-5-551-2009, 2009.

Kageyama, M., Paul, A., Roche, D. M., and Meerbeeck, C. J. V.: Modelling glacial climatic millennial-scale variability related to changes in the Atlantic meridional overturning circulation: a review, Quaternary Sci. Rev., 29, 2931–2956, doi:10.1016/j.quascirev.2010.05.029, 2010.

Kageyama, M., Merkel, U., Otto-Bliesner, B., Prange, M., Abe-Ouchi, A., Lohmann, G., Ohgaito, R., Roche, D. M., Sin-garayer, J., Swingedouw, D., and X Zhang: Climatic impacts of fresh water hosing under Last Glacial Maximum conditions: a multi-model study, Clim. Past, 9, 935–953, doi:10.5194/cp-9-935-2013, 2013.

Kanner, L. C., Burns, S. J., Cheng, H., and Edwards, R. L.: High-latitude forcing of the South American summer monsoon during the last glacial, Science, 335, 570–573, doi:10.1126/science.1213397, 2012.

Kanner, L. C., Buenning, N. H., Stott, L. D., Timmermann, A., and Noone, D.: The role of soil processes in ¹⁸O terrestrial climate proxies, Global Biogeochem. Cy., 28, 239–252, doi:10.1002/2013GB004742, 2014.

Kawamura, K., Parrenin, F., Lisiecki, L., Uemura, R., Vimeux, F., Severinghaus, J. P., Hutterli, M. A., Nakazawa, T., Aoki, S., Jouzel, J., Raymo, M. E., Matsumoto, K., Nakata, H., Motoyama, H., Fujita, S., Goto-Azuma, K., Fujii, Y., and Watan-
abe, O.: Northern hemisphere forcing of climatic cycles in Antarctica over the past 360,000 years, Nature, 448, 912–U4, doi:10.1038/nature06015, 2007.

Keeling, R. F.: Development of an interferometric oxygen analyzer for precise measurement of the atmospheric O₂ mole fraction, PhD thesis, Harvard University, 190 pp., 1988.

Kindler, P., Guillevic, M., Baumgartner, M., Schwander, J., Landais, A., and Leuenberger, M.: Temperature reconstruction from 10 to 120 kyr b2k from the NGRIP ice core, Clim. Past, 10, 887–902, doi:10.5194/cp-10-887-2014, 2014.

Kleppin, H., Jochum, M., Otto-Bliesner, B., Shields, C. A., and Yeager, S.: Stochastic atmospheric forcing as trigger for sudden Greenland warmings, J. Climate, 28, 19, doi:10.1175/JCLI-D-14-00728.1, 2015.

Kohfeld, K. E., Quere, C. L., Harrison, S. P., and Anderson, R. F.: Role of marine biology in glacial-interglacial CO₂ cycles, Science, 308, 74–78, doi:10.1126/science.1105375, 2005.

Krebs, U. and Timmermann, A.: Tropical air-sea interactions accelerate the recovery of the Atlantic meridional overturning circulation after a major shutdown, J. Climate, 20, 4940–4956, 2007.

Krinner, G., Viovy, N., de Noblet-Ducoudre, N., Ogee, J., Polcher, J., Friedlingstein, P., Ciais, P., Sitch, S., and Pren'tice, I. C.: A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system, Global Biogeochem. Cy., 19, 1, doi:10.1029/2003GB002199, 2005.

Landais, A., Barnola, J. M., Masson-Delmotte, V., Jouzel, J., Chappellaz, J., Caillon, N., Huber, C., Leuenberger, M., and Johnsen, S. J.: A continuous record of temperature evolution over a sequence of Dansgaard-Oeschger events during Marine Isotopic Stage 4 (76 to 62 kyr BP), Geophys. Res. Lett., 31, L22211, doi:10.1029/2004GL021193, 2004.

Landais, A., Lathiere, J., Barkan, E., and Luz, B.: Reconsidering the change in global biosphere productivity between the Last Glacial Maximum and present day from the triple oxygen isotopic composition of air trapped in ice cores, Global Biogeochem. Cy., 21, GB1025, doi:10.1029/2006GB002739, 2007a.

Landais, A., Masson-Delmotte, V., Nebout, N. C., Jouzel, J., Blunier, T., Leuenberger, M., Dahl-Jensen, D., and Johnsen, S.: Milenial scale variations of the isotopic composition of atmospheric oxygen over marine isotopic stage 4, Earth Planet. Sc. Lett., 258, 101–113, doi:10.1016/j.epsl.2007.03.027, 2007b.

Landais, A., Dreyfus, G., Capron, E., Masson-Delmotte, V., Sanchez-Goñi, M. F., Desprat, S., Hoffmann, G., Jouzel, J., Leuenberger, M., and Johnsen, S.: What drives the millennial and orbital variations of δ¹⁸O_atm?, Quaternary Sci. Rev., 29, 235–246, doi:10.1016/j.quascirev.2009.07.005, 2010.

Landais, A., Dreyfus, G., Capron, E., Pol, K., Loutre, M. F., Raynaud, D., Lipenkov, V. Y., Arnaud, L., Masson-Delmotte, V., Paillard, D., Jouzel, J., and Leuenberger, M.: Towards orbital dating of the EPICA Dome C ice core using δO₂/N₂, Clim. Past, 8, 191–203, doi:10.5194/cp-8-191-2012, 2012.

LeGrande, A. N. and Schmidt, G. A.: Ensemble, water isotope-enabled, coupled general circulation modeling insights into the 8.2 ka event, Paleooceanography, 23, PA3207, doi:10.1029/2008PA001610, 2008.

Lemieux-Dudon, B., Bazin, L., Landais, A., Toyé Mahamadou Kele, H., Guillevic, M., Kindler, P., Parrenin, F., and Martinerie, P.: Implementation of counted layers for coherent ice core chronology, Clim. Past, 11, 959–978, doi:10.5194/cp-11-959-2015, 2015.

Lewis, S. C., LeGrande, A. N., Kelley, M., and Schmidt, G. A.: Water vapour source impacts on oxygen isotope variability in tropical precipitation during Heinrich events, Clim. Past, 6, 325–343, doi:10.5194/cp-6-325-2010, 2010.

Li, C., Battisti, D. S., Schrag, D. P., and Tziperman, E.: Abrupt climate shifts in Greenland due to displacements of the sea ice edge, Geophys. Res. Lett., 32, L19702, 10.1029/2005GL023492, 2005.

Li, C., Battisti, D. S., and Bitz, C. M.: Can North Atlantic sea ice anomalies account for Dansgaard-Oeschger climate signals?, J. Climate, 23, 5457–5475, 10.1175/2010JCLI3409.1, 2010.

Liu, Z., Otto-Bliesner, B. L., He, F., Brady, E. C., Tomas, R., Clark, P. U., Carlson, A. E., Lynch-Stieglitz, J., Curry, W., Brook, E., Erickson, D., Jacob, R., Kutzbach, J., and Cheng, J.: Transient Simulation of Last Deglaciation with a New Mechanism for Bølling-Allerød Warming, Science, 325, 310–314, doi:10.1126/science.1171041, 2009.

Lloyd, J. and Farquhar, G.: ¹³C discrimination during CO₂ assimilation by the terrestrial biosphere, Oecologia, 99, 201–215, doi:10.1007/BF00627732, 1994.

Luz, B. and Barkan, E.: Proper estimation of marine gross O₂ production with ¹⁷O/¹⁶O and ¹⁸O/¹⁶O ratios of dissolved O₂, Geophys. Res. Lett., 38, L19606, 10.1029/2011GL049138, 2011a.

Luz, B. and Barkan, E.: The isotopic composition of atmospheric oxygen, Global Biogeochem. Cy., 25, GB3001, doi:10.1029/2010GB003883, 2011b.

Luz, B., Barkan, E., Bender, M., Thiemens, M., and Boering, K.: Triple-isotope composition of atmospheric oxygen as a tracer of biosphere productivity, Nature, 400, 547–550, doi:10.1038/22987, 1999.

Luz, B., Barkan, E., Yam, R., and Shemesh, A.: Fractionation of oxygen and hydrogen isotopes in evaporating water, Geochim. Cosmochim. Ac., 73, 6697–6703, doi:10.1016/j.gca.2009.08.008, 2009.

Luz, B., Barkan, E., and Severinghaus, J.: 5.14 – The stable isotopic composition of atmospheric O₂, in: Treatise on Geochemistry, edited by: Holland, H. D. and Turekian, K. K., Elsevier, Oxford, 2nd edn., 363–383, doi:10.1016/B978-0-08-095975-7.00419-8, 2014.

Marcott, S. A., Clark, P. U., Padman, L., Klinkhammer, G. P., Springer, S. R., Liu, Z., Otto-Bliesner, B. L., Carlson, A. E., Ungerer, A., Padman, J., He, F., Cheng, J., and Schmittner, A.: Ice shelf collapse from subsurface warming as a trigger for Heinrich events, P. Natl. Acad. Sci. USA, 108, 13415–13419, doi:., 2011.

Majoube, M.: Fractionnement en oxygene-18 et en deuterium entre l'eau et sa vapeur, J. Chim. Phys, 68, 1423–1436, 1971.

Mariotti, V., Bopp, L., Tagliabue, A., Kageyama, M., and Swingedouw, D.: Marine productivity response to Heinrich events: a model-data comparison, Clim. Past, 8, 1581–1598, doi:10.5194/cp-8-1581-2012, 2012.

Marti, O., Braconnot, P., Dufresne, J.-L., Bellier, J., Benshila, R., Bony, S., Brockmann, P., Cadule, P., Caubel, A., Codron, F., de Noblet, N., Denvil, S., Fairhead, L., Fichefet, T., Foujols, M.-A., Friedlingstein, P., Goosse, H., Grandpeix, J.-Y., Guilyardi, E., Hourdin, F., Idelkadi, A., Kageyama, M., Krinner, G., Lžvy, C., Madec, G., Mignot, J., Musat, I., Swingedouw, D., and Ta-landier, C.: Key features of the IPSL ocean atmosphere model
and its sensitivity to atmospheric resolution, Clim. Dynam., 34, 1–26, 2010.

McGuire, K., DeWalle, D., and Gburek, W.: Evaluation of mean residence time in subsurface waters using oxygen-18 fluctuations during drought conditions in the mid-Appalachians, J. Hydrol., 261, 132–149, doi:10.1016/S0022-1694(02)00006-9, 2002.

Menviel, L., Timmermann, A., Mouchet, A., and Timm, O.: Meridional reorganizations of marine and terrestrial productivity during Heinrich events, Paleooceanography, 23, PA1203, doi:10.1029/2007PA001445, 2008.

Menviel, L., Timmermann, A., Friedrich, T., and England, M. H.: Hindcasting the continuum of Dansgaard-Oeschger variability: mechanisms, patterns and timing, Clim. Past, 10, 63–77, doi:10.5194/cp-10-63-2014, 2014.

Merlivat, L.: Molecular diffusivities of H₂¹⁶O, HD¹⁶O, and H₂¹⁸O in gases, J. Chem. Phys., 69, 2864–2871, doi:10.1063/1.436884, 1978.

Mohtadi, M., Oppo, D. W., Steinke, S., Stuut, J.-B. W., De Pol-Holz, R., Hebbeln, D., and Lückge, A.: Glacial to Holocene swings of the Australian-Indonesian monsoon, Nat. Geosci., 4, 540–544, 2011.

Monnin, E., Indermühle, A., Dällenbach, A., Flückiger, J., Stauffer, B., Stocker, T. F., Raynaud, D., and Barnola, J.-M.: Atmospheric CO₂ Concentrations over the Last Glacial Termination, Science, 291, 112–114, doi:10.1126/science.291.5501.112, 2001.

NGRIP members, Andersen, K., Azuma, N., Barnola, J., Bigler, M., Biscaye, P., Caillon, N., Chappellaz, J., Clausen, H., DahlJensen, D., Fischer, H., Fluckiger, J., Fritzsche, D., Fujii, Y., Goto-Azuma, K., Gronvold, K., Gundestrup, N., Hansson, M., Huber, C., Hvidberg, C., Johnsen, S., Jonsell, U., Jouzel, J., Kipfstuhl, S., Landais, A., Leuenberger, M., Lorrain, R., Masson-Delmotte, V., Miller, H., Motoyama, H., Narita, H., Popp, T., Rasmussen, S., Raynaud, D., Rothlisberger, R., Ruth, U., Samyn, D., Schwander, J., Shoji, H., Siggard-Andersen, M., Steffensen, J., Stocker, T., Sveinbjornsdottir, A., Svensson, A., Takata, M., Tison, J., Thorsteinsson, T., Watanabe, O., Wilhelms, F., and White, J. W. C.: High-resolution record of Northern Hemisphere climate extending into the last interglacial period, Nature, 431, 147–151, doi:10.1038/nature02805, 2004.

Otto-Bliesner, B. L. and Brady, E. C.: The sensitivity of the climate response to the magnitude and location of freshwater forcing: last glacial maximum experiments, Quaternary Sci. Rev., 29, 56–73, doi:10.1016/j.quascirev.2009.07.004, 2010.

Pausata, F. S., Battisti, D. S., Nisancioglu, K. H., and Bitz, C. M.: Chinese stalagmite [delta]¹⁸O controlled by changes in the Indian monsoon during a simulated Heinrich event, Nat. Geosci., 4, 474–480, 2011.

Peltier, W.: Global glacial isostasy and the surface of the ice-age Earth: the ICE-5G (VM2) model and GRACE, Annu. Rev. Earth Pl. Sc., 32, 111–149, 2004.

Peterson, L. C. and Haug, G. H.: Variability in the mean latitude of the Atlantic intertropical convergence zone as recorded by riverine input of sediments to the Cariaco Basin (Venezuela), Palaeogeogr. Palaeocl., 234, 97–113, doi:10.1016/j.palaeo.2005.10.021, 2006.

Prentice, I. C., Jolly, D., and Afanaseva, N. B., Ager, T. A., Anderson, K., Anderson, P. M., Andrieu, V., Andreev, A. A., Ballouche, A., Bartlein, P. J., de Beaulieu, J. L., Bengo, M., Berezina, N. A., Bezusko, L. G., Bezusko, T. V., Bigelow, N. H., Blyakharchuk, T. A., Bolikhovskaya, N. S., Bonnefille, R., Bottema, S., Bržnac, P., Brubaker, L. B., Buchet, G., Burney, D., Bykova, G. V., Cheddadi, R., Chen, X., Chernavskaya, M. M., Chernova, G. M., Cwynar, L. C., Dorofeyuk, N. I., Dirksen, V. G., Edorh, T., Edwards, M. E., Eisner, W. R., Elenga, H., Elina, G. A., Elmoutaki, S., Filimonova, L. V., Glebov, F. Z., Guiot, J., Gunova, V. S., Hamilton, A. C., Han, H., Harrison, S. P., Hu, F.-S., Huang, C., Huntley, B., Jolly, D., Jonson, H., Ke, M., Khomutova, V. I., Kong Z., Kvavadze, E. V., Laarif, F., Lamb, H. E., Lžzine, A.-M., Li, S., Li, W., Liew, P., Liu, G., Liu, J., Liu, Q., Liu, K.-B., Lozhkin, A. V., Maley, J., Marchant, R., Mbenza, M., MacDonald, G. M., Miyoshi, N., Mock, C. J., Morita, Y., Newby, P., Ni, J., Osipova, I. R., Panova, N. K., Perez-Obiol, R., Peyron, O., Prennice, I. C., Qiu, W., Reille, M., Ren, G., Reynaud-Farrera, I., Richard, P. J. H., Riotlet, G., Ritchie, J. C., Roche, E., Saarse, L., Scott, L., Sevastyanov, D. V., Sher, A. V., Song, C., Spear, R. W., Ssemmanda, I., Straka, H., Sugita, S., Sun, X., Takahara, H., Tang, L., Tarasov, P. E., Taylor, D., Thompson, R. S., Uchiyama, T., Van Campo, E., Vilimumbalo, S., Vincens, A., Volkova, V. S., Waller, M., Webb, T., Williams, J. W., Xia, Y., Xu, Q., Yan, S., Yang, X., Yu, G., Zernitskaya, V. P., Zhao, J., and Zheng, Z.: Mid-Holocene and glacial-maximum vegetation geography of the northern continents and Africa, J. Biogeogr., 27, 507–519, 2000.

Raich, J. W. and Potter, C. S.: Global patterns of carbon dioxide emissions from soils, Global Biogeochem. Cy., 9, 23–36, doi:10.1029/94GB02723, 1995.

Rasmussen, S. O., Abbott, P. M., Blunier, T., Bourne, A. J., Brook, E., Buchardt, S. L., Buizert, C., Chappellaz, J., Clausen, H. B., Cook, E., Dahl-Jensen, D., Davies, S. M., Guillevic, M., Kipfstuhl, S., Laepple, T., Seierstad, I. K., Severinghaus, J. P., Steffensen, J. P., Stowasser, C., Svensson, A., Vallelonga, P., Vinther, B. M., Wilhelms, F., and Winstrup, M.: A first chronology for the North Greenland Eemian Ice Drilling (NEEM) ice core, Clim. Past, 9, 2713–2730, doi:10.5194/cp-9-2713-2013, 2013.

Rhodes, R. H., Brook, E. J., Chiang, J. C. H., Blunier, T., Maselli, O. J., McConnell, J. R., Romanini, D., and Severinghaus, J. P.: Enhanced tropical methane production in response to iceberg discharge in the North Atlantic, Science, 348, 1016–1019, doi:10.1126/science.1262005, 2015.

Risi, C., Bony, S., Vimeux, F., and Jouzel, J.: Water-stable isotopes in the LMDZ4 general circulation model: model evaluation for present-day and past climates and applications to climatic interpretations of tropical isotopic records, J. Geophys. Res.-Atmos., 115, D12, doi:10.1029/2009JD013255, 2010.

Roche, D., Wiersma, A., and Renssen, H.: A systematic study of the impact of freshwater pulses with respect to different geographical locations, Clim. Dynam., 34, 997–1013, doi:10.1007/s00382-009-0578-8, 2010.

Ruddiman, W. F.: Late Quaternary deposition of ice-rafted sand in the subpolar North Atlantic (lat 40 to 65 N), Geol. Soc. Am. Bull., 88, 1813–1827, 1977.

Sanchez Goni, M. F. and Harrison, S. P.: Millennial-scale climate variability and vegetation changes during the last glacial: concepts and terminology, Quaternary Sci. Rev., 29, 2823–2827, doi:10.1016/j.quascirev.2009.11.014, 2010.
Schilt, A., Brook, E. J., Bauska, T. K., Bagggenstos, D., Fischer, H., Joos, F., Petrenko, V. V., Schaefer, H., Schmitt, J., Severinghaus, J. P., Spahni, R., Stocker, T. F.: Isotopic constraints on marine and terrestrial N2O emissions during the last deglaciation, Nature, 516, 234–237, doi:10.1038/nature13971, 2014.

Schlesinger, W. and Andrews, J.: Soil respiration and the global carbon cycle, Biogeochemistry, 48, 7–20, doi:10.1023/A:1006247623877, 2000.

Schmittner, A.: Decline of the marine ecosystem caused by a reduction in the Atlantic overturning circulation, Nature, 434, 628–633, doi:10.1038/nature03476, 2005.

Schmittner, A. and Galbraith, E. D.: Glacial greenhouse-gas fluctuations controlled by ocean circulation changes, Nature, 456, 373–376, 2008.

Severinghaus, J. P.: Studies of the terrestrial O₂ and carbon cycles in sand dune gases and in biosphere 2, Tech. rep., Oak Ridge Associated Universities, Inc., TN, 159 pp., 1995.

Severinghaus, J. P., Beaudette, R., Headly, M. A., Taylor, K., and Brook, E. J.: Oxygen-18 of O₂ records the impact of abrupt climate change on the terrestrial biosphere, Science, 324, 1431–1434, doi:10.1126/science.1169473, 2009.

Stocker, B. D., Roth, R., Joos, F., Spahni, R., Steinacher, M., Zaehle, S., Bouwman, L., Zu, R. and Prentice, I. C.: Multiple greenhouse-gas feedbacks from the land biosphere under future climate change scenarios, Nature, 3, 666–672, doi:10.1038/nclimate1864, 2013.

Stouffer, R. J., Yin, J., Gregory, J., Dixon, K., Spelman, M., Hurlin, W., Weaver, A., Eby, M., Flato, G., Hasumi, H., Hu, A., Jungclaus, J. H., Kamenkovich, I. V., Levermann, A., Montoya, M., Murakami, S., Nawrath, S., Oka, A., Peltier, W. R., Robitaille, D. Y., Sokolov, A., Vettoretti, G., and Weber, S. L.: Investigating the causes of the response of the thermohaline circulation to past and future climate changes, J. Climate, 19, 1365–1387, 10.1175/JCLI3689.1, 2006.

Swingedouw, D., Mignot, J., Braconnot, P., Mosquet, E., Kageyama, M., and Alkama, R.: Impact of freshwater release in the North Atlantic under different climate conditions in an OAGCM, J. Climate, 22, 6377–6403, 2009.

Tagliabue, A., Bopp, L., and Aumont, O.: Evaluating the importance of atmospheric and sedimentary iron sources to Southern Ocean biogeochemistry, Geophys. Res. Lett., 36, 13, doi:10.1029/2009GL038914, 2009.

Takahara, H., Igarashi, Y., Hayashi, R., Kumon, F., Liew, P.-M., Yamamoto, M., Kawai, S., Oba, T., and Irino, T.: Millennial-scale variability in vegetation records from the East Asian islands: Taiwan, Japan and Sakhalin, Quaternary Sci. Rev., 29, 2900–2917, doi:10.1016/j.quascirev.2009.11.026, 2010.

Von Caemmerer, S.: Biochemical Models of Leaf Photosynthesis, 2, Csiro publishing, 528–528, 2000.

Waelbroeck, C., Labeyrie, L., Michel, E., Duplessy, J., McManus, J., Lambeck, K., Balbon, E., and Labracherie, M.: Sea-level and deep water temperature changes derived from benthic foraminifera isotopic records, Quaternary Sci. Rev., 21, 295–305, doi:10.1016/S0277-3791(01)00101-9, 2002.

Walker, C. D. and Brunel, J.-P.: Examining evapotranspiration in a semi-arid region using stable isotopes of hydrogen and oxygen, J. Hydrol., 118, 55–75, doi:10.1016/0022-1694(90)90250-2, 1990.

Walker, C. D., Leaney, F. W., Dighton, J. C., and Allison, G. B.: The influence of transpiration on the equilibration of leaf water with atmospheric water vapour, Plant Cell Environ., 12, 221–234, doi:10.1111/j.1365-3040.1989.tb01937.x, 1989.

Wang, L., Good, S. P., Caylor, K. K., and Cernusak, L. A.: Direct quantification of leaf transpiration isotopic composition, Agr. Forest Meteorol., 154–155, 127–135, doi:10.1016/j.agrformet.2011.10.018, 2012.

Wang, Y. J., Cheng, H., Edwards, R. L., An, Z. S., Wu, J. Y., Shen, C.-C., and Dorale, J. A.: A high-resolution absolute-dated late pleistocene monsoon record from Hulu Cave, China, Science, 294, 2345–2348, doi:10.1126/science.1064618, 2001.

Wang, Y. J., Cheng, H., Edwards, R. L., He, Y., Kong, X., An, Z., Wu, J., Kelly, M. J., Dykoski, C. A., and Li, X.: The holocene asian monsoon: links to solar changes and North Atlantic climate, Science, 308, 854–857, doi:10.1126/science.1106296, 2005.

Wang, Y. J., Cheng, H., Edwards, R. L., Kong, X., Shao, X., Chen, S., Wu, J., Jiang, X., Wang, X., and An, Z.: Millennial- and orbital-scale changes in the East Asian monsoon over the past 224,000 years, Nature, 451, 1090–1093, doi:10.1038/nature06692, 2008.

Weldeab, S.: Bipolar modulation of millennial-scale West African monsoon variability during the last glacial (75,000–25,000 years ago), Quaternary Sci. Rev., 40, 21–29, doi:10.1016/j.quascirev.2012.02.014, 2012.

Welp, L. R., Keeling, R. F., Meijer, H., Bollenbacher, A., F., Piper, S. C., Yoshimura, K., Francey, R. J., Allison, C. E., and Wahlen, M.: Interannual variability in the oxygen isotopes of atmospheric CO₂ driven by El Nino, Nature, 477, 579–582, doi:10.1038/nature10421, 2011.

Woillez, M.-N.: Modellisation des variations rapides du système atmospherique-ocean-vegetation-cryosphere en climats glaciaires, PhD thesis, Paris Versailles St-Quentin, 256 pp., 2012.

Woillez, M.-N., Kageyama, M., Krinner, G., de Noblet-Ducoudré, N., Viovy, N., and Mancip, M.: Impact of CO₂ and climate on the Last Glacial Maximum vegetation: results from the ORCHIDEE/IPSL models, Clim. Past, 7, 557–577, doi:10.5194/cp-7-557-2011, 2011.

Woillez, M.-N., Kageyama, M., Combourieu-Nebout, N., and Krinner, G.: Simulating the vegetation response in western Europe to abrupt climate changes under glacial background conditions, Biogeosciences, 10, 1561–1582, doi:10.5194/bg-10-1561-2013, 2013.

Yakir, D., D., DeNiro, M. J., and Gat, J. R.: Natural deuterium and oxygen-18 enrichment in leaf water of cotton plants grown under wet and dry conditions: evidence for water compartmentation and its dynamics, Plant Cell Environ., 13, 49–56, doi:10.1111/j.1365-3040.1990.tb01298.x, 1990.