Climate of the Past Reconstruction of a continuous high-resolution CO 2 record over the past 20 million years

The gradual cooling of the climate during the Cenozoic has generally been attributed to a decrease in CO2 concentration in the atmosphere. The lack of transient climate models and, in particular, the lack of high-resolution proxy records of CO2, beyond the ice-core record prohibit, however, a full understanding of, for example, the inception of the Northern Hemisphere glaciation and mid-Pleistocene transition. Here we elaborate on an inverse modelling technique to reconstruct a continuous CO2 series over the past 20 million year (Myr), by decomposing the global deep-sea benthic δ18O record into a mutually consistent temperature and sea level record, using a set of 1-D models of the major Northern and Southern Hemisphere ice sheets. We subsequently compared the modelled temperature record with ice core and proxy-derived CO2 data to create a continuous CO2 reconstruction over the past 20Myr. Results show a gradual decline from 450 ppmv around 15Myr ago to 225 ppmv for mean conditions of the glacial-interglacial cycles of the last 1Myr, coinciding with a gradual cooling of the global surface temperature of 10K. Between 13 to 3Myr ago there is no long-term sea level variation caused by ice-volume changes. We find no evidence of change in the long-term relation between temperature change and CO2, other than the effect following the saturation of the absorption bands for CO2. The reconstructed CO2 record shows that the Northern Hemisphere glaciation starts once the long-term average CO2 concentration drops below 265 ppmv after a period of strong decrease in CO2. Finally, only a small long-term decline of Correspondence to: R. S. W. van de Wal (r.s.w.vandewal@uu.nl) 23 ppmv is found during the mid-Pleistocene transition, constraining theories on this major transition in the climate system. The approach is not accurate enough to revise current ideas about climate sensitivity.


Introduction
The overall climate cooling reconstructed for the past 20 Myr has generally been attributed to a change in CO 2 concentration in the atmosphere (Zachos et al., 2008;Ruddiman, 2003), although the amount of CO 2 decreases and the amplitude of subsequent cooling are discussed widely (Jansen et al., 2007).Since data and modelling studies covering this time period are poorly integrated, our understanding of the inception of ice ages in the Northern Hemisphere (NH) (Raymo 1994), as well as the mechanisms causing the transition from 41 000-yr to 100 000-yr dominated climate cycles (Tziperman and Gildor, 2003;Clark et al., 2006;Huybers, 2007;Bintanja and Van de Wal, 2008), that occurred without apparent changes in the insolation forcing (Hays et al., 1976;Imbrie and Imbrie, 1980) is still incomplete.Current difficulties in assessing the role of CO 2 on the long timescales are the lack of reliable CO 2 data from the pre ice-core record (Ruddiman, 2010), and the limited data of sea level (Miller et al., 2005;Müller et al., 2008) and temperature (De Boer et al., 2010).Our current knowledge of long-term climate variability builds on the Milankovitch theory of solar-insolation variability (Milankovitch 1941), including scenarios that rely on highly parameterized nonlinear response mechanisms to the insolation forcing.Recent developments in the interpretation of marine δ 18 O records and new CO 2 proxies allow R. S. W. van de Wal et al.: Reconstruction of a continuous high-resolution CO 2 record us to reassess this understanding and to present a global overview of temperature, sea level and CO 2 changes over time.
We build on a model set-up that aims to integrate the main climate variables, temperature and sea level, of the glacialinterglacial variability over the past 20 Myr (Bintanja et al., 2005a, b;Bintanja and van de Wal, 2008;De Boer et al., 2010, 2011) to reconstruct changes in past CO 2 concentrations.The model takes advantage of the mass conservation of δ 18 O on a global scale, while an inverse routine guarantees that changes in ocean δ 18 O caused by both land ice sheet growth (sea-level change) and deep ocean temperature change are in agreement with marine benthic δ 18 O reconstructions.Output of the model is ice volume and the change in NH air temperature (T NH ) relative to the present day, calculated from the deep ocean temperature, which gives an estimate of the air temperatures at sea level in areas where the NH ice sheets developed (40 • -80 • N).Here we will use the T NH derived by De Boer et al. (2010) to reconstruct a continuous CO 2 record over the past 20 Myr.
For this purpose, we compared the T NH record, with existing proxies for CO 2 , first for the past 800 000 yr using primarily ice-core CO 2 data (Petit et al., 1999;Siegenthaler et al., 2005;Lüthi et al., 2008), and secondly with geochemical and stomata-derived CO 2 proxy data for the older time interval (Tripati et al., 2009;Kürschner et al., 1996Kürschner et al., , 2008;;Retallack 2009;Pearson and Palmer, 2000;Hönisch et al., 2009;Pagani et al., 2005Pagani et al., , 2009;;Seki et al., 2010).Through this comparison, a regression coefficient between temperature and CO 2 could be determined that allowed us to reconstruct a global mutually self-consistent and continuous overview of temperature, sea level and CO 2 over the past 20 Myr.Eventually the paper is concluded by a discussion on climate sensitivity following from our results.

