1 New atmospheric carbon isotopic measurements constrain the CO2 rise during the last 2 deglaciation 3 4 Anna Lourantou1, Jošt V. Lavrič1†, Peter Köhler2, Jean-Marc Barnola1+, Didier Paillard3, 5 Elisabeth Michel3, Dominique Raynaud1 and Jérôme Chappellaz1* 6 7 1 8 Joseph Fourier- Grenoble), St Martin d’Hères, France 9 2 Alfred Wegener Institute for Polar and Marine Research, Bremerhaven, Germany 10 3 Laboratoire des Sciences du Climat et de l‘Environnement (IPSL/CEA, CNRS, Université 11 Versailles-St Quentin), Gif-sur-Yvette, France 12 † Now at 3 13 + Deceased Laboratoire de Glaciologie et Géophysique de l'Environnement (LGGE, CNRS, Université 14 15 * Corresponding author: chappellaz@lgge.obs.ujf-grenoble.fr 16 17 Abstract 18 The causes of the ~80 ppmv increase of atmospheric carbon dioxide (CO2) during the 19 last glacial-interglacial climatic transition remain debated. We analyzed the parallel evolution 20 of CO2 and its stable carbon isotopic ratio (δ13CO2) in the EPICA Dome C ice core to bring 21 additional constraints. Agreeing well but largely improving the Taylor Dome ice core record 22 of lower resolution, our δ13CO2 record is characterized by a W-shape, with two negative 23 δ13CO2 excursions of 0.5‰ during Heinrich 1 and Younger Dryas events, bracketing a 24 positive δ13CO2 peak during the Bølling/Allerød warm period. The comparison with marine 25 records and the outputs of two C-cycle box models suggests that changes in Southern Ocean 26 ventilation drove most of the CO2 increase, with additional contributions from marine 27 productivity changes on the initial CO2 rise and δ13CO2 decline and from rapid vegetation 28 buildup during the CO2 plateau of the Bølling/Allerød. 29 30 1. Introduction 31 Atmospheric CO2 is the most important anthropogenic greenhouse gas and arguably the 32 largest contributor to the current global warming [IPCC, 2007]. The monitoring of its stable 33 carbon isotopic ratio (δ13CO2) evolution is useful for the identification of biogeochemical 34 processes driving the observed variations in CO2. Former studies [Friedli et al., 1986; Francey 35 et al., 1999] provided decisive evidence for the man-made origin of the CO2 rise during the 36 last 200 years, based on a ~1.5‰ decline of δ13CO2 to its modern value of –7.8‰. This 37 decrease is caused by the 13C-depleted signature of the two major anthropogenic CO2 sources, 38 fossil fuel burning and carbon release from deforestation, having δ13CO2 values of ~-30‰ and 39 -25‰, respectively. 40 In contrast, natural changes in CO2, such as the 80-ppmv rise over Termination I 41 (hereafter TI), i.e. the transition from the Last Glacial Maximum (LGM ~20 ky BP, Before 42 Present, the present being defined at 1950 Anno Domini AD; ky for kilo -103- years ) to the 43 Early Holocene (EH, ~10 ky BP), are still not well understood. Modeling studies attribute it to 44 various oceanic processes, but without consensus on their relative importance [Broecker & 45 Peng, 1986; Watson & Naveira Garabato, 2006]. Two major mechanisms in the ocean are 46 usually invoked to explain the CO2 glacial-interglacial (G-IG) changes: (i) a physical one, 47 mainly related to Southern (S.) Ocean ventilation changes eventually releasing during 48 Terminations old carbon stored in the deep ocean during the preceding glaciation 49 [Toggweiler, 1999] and (ii) a biological one, involving the efficiency of nutrient utilization by 50 phytoplankton in the Austral ocean, with decreased efficiency (and thus lower CO2 uptake) 51 when atmospheric dust fertilization gets reduced [Archer et al., 2000; Sigman & Boyle, 52 2000]. For more than three decades scientists tried to disentangle the relative role of these or 53 alternative processes, such as changes in oceanic pH and carbonate compensation [Archer et 54 al., 2000], in the evolution of atmospheric CO2. Currently, few models can reproduce the 55 observed amplitude in G-IG CO2 rise; they succeed only if all processes relevant on these 56 time scales are considered [Köhler et al., 2005a; Brovkin et al., 2007]. 57 To validate their hypothesis or to propose alternative ones, more observational 58 constraints are needed. Paleo-atmospheric δ13CO2 makes one of them and is central to our 59 study. So far, a unique record of atmospheric δ13CO2 through TI (including ~15 60 measurements) has been obtained from the Taylor Dome (TD) ice core [Smith et al., 1999 - 61 SM1999 thereafter-], filling the time jigsaw between LGM and EH first produced from the 62 Byrd core [Leuenberger et al., 1992]. Although of coarse resolution, the TD δ13CO2 record 63 has already been used to evaluate the output of several carbon (C) cycle models [Schulz et al., 64 2001; Brovkin et al., 2002; Köhler et al., 2005a; Obata, 2007]. For instance [Obata, 2007], 65 using a coupled climate- C cycle model, simulates a decrease in net primary productivity and 66 soil respiration during the Younger Dryas, in agreement with a combined increase of 67 atmospheric CO2 and minimum of δ13CO2 observed in the TD ice core at that time. [Brovkin 68 et al., 2002] emphasize the role of the G-IG reduced biological pump to explain the 69 simultaneous CO2 increase and δ13CO2 decrease suggested by the TD data during the early 70 part of the Termination. 71 In this study we: (1) present a new highly-resolved record of CO2 and δ13CO2 across TI 72 from the EPICA Dome C (EDC) ice core, (2) compare it to existing ice core data (CO2 from 73 EDC [Monnin et al., 2001] and δ13CO2 from TD [SM1999]), (3) propose a qualitative scenario 74 on the causes of the deglacial CO2 rise, based on a comparison with other proxies, and (4) test 75 this scenario with two C-cycle box models [Köhler et al., 2005a; Paillard et al., 1993]. 76 2. Method 77 A detailed description of the experimental method will be provided in [Lavrič et al., in 78 prep.]. In short, 40-50 g of ice are cut in a cold room, removing about 3 mm of the original 79 sample surface in order to avoid artefacts due to gas diffusion at the atmosphere/ ice interface 80 [Bereiter et al., 2009]. The sample is then sealed in a stainless steel ball mill, evacuated and 81 crushed to fine powder. The gas liberated from the bubbles is expanded over a -80°C ethanol/ 82 liquid nitrogen (LN) water trap onto an evacuated 10 cm3 sample loop. From there it is 83 flushed by an ultra pure helium stream through a partially-heated glass trap where the CO2 is 84 frozen out at LN temperature (-196°C). The trapped CO2 is then transferred into another ultra 85 pure helium stream of lower flow rate, to be cryofocused on a small volume uncoated glass 86 capillary tubing at LN temperature. The subsequent warming of the capillary allows the gas 87 transfer with ultrapure helium into a gas chromatograph to separate the CO2 from residual 88 impurities (e.g. N2O having the same mass over charge ratio as CO2, [Ferretti et al., 2000]), its 89 subsequent passage through an open split system to be finally directed to the isotope ratio 90 mass spectrometer (IRMS, Finnigan MAT 252). 91 92 2.1. Signal determination and correction 93 2.1.1. Standard gases 94 The CO2 mixing ratio in the ice samples is deduced from a linear regression between the 95 varying pressure of several external standard gas injections and the corresponding CO2 peak 96 amplitude measured by the IRMS. The external standard gas has been prepared at CSIRO 97 (Australia) and contains CO2 = 260.3 ± 0.2 ppmv in dry air, with a δ13CO2= –6.40±0.03‰ 98 versus the international standard Vienna Pee-Dee Belemnite, VPDB (δ13CO2 is reported in 99 standard δ notation as the per mil (‰) difference between the stable carbon isotope 100 composition of the sample and VPDB; δ13C = [(13C/12C)sample / (13C/12C)VPDB] –1). It is pre- 101 concentrated and transferred throughout the system similarly as ice-core gas samples. Each 102 sample or external standard introduction in the IRMS is bracketed with injections of a pure 103 CO2 standard reference gas (internal standard, ATMO MESSER, δ13C = -6.5±0.1‰ versus 104 VPDB) through another open split, to calibrate the IRMS and to correct for instrumental drift 105 at the scale of a few minutes. Each spectrogram contains the sample/external standard peak, 106 juxtaposed with peaks eluted from the internal standard gas. The mass over charge (m/z) 44 107 peak height of the internal standard injected with each gas sample is fitted as closely as 108 possible to the expected CO2 peak height from the ice-core gas sample or CSIRO standard, in 109 order to avoid linearity corrections due to the IRMS response. The amount of the external 110 standard gas processed before each ice-core gas sample expansion is also adjusted to the 111 expected gas sample peak height for the same reason. 112 During the experimental protocol, the CSIRO external standard gas is processed seven 113 times before, during and after the ice core gas sample measurement. The latter is usually 114 processed several times, with three consecutive expansions of the same sample gas stored in 115 the extraction container. Thus each data point corresponds to the average value of three 116 replicate measurements of the same extracted gas. The pooled standard deviation on these 117 replicates is 0.98 ppmv for CO2 and 0.098‰ for δ13CO2, while the pooled standard deviation 118 on the routine daily processing of the CSIRO external standard gas is 0.90 ppmv for CO2 and 119 0.15‰ for δ13CO2. The last number does not directly translate to ice core measurements, as it 120 integrates the large daily range of standard gas amount processed through the system and thus 121 the non-linearity of the IRMS response, whereas each ice core gas sample is measured for 122 13 CO2 against a single standard gas peak having a comparable CO2 amplitude. 123 On a daily basis, a correction is applied on the carbon isotopic ratios obtained on ice 124 samples, based on the deviation observed between the external air standard measurements and 125 the attributed CSIRO value. The correction relies on the seven external air standard injections 126 processed before the ice sample, in-between the three expansions of the ice sample, and after 127 the ice sample. On average, a systematic deviation of –0.30‰ from the attributed CSIRO 128 value was observed over the whole EDC measurement period, without any systematic trend 129 from day to day [Lavrič et al., in prep.]. 130 131 132 133 134 2.1.2. Blank tests Three different “blank tests” were conducted throughout the sampling period, with the following differences compared with the procedure described above: (I) No gas introduced in the sample loop. Results show a very low blank (residual 135 traces of CO2 in the transfer lines and carrier gas): we obtain on average (n=35) a CO2 136 amplitude equivalent to 0.33-1.7% of the external standard gas peak heights. 137 (II) A known quantity of external standard gas is introduced in an empty ice mill and 138 then processed to evaluate possible fractionations when expanding a known gas from the cold 139 mill to the sample loop. 140 (III) A known quantity of external standard gas is introduced in the ice mill together 141 with artificial bubble-free ice and then processed after crushing, to reproduce conditions 142 similar to those of a real ice core sample. 143 Results of the last two blank tests are shown in Table 1. 144 CO2 results of the two blank tests are identical to the external standard gas value within 145 the analytical uncertainty (Table 1). The same applies for δ13CO2 in test (II). On the other 146 hand, test (III) with bubble-free ice give an average δ13CO2 depleted by ~0.3‰ compared to 147 the CSIRO value. This may arise from a small fractionation taking place when a gas sample 148 including a small amount of water vapour (vapour pressure at -60°C, i.e. the temperature in 149 the container) is transferred into the vacuum line. We decided not to apply such correction to 150 our measurements, due to insufficient statistics. The absolute values presented here should 151 thus be considered with caution, until we obtain good statistics on applying our system for 152 instance on numerous samples of pre-industrial and industrial ice. The δ13CO2 signal for the 153 past 1000 y is well established [Francey et al., 1999]; numerical deviations obtained with our 154 system would confirm or infirm the need for such blank correction. If any, such small 155 possible bias does not affect the relative δ13CO2 changes observed throughout Termination I. 156 Our results can thus safely be compared one to the other and discussed within the 157 experimental uncertainty range, being on average of 0.1‰. 158 159 2.2. Corrections due to diffusion processes in the firn column 160 Gas molecules in interstitial firn air mostly fractionate by molecular diffusion, in 161 addition to gravitational settling. The latter provokes a preferential accumulation of heavier 162 molecules (for the case of gases) or isotopologues (for the case of isotopes) at the bottom of 163 the firn column compared with the atmosphere [Craig et al., 1988; Schwander et al., 1993]. 164 The fractionation is proportional to the mass difference between the involved gases; the one 165 between 13CO2 and 12CO2, is identical to 15N versus 14N of N2. Therefore we use δ15N of N2 166 data from the EDC core, or modelled δ15Ν of N2 from an empirical relationship with δD in the 167 ice [both provided by G. Dreyfus, pers. communication] to correct δ13CO2 for gravitational 168 fractionation. The CO2 mixing ratio was also corrected for gravitational fractionation, 169 following [Etheridge et al., 1996]. 170 Using measured or modelled δ15Ν of N2 changes the correction by a maximum of 0.03- 171 0.04‰. We finally used the modelled δ15Ν of N2, due to the limited depth coverage of the 172 measured δ15Ν of N2 data. For CO2, the gravitational correction varies from -1.16 to -2.20 173 ppmv, while for δ13CO2; it amounts between -0.41‰ (glacial ice) and -0.55‰ (Holocene ice). 174 Note that such correction was not applied to the previous EDC CO2 record [Monnin et al., 175 2001]. 176 The difference of diffusion coefficient in air between 12 CO2 and 13 CO2 generates 177 changes in the δ13CO2 signal in firn air and trapped bubbles due to molecular diffusion, 178 whenever CO2 varies in the atmosphere, even when atmospheric δ13CO2 remains unchanged. 179 The magnitude of this effect can be calculated with firn air diffusion models [Trudinger et al., 180 1997]. Under present-day conditions when CO2 increases by about 2 ppmv/y, the diffusion 181 correction on firn air and trapped bubbles composition amounts to about 0.10‰ on a 70-m 182 thick firn column [Trudinger et al., 1997]. Since the correction is at first order proportional to 183 the CO2 rate of change, and as the largest observed CO2 rate of change during TI is about 20 184 times smaller than the present-day increasing rate [Joos & Spahni, 2008], the molecular 185 diffusion correction would amount to less than 0.01‰ on the EDC δ13CO2 profile, and is thus 186 neglected here. 187 A final possible correction on gas mixture measured in air bubbles is related to thermal 188 fractionation [Severinghaus et al., 2001; Grachev and Severinghaus, 2003]. As surface 189 temperature changes at EDC were too slow to generate large thermal gradients and gas 190 fractionation, and as no thermal anomaly was detected in the measured δ15Ν of N2 at EDC, no 191 thermal correction was applied to the measured δ13CO2. 192 193 2.3. Reliability of the record 194 Greenland ice has been found to include in situ produced CO2, involving either 195 carbonate/acid reaction or oxidation of organic compounds [Anklin et al., 1995; Tschumi & 196 Stauffer, 2000; Ahn et al., 2004]. No such artifact has been observed so far in Antarctic ice, 197 probably due to the much lower impurity content compared with Greenland ice. 198 All samples measured here originate from the EDC ice core drilled at Concordia Station 199 in Antarctica (75°06’S, 123°21’E; 3233m. above sea level) during the field season 1997-98. 200 Experimental or chemical artifacts affecting CO2 and/or δ13CO2 can be detected when the 201 scatter of duplicates exceeds 3σ of the external precision of the analytical technique. None of 202 the investigated depth levels show such anomaly, thus indicating that the signal can be 203 interpreted within the experimental uncertainty limits. On the other hand, one of the bag 204 sections (dated at 12.6 ky in the EDC3_gas_a scale [Loulergue et al., 2007] cf. next section) 205 provided reproducible mixing and isotopic ratios on duplicate measurements, but its average 206 δ13CO2 differed from neighboring bags (including trapped gas younger or older by less than 207 100 y) by more than 0.2‰. We hypothesize that the corresponding core section has been 208 affected by anomalous storage and local transportation conditions (exposure to warm 209 temperatures), leading to a suspicious result. We thus discard it in the following discussion. 210 211 2.4. Age scale 212 All EDC records are officially dated on the EDC3beta6 [Parrenin et al., 2007] and 213 EDC3_gas_a [Loulergue et al., 2007] age scales for ice and gas data, respectively. However, 214 in order to compare our EDC data with data from other cores (of marine or polar origin) and 215 with model simulations constrained by other datasets, we synchronised both EDC and TD, 216 using CH4 as a time marker, to the newest Greenland chronology GICC05 [Rasmussen et al., 217 2006], using the Analyseries software [Paillard et al., 1996]. The tie-points for each core are 218 presented in Table 2. The synchronised TD chronology is less constrained than the EDC one, 219 due to the poorer time resolution of the TD CH4 record [Köhler et al., 2005a]. The EDC ice 220 chronology (for e.g. δD in Fig. 1a) is obtained by combining the CH4 gas age fit on the 221 GICC05 time scale and the Δage calculated with the EDC3beta6 chronology. 222 223 3. Results 224 Sixty three samples were measured from 50 different depth intervals (345 to 580 m of 225 depth), covering the time period from 9 to 22 ky BP. This provides a mean time resolution of 226 220 y through the transition, whereas the previous published TD record offered a mean time 227 resolution of only ~1000 y. Duplicate analyses of thirteen samples cut on the same ice bags 228 yielded a reproducibility (1σ) of 0.99 ppmv and 0.1‰, respectively. The good correspondence 229 between the reproducibility of CSIRO external standard measurements and of duplicate 230 measurements of neighboring ice samples gives confidence in our main ice core signal 231 structure. Measurements were performed exclusively on clathrate-free ice samples, at depths 232 shallower than 600 m. 233 234 3.1. Comparison with previous datasets 235 The new CO2 and δ13CO2 datasets are plotted together with previously published data 236 (CO2, δD and CH4) from EDC [Monnin et al., 2001; Jouzel et al., 2007; Loulergue et al., 237 2008] and TD [SM1999; Brook et al., 2000], as well as the δ18O data from NGRIP core 238 [NGRIP Members, 2004] in Fig. 1. The agreement between the detailed trends of both CO2 239 records from the same EDC core [Monnin et al., 2001] is remarkable (R2 = 0.996, Fig. 1d). 240 Minor differences in the absolute values result from the use of different CO2 international 241 scales (SIO for the data of [Monnin et al., 2001], CSIRO in this study) and from the 242 gravitational correction only applied to our dataset. The high temporal resolution allows the 243 division of TI into four sub-periods (SP-I to SP-IV) as initiated by [Monnin et al., 2001], 244 characterized by different rates of CO2 change. With 40 measurements throughout TI, the data 245 resolution is improved by more than a factor of two compared with SM1999 (Fig. 1d;e). 246 Overall, the EDC and TD δ13CO2 show similar mean values and trends in the course of TI, 247 with 75% of the TD data falling within the 1σ EDC uncertainty (taking into account dating 248 errors in the comparison). On the other hand, the TD CO2 data are more scattered than the 249 EDC ones. 250 Both EDC and TD δ13CO2 records reveal a W-shape through TI, much more obvious in 251 this new EDC record, with maximum amplitude of contiguous change of ~0.5‰, and a full 252 δ13CO2 range of 0.7‰. The better time resolution of the EDC profile reveals a more 253 structured signal than the TD one within the ~0.1‰ experimental uncertainty, depicting 254 notably faster transitions. This permits for the first time a detailed comparison of the isotopic 255 signal with the changes in the CO2 slope, within an uncertainty range comparable to the TD 256 dataset (given as ±0.085‰ by SM1999). The latter value is probably a low estimate, as the 257 atmospheric N2O trend, needed to apply a correction on the TD δ13CO2 measurements, was 258 considered linear through the deglaciation, whereas the real N2O signal reconstructed since 259 shows a much different structure [Flückiger et al., 1999]. We remind that in our case no such 260 N2O correction is needed (cf. methods section). 261 262 3.2. CO2 and δ13CO2 trends throughout TI 263 Fig. 1d;e reveal a much different behavior between CO2 and δ13CO2: while CO2 mostly 264 shows linear trends within each sub-period (SP), δ13CO2 exposes a more dynamic pattern 265 during the SPs II to IV, with spikes and troughs superimposed on relatively stable boundary 266 values. 267 LGM δ13CO2 also shows a large variability whereas CO2 bears little changes, a feature 268 already observed with similar amplitude in previous datasets [Leuenberger et al., 1992; 269 SM1999]. Part of the LGM δ13CO2 variability parallels very small fluctuations in the CO2 270 rate of change observed in the [Monnin et al., 2001] dataset. Between 22 and 17.6 ky BP, we 271 obtain an average of 188±1 ppmv for CO2 and –6.6±0.1‰ for δ13CO2 (n=10). 272 The evolution of both CO2 and δ13CO2, with respect to Northern and Southern 273 Hemisphere (hereafter NH and SH, respectively) climatic events, can be summarized as 274 follows: 275 276 - Subsequent to the late LGM (22-17.6 ky BP), the early part of TI (SP-I, from 17.6 to 16.2 ky BP) is associated with a 25-ppmv rise of CO2 and a 0.3‰ fall of δ13CO2. 277 - SP-II (16.2 to 14.7 ky BP), during which the Heinrich 1 (H1) event ends in the NH (as 278 deduced from ice-rafted debris in the N. Atlantic [Hemming, 2004] and also seen in NGRIP 279 temperature data in Fig. 1b), reveals a two-step CO2 rise; the first occurs until 15 ky with a 280 progressive 14-ppmv increase and the second with a 12-ppmv rise within only 300 y. 281 Meanwhile, δ13CO2 experiences an oscillation of ~0.2‰ amplitude and reaches a minimum of 282 –7.0±0.1‰ at about 15.5 ky BP, followed by a return to heavier values of ~ -6.8‰. A small 283 δ13CO2 peak also takes place at the start of SP-II, which coincides with a slightly smaller rate 284 of CO2 increase in the detailed Monnin et al. (2001) record. In a recent study, [Barker et al., 285 2009] introduced the notion of “Heinrich Stadial 1” to characterize oceanic conditions during 286 the first two SP; we will refer to this notion in the following. 287 - SP-III (from 14.7 to 12.8 ky BP), coincident with the Antarctic Cold Reversal (ACR) in 288 the SH and the Bølling/Allerød (B/A) warm event in the NH, is marked by a progressive 3- 289 ppmv decrease of CO2, while a positive excursion culminating at –6.5±0.1‰ during the mid- 290 SP-III (~14.1 ky BP) is observed for δ13CO2. 291 - SP-IV (between 12.8 and 11.6 ky BP), during which the Younger Dryas (YD) cold 292 event in the NH and the post-ACR warming in the SH took place, reveals similar patterns for 293 both CO2 and δ13CO2 as for SP-II. Thus, a progressive 13-ppmv CO2 increase is observed 294 until 12 ky, while a more abrupt rise of 10 ppmv is seen during the last 300 y. δ13CO2 295 experiences a negative excursion of more than 0.2‰ amplitude, down to –7.0±0.1‰ (n=6). 296 - The EH (11.6 to 9 ky BP) δ13CO2 mean level is more 13 C-enriched than during SP-IV 13 297 and amounts to -6.8±0.1‰ (n=14). It also seems by 0.2±0.2‰ more depleted in C than at 298 the end of LGM. In contrast, previous studies concluded to more enriched δ13CO2 values 299 during the Holocene (by 0.16‰ in SM1999 to 0.2±0.2‰ [Leuenberger et al., 1992]) than at 300 the LGM. They were based on measurements performed on older LGM ice (Fig. 1d), while 301 Holocene data covered a different time window than considered here (from 9 to 7 ky BP, 302 GICC05 age scale, Fig. 1d;e). In addition, the Holocene δ13CO2 level may be subject to 303 significant fluctuations, as pointed out by SM1999. 304 Our measurements show that, although δ13CO2 starts to decrease in parallel with the 305 early CO2 increase (a trend not captured in the less resolved SM1999 signal), its rapid drop 306 takes place ~1 ky later, when CO2 has already increased by more than 10 ppmv. On the other 307 hand, the EH δ13CO2 rise appears more modest in our dataset than in the TD record. 308 309 4. Discussion 310 Despite the small size of the δ13CO2 signal to be deciphered and the relatively small 311 signal to noise ratio, some clear conclusions can now be drawn on its evolution during TI. Our 312 more detailed EDC δ13CO2 signal compared to the TD one supports some of the earlier 313 conclusions drawn by SM1999. It also sheds some new light on the C-cycle dynamics during 314 the last deglaciation. The EDC record highlights an overall W-shape of atmospheric δ13CO2 315 first broadly depicted by the SM1999 record throughout TI. It differs from SM1999 on the 316 following patterns: 317 1) the two well-resolved minima taking place at times of steadily and important rises of 318 CO2 levels (late part of H1, and YD) reach comparable δ13CO2 levels, around -7.0‰, 319 2) the CO2 plateau accompanying the ACR goes together with a δ13CO2 peak, 320 3) the average δ13CO2 during the EH seems slightly more 13C-depleted than at the end of 321 LGM, 322 4) SM1999 used a plot of δ13CO2 as a function of the inverse of CO2 (a so-called 323 “Keeling plot”, i.e. a mixing diagram where the y-intercept should provide the isotopic 324 composition of the added CO2 in the atmosphere), taking all T-I data together to discuss the 325 possible cause of the CO2 increase. The improved time resolution of our dataset permits us to 326 sub-divide T-I with distinct intercepts through time. A y-intercept of -6‰ is obtained, similar 327 to all deglacial data of SM1999, but only through SP-II and SP-III data (not shown). On the 328 other hand, the two periods when CO2 largely increases and δ13CO2 simultaneously decreases 329 in our record (SP-I and SP-IV) reveal another “Keeling plot” type of isotopic signature for the 330 additional CO2, similar for both sub-periods: ~-11‰ (Fig. 4). The main conclusion of 331 SM1999 that the C-cycle behaved in a dual mode depending on the speed of climatic changes 332 i.e. a slow mode taking place during the EH and LGM, and a fast mode during TI, is thus not 333 supported by the new Keeling plot (see supplementary material). 334 The similar “Keeling plot” signature of SP-I and SP-IV suggests at first hand that the 335 two main steps of atmospheric CO2 increase during TI involved similar C-cycle mechanisms. 336 But their common y-intercept cannot be directly interpreted as the isotopic signature of such 337 mechanisms. Keeling plots work well only in an atmosphere-biosphere two-reservoir system 338 experiencing fast exchanges [e.g. Pataki et al., 2003]. On time scales of centuries to millennia 339 such as during TI, the isotopic buffering effect of the ocean (air/sea exchanges, carbonate 340 system) modifies the y-intercept in a three-reservoir model, as shown for instance by [Köhler 341 et al., 2006a] using the pre-industrial to industrial CO2 increase and δ13CO2 decrease as a case 342 study. Therefore, other approaches are required to extract possible scenarios out of our new 343 dataset, relevant to carbon exchanges between the atmosphere, ocean and biosphere during 344 TI. We use two of them here: a comparison with proxy records relevant to C-cycle processes, 345 and simulations of CO2 and δ13CO2 with two C-cycle box models. 346 347 4.1. Comparison with other C-cycle proxy records 348 The good correlation between CO2 and Antarctic deuterium throughout TI (Fig. 1a;d), 349 already noticed in numerous works [e.g. Monnin et al., 2001; Bianchi and Gersonde, 2004], 350 points towards a leading role of the S. Ocean to drive the corresponding CO2 evolution. As 351 pointed out in the introduction, two types of S. Ocean processes, a biological and a physical 352 one, can be evoked. 353 According to the first one, the S. Ocean during the LGM experienced a higher 354 productivity due to higher atmospheric dust fluxes bringing more iron, a limiting 355 micronutrient in high nutrient low chlorophyll (HNLC) regions [Martin, 1990]. As recorded 356 in e.g. EDC ice [Lambert et al., 2008; Gaspari et al., 2006] and shown in Fig. 2a, the 357 atmospheric dust (and iron) flux considerably decreases between ~18 and 14.6 ky BP, 358 corresponding to the first half of the CO2 deglacial increase. This would imply a decreasing 359 biological pump in the S. Ocean exporting less carbon to the ocean interior and thus 360 increasing atmospheric CO2. As phytoplankton preferentially assimilates the lighter carbon 361 isotope (12C), a decreasing productivity would be accompanied by a decreasing atmospheric 362 δ13CO2 [Brovkin et al., 2002], in agreement with our record. During the second half of TI, the 363 low dust values encountered in EDC ice suggest that the biologically-mediated mechanism in 364 the S. Ocean did not influence the CO2 and δ13CO2 trends. 365 The physical mechanism involves the rate of vertical mixing of the S. Ocean: the cold 366 LGM was associated with increased sea ice extent (mostly in winter) and with considerable 367 stratification of the S. Ocean water column [Sigman and Boyle, 2000; Stephens and Keeling, 368 2000; Marchitto et al., 2007]. The deep S. Ocean thus held a large amount of CO2, due to 369 organic matter remineralization, with a strongly 370 decomposed organic matter. [Duplessy et al., 1988] showed that changes in Atlantic 371 circulation at the end of LGM might have transferred low-δ13C deep waters towards the ocean 372 surface, a phenomenon validated subsequently by [Curry and Oppo, 2005]. Overall, the 373 deglaciation, combining sea ice retreat, possible shifts of westerlies, and collapse of North 374 Atlantic Deep Water (NADW) formation during its early phase, would have generated a S. 13 C-depleted signature originating from 375 Ocean stratification breakdown and hence, the release of deep ocean 13C-depleted CO2 in the 376 atmosphere, leading to an atmospheric CO2 increase paralleled with a decrease of δ13CO2 [e.g. 377 Toggweiler et al., 2006; Menviel et al., 2008]. Such mechanism could have also acted during 378 the YD, with a pause in-between during the B/A, when the NADW was probably switched on 379 again [Knorr and Lohmann, 2003]. As shown in Fig. 2, several proxy records, matching the 380 general shape of our δ13CO2 record within their respective age model uncertainties, are in line 381 with this physical scenario: 382 1) During the Heinrich Stadial 1 and YD, NADW formation got weakened [Marchitto et 383 al., 1998], as evidenced by an increased 231Pa/230Th ratio in North Atlantic sediments (Fig. 2c) 384 [McManus et al., 2004]. This implies less NADW signal propagation towards the S. Ocean 385 waters. 386 2) NADW reduction is accompanied by a flushing of deep waters from the S. Ocean into 387 the Atlantic basin, thus equilibrating the water mass loss in the North Atlantic region. The 388 nutrient-enriched and 389 Intermediate Waters, AAIW) compared with deep waters from North Atlantic, is registered in 390 North Atlantic marine sediments through two negative δ13C-excursions during the Heinrich 391 Stadial 1 and YD [Rickaby & Elderfield, 2005] (Fig. 2d). 13 C-depleted signal of old Antarctic-dwelled waters (e.g. Antarctic 392 3) A high-resolution Δ14C record from the North Pacific [Marchitto et al., 2007], shown 393 in Fig. 2e, reveals two negative excursions of more than 200‰ during the Heinrich Stadial 1 394 and the YD. They are interpreted as two episodes of transfer of old AAIW (including aged C 395 of up to 4-5 ky) into intermediate waters of the North Pacific, associated with the sea-ice 396 retreat [Stephens and Keeling, 2000] and the S. Ocean stratification breakdown [Marchitto et 397 al., 2007; Schmittner et al., 2007]. This should be accompanied by the release of sequestered 398 and 13C-depleted deep oceanic carbon into the atmosphere. The scenario is corroborated at the 399 start of TI by a low resolution planktonic δ13C record in subantarctic and equatorial Pacific 400 [Ninnemann and Charles, 1997; Spero and Lea, 2002]. A recent high-resolution tropical 401 Pacific planktic δ13C record from [Stott et al., 2009], also reveals a nice W-trend throughout 402 the deglaciation, reinforcing the scenario of G-IG S. Ocean upwelling changes mentioned 403 above (Fig. 2f). 404 4) A high-resolution record of opal flux in the S. Ocean [Anderson et al., 2009], plotted 405 in Fig. 2g, shows an increase of upwelling strength in two steps, coincident with the Heinrich 406 Stadial 1 and YD, thus also pointing towards increased ventilation of deep S. Ocean waters as 407 a main trigger of the two steps in the deglacial CO2 increase. 408 The physical mechanism involving S. Ocean stratification breakdown in two episodes 409 during TI thus qualitatively matches the CO2 and δ13CO2 trends, and could also explain the 410 common “Keeling plot” isotopic signature of the added carbon during the two episodes. Aside 411 from the biological pump and ocean circulation hypotheses, a possibly straightforward 412 explanation of the co-evolution between CO2 and δ13CO2 concerns changes in Sea Surface 413 Temperature (SST) during TI: due to isotopic fractionation during air/sea exchanges, a 414 warmer ocean will leave a 415 portion of their δ13CO2 signal to this mechanism, by considering a globally averaged SST 416 increase of 5°C between the LGM and EH. On the other hand [Brovkin et al., 2002], using the 417 climate model CLIMBER-2, pointed out that the parallel changes of CO2, alkalinity and 418 bicarbonate ion concentration significantly affect the isotopic fractionation during air/sea 419 exchanges, thus reducing the atmospheric δ13CO2 imprint of SST changes. Moreover, the 420 rapid δ13CO2 changes observed in our record may be difficult to reconcile with the speed of 421 SST changes in areas of deep water formation. 13 C-enriched signal in the atmosphere. SM1999 assigned a large 422 SP-III encounters terrestrial carbon buildup in vegetation, soils and peat deposits 423 [MacDonald et al., 2006] which could have contributed to the small CO2 decrease and to a 424 positive δ13CO2 anomaly in the atmosphere, as the biosphere preferentially assimilates 12 C. 425 This is qualitatively corroborated by the CH4 evolution (Fig. 1c), pointing towards a switch- 426 on of boreal wetland CH4 emissions (requiring a concomitant intensification of the terrestrial 427 C-cycle) at that time [Fischer et al., 2008]. Alternative scenarios attribute even more control 428 of the δ13CO2 variability by terrestrial biosphere carbon uptake and release as a consequence 429 of abrupt temperature changes in the NH caused by the AMOC shutdown during H1 and YD 430 (e.g. Köhler et al., 2005b). For the deglacial CO2 rise, the contribution from progressively 431 flooded continental shelves might also need some consideration [Montenegro et al., 2006]; 432 however this scenario is challenged by the lag of the sea level increase with respect to CO2 433 [Pépin et al., 2001 434 In summary, the qualitative comparison of the EDC CO2 and δ13CO2 records with C- 435 cycle proxies suggest a dominant role of increased overturning in the S. Ocean (as mainly 436 evidenced by Δ14C, opal flux and δ13C records) during SP-I and SP-II to explain the two main 437 steps of CO2 increase, with an additional contribution of reduced biological pump during SP- 438 I. The two mechanisms would have stalled during SP-III when NADW became stronger, and 439 would also have been counterbalanced by terrestrial carbon buildup. To go further into a 440 quantitative evaluation of mechanisms able to explain the CO2 and δ13CO2 signals, modeling 441 is required. In the following, we explore the problem using two C-cycle box models. 442 443 4.2. BICYCLE model runs 444 We employed the BICYCLE model [Köhler et al., 2005a], a coupled atmosphere/ 445 ocean/ sediment/ biosphere C-cycle box model, run in a transient mode and forced with 446 various time-dependent paleoclimatic data over TI. It consists of a single atmospheric box 447 interacting with a 10-box ocean reservoir and the terrestrial biosphere, which is sub-divided 448 into 7 compartments [Köhler et al., 2005a]. The ocean further communicates with a sediment 449 reservoir. Mass balance equations are solved for the carbon stocks of the biospheric 450 compartments, for DIC, TAlk, PO4 and O2 in the 10 oceanic reservoirs, for CO2 in the 451 atmosphere and for the carbon isotopes in all reservoirs. 452 BICYCLE is the only C-cycle model we are aware of which was run in transient mode 453 over TI. It was also used for the interpretation of atmospheric carbon records and deep ocean 454 δ13C data over TI and much longer time scales of up to 2 My [Köhler et al., 2005a; Köhler et 455 al., 2006a; Köhler et al., 2006b; Köhler et al., 2006c; Köhler & Bintanja, 2008; Köhler et al., 456 subm to Paleoceanography]. Since its application over TI [Köhler et al., 2005a], 457 improvements were performed in the parameterization of ocean circulation and sediment- 458 ocean interaction, we thus use new simulation results, instead of the model output published 459 in [Köhler et al., 2005a]. 460 461 462 463 4.2.1. Parameterizations The main model parameterizations, based on data obtained from ice core or marine cores (Fig. 3) are the following: 464 1) Sea level rises by ~110m between 22 and 8 ky BP based on reconstructions of coral 465 reef terraces [Fairbanks, 1989] (Fig. 3A). This leads to changes in the salinity, in the 466 concentrations of all oceanic tracers and in the volumes of the oceanic boxes. 467 2) Temperature of all oceanic boxes is prescribed for present day from [Levitus & Boyer, 468 1994]. It changes over time according to oceanic proxy evidences for equatorial SST [Visser 469 et al., 2003] and deep ocean temperature [Labeyrie et al., 1987] (Fig. 3C). At high latitudes, it 470 is represented by ice core isotopic profiles (North Atlantic and North Pacific: δ18O on 471 GICC05 age scale from NorthGRIP [NGRIP Members, 2004; Andersen et al., 2007], Fig. 3B; 472 S. Ocean: δD (corrected for the effect of sea level rise) from EDC [Parrenin et al., 2007; 473 Jouzel et al., 2001] synchronised to GICC05 [EPICA Community Members, 2006; Andersen 474 et al., 2007], as shown in Fig 3D. Both ice core records are scaled to provide a SST ΔT of 4 K 475 between the minimum glacial values and the present-day. 476 3) Marine productivity in the S. Ocean is scaled (if allowed by macro-nutrient availability) 477 on dust input to the S. Ocean as approximated by the non sea-salt-dust record measured in 478 EDC [Roethlisberger et al., 2002] (Fig 3E). 479 4) Ocean circulation between the 10 boxes for present conditions is parameterized with 480 data from the World Ocean Circulation Experiment WOCE [Ganachaud & Wunsch, 2000] 481 (Fig 3F). Compared to the initial BICYCLE runs over TI [Köhler et al., 2005a], it was slightly 482 modified to get a better agreement between simulated and reconstructed oceanic 13C [Köhler 483 et al., subm. to Paleoceanography]. About 30% of the upwelled waters in the S. Ocean are 484 immediately redistributed to the intermediate equatorial Atlantic Ocean to account for the 485 effect that, in the natural carbon cycle, upwelling waters in the S. Ocean (which are then 486 flowing as water masses of intermediate depth to the north) are still enriched in DIC (Gruber 487 et al., 2009). 488 Three major ocean currents are parameterized as follows (Fig. 3F): 489 - The strength of NADW (i.e. of its overturning, cf. [Köhler et al., 2005]) is assumed to 490 be about 40% weaker [Meissner et al., 2003] during the LGM than at present day (10 versus 491 16 Sv), 2 Sv, 13 Sv and 11 Sv during the H1, B/A and YD, respectively [McManus et al., 492 2004]. 493 - Antarctic-dwelled waters (e.g. Antarctic Bottom Waters, AABW), penetrating both 494 the deep Atlantic (AABW_A) and deep Pacific (AABW_P), is strengthened when NADW 495 weakens (i.e. during H1 and YD) and vice versa (for B/A), considering the north-south 496 opposite trend during these abrupt climatic changes [Broecker, 1998; Rickaby & Elderfield, 497 2005; Kissel et al., 2008]. The EH and LGM AABW overall flux is set at of 15 Sv (6 SV for 498 the Atlantic branch and 9 Sv for the Pacific one). During H1 and YD (B/A), each section is 499 strengthened (weakened) by 3 Sv. 500 - Vertical mixing in the S. Ocean (SOX) is set to 0 Sv during LGM in accordance with 501 proxy evidence (e.g. δ13C data from [Hodell et al., 2003; Spero & Lea, 2002]). Just after the 502 dust proxy (nss-Ca2+) drop during SP-I, it is set to 15 Sv and maintained throughout B/A to be 503 finally increased by another 5 Sv at the end of YD [Köhler et al., 2005]. 504 5) Changes in the terrestrial biosphere carbon pool are assumed to be primarily 505 temperature-dependent and made proportional to ¾ of the NorthGRIP δ18O and ¼ of the EDC 506 δD changes (taken as temperature proxies), reflecting the latitudinal distribution of vegetated 507 land. The G-IG amplitude of land temperature change is considered as 8 K in the North and 5 508 K in the South. Net primary productivity is also parameterized on the modelled atmospheric 509 CO2 values to take into account CO2 fertilisation. More details on the terrestrial biosphere 510 module can be found in [Köhler & Fischer, 2004]. 511 6) The additional effect of carbonate compensation [Archer & Maier-Reimer, 1994] to all 512 temporal changing processes listed above, is considered with a relaxation approach bringing 513 the deep ocean carbonate ion concentration back to initial values. 514 515 All ice core records (isotopic temperature proxies and nss-dust) are implemented as 500-y running means in the different parameterizations. 516 517 4.2.2. BICYCLE model output 518 The left panel of Fig. 5 illustrates the imprint Δ (with respect to the EH value) of major 519 processes simulated with BICYCLE, on atmospheric (a) CO2 and (b) δ13CO2. The reduction 520 of S. Ocean biological productivity (due to the onset of Fe-limitation in HNLC [Martin, 521 1990]), as well as the S. Ocean stratification breakdown [Spero & Lea, 2002] (associated with 522 sea-ice retreat and decreasing salinity [Watson & Naveira Garabato, 2006; Stephens & 523 Keeling, 2000]) are the main processes at work in the BICYCLE simulation at the TI 524 inception, provoking a 15 and 22-ppmv CO2-rise and a 0.20 and 0.32 ‰ δ13CO2 decline, 525 respectively. 526 During the NH cold events (H1 and YD), NADW weakens [McManus et al., 2004], 527 dampening the CO2 increase and δ13CO2 decrease related to AABW enhancement (+4.5 528 ppmv; -0.04‰) [Rickaby & Elderfield, 2005]; NADW (AABW) strengthening (weakening) at 529 the end of SP-II and SP-IV, combined with stronger S. Ocean water mixing at the end of YD, 530 lead to a CO2 out-gassing of 10 and 7 ppmv, respectively (see supplementary material). 531 Sea level rise [Fairbanks, 1989] processes do not leave an important imprint on δ13CO2 532 within the SPs, although they significantly affect CO2 during the ACR, by provoking a 3- 533 ppmv reduction (see supplementary material). In contrast, vegetation growth, lagging S. 534 Ocean warming [Hughen et al., 2004] and forced by CO2 fertilization and NH warming, starts 535 affecting δ13CO2 during SP-II and becomes a major driver of this signal during SP-III and SP- 536 IV (green line of Fig. 5a;b). The rise-and-fall of total biospheric carbon by 200PgC during 537 SP-III and SP-IV respectively [Köhler & Fischer, 2004; Köhler et al., 2005a], lead to a 15- 538 ppmv decrease and a 17-ppmv rise of CO2, also causing a +0.35‰ and -0.40‰ δ13CO2 539 anomaly. This is in corroboration with previous numerical studies: for e.g. [Scholze et al., 540 2003] assuming a total 180PgC decline of terrestrial carbon pools during YD, result in an 541 atmospheric CO2 rise of 30 ppmv due to land cooling and precipitation decline, both 542 following the overturning circulation reduction. Köhler et al., 2005b find, as a consequence of 543 the AMOC shutdown and the accompanying northern hemispheric cooling, a total decline of 544 terrestrial carbon pools of up to 140 PgC, resulting in peak-to-peak changes in atmospheric 545 CO2 and δ13CO2 of 13 ppmv and 0.25‰, respectively. Recently, [Brovkin et al., 2007] used a 546 model of intermediate complexity to evaluate the shared contributions of different C-cycle 547 mechanisms on CO2 and δ13CO2 G-IG changes. They conclude to relative imprints of ocean 548 circulation, SST, land and marine productivity changes on δ13CO2 very close to the BICYCLE 549 results. 550 4.2.3. Data / model comparison 551 The integrated signal from all the processes simulated with BICYCLE is compared with 552 our data in Fig. 5c for CO2 and Fig. 5d for δ13CO2. Throughout TI, BICYCLE produces an 553 increasing CO2 trend, interrupted by phases of slower rate of increase or of stabilization, and 554 accompanied by a marked δ13CO2 W-shaped trend. The lowest δ13CO2 is simulated during 555 SP-II and SP-IV, with similar values around -7.0‰. A δ13CO2 peak at the start of SP-III 556 reaches ~-6.5‰. BICYCLE thus captures the main features of the EDC records. A direct 557 comparison of Keeling plots, obtained with the data and the simulations is provided in Fig. 4. 558 Similar y-intercepts are obtained for the two main periods of abrupt CO2 rise, suggesting that 559 the sequence and amplitude of the involved processes are well captured by the model 560 configuration. The Keeling plot comparison also highlights the limit of such plot, as the 561 similar y-intercepts generated by BICYCLE for the two periods come from a different 562 combination of C cycle mechanisms at work. 563 On the other hand, the timing of changes (for both CO2 and δ13CO2) can differ between 564 observations and model outputs. The EDC δ13CO2 peak of SP-III and the minimum of SP-IV 565 appear earlier in the BICYCLE simulation. Both δ13CO2 features mainly result from the 566 terrestrial component in BICYCLE, itself mainly parameterized on NH temperature. As both 567 EDC δ13CO2 and the biosphere imprint in BICYCLE are on a common time scale (GICC05), 568 the shift cannot be attributed to dating errors. One explanation lies in the time response of 569 biospheric components to climate change being possibly underestimated in BICYCLE (a lag 570 of ~400y, also found by [Scholze et al., 2003]). 571 Another data/model difference appears during SP-IV, when BICYCLE simulates a CO2 572 plateau and a large δ13CO2 increase whereas the EDC data reveal a steady increase and a 573 minimum, respectively. Such BICYCLE output clearly depends on how well the timing and 574 amplitude of SST and ocean circulation changes are parameterized during the YD and EH. As 575 pointed out in the previous section, C-cycle proxy data suggest that the S. Ocean mixing 576 should have driven the YD CO2 increase (it is parameterized as constant in the model), 577 whereas the BICYCLE simulation gives more weight to terrestrial carbon release and SST 578 increase. 579 The box model bears other uncertainties and potential biases, such as: (i) the relative 580 dating of the various input signals and their synchronization with ice cores, (ii) the pertinence 581 of the proxies used for each process (e.g. the magnitude of oceanic fluxes throughout TI), (iii) 582 the coarse spatial resolution of low-latitudes. Still, it shows that the general shape of the EDC 583 CO2 and δ13CO2 signals can be reproduced with a reasonable temporal sequence of C-cycle 584 mechanisms. It supports a scenario where S. Ocean stratification breakdown and decrease of 585 marine productivity jointly explain the early half of the CO2 and δ13CO2 signals, with the 586 terrestrial biosphere intervening in the shape of both signals during the B/A. The conclusion 587 to be drawn for the YD episode is less clear, as the box model produces a δ13CO2 minimum 588 but fails to simulate the parallel CO2 increase. 589 A crucial point when comparing temporally highly-resolved atmospheric records 590 derived from ice cores with transient model simulations is that the gas records are smoothed 591 by gas diffusion in the firn and by progressive bubble close-off. Therefore, they do not 592 represent one single point in time, but are averaged over decades to centuries, mainly 593 depending on accumulation rate and temperature. A gas diffusion and enclosure model 594 [Spahni et al., 2003] was used earlier to calculate the age distribution for CO2 in EDC and the 595 attenuation of atmospheric signals during the enclosure process. It has been calculated that the 596 gas records represent averages between 213 (preindustrial) and 590 (LGM) y with a 597 lognormal-shapelike age distribution [Joos and Spahni, 2008]. As a consequence, BICYCLE 598 model simulations, which should represent atmospheric records before gas enclosure, might 599 not be directly comparable with ice core records, especially for low accumulation sites and 600 fast processes, because original atmospheric amplitudes are attenuated during the enclosure 601 process [Köhler et al, submitted]. A solution would be to proceed to similar measurements for 602 the same time interval on a core with larger accumulation rate. 603 604 4.3. BOXKIT model 605 We also applied BOXKIT, a conceptual ocean/atmosphere model run under equilibrium 606 states [Paillard et al., 1993]. The ocean is splitted in 10 boxes, 5 for the surface, 2 for 607 intermediate waters and 3 for the deep ocean. BOXKIT includes a single atmospheric box, but 608 no terrestrial biosphere. The same forcings as BICYCLE are applied for 6 ‘snapshots’ over 609 TI. 610 Similar overall trends as those of BICYCLE, both for CO2 and δ13CO2 are obtained (red 611 squares in Fig. 5c;d). As BOXKIT provides easier tuning than BICYCLE to carry on 612 sensitivity tests, we used it to evaluate the output sensitivity to low latitudes SST. Increasing 613 tropical SSTs by 3°C (instead of 0.5°C as done for BICYCLE forcings) for the SP-III 614 simulation, concomitantly with NH warming as is seen in N. Atlantic sediment data from e.g. 615 [Lea et al., 2003], leads to a δ13CO2 increase by ~0.2‰, more in line with the EDC data. This 616 example shows the non-uniqueness of solutions when interpreting the C-cycle data with box 617 models. 618 619 5. Conclusions 620 Our new record of δ13CO2 from the EDC ice core over the last deglaciation reveals 621 sharp fluctuations mostly associated with variations in the CO2 rate of change. A comparison 622 with other CO2 and δ13CO2 ice core data gives confidence in the validity of this new dataset. 623 In addition, consistent standard deviations are observed between different statistical 624 approaches of the experimental system. The general shape of the deglacial δ13CO2 signal can 625 be summarized as a “W”, with two minima accompanying the two major steps of CO2 626 increase, and a peak when CO2 gets stabilized or slightly decreasing. 627 The comparison with C-cycle related proxies highlights similarities with marine signals 628 associated with the strength of S. Ocean ventilation and upwelling, suggesting that this 629 physical mechanism would be the main driver of the deglacial CO2 increase. 630 Two C-cycle box models (BICYCLE and BOXKIT), run under the same input 631 parameters support the dominant role of S. Ocean physical processes and add the marine 632 productivity decline during the early part of the deglaciation as another mechanism 633 contributing to the δ13CO2 decrease and CO2 increase. The BICYCLE model supports an 634 additional role of terrestrial carbon buildup to explain the CO2 plateau and δ13CO2 peak 635 paralleling the ACR. It also simulates an early YD δ13CO2 minimum followed by an increase 636 to EH values, attributed to terrestrial carbon and SST decrease and subsequent increase, an 637 explanation conflicting with C-cycle proxy data which suggest a dominant role of 638 strengthening S. Ocean ventilation. The failure of BICYCLE to simulate a parallel CO2 639 increase shows the limit of this modeling exercise, which crucially depends on assumptions 640 regarding SO upwelling changes. 641 More sophisticated approaches using coupled carbon-climate Earth system models will 642 be needed in the future to better disentangle the contribution of each process, with their direct 643 parameterizations in the models instead of the use of proxies. Our detailed EDC profile 644 clearly highlights the need for fine time resolution in producing future δ13CO2 records 645 throughout major climatic events. 646 647 648 Acknowledgments 649 This work is a contribution to the European Project for Ice Coring in Antarctica 650 (EPICA), a joint ESF (European Science Foundation)/ EC scientific program, funded by the 651 European Commission and by national contributions from Belgium, Denmark, France, 652 Germany, Italy, the Netherlands, Norway, Sweden, Switzerland and the United Kingdom. The 653 main logistic support was provided by IPEV and PNRA (at Dome C). AL was funded by the 654 European Research Training and Mobility Network GREENCYCLES. Additional funding 655 support was provided by the QUEST-INSU project DESIRE, the FP6 STREP EPICA-MIS, 656 and the French ANR PICC (ANR-05-BLAN-0312-01). Long-term support for the mass 657 spectrometry work at LGGE was provided by the Fondation de France and the Balzan Price. 658 Discussions with G. Dreyfus, H. Schaefer and G. Delaygue were very much appreciated. We 659 particularly thank C. Lorius for his confidence in our earlier work and five anonymous 660 reviewers for fruitful comments on previous versions of this manuscript. This is EPICA 661 publication no XXX. 662 663 Reference list 664 Ahn J., Wahlen M., Deck B. L., Brook E. J., Mayewski P. A., Taylor K. C. and White J. W. 665 C. (2004): A record of atmospheric CO2 during the last 40,000 years from the Siple Dome, 666 Antarctica 667 10.1029/2003JD004415. 668 Andersen K. K., Svensson A., Johnsen S. J., Rasmussen S. O., Bigler M., Röthlisberger R., 669 Ruth U., Siggaard-Andersen M.-L., Steffensen J. P., Dahl-Jensen D., Vinther B. M. and 670 Clausen H. B. (2007): The Greenland Ice Core Chronology 2005, 15–42 ka. Part 1: 671 constructing the time scale. Quaternary Science Reviews 25(23-24), pp. 3246-3257. ice core. Journal of Geophysical Research 109, D13305, DOI: 672 Anderson R. F., Ali S., Bradtmiller L. I., Nielsen S. H. H., Fleisher M. Q., Anderson B. E. and 673 Burckle L. H. (2009): Wind-driven upwelling in the Southern Ocean and the deglacial rise in 674 atmospheric CO2. Science 323(5920), pp. 1443-1448. 675 Anklin M., Barnola J.-M., Schwander J., Stauffer B. and Raynaud D. (1995): Processes 676 affecting the CO2 concentrations measured in Greenland ice. Tellus 47B, pp. 461-470. 677 Archer D., and E. Maier-Reimer (1994): Effect of deep-sea sedimentary calcite preservation 678 on atmospheric CO2 concentration. Nature 367(6460), pp. 260 - 263. 679 Archer D., Winguth A., Lea D. and Mahowald N. (2000): What caused the glacial/interglacial 680 atmospheric pCO2 cycles? Reviews of Geophysics 38(2), pp. 159-190. 681 Barker S., Diz P., Vautravers M. J., Pike J., Knorr G., Hall I. R. and Broecker W. S. (2009): 682 Interhemispheric Atlantic seesaw response during the last deglaciation. Nature 457(7233), pp. 683 1097-1102. 684 Bereiter B., Schwander J., Lüthi D. and Stocker T. F. (2009): Change in CO2 concentration 685 and O2/N2 ratio in ice cores due to molecular diffusion Geophysical Research Letters 36, 686 L05703, DOI: 10.1029/2008GL036737. 687 Bianchi C., and R. Gersonde (2004): Climate evolution at the last deglaciation: the role of the 688 Southern Ocean. Earth and Planetary Science Letters, 228(3-4), pp. 407-424. 689 Broecker W. S., and T. H. Peng (1986): Carbon cycle: 1985 glacial to interglacial changes in 690 the operation of the global carbon cycle. Radiocarbon, 28(2A), pp. 309-327. 691 Broecker W. S. (1998): Paleocean circulation during the last deglaciation: A bipolar seesaw? 692 Paleoceanography 13(2), pp. 119-121. 693 Brook E. J., Harder S., Severinghaus J., Steig E. J. and Sucher C. M. (2000): On the origin 694 and timing of rapid changes in atmospheric methane during the last glacial period. Global 695 Biogeochemical Cycles 14(2), pp. 559–572. 696 Brovkin V., Hofmann M., Bendtsen J. and Ganopolski A. (2002): Ocean biology could 697 control atmospheric δ13C during glacial-interglacial cycle. Geochem., Geophys., Geosyst. 698 3(5), DOI: 10.1029/2001GC000270. 699 Brovkin V., Ganopolski A., Archer D. and Rahmstorf S. (2007): Lowering of glacial 700 atmospheric CO2 in response to changes in oceanic circulation and marine biogeochemistry. 701 Paleoceanography 22, PA4202, DOI: 10.1029/2006PA001380. 702 Craig H., Horibe Y. and Sowers T. (1988): Gravitational separation of gases and isotopes in 703 polar ice caps. Science 242(4886), pp. 1675 - 1678. 704 Curry W. B. and Oppo D. W. (2005): Glacial water mass geometry and the distribution of 705 δ13C of ΣCO2 in the western Atlantic Ocean. Paleoceanography 20, PA1017, DOI: 706 10.1029/2004PA001021. 707 Duplessy J. C., Shackleton N. J., Fairbanks R. G., Labeyrie L., Oppo D. and Kallel N. (1988): 708 Deepwater source variations during the last climatic cycle and their impact on the global 709 deepwater circulation. Paleoceanography 3(3), pp. 343-360. 710 EPICA Community Members (2006): One-to-one coupling of glacial climate variability in 711 Greenland and Antarctica. Nature, 444(7116), pp. 195-198. 712 Etheridge D. M., Steele L. P., Langenfelds R. L., Francey R. J., Barnola J.-M. and Morgan V. 713 I. (1996): Natural and anthropogenic changes in atmospheric CO2 over the last 1000 years 714 from air in Antarctic ice and firn. Journal of Geophysical Research 101(D2), pp. 4115–4128. 715 Fairbanks R. G. (1989): A 17,000-year glacio-eustatic sea level record: influence of glacial 716 melting rates on the Younger Dryas event and deep-ocean circulation. Nature, 342(6250), pp. 717 637 – 642. 718 Ferretti D. F., Lowe D. C., Martin R. J. and Brailsford G. W. (2000): A new gas 719 chromatograph-isotope ratio mass spectrometry technique for high-precision, N2O-free 720 analysis of δ13C and δ18O in atmospheric CO2 from small air samples. Journal of Geophysical 721 Research 105(D5), pp. 6709-6718. 722 Fischer H., Behrens M., Bock M., Richter U., Schmitt J., Loulergue L., Chappellaz J., Spahni 723 R., Blunier T., Leuenberger M. and Stocker T. F. (2008): Changing boreal methane sources 724 and constant biomass burning during the last termination. Nature 452(7189), pp. 864-867. 725 Flückiger J., Dällenbach A., Blunier T., Stauffer B., Stocker T. F., Raynaud D. and Barnola 726 J.-M. (1999): Variations in atmospheric N2O concentration during abrupt climatic changes. 727 Science 285(5425), pp. 227-230. 728 Francey R. J., Allison C. E., Etheridge D. M., Trudinger C. M., Enting I. G., Leuenberger M., 729 Langenfelds R. L., Michel E. and Steele L. P. (1999): A 1000-year high precision record of 730 δ13C in atmospheric CO2. Tellus B 51(2), pp. 170–193. 731 Friedli H., Lötscher H., Oeschger H., Siegenthaler U. and Stauffer B. (1986): Ice core record 732 of the 13C/12C ratio of atmospheric CO2 in the past two centuries. Nature 324(6094), pp. 237- 733 238. 734 Ganachaud A. and Wunsch C. (2000): Improved estimates of global ocean circulation, heat 735 transport and mixing from hydrographic data. Nature 408(6811), pp. 453-457. 736 Gaspari V., Barbante C., Cozzi G., Cescon P., Boutron C. F., Gabrielli P., Capodaglio G., 737 Ferrari C., Petit J. R. and Delmonte B. (2006): Atmospheric iron fluxes over the last 738 deglaciation: Climatic implications. Geophysical Research Letters 33, L03704, DOI: 739 10.1029/2005GL024352. 740 Grachev A. M. and Severinghaus J. P. (2003): Laboratory determination of thermal diffusion 741 constants for 742 magnitudes of abrupt climate changes using the ice core fossil–air paleothermometer. 743 Geochimica et Cosmochimica Acta 67(3), pp. 345-360. 29 N2/28N2 in air at temperatures from −60 to 0 °C for reconstruction of 744 Gruber N., Gloor M., Mikaloff Fletcher S. E., Doney S. C., Dutkiewich S., Follows M. J., 745 Gerber M., Jacobson A. R., Joos F., Lindsay K., Menemenlis D., Mouchet A., Müller S. A., 746 Sarmiento J. L., Takahashi T. (2009): Oceanic sources, sinks and transport of atmospheric 747 CO2. Global Biogeochemical Cycles 23, GB1005, DOI: 10.1029/2008GB003349. 748 Hemming S. R. (2004): Heinrich events: Massive late Pleistocene detritus layers of the North 749 Atlantic and their global climate imprint. Reviews of Geophysics 42, RG1005, DOI: 750 10.1029/2003RG000128. 751 Hodell D. A., Venz K. A., Charles C. D. and Ninnemann U. S. (2003): Pleistocene vertical 752 carbon isotope and carbonate gradients in the South Atlantic sector of the Southern Ocean. 753 Geochemistry, Geophysics, Geosystems 4(1), CiteID 1004, DOI: 10.1029/2002GC000367. 754 Hughen K. A., Eglinton T. I., Xu L. and Makou M. (2004): Abrupt tropical vegetation 755 response to rapid climate changes. Science 304(5679), pp. 1955-1959. 756 IPCC (2007): Climate Change 2007: The Physical Science Basis. Contribution of Working 757 Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change 758 [Solomon, S., D. Qin, M. Manning, Z. Chen, M. Marquis, K.B. Averyt, M. Tignor and H.L. 759 Miller (eds.)]. Cambridge University Press, Cambridge, United Kingdom and New York, NY, 760 USA, 996 pp. 761 Joos F. and Spahni R. (2008): Rates of change in natural and anthropogenic radiative forcing 762 over the past 20,000 years. PNAS 105(5), pp. 1425-1430. 763 Jouzel J., Masson V., Cattani O., Falourd S., Stievenard M., Stenni B., Longinelli A., Johnsen 764 S. J., Steffenssen J. P., Petit J. R., Schwander J., Souchez R. and Barkov N. I. (2001): A New 765 27 Ky high resolution east Antarctic climate record. Geophysical Research Letters 28(16), pp. 766 3199-3202. 767 Jouzel J., Masson-Delmotte V., Cattani O., Dreyfus G., Falourd S., Hoffmann G., Minster B., 768 Nouet J., Barnola J. M., Chappellaz J., Fischer H., Gallet J. C., Johnsen S., Leuenberger M., 769 Loulergue L., Luethi D., Oerter H., Parrenin F., Raisbeck G., Raynaud D., Schilt A., 770 Schwander J., Selmo E., Souchez R., Spahni R., Stauffer B., Steffensen J. P., Stenni B., 771 Stocker T. F., Tison J. L., Werner M. and Wolff E. W. (2007): Orbital and Millennial 772 Antarctic Climate Variability over the Past 800,000 Years. Science 317(5839), pp. 793 - 796. 773 Kissel C., Laj C., Piotrowski A. M., Goldstein S. L. and Hemming S. R. (2008): Millennial- 774 scale propagation of Atlantic deep waters to the glacial Southern Ocean. Paleoceanography 775 23, PA2102, DOI: 10.1029/2008PA001624. 776 Knorr G. and Lohmann, G. (2003): Southern Ocean Origin for Resumption of Atlantic 777 Thermohaline Circulation during Deglaciation, Nature 424, pp. 532-536. 778 Köhler P. and Fischer H. (2004): Simulating changes in the terrestrial biosphere during the 779 last glacial/interglacial transition. Global and Planetary Change 43(1-2), pp. 33-55. 780 Köhler P., Fischer H., Munhoven G. and Zeebe R. E. (2005a): Quantitative interpretation of 781 atmospheric carbon records over the last glacial termination. Global Biogeochemical Cycles 782 19, GB4020, DOI: 10.1029/2004GB002345. 783 Köhler P., Joos, F., Gerber, S., Knutti, R.(2005b): Simulated changes in vegetation 784 distribution, land carbon storage, and atmospheric CO2 in response to a collapse of the North 785 Atlantic thermohaline circulation, Climate Dynamics, 25, pp. 689-708, DOI:10.1007/s00382- 786 005-0058-8. 787 Köhler P., Fischer H., Schmitt J. and Munhoven G. (2006a): On the application and 788 interpretation of Keeling plots in paleo climate research – deciphering δ13C of atmospheric 789 CO2 measured in ice cores. Biogeosciences 3, pp. 539-556. 790 Köhler P., Muscheler R. and Fischer H. (2006b): A model-based interpretation of low- 791 frequency changes in the carbon cycle during the last 120,000 years and its implications for 792 the reconstruction of atmospheric Δ14C. Geochemistry Geophysics Geosystems 7, Q11N06, 793 DOI: 10.1029/2005GC001228. 794 Köhler P. and Fischer H. (2006c): Simulating low frequency changes in atmospheric CO2 795 during the last 740,000 years. Climate of the Past 2(2), pp. 57-78. 796 Köhler P. and Bintanja R. (2008): The carbon cycle during the Mid Pleistocene transition: the 797 Southern Ocean decoupling hypothesis. Clim. Past 4, pp. 311-332. 798 Köhler P., Fischer H. and Schmitt J. (submitted): Atmospheric δ13CO2 and its relation to 799 pCO2 and deep ocean δ13C during the last Pleistocene. Paleoceanography, DOI: 800 10.1029/2008PA001703. 801 Labeyrie L. D., Duplessy J. C. and Blanc P. L. (1987): Variations in mode of formation and 802 temperature of oceanic deep waters over the past 125,000 years. Nature 327(6122), pp. 477- 803 482. 804 Lambert F., Delmonte B., Petit J. R., Bigler M., Kaufmann P. R., Hutterli M. A., Stocker T. 805 F., Ruth U., Steffensen J. P. and Maggi V. (2008): Dust-climate couplings over the past 806 800,000 years from the EPICA Dome C ice core. Nature 452(7187), pp. 616-619. 807 Lavrič J. V., Lourantou A., Barnola J.-M., Michel E., Raynaud D. and Chappellaz J (in prep.): 808 Measurement of carbon isotope composition and mixing ratio of CO2 in ancient air from ice 809 core samples. 810 Lea D. W., Pak D. K., Peterson L. C. and Hughen K. A. (2003): Synchroneity of tropical and 811 high-latitude Atlantic temperatures over the last glacial Termination. Science 301(5638), pp. 812 1361 - 1364. 813 Leuenberger M., Siegenthaler U. and Langway C. C. (1992): Carbon isotope composition of 814 atmospheric CO2 during the last ice age from an Antarctic ice core. Nature 357(6378), pp. 815 488-490. 816 Levitus and Boyer. (1994): World Ocean Atlas Vol 4: Temperature. 817 Loulergue L., Parrenin F., Blunier T., Barnola J.-M., Spahni R., Schilt A., Raisbeck G. and 818 Chappellaz J. (2007): New constraints on the gas age-ice age difference along the EPICA ice 819 cores, 0–50 kyr. Clim. Past 3(3), pp. 527-540. 820 Loulergue L., Schilt A., Spahni R., Masson-Delmotte V., Blunier T., Lemieux B., Barnola J.- 821 M., Raynaud D., Stocker T. F. and Chappellaz J. (2008): Orbital and millennial-scale features 822 of atmospheric CH4 over the past 800,000 years. Nature 453(7193), pp. 383-386. 823 MacDonald G. M., Beilman D. W., Kremenetski K. V., Sheng Y., Smith L. C. and Velichko 824 A. A. (2006): Rapid early development of circumarctic peatlands and atmospheric CH4 and 825 CO2 variations. Science 314(5797), pp. 285-288. 826 Marchitto J. T. M., Curry W. B. and Oppo D. W. (1998): Millennial-scale changes in North 827 Atlantic circulation since the last Glaciation. Nature 393(6685), pp. 557-561. 828 Marchitto T. M., Lehman S. J., Ortiz J. D., Flückiger J. and Geen A. v. (2007): Marine 829 radiocarbon evidence for the mechanism of deglacial atmospheric CO2 rise. Science 830 316(5830), pp. 1456 - 1459. 831 Martin J. H. (1990): Glacial-Interglacial CO2 Change: The iron hypothesis Paleoceanography 832 5(1), 1-13. 833 McManus J. F., Francois R., Gherardi J.-M., Keigwin L. D. and Brown-Leger S. (2004): 834 Collapse and rapid resumption of Atlantic meridional circulation linked to deglacial climate 835 changes. Nature 428(6985), pp. 834-837. 836 Meissner K. J., Schmittner A., Weaver A. J. and Adkins J. F. (2003): Ventilation of the North 837 Atlantic Ocean during the Last Glacial Maximum: A comparison between simulated and 838 observed radiocarbon ages. Paleoceanography 18(2), ID1023, DOI: 10.1029/2002PA000762. 839 Menviel L., Timmermann A., Mouchet A. and Timm O. (2008): Climate and marine carbon 840 cycle response to changes in the strength of the Southern Hemispheric westerlies. 841 Paleoceanography 23, PA4201, DOI: 10.1029/2008PA001604. 842 Monnin E., Indermühle A., Dällenbach A., Flückiger J., Stauffer B., Stocker T. F., Raynaud 843 D. and Barnola J.-M. (2001): Atmospheric CO2 concentrations over the last glacial 844 Termination. Science 291(5501), pp. 112 - 114. 845 Montenegro A., Eby M., Kaplan J. O., Meissner K. J. and Weaver A. J. (2006): Carbon 846 storage on exposed continental shelves during the glacial-interglacial transition. Geophysical 847 Research Letters 33, L08703, DOI: 10.1029/2005GL025480. 848 NGRIP Members (2004): High-resolution record of Northern Hemisphere climate extending 849 into the last interglacial period. Nature 431(7005), pp. 147-151. 850 Ninnemann U. S. and Charles C. D. (1997): Regional differences in Quaternary subantarctic 851 nutrient cycling: link to intermediate and deep water ventilation. Paleoceanography 12(4), pp. 852 560–567. 853 Obata A. (2007): Climate carbon cycle model response to freshwater discharge into the North 854 Atlantic. Journal of Climate 20(24), pp. 5962-5976. 855 Paillard D., Ghil M. and Treut H. L. (1993): Dissolved organic mater and the glacial- 856 interglacial pCO2 problem. Global Biogeochemical Cycles 7(4), pp. 901-914. 857 Paillard D., Labeyrie L. and Yiou P. (1996): Macintosh program performs time-series 858 analysis. EOS, transactions American Geophysical Union 77(39), p. 379. 859 Parrenin F., Barnola J.-M., Beer J., Blunier T., Castellano E., Chappellaz J., Dreyfus G., 860 Fischer H., Fujita S., Jouzel J., Kawamura K., Lemieux-Dudon B., Loulergue L., Masson- 861 Delmotte V., Narcisi B., Petit J.-R., Raisbeck G., Raynaud D., Ruth U., Schwander J., Severi 862 M., Spahni R., Steffensen J. P., Svensson A., Udisti R., Waelbroeck C. and Wolff E. (2007): 863 The EDC3 chronology for the EPICA Dome C ice core. Clim. Past 3(3), pp. 485-497. 864 Pataki D.E., Ehleringer J.R., Flanagan L.B., Yakir D., Bowling D.R., Still C.J., Buchmann N., 865 Kaplan J.O. and Berry J.A. (2003): The application and interpretation of Keeling plots in 866 terrestrial carbon cycle research. Global Biogeochemical Cycles, 17(1), 1022, DOI: 867 10.1029/2001GB001850. 868 Pépin L., Raynaud D., Barnola J.-M. and Loutre M. F. (2001): Hemispheric roles of climate 869 forcings during glacial-interglacial transitions as deduced from the Vostok record and LLN- 870 2D model experiments. Journal of Geophysical Research 106, D23, pp. 31885–31892. 871 Rasmussen S. O., Andersen K. K., Svensson A. M., Steffensen J. P., Vinther B. M., Clausen 872 H. B., Siggaard-Andersen M.-L., Johnsen S. J., Larsen L. B., Dahl-Jensen D., Bigler M., 873 Röthlisberger R., Fischer H., Goto-Azuma K., Hansson M. E. and Ruth U. (2006): A new 874 Greenland ice core chronology for the last glacial termination. Journal of Geophysical 875 Research 111, D06102, DOI: 10.1029/2005JD006079. 876 Rickaby R. E. M., and H. Elderfield (2005): Evidence from the high-latitude North Atlantic 877 for variations in Antarctic Intermediate water flow during the last deglaciation. Geochemistry 878 Geophysics Geosystems 6, Q05001, DOI: 10.1029/2004GC000858. 879 Roethlisberger R., Mulvaney R., Wolff E. W., Hutterli M. A., Bigler M., Sommer S. and 880 Jouzel J. (2002): Dust and sea salt variability in central East Antarctica (Dome C) over the 881 last 45 kyrs and its implications for southern high-latitude climate. Geophysical Research 882 Letters 29(20) pp. 24-1, Cite ID 1963, DOI: 10.1029/2002GL015186. 883 Schmittner A., Brook E. and Ahn J. (2007): Impact of the Ocean’s Overturning Circulation on 884 Atmospheric CO2. AGU Geophysical Monograph Series 173, pp. 209-246. 885 Scholze M., Knorr W. and Heimann M. (2003): Modelling terrestrial vegetation dynamics and 886 carbon cycling for an abrupt climatic change event. The Holocene 13(3), pp. 327-333. 887 Schulz M., Seidov D., Sarnthein M. and Stattegger K. (2001): Modeling ocean-atmosphere 888 carbon budgets during the Last Glacial Maximum - Heinrich 1 Meltwater Event - Bølling 889 transition. Int. Journ. Earth Sciences 90, pp. 412-425. 890 Schwander J., Barnola J.-M., Andrié C., Leuenberger M., Ludin A., Raynaud D. and Stauffer 891 B. (1993): The age of the air in the firn and the ice at Summit, Greenland. Journal of 892 Geophysical Research 98(D2), pp. 2831–2838. 893 Severinghaus J. P., Grachev A. and Battle M. (2001): Thermal fractionation of air in polar firn 894 by seasonal temperature gradients. Geochemistry, Geophysics, Geosystems 2(7) 1048, DOI: 895 10.1029/2000GC000146. 896 Sigman D., and E. Boyle (2000): Glacial/interglacial variations in atmospheric carbon 897 dioxide. Nature 407(6806), pp. 859-869. 898 Smith H. J., Fischer H., Wahlen M., Mastroianni D. and Deck B. (1999): Dual modes of the 899 carbon cycle since the Last Glacial Maximum. Nature 400(6741), pp. 248-250. 900 Spahni R., Schwander J., Flückiger J., Stauffer B., Chappellaz J. and Raynaud D. (2003): The 901 attenuation of fast atmospheric CH4 variations recorded in polar ice cores. Geophysical 902 Research Letters 30(11), 1571, DOI: 10.1029/2003GL017093. 903 Spero H. J. and D. W. Lea (2002): The cause of carbon isotope minimum events on glacial 904 terminations. Science, 296(5567), pp. 522-525. 905 Stephens B. B., and R. F. Keeling (2000): The influence of Antarctic sea ice on glacial- 906 interglacial CO2 variations. Nature, 404(6774), pp. 171-174. 907 Stott L., J. Southon, A. Timmermann, and A. Koutavas (2009): Radiocarbon age anomaly at 908 intermediate water depth in the Pacific Ocean during the last deglaciation, Paleoceanography, 909 24, PA2223, DOI: 10.1029/2008PA001690 910 Toggweiler J. R. (1999): Variation of atmospheric CO2 by ventilation of the ocean's deepest 911 water. Paleoceanography 14(5), pp. 571-588. 912 Toggweiler J. R., Russell J. L. and Carson S. R. (2006): Midlatitude westerlies, atmospheric 913 CO2, and climate change during the ice ages. Paleoceanography 21, PA2005, DOI: 914 10.1029/2005PA001154. 915 Trudinger C. M., Enting L. G., Etheridge D. M., Francey R. J., Levchenko V. A., Steele L. P., 916 Raynaud D. and Arnaud L. (1997): Modeling air movement and bubble trapping in firn. 917 Journal of Geophysical Research 102(D6), pp. 6747–6763. 918 Tschumi J. and Stauffer B. (2000): Reconstructing past atmospheric CO2 concentration based 919 on ice-core analyses: open questions due to in situ production of CO2 in the ice. Journal of 920 Glaciology 46(152), pp. 45-53. 921 Visser K., Thunell R. and Stott L. (2003): Magnitude and timing of temperature change in the 922 Indo-Pacific warm pool during deglaciation. Nature 421(6919), pp. 152-155. 923 Watson A. J., and Naveira Garabato A. C. (2006): The role of Southern Ocean mixing and 924 upwelling in glacial-interglacial atmospheric CO2 change. Tellus B, 58(1), pp. 73-87. 925 926 Figure captions 927 928 Figure 1. 929 CO2 and δ13CO2 evolution during the last deglaciation from the EPICA Dome C (EDC) ice 930 core, superimposed with other ice core data: (a) δD of ice in EDC (grey line, [Jouzel et al., 931 2007]), averaged over 500y (black line, [EPICA, 2006]) (b) δ18O of NGRIP ice [NGRIP, 932 2004], with a running average over 500y (dark grey line), (c) atmospheric CH4 mixing ratio 933 (red line and triangles: EDC [Loulergue et al., 2008]; green dots: TD [Brook et al., 2000]), (d) 934 atmospheric CO2 mixing ratio (red line and dots: EDC [Monnin et al., 2001]; blue line and 935 diamonds: this study; green dots: TD [Smith et al., 1999]) and (e) δ13CO2 data (blue line and 936 diamonds: this study; green dots: TD [Smith et al., 1999]). When duplicate measurements 937 were performed, the line runs through the mean. The dotted blue lines in (e) correspond to the 938 1σ (0.1 ‰ average) uncertainty envelope. The two blue open diamonds indicate a suspicious 939 result that we discarded in the discussion. 940 All gas records are plotted versus the Greenland GICC05 age scale. The upper x-axis 941 represents the EDC depth for the gas records. δD is plotted on a chronology combining the 942 CH4 fit to GICC05 and the EDC3 Δage [EPICA, 2006]. The vertical dotted lines correspond 943 to boundaries between different CO2 rates of change during the deglaciation, as defined by 944 [Monnin et al., 2001], adapted to the new age scale. The time periods in-between are noted 945 SP-I to SP-IV. YD: Younger Dryas; B/A: Bølling/Allerød; ACR: Antarctica Cold Reversal; 946 H1: Heinrich 1. 947 948 Figure 2. 949 Comparison of atm. δ13CO2 data with C-cycle related tracers during TI: (a) dust concentration 950 in the EDC core, taken as proxy of the biological pump G-IG patterns in the S. ocean 951 [Lambert et al., 2008]; (b) atmospheric δ13CO2, from this study; (c) 952 subtropical North Atlantic, a tracer of North Atlantic Deep Waters formation strength 953 [McManus et al., 2004]; (d) benthic δ13C in North Atlantic intermediate waters, reflecting the 954 relative contribution between NADW and Antarctic intermediate waters throughout TI in the 955 N. Atlantic basin [Rickaby & Elderfield, 2005]; (e) Δ14C data of intermediate waters in North 956 Pacific, a proxy for S. Ocean overturning strength [Marchitto et al., 2007]; (f) planktonic δ13C 957 data from the western tropical Pacific [Stott et al., 2009], also depicting changes in S. Ocean 958 G-IG overturning changes; (g) opal flux data from the Atlantic sector of the S. Ocean, a proxy 959 for S. Ocean upwelling [Anderson et al., 2009]. Atmospheric δ13CO2 is plotted versus the 960 Greenland GICC05 age scale, while the dust data are presented on the GICC05 ice scale. The 961 oceanic proxies are on their original time scale. Shaded parts represent the cold periods of the 962 North Hemisphere, as deduced from the individual time scales for each proxy. 963 964 231 Pa/230Th in the 965 Figure 3. 966 Proxy data sets used as BICYCLE input parameterizations. Shadings highlight the definition 967 of sub-periods given in the main text. 968 A: Coral reef terraces as indicator for sea level rise [Fairbanks, 1989]. 969 B: NorthGRIP δ18O as northern high latitude temperature proxy [NGRIP Members, 2004]. 970 C: Changes in equatoral SST [Visser et al., 2003] and deep ocean temperature in different 971 oceanic compartments [Labeyrie et al., 1987]. 972 D: EDC δD as southern high latitude temperature proxy [Jouzel et al., 2001]. 973 E: EDC non-seasalt (nss) dust as proxy of aeolian iron input into the S. Ocean [Roethlisberger 974 et al., 2002]. 975 F: Assumed changes in strengths of the main oceanic currents. 976 For B, D and E cases, data show large short term fluctuations; therefore a 500-y running mean 977 is used in the simulations. 978 All ice core records (B, D, E) are plotted versus the GICC05 age scale. 979 980 Figure 4 981 Mixing diagram depicting the relationship between atmospheric δ13CO2 and the inverse of 982 CO2 (Keeling plot). The new data are shown as open circles with different colors, 983 corresponding to the different sub-periods defined in the main text. BICYCLE model ouput is 984 represented by open black diamonds. Data points used for calculating the regression lines 985 corresponding to the two periods of abrupt CO2 rise and δ13CO2 decline are filled with light 986 blue (for the first 987 represent model results used for plotting the corresponding regression black lines. The y- 988 intercept values are shown next to the regression lines, together with the number data points. 13 C dip) and dark blue (for the second 13 C dip). Filled black diamonds 989 The y-intercepts of the two rapid δ13CO2 decreases give reasonably consistent values of ~- 990 11‰, comparable with the model results. 991 992 Figure 5. 993 Comparison between EDC CO2 and δ13CO2 data and box-model simulations. 994 Left panels: 995 Imprint of individual major C-cycle processes on atmospheric (a) CO2 and (b) δ13CO2, 996 simulated with the BICYCLE model. All curves express an anomaly ΔpCO2 and Δδ13CO2 997 versus a reference corresponding to boundary EH conditions. The following processes are 998 shown at this point: 999 (1) S. Ocean mixing; (2) marine productivity; (3) ocean temperature and (4) terrestrial 1000 biosphere. 1001 Right panels: 1002 Superposition of the BICYCLE simulation integrating all individual processes of the left 1003 panels (grey line) with our data (deep blue line and diamonds). The equilibrium-state 1004 BOXKIT model outputs, using similar boundary conditions as BICYCLE for each time 1005 period (red triangles) are also plotted for (c) CO2 and (d) δ13CO2. Red squares correspond to 1006 BOXKIT simulations using higher equatorial SST magnitudes. All series are plotted versus 1007 GICC05 age scale. 1008 1009 Tables 1010 1011 Table 1: Blank tests results of the experimental setup on standard gas, with their 1σ standard 1012 deviation and the number of tests. 1013 Test II III CO2 (ppmv) 261.1 ± 1.2 261.4 ± 1.8 δ13CO2 (‰) -6.4 ± 0.1 -6.7 ± 0.1 tests number (n) 5 14 1014 1015 1016 Table 2: Tie-points between the EDC3_gas_a [Loulergue et al., 2007] and GICC05 time 1017 scales [EPICA Community Members, 2006; Andersen et al., 2007; NGRIP Members, 2004] 1018 using the ANALYSERIES software [Paillard et al., 1996] 1019 EDC3_gas-a age (y BP) 7890 11330 11920 12340 13070 14010 15870 17790 19690 21220 GICC05 gas age (y BP) 8240 11680 12330 12790 13600 14640 16200 17800 19670 21100 CH4 value (ppbv) 590 560 460 540 670 570 470 370 350 350 1020 event description CH4 minimum during Holocene CH4 mid-rise / ending of YD CH4 YD minimum CH4 mid decrease / ending of B/A B/A CH4 peak CH4 mid-rise / towards B/A CH4 peak CH4 drop CH4 drop CH4 drop during LGM 1021 Tie-points between TD and GICC05 age scales, using the same software of [Paillard et al., 1022 1996]. The TD-core was initially plotted versus GISP2 age scale [Brook et al., 2000]. GISP2 1023 is almost synchronous to GICC05; still, for the LGM time-period, GISP2 had to be rescaled: 1024 TD gas age GICC05 gas CH4 value (ppbv) (y BP) age (y BP) 8300 8290 570 11690 11660 660 11890 11860 430 12910 12830 600 13570 13430 670 14880 14790 510 16770 16200 500 26470 22810 420 1025 1026 event description CH4 minimum during Holocene CH4 peak after YD CH4 YD minimum CH4 mid decrease / ending of B/A B/A CH4 peak Just before the B/A CH4 rise CH4 peak before B/A CH4 peak during LGM sea level [m] 0 EH IV III II I LGM -20 reef terraces -40 -60 -80 -100 -120 -140 NorthGRIP A o D [ /oo] C -380 -390 -400 -410 -420 -430 -440 -450 EQ SST Deep Atlantic Deep Southern Ocean Deep Pacific D E Flux [Sv] 20 EDC EDC 0 10 20 30 40 50 60 NADW AABW_A AABW_P SOX 15 10 5 0 F 8 10 12 14 16 18 20 22 Age [kyr BP] [ppb] -44 nss-Ca 1 0 -1 -2 -3 -4 -5 -42 o -40 B T [ C] -38 2+ o -36 18 O [ /oo] -34 1 Supporting non-print material 2 3 Keeling plot as a function of each sub-period 4 We provide a Keeling plot where our data are distinguished as a function of each sub- 5 period (open circles and continuous lines, Fig. S1). We find similar regressions for each sub- 6 period between our data and Smith et al., 1999 data (crosses and dotted lines, Fig. S1). The y- 7 intercepts are of little use in the case of SP-II, SP-III and SP-IV, due to the very large slope. 8 This highlights again the limit of Keeling plots when used in such context. 9 10 BICYCLE model updates 11 The overall model configuration is seen in fig. S2a, while fig. S2b illustrates in detail 12 the oceanic boxes interactions. Here we use the model version, where terrestrial net primary 13 productivity is more influenced by climate (notably temperature) change than by CO2 14 fertilization. This version is labeled “TB2” in former applications [Köhler & Fischer, 2004] 15 (from now on called “GBC2005-version”). 16 Model outputs are given in atmospheric partial pressure (pCO2) in µatm units, which, 17 only in dry air and at standard pressure conditions are identical to ice core nomenclature 18 (ppmv); we assume equality between the two, neglecting a relatively constant offset between 19 both quantities of a few ppmv. 20 Here an update of the “GBC2005” version is presented, named (“GBC2009”). The 21 modifications mainly concern oceanic processes 22 (a) Carbonate compensation is represented by a relaxation function, which brings deep ocean 23 carbonate ion concentration, CO32-, back to initial values after every perturbation. This 24 relaxation operates with a time delay (e-folding time) of τ = 1.5 ky to account for the 25 relatively slow processes in the sediments. This value of τ was chosen based on reconstructed 26 deep ocean carbonate ion dynamics [Marchitto et al., 2005]. Details of the approach are 27 described in [Köhler & Fischer, 2006]. 28 (b) In terms of ocean circulation changes, we assumed a reduction of NADW by 8 Sv during 29 H1, instead of the complete shutdown initially proposed [Köhler et al., 2005]. We also 30 changed AABW in antiphase with NH deep waters, in order to better represent the bipolar 31 seesaw in the simulations. Furthermore, 30% of upwelling fluxes in the S. Ocean are directly 32 redistributed to the Atlantic intermediate waters. 33 Fig. S3 illustrates the new GBC2009 integrated result against our data and the older 34 GBC2005 model simulations, (a) for CO2 and (b) for δ13CO2. This new GBC2009 version of 35 BICYCLE leads to the following improvements: 36 1. LGM and EH boundary values are more consistent with our data 37 2. The timing and trend of the early deglaciation are better reconstructed, due to updated 38 runs using the most recent age scale GICC05. Still, inconsistencies exist for the last two SPs 39 of TI, as for the case of TB2; the model trends lead the data. 40 3. The magnitudes of both CO2 and δ13CO2 changes throughout the different SPs are 41 more consistent with our data. 42 The largest discrepancy between our data and BICYCLE new version lies in the EH. It could 43 be explained by inadequate ocean circulation parameterisation: for instance AABW gets 44 stable in BICYCLE simulations at the end of YD, whereas it should have been enhanced at 45 that time, as for NADW. 46 Fig. S4 finally shows the contribution of the three less important forcing factors towards the 47 atmospheric CO2 and δ13CO2, and are already been commented in the main text. 48 49 50 51 Figure captions 52 53 Figure S1 54 Mixing diagram depicting the relationship between atmospheric δ13CO2 and the inverse of 55 CO2 (Keeling plot). The new data are shown as open circles with different colors, 56 corresponding to the different sub-periods defined in the main text. The data from Smith et 57 al., 1999 are provided as crosses with a color coding similar to our data, for each sub-period 58 they belong to. Regression lines for both datasets are plotted as well for each sub-period 59 (continuous lines and dashed lines respectively). 60 61 Figure S2 62 Sketch of the “Box model of the Isotopic Carbon cYCLE” (BICYCLE). 63 (a) Overall model setup. 64 Compartments in the terrestrial biosphere distinguish different primary production schemes 65 for grasses (C3, C4), non woody (NW) and woody (W) biomass of tree, detritus (D) and fast 66 and slow (FS, SS) decomposing soils. 67 (b) Close-up on the definition of ocean boxes and the circulation scheme, fluxes quantified for 68 the pre-industrial period (PRE). 69 Fluxes are given in Sverdrup (1Sv = 106 m3/s) 70 71 Figure S3 72 Comparison of both BICYCLE model versions (black lines) to our dataset (blue diamonds). 73 The dotted line corresponds to the initial BICYCLE (GBC2005) simulation throughout TI, as 74 described in the main manuscript. The straight line corresponds to our modified BICYCLE 75 GBC2009-configuration. (a) gives the result on CO2 and (b) for δ13CO2. Both plots are on the 76 GICC05 age scale. 77 78 Figure S4 79 Impact of less important forcing factors on the atmospheric (a) CO2 and (b) δ13CO2 signal: (5) 80 sea level rise; (6) sea ice retreat and (7) Northern (NADW) and Southern (AABW)- sourced 81 deep water fluxes changes, throughout TI. All curves express an anomaly ΔpCO2 and 82 Δδ13CO2 versus a reference representing the boundary EH conditions. 83 84 Figures 85 Figure S1 86 87 -6.2 δ13CO2 (‰) -6.4 -6.6 -6.8 LGM SP-I SP-II SP-III SP-IV EH -7 -7.2 0.0036 88 89 90 91 92 93 94 95 96 97 98 99 0.004 0.0044 0.0048 1/CO2 (ppmv ) -1 0.0052 0.0056 Sm99-LGM Sm99-I Sm99-II Sm99-III Sm99-IV Sm99-EH 100 101 102 103 104 105 106 107 108 109 Figure S2 110 Figure S3 111 112 IV III II I -6.2 b -6.3 -6.4 -6.6 -6.7 δ13C (‰) -6.5 -6.8 -6.9 -7 a 260 240 220 200 180 EDC - points GBC2005 run GBC2009 run 8 113 114 115 116 117 118 119 120 121 10 12 160 14 16 Age (ky BP) 18 20 22 pCO2 or CO2 (µatm or ppmv) -7.1 280 122 Figure S4 123 124 0.5 0.4 b IV III II I 0.3 Δδ13C (‰) 0.2 6 7 5 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 20 a 5 10 6 ΔpCO2 (µatm) 0 -10 -20 7 -30 -40 Sea level Sea ice NADW+AABW -50 -60 8 10 12 14 16 18 20 22 Age (ky BP) 125 126 127 Supplementary References 128 Köhler P. and Fischer H. (2004): Simulating changes in the terrestrial biosphere during the 129 last glacial/interglacial transition. Global and Planetary Change 43(1-2), pp. 33-55. 130 Köhler P., Fischer H., Munhoven G. and Zeebe R. E. (2005): Quantitative interpretation of 131 atmospheric carbon records over the last glacial termination. Global Biogeochemical Cycles 132 19, GB4020, DOI: 10.1029/2004GB002345. 133 Köhler P. and Fischer H. (2006): Simulating low frequency changes in atmospheric CO2 134 during the last 740,000 years. Climate of the Past 2(2), pp. 57-78. 135 Marchitto T. M., Lynch-Stieglitz J. and Hemming S. R. (2005): Deep Pacific CaCO3 136 compensation and glacial–interglacial atmospheric CO2. Earth and Planetary Science Letters 137 231(3-4), pp. 317-336. 138 Smith H. J., Fischer H., Wahlen M., Mastroianni D. and Deck B. (1999): Dual modes of the 139 carbon cycle since the Last Glacial Maximum. Nature 400(6741), pp. 248-250. 140 View publication stats