Inverse δ 18 O modelling approach
The inverse modelling approach enables the deep-sea benthic δ 18 O record to be decomposed in a temperature and ice-volume component, using an ice-sheet model.Key processes in the ice-sheet model are a variable isotopic sensitivity and isotopic lapse rate, the mass balance height feedback, the mass balance albedo feedback and the adjustment of the underlying bedrock.In the early stages, the applied model was forced with sea-level reconstructions over the past 120 000 yr (Bintanja et al., 2005a).Through the inverse routine, temperature was adjusted so that modelled ice volume matched the observations.This constraint ensured that sea level and temperature are mutually consistent.In addition, it allowed a quantification of model errors, and errors arising from the uncertainty in the sea-level observations or reconstructions.Results have been compared favourably with sea-level (Rohling et al., 2009;Lambeck and Chapell, 2001) and temperature data (Lear, 2000).The paper by Bintanja et al. (2005a) yielded a T NH reconstruction, which is highly coherent with the classical Vostok temperature record (Petit et al., 1999).Nevertheless, an obvious limitation of this work was that global sea-level observations are limited to the last 0.5 Myr.
In subsequent studies, we, therefore, started to incorporate the marine benthic δ 18 O record as forcing (e.g.Bintanja et al., 2005b).This was achieved by taking advantage of mass conservation of δ 18 O on the global scale.First, it was applied to calculate temperature and sea level over the past million years (Bintanja et al., 2005b) and later to explore the mechanisms of the Mid-Pleistocene Transition (Bintanja and Van de Wal, 2008), both focusing on the climate in the NH, as only the Eurasian and North American ice sheet complexes were modelled explicitly.Their temperature reconstruction was compared to alkenone-derived equatorial temperatures (Lawrence et al., 2006) for the last 3 Myr, indicating similar strength for most of the glacials.In order to use the benthic δ 18 O record as forcing, a simple parameterization was used to separate deep-water temperature from icevolume changes, based on a linear relation between deepwater temperature and δ 18 O (Duplessy et al., 2002), and an idealised climate model (Bintanja and Oerlemans, 1996) relating changes in deep-ocean temperature to atmospheric temperature changes (see Bintanja et al., 2005bBintanja et al., , 2008 for details) for details).
The last step in the model sequence until now is the explicit inclusion of ice sheets in the Southern Hemisphere (SH), allowing a longer time span to be covered, since for warmer conditions ice-volume changes are dominated by changes in the SH (De Boer et al., 2010).This has been done at the expense of the complexity of the icesheet models used.Hence, five 1-D ice sheets, rather than the two 3-D ice-sheet models used previously (Bintanja et al., 2005b;Bintanja and Van de Wal, 2008) were used to reconstruct temperature and sea level over 20 Myr (De Boer at al., 2010).The five 1-D ice-sheet models simulate ice flow over a cone shaped continent.(De Boer et al., 2010).They represent glaciation in Eurasia, North America, Greenland and East and West Antarctica, further abbreviated to NAIS, EAS, GrIS, EAIS, WAIS, respectively, where each has a different geometry, mass balance forcing and isotopic content.Differences between NH and SH temperatures are constant in time and calibrated to match a strong ice volume increase in Antarctica at the Eocene-Oligocene transition (De Boer et al., 2010).
Results of these late Cenozoic runs are compared to sea level by De Boer et al. (2010) with sea level estimates by Miller et al. (2005) and Müller et al. (2008), suggesting considerable differences between 10 and 35 Myr ago.Over the past 10 Myr similarities are high, as the Miller record is derived from the δ 18 O by scaling.In terms of change in seawater isotope per m sea level change (δw) our results are in good agreement to commonly used values in the literature (De Boer et al., 2011).However, it should be noted that our modelling approach allows for temporal variations in δw, which is usually not the case.Results of the deep water temperature change since the Miocene by De Boer et al. (2011) are comparable to the results by Lear et al. (2000) based on the Mg/Ca proxy if we use the paleo temperature equation from Shackleton (1974) with a δw present-day value of −0.28 per mille (VPDB).Data by Lear et al. (2010) indicate, however, considerably higher values for the Middle Miocene.
In this model, T NH , determines (1) the growth and retreat of the NAIS, EAS and GrIS, (2) changes in deep-water temperature from which it has been derived, and (3) the SH air temperatures which determine the growth and retreat of the EAIS and WAIS.Hence, besides a simple parameterization to relate deep-water temperature to atmospheric temperature (Bintanja et al., 2005b), we included a constant value to set the temperature difference between the NH and SH.Both parameterizations contribute to the uncertainty of the model as will be explained later.The conceptual approach used here was developed for orbital timescales.Thus, the antiphase dynamics of temperature in the northern and southern high latitudes as observed for the bipolar seesaw (e.g.Barker et al., 2009) is not embedded here, neither are Dansgaard/Oeschger events resolved.Far more detailed models are requested to address these temporal patterns, but the lack of high temporal resolution is not believed to affect the general temporal trends of the orbital timescales addressed here.
In the ice-sheet model, isotopic content and ice volume are calculated with a time step of 1 month and are implemented every 100 yr in the ocean isotope module.Every 100 yr, the modelled benthic isotope is evaluated and forwarded to calculate the temperature anomaly for the next time step (Bintanja et al., 2005b).As forcing we use the stacked benthic δ 18 O record of Zachos et al. (2008), which is smoothed and interpolated to obtain a continuous record with a resolution of 100 yr.This implies that the timescale of the reconstruction is implicitly determined by the benthic record.The methodology ensures that the phasing between temperature and sea level is consistent with respect to the benthic δ 18 O data.Further details and a more thorough model description are presented by De Boer et al. (2010).

Results in terms of sea level and temperature trends
Our model-based deconvolution shows a long-term decrease in T NH by 27 K since the Miocene (about 10 K in the global surface temperature) with superimposed orbitally forced changes (Fig. 1b).Eustatic sea level, more strictly sea level from ice-volume changes only, gradually falls, but is roughly constant from 13 Ma (+15 m) to 3 Ma (+5 m) as the ice sheets in the SH are full grown and major ice sheets in the NH are not yet developed (Fig. 1c).Moreover, the deviation of the sea level changes from the 400 kyr running mean revealed only low amplitude sea level changes of 10 m during this time period, whereas it fluctuated to ±20 m prior to 13 Ma For CO 2 the error bar is calculated as 400-kyr running mean, for the other records it is the standard deviation on the 0.1 kyr value as used in the model.Data are available as a Supplement to this paper.and up to ±66 m after 3 Ma.Maximum sea level high-stand of +55 m occurred around 15 Ma, probably caused by a reduced EAIS (De Boer et al., 2010).
Figure 2 shows that there is not a unique solution for sea level given a certain temperature.This results from the different timescales in the coupled system of ice sheets, changing deep-water temperatures, surface temperatures, bedrock adjustment, and forcing and feedbacks of the mass balance height and albedo-temperature feedback.Obviously, sea level rises on average with temperature as illustrated by the thick lines in Fig. 2a.On average the sea level change is  (Bintanja and Van de Wal, 2008;De Boer et al., 2010).The more sophisticated 3-D results are validated by observation of sea level (Lambeck and Chapell, 2001;Rohling et al., 2009).The sea level sensitivity to temperature change is approximately similar for warm (6.5 m K −1 ) and cold (7.8 m K −1 ) climates as indicated by the thick line in Fig. 2a.In addition, Fig. 2b shows the volume change for the individual ice sheets as a function of temperature leading by summation to the complex pattern in Fig. 2a.Also on the level of an individual ice sheet, transient effects impede a simple and unique solution between temperature and sea level, which implies that inverting climate information from sea level, has to be considered with great care.
In contrast to the sea level record, temperature shows a more gradual decline from the Miocene maximum around 15 Myr ago to the start of the major glaciation in the NH around 3 Ma.The gradual increase in the benthic δ 18 O record leads to a long-term cooling of the climate between 13 and 3 Ma.The amplitude of temperature and sea level variability both increase once the major ice sheets develop in the NH around 3 Myr ago (Fig. 1).
Many tests have been performed with the model to assess the uncertainties in the input and model parameters on sea level and temperature results.The most important tests allow us to estimate the uncertainty range displayed in Fig. 1.For the δ 18 O input we defined an uncertainty of 0.16 per mille, which is derived from the root-mean-squared difference between the smoothed marine record and the actual data points.The key model parameters contributing to the uncertainty are (1) the deep-water to surface-air temperature coefficient (range 0.15 to 0.25), (2) the temperature difference between T NH and the temperatures around Antarctica (range: 6-14 K for EAIS, range: 2-10 K for WAIS), and (3) the isotopic content of the ice sheets (range from −43, −32, −28 to −55, −42, −36 per mille, respectively, for EAIS, WAIS, GrIS; see De Boer et al., 2010 for details).For the three model parameters, maximum and minimum values are used to test the effect on modelled temperature and sea level.The resulting standard deviation varies over time, but is on average 1.9 K for temperature and 6.2 m for sea level over the past 20 Myr.In order to interpret the results, one has to bear in mind that the reconstructed temperatures are strictly only valid in the continental areas of the NH where ice sheets develop (about 40 • -80 • N), implying that they are, therefore, not a representative for the entire globe (T g ).A relationship between T NH and T g is developed later on.
2009), a combination of alkenones and δ 11 B (Seki et al., 2010), and ice cores (Petit et al., 1999;Siegenthaler et al., 2005;Lüthi et al., 2008), all shown in Fig. 3.We did not consider the paleosol record presented by Beerling and Royer (2011) because of their large uncertainty.All data points are representative for different discrete time intervals, with obviously a bias towards the more modern data points and each having its advantages and drawbacks.For example, the boron isotope derived estimates of the CO 2 concentration are based on the fact that higher atmospheric concentrations lead to more dissolved CO 2 in the surface ocean, which causes a reduction in the pH of the ocean.As the pH can be derived from measurements of the δ 11 B of calcium carbonate (Pearson and Palmer, 2000), CO 2 can be calculated provided that another parameter of the marine carbonate system (e.g.alkalinity) is known (Zeebe and Wolf-Gladrow,  Hönisch et al., 2009), B/Ca (Tripati et al., 2009), alkenones+ δ 11 B s (Seki et al., 2010), stomata data (Kürschner et al., 1996(Kürschner et al., , 2008;;Retallack 2009), and the ice-core record (Petit et al., 1999;Siegenthaler et al., 2005;Lüthi et al., 2008) are used for further analysis.
For reasons of transparency CO 2 is plotted in ppmv.If CO 2 would be plotted as ln(CO 2 /CO 2 , ref ) a similar picture emerges.The latter is physically more consistent as it takes the saturation of the absorption bands into account.

2001
).The method is expensive, time consuming and only well-preserved foraminiferal specimen are suitable for the analysis, resulting in only low-resolution records up to now.Ice cores provide the most robust and high-resolution CO 2 archive as they directly preserve the atmospheric concentrations, but only for the past 800 000 yr (Lüthi et al., 2008).
Here, we accept all data as they are published without any further correction.The general picture is that the scatter in the different approaches is large, but there is a tendency for higher CO 2 values in the early Cenozoic (Ruddiman, 2003;Zachos et al., 2008), with ambiguous results for the last 20 Myr.Moreover, none of the proxies has a continuous record for the last 20 Myr (Fig. 3).For this reason there is a need to compile all available records in a consistent manner.The decomposition of the marine benthic δ 18 O record offers a framework to do so.We use the modelled temperature as a tool to select mutually consistent CO 2 records by assuming that there is a relation between CO 2 and T NH , which is comparable to the relation found in ice cores.In fact, this is justified as several independent proxies do show a similar linear relation (Fig. 4).Different methodologies may explain why the δ 11 B h (δ 11 B from Hönisch et al., 2009) is more consistent with the ice cores CO 2 data then the δ 11 B p (δ 11 B from Pearson and Palmer, 2000).Hönisch et al. (2009) selected samples R. S. W. van de Wal et al.: Reconstruction of a continuous high-resolution CO 2 record around glacial and interglacial extremes, which was not done by Pearson and Palmer (2000).In addition, it has been argued that the Pearson and Palmer (2000) data need to be rejected for reasons related to diagenesis, use of incorrect fractionation factors, and poor modelling of seawater alkalinity and δ 11 B (Foster et al., 2006).The comparison in Fig. 4 reveals that the CO 2 estimates derived from the ice cores, B/Ca, stomata, δ 11 B h and the combination of alkenones and δ 11 B s (δ 11 B from Seki et al., 2010) are mutually consistent, because they reveal similar slopes, whereas the δ 11 B p , and alkenones-derived CO 2 estimates do not show a consistency with the ice-core record.
Hence, we selected only the consistent records to derive an empirical relationship between temperature and CO 2 .This relation is used to calculate CO 2 from the reconstructed T NH in order to generate a continuous CO 2 proxy series that is mutually consistent with the benthic δ 18 O record.The application of the correlation between CO 2 and temperature implies, however, that the regression needs to cover the temperature range as shown in Fig. 1 without having too much bias to the data-rich cold climate state.For this reason, we binned the CO 2 observations in intervals of 1 K NH temperature change, for which results are shown in Fig. 5.The temperature records are running averages over 2000 yr, in order to prevent outliers due to a mismatch in dating of the CO 2 proxy and the benthic record.Furthermore, several tests have been performed to weigh the different accepted CO 2 proxies, by uncertainty in modelled temperature and measured CO 2 .In addition, we tested the effect of the binning size and averaging period, which contribute to the uncertainty in the reconstructed CO 2 .Omitting one of the proxies effects the end result by at most 6 % (B/Ca data) in the slope between ln(CO 2 /CO 2,ref ).Based on all these tests, we eventually estimated an uncertainty of 10 % in the slope between ln(CO 2 /CO 2,ref) and T NH around a central value of 39 K.A log-linear regression between T NH and CO 2 is used because of the saturation of the absorption bands for CO 2 (Myhre et al., 1998).Accordingly, the CO 2 record as presented in Fig. 1 has an uncertainty of 20 ppmv for cold climates and up to 45 ppmv for warm climates.The larger uncertainty for warmer climates is due to the logarithmic relation between CO 2 and temperature.
Over the past 800 kyr the reconstructed CO 2 record is in good agreement with the ice-core record, (Fig. 6c), which is, however, input to the reconstruction and, therefore, not an independent result.Over the mid-Pleistocene transition (defined here from 1.5 to 0.5 Myr), our results indicate a gradual decline of about 23 ppmv, and at the same time an increase in the amplitude.Carbon-cycle simulation results over the last 2 Myr across the Mid-Pleistocene Transition (Köhler and Bintanja, 2008) support the change in amplitude, but suggest stable glacial CO 2 values and reduced interglacial CO 2 .It is also unclear why the combined δ 11 B and Alkenone record is higher than our reconstruction for the last 1.5 Myr.Around 10 Myr ago the B/Ca data measured on planktonic foraminifers indicate much lower CO 2 concentrations, in fact, more in line with the GEOCARB (carbon-cycle model, Berner, 1994) estimates (Fig. 6a).Ultimately this implies an inconsistency between benthic δ 18 O reconstructions and B/Ca.The difference is too large to be attributed to model uncertainties.

Long-term knowledge on climate sensitivity
Since we now obtained a continuous record of T NH and CO 2 , we can address the long-term climate sensitivity in more detail.There are various ways to define climate sensitivity.
Here we define climate sensitivity (S) as the functional dependency of changes in global surface temperature (T g ) on CO 2 , thus, T g = f (CO 2 ).It is calculated from the radiative forcing (R) caused by changes in CO 2 , other greenhouse gases, and various fast and slow feedbacks (f ), which will be specified below.A general formulation for the global temperature is: In this general setting, changes in CO 2 might be the cause for climate change, thus, representing the forcing term R or a feedback.A functional relationship between the global temperature and CO 2 can subsequently be developed, assuming that CO 2 is causing the radiative imbalance, R = f (CO 2 ), while the initial perturbation of this imbalance might be caused by other processes.This by no means implies that we believe that changes in CO 2 were always the driver for climate change over the last 20 Myr, but it is used to derive a functional relationship between T g and CO 2 .The opposite procedure (forcing by other processes and feedbacks by CO 2 ) is certainly a valid option, but, for reasons of simplicity, here we follow only one of the two possible calculations.
From radiative transfer theory, we know that due to the saturation of the absorption bands, a logarithmic relationship has to be applied for the radiative forcing of CO 2 : where R is the radiative forcing in W m −2 and β is estimated to be 5.35 W m −2 (Myhre et al., 1998).This implies a radiative forcing of −2.4 W m −2 for the observed changes in CO 2 from LGM to present day, and +3.7 W m −2 for a doubling of CO 2, with CO 2 , ref = 278 ppmv, the pre-industrial level.Non-CO 2 greenhouse gases like CH 4 and N 2 O enhance this direct radiative forcing of CO 2 .Hence, for the last 800 kyr this enhancement was about 30 % (Köhler et al., 2010), which is represented by γ = 1.3.The sensitivity S of the climate system to external forcing is typically described by the Charney sensitivity S c (Charney et al., 1979), which includes the fast feedbacks of the system (water vapour, lapse rate, albedo, snow and sea ice, clouds).It is the quantity usually calculated by coupled ocean-atmosphere models.Here, we use a sensitivity S p derived from paleo data of 0.72 K/(W/m 2 ) (Köhler et al., 2010).It is based on a LGM cooling of T g,LGM = −5.8K (Schneider von Deimling et al., 2006) and a total radiative forcing R LGM = −9.5 W m −2 (Köhler et al., 2010).
The total forcing of the system (R ) includes the forcing R caused by all greenhouse gases, which is amplified by a feedback factor f consisting of the slow feedbacks not included in S p .It represents the feedbacks from albedo changes caused by land ice, vegetation and dust.
For the last 800 kyr, a value for f = 0.71 and γ = 1.3 is derived from proxy-based evidence (Köhler et al., 2010).For T NH we obtain a value of −15.8 K averaged over the period 23 to 19 kyr BP (LGM), which is 2.7 times larger than the global temperature change of −5.8 K.This is in line with estimates of the polar amplification factor (α) from a GCM model of 2.5-3 by Singarayer and Valdes (2010).The final expression for the change of T NH can now be written as: Calculation of C (α = 2.7, β = 5.35, γ = 1.3, S p = 0.72, f = 0.71, CO 2,ref = 278) results in an indicative value of 47.0 K for cold conditions (i.e.past 800 000 yr) as f is based on LGM proxy data.Where it might be noted that application of CO 2 , ref = 278 ppmv implies that T NH is expressed relative to pre-industrial levels.This calculated value for C is, however, considerably larger than the 39 ± 3.9 K we found for the slope between ln(CO 2 /CO 2,ref) and T NH over the past 20 Myr (Fig. 5).Remarkably, Fig. 5 does not indicate that C depends on the CO 2 concentration itself.One might expect f to be smaller for warmer conditions, but this is apparently compensated by a change in one of the other parameters.From Fig. 5, it is clear that the scatter for warmer conditions is large and we have to await more proxy data for warmer conditions to see whether this log-linear relation holds.
A source of uncertainty is the value for S p , which is derived from LGM conditions.Hargreaves et al. (2007) argued that this value is 15 % smaller than the value for 2 × CO 2 , close to our derived values for the early Miocene climatic optimum.A similar change in the sensitivity implies that C would decrease to a value of 40.0 K, which is well within the uncertainty of our estimated value for C of 39 ± 3.9 K (Fig. 5).
Another source of uncertainty in this analysis is probably the assumption that the polar amplification factor, α, is constant over time.Theories and observations on much warmer climate states suggest a decrease in the meridional temperature gradient (e.g.Huber and Cabellero, 2011) implying a decrease in α.The fact that C does not vary with increased temperature suggests that a possible change in α or the longterm feedbacks compensate one another.The applied method does not allow separating α from the long-term feedbacks or the GHG forcing.In theory, if α were much smaller for warmer climate conditions and f smaller, it would imply that considerably higher CO 2 concentrations in the past are necessary to explain the temperatures derived from the benthic δ 18 O record.However, the stomata-derived CO 2 data, the GeoCarb data (Berner 1994), and the B/Ca data do not indicate much higher CO 2 concentrations, at least not over the period considered here, suggesting that this theoretical example is not likely.
Too little information is available to attribute individual changes in the parameters over 20 Myr.But given the fact that the fitted value of C based on the presented data in this paper, and the estimated value of C based on our knowledge of the system (Köhler et al., 2010) are close to each other, implies that the combined effect of the key processes affecting benthic δ 18 O records, temperature and CO 2 are incorporated sufficiently accurately for at least the period that there is ice on Earth.The implication of Eqs. ( 4) and ( 5) is that systematic errors in the various components can cancel each other very easily, complicating the interpretation of the sensitivity.A 20 % increase in α can be compensated by a 20 % decrease in β, γ , S c or equally large increase of (1 − f ).As we determine C directly from the model inversion, one cannot expect improved insights in the different values for α, β, γ , S c or f from our methodology.The fact that the values used in the literature eventually lead to a similar value for C merely indicates that there is no reason to adjust literature values for climate sensitivity based on the inversion used here.
To set our approach into context with existing calculations of climate sensitivity one might simply calculate the pure Charney climate sensitivity out of Eqs. ( 1), ( 4) and ( 5) by choosing f = 0 (no slow feedbacks of ice sheets and vegetation) and γ = 1 (no non-CO 2 GHG), which leads to C = 7.2 K and T g = 2.7 K for a CO 2 doubling.This is close to the original value and well within the uncertainties of 3 ± 1.5 K (Charney et al., 1979).
For the interpretation of the climate sensitivity values one needs to be careful.In order to calculate temperature changes from CO 2 changes as presented in the literature often different feedbacks are included.Here, we use the feedback factor obtained by Köhler et al. (2010) based on paleo proxy evidence for the LGM to show that the value of C as obtained from our analysis are in agreement with estimates of α, β, γ and S p and the 15 % reduction in S p depending on the climate state as proposed by Hargreaves et al. (2007).Adopting a lower value for f as maybe deduced from Hansen et al. (2008) would indicate that our value for C would be too large.Of course this can be compensated by assuming a larger value for α, implying an even larger meridional temperature gradient, but that seems unlikely.In other words, the temperatures as we find them, based on the benthic records, support the feedback factor as derived by Köhler et al. (2010).If we extract this value (f = 0.71) from C to consider the short-term climate sensitivity and facilitate a more realistic comparison with Hansen et al. (2008), we yield a value for the sensitivity, which is smaller than the sensitivity by Hansen et al. (2008).Furthermore, Hansen et al. (2008) uses in their calculations a global LGM cooling of T g = −5 K, while we refer to the temperature anomaly of −5.8 K, which is our understanding best supported by proxy data.This readily explains a similar larger difference of 16 % in the projected temperature changes of the two approaches calculated for future climate with doubled CO 2 .

Discussion and conclusion
Accepting the CO 2 concentration as presented in Fig. 1 with all its caveats and limitations, completes the picture of the key climate variables over the last 20 Myr.The figure shows a gradual decline from about 450 ppmv around 15 Myr ago to a mean level during the last 1 Myr of 225 ppmv or a decrease of 225 ppmv.This is about 2.2 times the increase in CO 2 concentration over the last century as well as 2.2 times the range in the ice-core record over the past 800 kyrs.If we would have used only the ice-core record, we would have obtained Middle Miocene values 300 ppmv above present-day level and the sensitivity would not agree with the analyses presented in the previous paragraph as the sensitivity (C) would decrease to a value as low as 28.5 K. Hence, the application of the inverse model and the stacked binning procedure is crucial for the results.
It should be realised that the results are as good as the input is for the model.We use the widely used benthic record by Zachos et al. (2008), which has δ 18 O values during the Middle Miocene, which are comparable to their values at the E-O transition.As a consequence, our Middle Miocene ice volume is small.A test with the stacked record by Cramer et al. (2009), which addresses in more detail interbasin changes of δ 18 O records, indicate a larger ice volume for the Middle Miocene.As a consequence temperature and CO 2 changes are differently, yielding Middle Miocene CO 2 values of 410 ppm.
The question remains of course what causes these subtle changes in the carbon cycle on the long timescale.In order to answer this question, much higher resolution and accuracy of CO 2 records are necessary.The large sensitivity implies that, in contrast to earlier conclusions (Hönisch et al., 2009), subtle changes in CO 2 (possibly internal), may have caused the MPT, when dominant 41-kyr glacial cycles evolved into a dominant 100-kyr rhythm (Van de Wal and Bintanja, 2009).Our results indicate an average change of only 23 ppmv between 1.5 Ma and 0.5 Ma, and also an increasing amplitude.This result seems to be more in line with a recent estimate by Lisiecki (2010) based on marine δ 13 C measurements and the δ 11 B data by Hönisch et al. (2009) than with the B/Ca derived CO 2 data by Tripati (2009), which indicates a larger change in CO 2 over the MPT.However, the trend in CO 2 over time is too small given the accuracy of the applied methods to draw firm conclusions on this point.
With respect to the inception of the Northern Hemisphere ice around 2.7 Myr ago, our results indicate that the trend in CO 2 before the inception is higher than the average rate of change (see Fig. 1d), and that the inception takes place once the long-term average concentration drops below 265 ± 20 ppmv (Fig. 6b).So for this climate transition a change in CO 2 seems to be more important than for the mid-Pleistocene transition.
In conclusion, the self-consistency of our approach should enable researchers from various disciplines to identify more easily whether new CO 2 proxies are in line with the reconstructed temperature change derived from the δ 18 O record and the ice core derived CO 2 , δ 11 B h , δ 11 B s + alkenones, B/Ca and stomata.It is tempting to use the newly reconstructed CO 2 to address the long-term climate sensitivity, but the approach is not accurate enough to revise current ideas about climate sensitivity because no clear distinction can be made between the various feedbacks and possibilities for compensating mechanisms.Various geological processes important during the last 20 Myr such as mountain uplift (e.g.Foster et al., 2010) and changes in the gateways are not considered here.However, for global climate changes, CO 2 induced changes dominate, as shown by Henrot et al. (2010) who based their argument on a model of intermediate complexity that geological processes like mountain building and changes in ocean gateways, are of secondary importance for global temperature and can not explain the proxy reconstructions of the change in temperature within their modelling framework.Based on the analyses of climate sensitivity, we argue that the applied method, in combination with the available data, is not accurate enough to revise ideas on climate sensitivity.Compensating effects prevent the unravelling of the different contributions to the log-linear relation between CO 2 and temperature.
As a final remark, we stress that the relation between CO 2 and the Northern Hemisphere temperature may change, for instance, by a changing strength of feedbacks or a different relation between global and Northern Hemisphere temperature changes, and as such it complicates the interpretation of paleo data as analogue for present day conditions.Paleo data provide the range of natural fluctuations, but the rate of change of key variables is shown to be depending on the state of the system (Köhler et al., 2010), the timescale of interest and the processes at stake, which are not necessarily similar in the past as for present day climate change.Paleo data should be interpreted in the context of the conditions and forcing prevailing at that time.

Fig. 1 .
Fig. 1.Records of key climate variables over the last 20 Myr.Forcing of the model are changes in the stacked benthic δ 18 O record with respect to preindustrial times (a, dark blue, Zachos et al., 2008).Output is a consistent record for the Northern Hemisphere temperature change with respect to pre-industrial conditions (b, green) and sea level (c, light blue).The reconstructed CO 2 record (d, orange) is obtained by inverting the relation between NH temperatures and CO 2 data.The δ 18 O curve is smoothed in order to clarify the gradual decrease over time.All data are available every 0.1 kyr.The thick lines represent 400-kyr running mean.Grey error bars indicate the standard deviation of model input and output.For CO 2 the error bar is calculated as 400-kyr running mean, for the other records it is the standard deviation on the 0.1 kyr value as used in the model.Data are available as a Supplement to this paper.

Fig. 2 .
Fig. 2. (a) Sea level change is shown as a function of the reconstructed temperature for a set of 3-D NH ice sheets (blue) and for a set of five 1-D ice-sheet models (red) (Bintanja and Van de Wal, 2008; De Boer et al., 2010).The more sophisticated 3-D results are validated by observation of sea level (Lambeck and Chapell, 2001; Rohling et al., 2009).The 1-D results are in line for the colder climate condition with the 3-D results.The warm temperatures in combination with the sea level change resemble the melt of SH ice.The thick lines show the mean trends, emphasizing the low gradient for the present-day climate centred around zero.(b) The response of the individual ice sheets.Note the strong transient and nonlinear response for each ice sheet.
Fig. 2. (a) Sea level change is shown as a function of the reconstructed temperature for a set of 3-D NH ice sheets (blue) and for a set of five 1-D ice-sheet models (red) (Bintanja and Van de Wal, 2008; De Boer et al., 2010).The more sophisticated 3-D results are validated by observation of sea level (Lambeck and Chapell, 2001; Rohling et al., 2009).The 1-D results are in line for the colder climate condition with the 3-D results.The warm temperatures in combination with the sea level change resemble the melt of SH ice.The thick lines show the mean trends, emphasizing the low gradient for the present-day climate centred around zero.(b) The response of the individual ice sheets.Note the strong transient and nonlinear response for each ice sheet.

Fig. 4 .
Fig. 4.Scatter plot of the different CO 2 proxies as a function of the reconstructed temperature, which is derived from the benthic δ 18 O record as shown in Fig.1.Only records with filled symbols δ 11 B h(Hönisch et al., 2009), B/Ca(Tripati et al., 2009), alkenones+ δ 11 B s(Seki et al., 2010), stomata data(Kürschner et al., 1996(Kürschner et al., , 2008;;Retallack 2009), and the ice-core record(Petit et al., 1999;Siegenthaler et al., 2005;Lüthi et al., 2008) are used for further analysis.For reasons of transparency CO 2 is plotted in ppmv.If CO 2 would be plotted as ln(CO 2 /CO 2 , ref ) a similar picture emerges.The latter is physically more consistent as it takes the saturation of the absorption bands into account.

Fig. 5 .
Fig. 5.The selected (n = 1302) proxy CO 2 data (red dots) binned in intervals of 1 K NH temperature change.The error bars represent one standard deviation variability of the data in the selected temperature interval.The additional lines show the range in C values from different weighing tests, blue C +10 %, red C −10 %.The slope of the regression line corresponding to C is 39 ± 3.9 K.

Fig. 6 .
Fig. 6.Comparison of reconstructed CO 2 record with C = 39 K, with proxy records (symbols as in Fig. 3).Panel (a) for the full 20 Myr period, (b) for the period around the Northern Hemisphere glacial inception and (c) for the mid-Pleistocene transition.Note that the vertical scale is different for the different panels.The Horizontal bar in panel (b) indicates the onset of major glaciation in the Northern Hemisphere.The crosses with open squares are the GEOCARB data.
. W. van de Wal et al.: Reconstruction of a continuous high-resolution CO 2 record