Palaeogeographic controls on climate and proxy interpretation

. During the period from approximately 150 to 35 million years ago, the Cretaceous–Paleocene–Eocene (CPE), the Earth was in a “greenhouse” state with little or no ice at either pole. It was also a period of considerable global change, from the warmest periods of the mid-Cretaceous, to the threshold of icehouse conditions at the end of the Eocene. However, the relative contribution of palaeogeographic change, solar change, and carbon cycle change to these climatic variations is unknown. Here, making use of recent advances in computing power, and a set of unique palaeogeographic maps, we carry out an ensemble of 19 General Circulation Model simulations covering this period, one simulation per stratigraphic stage. By main-taining atmospheric CO 2 concentration constant across the

plications for the interpretation of single-site palaeo proxy records. In particular, our results allow the non-CO 2 (i.e. palaeogeographic and solar constant) components of proxy records to be removed, leaving a more global component associated with carbon cycle change. This "adjustment factor" is used to adjust sea surface temperatures, as the deep ocean is not fully equilibrated in the model. The adjustment factor is illustrated for seven key sites in the CPE, and applied to proxy data from Falkland Plateau, and we provide data so that similar adjustments can be made to any site and for any time period within the CPE. Ultimately, this will enable isolation of the CO 2 -forced climate signal to be extracted from multiple proxy records from around the globe, allowing an evaluation of the regional signals and extent of polar amplification in response to CO 2 changes during the CPE. Finally, regions where the adjustment factor is constant throughout the CPE could indicate places where future proxies could be targeted in order to reconstruct the purest CO 2 -induced temperature change, where the complicating contributions of other processes are minimised. Therefore, combined with other considerations, this work could provide useful information for supporting targets for drilling localities and outcrop studies.

Introduction
Over the last 150 million years, the climate of the Earth has experienced change across a broad range of timescales, from geological (10's of millions of years), to orbital 10's-100's of thousands of years), to millenial, to decadal.
Variability on timescales from geological to orbital has been characterised by measurements of the isotopic composition of oxygen in the calcium carbonate of benthic foraminifera (δ 18 O benthic , e.g. Friedrich et al. (2012); Cramer et al. (2009) , Fig. 1); this indicates, in the broadest sense, a long-term cooling (and increasing ice volume) trend from the mid-Cretaceous (∼ 100 million years ago, 100 Ma) to the modern. Imprinted on this general cooling are several shorter geological-timescale variations such as cooling through the Paleocene (65 to 55 Ma), and sustained warmings in the early Eocene (55 to 50 Ma) and middle Miocene (17 to 15 Ma). The δ 18 O benthic record also shows evidence for orbital scale variability in icehouse periods (for example Quaternary glacial-interglacial cycles, from around ∼ 2 Ma) and greenhouse periods (for example Paleogene hyperthermals, ∼ 55 Ma), and other events occurring on sub-geological timescales, such as the Eocene-Oligocene boundary (34 Ma) and Cretaceous Oceanic Anoxic Events (OAEs, ∼ 100 Ma).
A methodology for reconstructing the global-mean surface temperature of the past ∼ 65 million years has been developed by Hansen et al. (2013), making assumptions about the relationship between δ 18 O benthic and deep ocean temperature, and between deep ocean temperature and surface temperature (Fig. 1). The assumptions mean that the absolute values of temperature need to be treated with considerable caution, but the record indicates global mean surface temperatures of ∼ 25 • C in the Paleocene, peaking at over 28 • C in the early Eocene. From the early Eocene to the present day, there is a general cooling trend, with temperatures decreasing to ∼ 24 • C degrees by the late Eocene. Today, global mean surface temperature is close to 15 • C.
Although it has long been thought that greenhouse gas concentrations are the primary cause of Cretaceous and Eocene warmth (Barron and Washington, 1984;Barron et al., 1995), the reasons for variability within this period are currently largely unknown. Possible candidates for the forcing on the Earth system on these timescales include changes in solar forcing, direct tectonic forcing (i.e. changes in orography and bathymetry and continental position), and greenhouse gas forcing (most likely CO 2 ) through changes in the carbon cycle, or some combination of these.
The solar forcing is related to an increase in solar constant, resulting from an increasing luminosity of the sun over 10's of millions of years. This itself is due to continued nuclear fusion in the core of the Sun, converting hydrogen into helium, reducing the core's density (Fig. 2, Gough, 1981). As a result, the core contracts and temperature increases. The increase in luminosity is a consequence of this core temperature increase. According to the solar model of Gough (1981), the Ts estimate (Cramer, Hansen) Figure 1. Climate evolution over the last 150 million years, as expressed by the benthic oxygen isotope records of Friedrich et al. (2012) and Cramer et al. (2009) (coloured dots), and a surface temperature record (H13C09) produced by applying the methodologies of Hansen et al. (2013) to the Cramer et al. (2009) δ 18 O benthic data, and applying a 10-point running average (grey line).  Gough (1981) and used for the CPE simulations in this paper (green line), and from Caldeira and Kasting (1992) (orange line). Horizontal line shows the modern value. solar constant increases nearly linearly from ∼ 1348 W m −2 at 150 Ma to 1365 W m −2 for the modern.
Also over the timescales of tens of millions of years, plate tectonics has led to major changes in the position and configuration of the continents, and bathymetric and topographic depths and heights. These changes can have a direct effect on climate, for example through changing the global distribution of low-albedo (ocean) vs. high-albedo (land) surfaces (e.g. Barron et al., 1980), and/or changing ocean gateways leading to modified ocean circulation (e.g. Zhang et al., 2011), and/or topography modifying the area of land above the snowline (e.g. Foster et al., 2010) or atmospheric circulation (e.g. Ruddiman and Kutzbach, 1989). Palaeogeographical reconstructions for certain periods in Earth's history exist (e.g. Scotese, 2001), but not always at high temporal resolution, nor in Overlain is the CO 2 concentration in each simulation described in this paper (red dots). The open blue circles show CO 2 concentration in a set of model simulations of the Oligocene and Neogene, which are not discussed in this paper but are provided for reference. Horizontal black lines show 280, 560, and 1120 ppmv. a form which can readily be implemented into a climate model. In Sect. 2.1 we present a new set of palaeogeographical reconstructions which span the Cretaceous-Paleocene-Eocene, from approximately 150 to 35 Ma (Fig. S1 in the Supplement), and provide one of the reconstructions (Ypresian stage) in digital format in the Supplement.
The greenhouse gas forcing is indirectly related to tectonic changes, through changes in the balance of sources (e.g. volcanism) and/or sinks (e.g. weathering of silicate rocks, Raymo et al., 1988) of CO 2 , or other greenhouse gases, and/or changes to the sizes of the relevant reservoirs (e.g. due to changes in the residence time of carbon in the ocean due to changes in ocean circulation, Lauderdale et al., 2013, them-selves driven ultimately by tectonic changes). In addition, greenhouse gases will also likely be modified in response to climate changes caused by the solar or direct tectonic forcing. Efforts have been made to reconstruct the history of CO 2 over geological timescales (e.g. see compilation in Fig. 3, Honisch et al., 2012). However, CO 2 proxies are still associated with relatively large uncertainties, despite currently undergoing a period of rapid development (e.g. Zhang et al., 2013;Franks et al., 2014;Martinez-Boti et al., 2015).
Disentangling these various forcings on long-term climate evolution is a key challenge. Previous work has often been in a modelling framework, and has focused on either the role of palaeogeography across time (e.g. Donnadieu et al., 2006), or the role of CO 2 for a particular period (e.g. Caballero and Huber, 2013).
Given the uncertainties in CO 2 concentration and the carbon cycle, and to avoid complications of the feedbacks associated with continental ice, in this paper we focus on the direct role of palaeogeography and solar forcing on controlling the greenhouse climates from the earliest Cretaceous to the end of the Eocene (Cretaceous-Paleocene-Eocene, CPE). This also allows us to provide an "adjustment factor" for palaeo proxy records, which accounts for the non-CO 2 component of climate change (see Sect. 3.4).
The very first attempts to model time periods within the CPE were carried out in the laboratory, with rotating water tanks covered with moulded foam representing palaeogeography, and jets of compressed air simulating wind stress (Luyendyk et al., 1972). Several early numerical modelling studies also focused on the ocean circulation, and in particular the flow regime through the Tethys seaway (e.g. Barron and Peterson, 1990). The relative importance of palaeogeography vs. surface albedo vs. greenhouse gases in warm palaeoclimates was examined using what would now be considered low resolution GCMs (Barron and Washington, 1984;Barron et al., 1995) or using energy-balance models (Barron et al., 1980). These indicated that CO 2 was likely the primary driver of Cretaceous and Eocene warmth. The majority of modelling work since then has focused on the periods of maximum warmth, i.e. during the mid-Cretaceous (e.g. Sellwood et al., 1994;Poulsen et al., 2001Poulsen et al., , 2003Zhou et al., 2012) or Early Eocene (e.g. Huber and Caballero, 2011); see summary in Lunt et al. (2012), or transient hyperthermals such as the PETM (e.g. Winguth et al., 2010). Other work has examined other relatively warm periods such as the Late Cretaceous (Otto-Bliesner et al., 2002;Upchurch et al., 2015;Donnadieu et al., 2016). Several studies have addressed issues of model-data comparisons, including the interpretation of oxygen isotope proxies in both continental and oceanic proxies (e.g. Poulsen et al., 2007;Zhou et al., 2008), or uncertainties in Mg/Ca calibrations (e.g. Bice et al., 2003Bice et al., , 2006. Recently, it has been argued that if these uncertainties, and other issues such as seasonality of the proxies, are taken into account, then some models can simulate the climate of the early Eocene consistently with the data . Although no previous study has explored the role of varying palaeogeography throughout the CPE as we do here, several previous modelling studies are worth noting, which carried out sensitivity studies to palaeogeography with a more limited scope. Poulsen et al. (2001Poulsen et al. ( , 2003 carried out model simulations under two different Cretaceous palaeogeographies, representing conditions before and after the separation of the African and South American continents to form the Atlantic. They found that continental positions strongly influenced ocean circulation, in particular regions of deep water formation. Bice and Marotzke (2002) examined the role of ocean gateways in the Eocene, and found that the configuration of polar seaways affected the sensitivity of climate to hydrological forcing, through changes in ocean overturning. Spicer et al. (2008) used three Cretaceous palaeogeographies, and compared a number of model simulations with data from the Cretaceous Siberian continental interior, but the sensitivity studies were not consistent across the time periods. Donnadieu et al. (2006) also examined three palaeo-geographies through the Cretaceous, using the FOAM model coupled to a slab ocean. They focused on the influence of continentality on seasonality, but noted that changing palaeogeography alone could give a ∼ 4 • C global-mean warming at a constant CO 2 level.
Our work presented here builds on these and other previous studies, but represents an advance because (a) new palaeogeographic maps of this time period have become available, which improve on previous representations in terms of both accuracy and temporal resolution (see Sect. 2.1) and (b) increases in available computing power means that for the first time we can (i) spin up a large number (19) of simulations through this time period towards equilibrium, allowing unprecedented temporal resolution, and (ii) use more advanced models than have typically been used previously.
The two key questions which we address are -How and why has palaeogeography and solar output affected global mean, zonal mean, and local temperatures through the CPE?
-What are the implications of our results for interpreting proxy reconstructions on geological timescales?

Experimental design
A set of 19 simulations was carried out, one for each stratigraphic stage (henceforth "Stage") between the earliest Cretaceous (Berriasian, 146-140 Ma) and the latest Eocene (Priabonian, 37-34 Ma). This section describes the boundary conditions implemented, the model used, and the simulation design (see Table 1).

Palaeogeographies
The reconstruction of tectonics, structures, and depositional environments which underpin this study was created by Getech Plc using methods based on those of Markwick and Valdes (2004). The palaeo digital elevation models (henceforth "palaeogeographies") used as boundary conditions in the model for each Stage are informed by these reconstructions, which are in turn constrained by extensive geological databases. These data include published lithologic, tectonic and fossil studies, the lithologic databases of the Paleogeographic Atlas Project (University of Chicago), and deep sea (DSDP/ODP) data. They are extensively updated from the series described in Markwick (2007); critically, they include bathymetric information which is essential for running coupled atmosphere-ocean climate models, and which is absent from the Markwick (2007) maps. These data are also used to develop the plate model on which the palaeogeographies are built. The palaeogeographies were produced at an original resolution of 0.5 • × 0.5 • , and from these we generated model-resolution (3.75 • × 2.5 • ) land-sea mask, topography and bathymetry, and the sub-gridscale orographic variables required by the model. In order to maintain stability, the palaeogeographies were smoothed globally, with additional smoothing applied in the Arctic. Furthermore, some palaeogeographies required additional flattening of the Arctic and/or Antarctic bathymetry, and/or smoothing of the topography, and/or minor changes to the land-sea mask, at various phases of the spinup of the model simulations, in order to maintain model stability. Details are given in Table 1. The runoff routing is carried out by assuming that rivers run downhill at the resolution of the model orography (see Fig. S2). The palaeogeographies are proprietary and cannot all be distributed digitally, but figures showing the orography and bathymetry at the model resolution (of Phase 4, see Sect. 2.5), for each Stage discussed in this paper, are included in the Supplement, Fig. S1. In additon, for one stage (Ypresian, early Eocene), we provide a digital version of the palaeogeography at the model resolution in the Supplement. Also provided with the palaeogeographies are distributions of lakes (shown in the Supplement, Fig. S3). Finally, for the latest two Stages in the CPE (Bartonian and Priabonian), there are also very small ice caps prescribed on Antarctica, which for the purposes of this study we consider part of the palaeogeography, and which we assume have only a very small effect on the results (also shown in the Supplement, Fig. S3).

CO 2 forcing
Geological proxy data for atmospheric CO 2 on the timescale of the CPE have large uncertainties, but, in general, indicate a mean of between 2× and 4× pre-industrial (PI) CO 2 concentrations, i.e. 560-1120 ppmv (Honisch et al., 2012;Beerling and Royer, 2011;Royer et al., 2012). This paper focuses on the effects of changing palaeogeography and solar output. As such, we keep CO 2 constant at a prescribed value of 1120 ppm (4 × PI) for all simulations. Therefore, any changes in climate between different Stages are due to the palaeogeographical and solar constant changes alone. The value of 1120 ppmv is chosen to represent a reasonable estimate of atmospheric CO 2 for the duration of the CPE (see Fig. 3), and can be considered as also incorporating the contribution to radiative forcing form other greenhouse gases, such as methane, which may also have been elevated compared to modern . Future work will analyse the climates of the CPE at lower or higher atmospheric CO 2 , and compare the climate sensitivities of these different Stages.

Solar forcing
The insolation at the top of the atmosphere (Total Solar Irradiance, TSI) for each Stage was calculated following Gough (1981), and is shown in Fig. 2. The evolution of TSI is very Table 1. Summary of all model simulations. The age column shows the age of the middle of the respective Stage. The solar constants are calculated using these ages according to the formula described in Gough (1981). The smoothing indicates the changes that had to applied to each Stage to ensure stability. F = fourier filtering at high latitudes, Ar = flat Arctic ocean, An = flat polar Southern ocean, O1 = polar orographic smoothing, O2 = polar orographic smoothing and capping of polar topography, L = minor changes to land-sea mask. P3 and P4 indicate Phases 3 and 4 of the simulations.

Stage
Age ( similar to that from Caldeira and Kasting (1992), also shown in Fig. 2, illustrating that this forcing likely has relatively small uncertainty.

Model description
The simulations described in this paper are all carried out using the UK Met Office coupled ocean-atmosphere general circulation model HadCM3L version 4.5. HadCM3L has been used in several palaeoclimate studies for different geological periods including the early Eocene (e.g. Lunt et al., 2010;Loptson et al., 2014) and the late Miocene (e.g. Bradshaw et al., 2012). The resolution of the model is 3.75 • in longitude by 2.5 • in latitude, with 19 vertical levels in the atmosphere and 20 vertical levels the ocean. The HadCM3L model is very similar to HadCM3, a description of which can be found in Gordon et al. (2000), but HadCM3L has a lower horizontal resolution in the ocean (3.75 • × 2.5 • compared with 1.25 • × 1.25 • ). In addition, HadCM3L is coupled to the dynamic global vegetation model TRIFFID (Top-down Representation of Interactive Foliage and Flora Including Dynamics) (Cox et al., 2001) via the land surface scheme MOSES 2.1 (Cox et al., 1999). TRIFFID calculates the fraction of each gridcell occupied by each of five plant functional types: broadleaf trees, needleleaf trees, C 3 grasses, C 4 grasses and shrubs. Although TRIFFID simulates modern plant functional types, it has been argued that such a model can provide the first order signal from vegetation feedbacks through the last 250 million years (Donnadieu et al., 2009;Zhou et al., 2012).

Simulation description
All simulations have undergone the same spin-up procedure, totaling 1422 years, with the same initial conditions and boundary conditions, with the exception of solar constant and palaeogeography (as discussed in Sects. 2.3 and 2.1). The ocean is initialised as stationary, with a zonal mean temperature structure given by an idealised cosine function of latitude, ( 21−z 20 )( 22 cos(φ) + 10), where φ is latitude and z is the model vertical level from 1 at the ocean surface to 20 at a depth of ∼ 5200 m, and a constant salinity of 35 ppt. The atmosphere is initialised from an arbitrary atmospheric state from a previous preindustrial simulation. Land-surface initial conditions (e.g. soil moisture, soil temperature) are globally homogeneous.
The spinup procedure consists of four "Phases". The first 50 years (Phase 1) of each simulation are run with an atmospheric CO 2 concentration of 1 × PI, and with global vegetation fixed as bare soil. For the next 319 years (Phase 2), the atmospheric CO 2 is increased to 4 × PI, the TRIFFID vegetation component of the model is turned on, and the vertical structure of atmospheric ozone is diagnosed from the modelled troposphere height as opposed to from a prescribed field (in order to avoid a runaway greenhouse encountered in previous high-CO 2 simulations with fixed ozone distributions). Phase 3 consists of 53 years of simulation in which prescribed lakes and glaciers are added to the model (see Sect. 2.1). In the first three Phases, the ocean dynamics are simplified to enhance stability, by imposing a purely baroclinic ocean circulation in which the vertically integrated flow is zero. The final phase, Phase 4, consists of a final 1000 years in which both the baroclinic and barotropic ocean dynamics are turned on, giving a total of 1422 years of simulation. The barotropic streamfunction calculation requires islands to be defined manually -these are shown in Fig. S4. In addition, to maintain model stability, for some Stages additional smoothing and/or flattening of the topography was required at the beginning of Phase 4.
The details of the simulations are summarised in Table 1.

Time series
In order to assess the extent to which the models are spun up, we first examine the time series of evolution of the global ocean. Figure 4 shows the temporal evolution of ocean temperature, at three vertical levels (5, 670, and 2700 m), for each Stage, over the 1422 years of model simulations. As can be seen the simulations are not fully in equilibrium at any of the depths, but are approaching equilibrium in the upper and mid-ocean. In the deepest ocean, there are still significant trends and little sign of equilibrium. As such, for the rest of this study we focus on the surface and upper ocean climatologies. Analysis of the top of the atmosphere (TOA) energy fluxes indicates that for all Stages the system is within 1 W m −2 of radiation balance, but is losing energy at the end of the simulations, at a rate that varies from −0.8 (Berriasian) to −0.34 W m −2 (Campanian).
There are some discontinuities apparent in the timeseries at 670 m depth, at the beginning of Phase 4 after 422 years. This is because at this time several of the simulations had flat Arctic bathymetry imposed, to ensure stability as the barotropic component of circulation was turned on (see Sect. 2.5 and Table 1). This results in cold Arctic waters being removed at this depth, resulting in an apparent sudden warming in the global mean.  The open circles show model results for the Oligocene and Neogene, with lower CO 2 , which are not discussed in this paper but are provided for reference. The modelled trend (ignoring the outlying Berriasian stage) is shown as a solid line, and the expected trend assuming solar forcing only is shown as a pair of dashed lines which start/end at the beginning/end of the modelled trend. Figure 5 shows the global mean annual surface air temperatures (MATs) for the simulations described in this paper, superimposed onto a surface temperature record produced by applying the methodology of Hansen et al. (2013) to the δ 18 O benthic record of Cramer et al. (2009); henceforth "H13C09 record". Model results from the Oligocene (∼ 34 Ma) through to the Pleistocene (last 2 million years), in which CO 2 varies and in which large ice sheets are prescribed, are also shown on this figure. CO 2 is modified to the value shown in Figure 3 at the beginning of Phase 3, otherwise these additional simulations are identical to the CPE simulations. These are shown simply to illustrate that the model can reproduce the first-order response seen in the data during periods when atmospheric CO 2 is better constrained -in particular, the temperature in the Pleistocene model simulation is in good agreement with the Pleistocene glacials in the H13C09 record. However, these will not be discussed in this paper; instead, they will be presented in detail in future publications.

Global annual mean temperatures
It is clear that for the Paleocene and Eocene, there is much less variability and trend in the modelled simulations than suggested by the proxy surface temperature record. Even accounting for the considerable uncertainties in the proxy record, this implies that the majority of the variability in the proxy data is caused by forcings which are not included in the model simulations. The most likely missing forcing is greenhouse gas forcing, probably primarily CO 2 , caused by changes in the carbon cycle on multi-million year timescales.
The H13C09 record implies warmer temperatures than given by our results. Our latest Eocene simulation has a similar temperature to that of the earliest Oligocene in the H13C09 record. If both the model and proxy record are correct, then this implies that the earliest Oligocene had an atmospheric CO 2 concentration of about 1120 ppmv. However this is higher than the values recently reconstructed from CO 2 proxies (Pearson et al., 2009;Pagani et al., 2011), which imply Oligocene CO 2 concentrations closer to 600 ppmv. Indeed, in principle it would be possible to obtain a perfect match between the modelled and observed global means, by choosing an appropriate CO 2 level in the model, thereby generating a model-derived CO 2 record, for comparison with other proxy CO 2 records such as Beerling and Royer (2011). However, given the considerable uncertainties in the H13C09 record, we consider that this would be of little value. Instead, we await the development of more long-duration single-site SST proxy records, across a wide geographical range, and with full consideration of uncertainties, with which to compare our simulations.
The variability in global mean annual temperature between our CPE simulations is much less than the available climate data records imply, and the trends in our simulations do not match those in the data. However, there is some variability present in our CPE simulations. In particular there is a long-term warming trend through the CPE. The trend is 0.0043 • C per million years (correlation coefficient of 0.42). There is a maximum scatter around this trend of about ±0.5 • C, the warmest anomaly being the Berriasian at +0.74 • C, and the coldest anomaly being the Campanian at −0.52 • C. Excluding the outlying Berriasian stage gives a warming trend of 0.0068 • C per million years (correlation coefficient of 0.64, see Fig. 5, black solid line). If the solar forcing was the only forcing acting on the system, the expected temperature trend, δT solar /δt would be where δT /δF is the climate sensitivity of the model (KW −1 m 2 ), S 0 is the solar constant (W m −2 ), and α is the planetary albedo. Under early Eocene conditions with the HadCM3L model, Loptson et al. (2014) found a climate sensitivity to CO 2 doubling of 4.8 • C, which, assuming a forcing due to CO 2 doubling of 3.7 W m −2 , equates to a climate sensitivity, δT /δF , of 1.3 KW −1 m 2 ). Taking a typical value of α from our simulations of 0.27, and assuming this does not change greatly through the simulations, gives an expected warming trend of δT solar /δt = 0.0272 • C per million years (see Fig. 5, black dashed line), about 5 times greater than that of our simulations. We interpret this as implying that there is a long-term effect of palaeogeography which is opposing the trend expected by solar forcing alone; that is, the changing palaeogeography is resulting in a cooling trend of 0.0272-D. J. Lunt et al.: Palaeogeographic controls on the climate system 0.0068 = 0.0204 ∼ 0.02 • C million years −1 . Further sensitivity studies with constant S 0 across the simulations could aid investigation of this long-term trend, and the component due to palaeogeography vs. solar constant.
In the following sections we examine the differences between the simulations in more detail, and investigate the causes of the similarities and differences between the Stages.

Causes and mechanisms of temperature change
Although the variations in global mean temperature due to changes in palaeogeography over time are relatively small, there are substantial variations in regional temperatures. In order to focus on the influence of palaeogeography, and to minimise the effect of the change in solar constant through the simulations, this is best expressed in terms of the temperature anomaly in each Stage, relative to the previous Stage ( Fig. 6; see also Fig. S5 for the absolute temperature of each Stage, and Fig. S6 for the temperature anomaly of each Stage relative to the mean of all Stages). Looking Stage-to-Stage also minimises the effects of the changing land-sea contrast due to continental plate movements.
The response of the system to the palaeogeographic forcing is highly complex, mediated by positive and negative feedbacks. Here we use a number of approaches to investigate the causes of the regional and global differences. Figure 6 shows that the largest local changes are in general over regions which have experienced large local orographic change (e.g. associated with changes in the Western Cordillera range in North America), especially where this is also associated with lateral shifts of the mountains, which is expressed as regions of localised warming adjacent to localised cooling, for example in the Aptian-Albian (Fig. 6n). However, these local changes are not significant in terms of the global mean. Figure 7 shows how much of the global mean temperature change from Stage-to-Stage is due to changes over land, ocean, or regions which switch between land and ocean. It is clear that for the largest Stage-to-Stage transitions (for example Berriasian-Valanginian, ∼ 143-138 Ma; Campanian-Maastrichtian, ∼ 77-68 Ma), the ocean is the dominant contributer to the global mean temperature. On average over all the Stage-Stage transitions, the ocean contributes 0.22 • C, the land 0.10 • C, and transitions from land to ocean contribute 0.05 • C. In addition, the temperature change over land correlates well with the temperature change over ocean (correlation coefficient = 0.74, see Fig. 8a). However, it is unclear whether the ultimate cause of the changes relates to ocean processes (in response e.g. to bathymetry or gateway changes), or whether the ocean is amplifying changes that originate over land.
To investigate this further, we explore the relationship between possible drivers of climate change, and the response. Figure 8b shows the modelled global surface temperature as a function of the change in continental land area. There is a weak negative correlation (correlation coefficient of −0.47), implying some influence on global temperature from the relative albedo of land compared with ocean. Figure 8c shows the modelled land surface temperature as a function of the change in mean orography. There is a weak negative correlation (correlation coefficient of −0.49), implying some influence on land temperature from the local lapse-rate effect. The relationship between global temperature and orography is weaker (not shown, correlation coefficient of −0.37), implying that the lapse-rate effect primarily affects local continental temperatures.
We also carry out an energy balance analysis of the causes of temperature change between each Stage, following Heinemann (2009) and Lunt et al. (2012). This allows the global and zonal mean surface temperature change between Stages to be partitioned into contributions from changes in planetary albedo, emissivity, solar constant, and heat transport. We do not save clear-sky flux output from the model, so further partitioning into cloud albedo vs. surface albedo and cloud emissivity vs. greenhouse gas emissivity is not possible in this case. The global mean temperature change correlates very well with the contribution due to emissivity (Fig. 8d) and planetary albedo (Fig. 8e), implying that on a global scale, both emissivity feedbacks (due to water vapour and clouds interacting with long-wave radiation) and planetary albedo feedbacks (due to clouds and the surface interacting with short-wave radiation) play an important role in amplifying the underlying forcing related to the palaeogeographic changes.
It is instructive to focus on the largest transitions in the modelled record. From Fig. 7, we identify 3 such transitions, which all have a global mean temperature change of over 0.7 • C (for comparison with the fourth largest transition which is 0.57 • C): the Berriasian-Valanginian (∼ 142-138 Ma), the Aptian-Albian (∼ 119-106 Ma), and the Campanian-Maastrichtian (∼ 77-68 Ma).

Berriasian-Valanginian (Fig. 6r)
This transition is a cooling of 0.84 • C in the global mean. In the highest Arctic, there is a warming of more than 10 • C, due to a transition from a polar Arctic continent in the Berriasian, to open ocean in the Valanginian. The opposite effect occurs around the margins of Antarctica around 130-30 • W, where there is cooling associated with a transition from open ocean to land. Almost all of the continental changes can be linked directly to orographic changes, via the lapse-rate effect (see Fig. S7). An exception is in the subtropical ocean just south of North America. Here, there is a cooling whereas from the change in topography a warming is expected (see Fig. S7r).
There is a substantial cooling in the proto-Arctic Ocean. This can be linked to the formation of an island chain at the west end of the Arctic Ocean, which constricts the transport of relatively warm ocean waters into the Arctic, and therefore cools this region. The cooling is amplified by an expansion of Arctic sea ice in the Valanginian (Fig. S8). This cooling ef- fect appears to extend beyond the Arctic, and into the North Pacific. This is also related to a decrease in ocean overturning (see Fig. S9) and in the extent and magnitude of regions of deep water formation (see Fig. S10). As such, we attribute this global cooling transition primarily to the closure of the Pacific-Arctic gateway. The timeseries of SST change indicates that the Arctic itself becomes cooler almost immediately (in the first year of the simulation) in the Valanginian. This cool anomaly then spreads southwards, increasing in magnitude over several hundred years, and is amplified at Phase 4 when the barotropic circulation is initialised.
The energy balance analysis for the Berriasian-Valanginian transition (Fig. S11r) shows that, on a global scale, changes in emissivity contribute about 60 % of the cooling, and planetary albedo changes 40 %. As CO 2 is constant across the transition, the emissivity change is a longwave water vapour and cloud feedback effect. The Northern Hemisphere cooling between 50 and 70 • N is due to a combination of emissivity and heat transport changes, whereas between 70 and 80 • N, at the latitudes of the Arctic Ocean, planetary albedo and emissivity changes dominate. (Fig. 6n) This transition is a warming of 0.77 • C in the global mean. As for the previous transition, the continental changes are dominated by a lapse-rate effect, and correlate very closely with orographic change (Fig. S7). Note that some of the largest signals, for example in eastern North America, are essentially artefacts associated with the movement of plates, which manifest as apparent warm-cold dipoles as a mountain range shifts horizontally in the model reference frame, as opposed to true tectonic effects such as uplift. An exception to the strong correlation is in the region of the Andes in South America, where there is a warming, whereas the change in topography would be expected to generate a cooling.

Aptian-Albian
In the ocean, the warmest anomalies are in the northern Pacific, and in the equatorial region that lies between S.America/N.America and Africa/Europe. There is also warming over much of the tropical and subtropical Pacific. However, in the southern Pacific there is a cooling. This is associated with a significant change in ocean circulation. In the Aptian stage, there is a region of deep-water formation off the coast of Antarctica (Fig. S10o), which is associated with a deep overturning cell in the Pacific sector (Fig. S9o). This is much weaker or nonexistent in the Albian stage. This reduction in southern deep water formation reduces surface poleward warm water transport, leading to a reduction in south Pacific temperatures in the Albian compared with the Aptian. The opposite signal in the north Pacific is likely a bipolar seesaw type response, amplified by seaice feedbacks ( Fig. S8n and o).
Again, this is supported by the energy balance analysis (Fig. S11n), which shows a cooling contribution due to heat transport changes through most of the region 40 to 80 • S, and a warming contribution 50 to 75 • N. On a global scale, 10 % of the warming is due to the solar constant increase directly (the Albian and Aptian are relatively far apart in time compared to many other consecutive Stages), with emissivity and planetary albedo feedback contributing roughly equally. (Fig. 6h) This transition is a warming of 0.74 • C in the global mean. In the ocean, there is warming globally, with the exceptions of the NE Pacific and the southern Atlantic. The largest ocean warmings are in the Pacific sector of the Southern Ocean (associated with a transition from land to ocean) and in the Indian sector of the Southern Ocean. The continental temperatures largely follow topographic change (Fig. S7h), although there is significant warming in northern Africa and western Eurasia which does not appear to be associated with topography, and may instead be related to the adjacent ocean warming.

Campanian-Maastrichtian
This Indian sector warming appears to be associated with an increase in deep water formation off the Antarctic coast in this sector ( Fig. S10h and i), likely leading to an increase in poleward heat transport from equatorial regions. Although there is some change in the overturning associated with this, it is relatively muted on the global scale ( Fig. S9h and i). The reason for the change in ocean circulation is not clear, but it may be due to the northward migration of India, allowing greater transport towards Antarctica in the Maastrichtian stage.
The important role of ocean circulation changes in the Southern Hemisphere is highlighted in the energy balance analysis (Fig. S11h), which shows a significant contribution to the warming polewards of 50 • S due to heat trans- port change. Globally, planetary albedo changes contribute 60 % of the signal, emissivity changes 30 % and solar constant change less than 10 %.

Summary of mechanisms
It appears that the three largest climate transitions are associated with changes in ocean circulation, and driven by quite subtle changes in palaeogeography. Whether climate is ultimately driving ocean circulation, or vice versa, remains difficult to assess without further sensitivity studies. However, ocean circulation changes do seem to be key (for example, the fourth largest transition (Selandian-Thanetian), is also associated with a change in mixed-layer depth ( Fig. S10e and f), and overturning stream function ( Fig. S9e and f). It should be noted that the length of the model simulations means that the deep ocean is not in equilibrium. As such, the associated mechanisms should be regarded as hypotheses at this stage (see Sect. 4). All changes are ultimately tectonically driven, but strong planetary albedo and emissivity feedbacks amplify the initial forcing.

Implications for interpretation of palaeo proxy records
Many records of CPE climate change have been developed, using a variety of proxies, e.g. using planktic δ 18 O (e.g. Friedrich et al., 2008;Bice et al., 2003;MacLeod et al., 2013;Erbacher et al., 2011;Huber et al., 1995Huber et al., , 2002, or GDGTs/TEX 86 (e.g. Bornemann et al., 2008;Jenkyns et al., 2012;Forster et al., 2007;Linnert et al., 2014;Littler et al., 2011;Inglis et al., 2015). Often the variations seen in a longduration proxy record from a single site are interpreted as being related to global phenomena, and are often linked to hypothesised atmospheric greenhouse gas and carbon cycle change. However, a component of the variations will be a local signal due to palaeogeographic change, either directly (e.g. due to lapse rate changes for terrestrial sites) or indirectly (e.g. due to ocean or atmospheric circulation change related to palaeogeographic change). Similarly, some of the change in very long records will be due to solar constant change. In addition, any proxy record will experience a change due to the horizontal movement of a site due to the movement of the underlying plate, even if the background climate is constant. This component will be particularly large if the location moves significantly latitudinally. In many cases, it would be of interest to "adjust" a record for the temperature changes associated with the local palaeogeographic components, solar components, and plate movement components, in order to leave a component which is likely to have a more global significance, likely related to greenhouse gas changes through changes in the carbon cycle. The model simulations presented here can aid in this process, by generating an "adjustment factor" which can be applied to long-term proxy records through the CPE. We illustrate this process below.
The work presented so far has been in a fixed Eulerian longitude-latitude reference frame (which resulted in the artefacts discussed in Sect. 3.3.2), but in order to generate such an adjustment factor it is necessary to use a Lagrangian frame which can take into account the effects of rotating plates. Getech Plc have provided us with palaeolongitudes and palaeolatitudes for each model gridcell and each Stage, which allows us to ascertain the palaeolocation of any modern location, consistent with the palaeogeographies used in the climate model simulations. The palaeolocations of seven key sites (Blake Nose (e.g. Huber et al., 2002), Demerara Rise (e.g. Bornemann et al., 2008), Falkland Plateau (e.g. Huber et al., 1995), Walvis Ridge (e.g. Friedrich et al., 2009), Maud Rise (e.g. Barrera and Huber, 1990), Tanzania (e.g. MacLeod et al., 2013), and the Saxony Basin (e.g. Erbacher et al., 2011)), which have previously been used to reconstruct CPE SSTs from planktic δ 18 O, is shown in Fig. 9a. This shows that all the sites move a significant distance over the course of the CPE. Note that the "oldest" location is different for each site, as some sites exist on ocean crust which was not formed until after the earliest Cretaceous (e.g. Walvis Ridge). Figure 9b shows the modelled annual mean surface air temperature at each of these seven locations, as climate changes through the CPE and as the plates underlying each site move. These variations are large; the maximum temperature in the CPE minus minimum temperature in the CPE varies from 1.2 • C at Demerara Rise, to 9.3 • C at Saxony Basin. This is in the context of a global mean modelled climate which is only varying by a fraction of a degree over this interval (Fig. 5). These modelled temperature records are the "adjustment factor" we describe above.
Some of the temporal variations in the adjustment factor, in particular at Saxony Basin, Tanzania, and Falkland Plateau, are partly due to transitions from the site being oceanic to continental ( Fig. 9b; filled circles compared with open triangles). For coastal sites, such as Tanzania, the transitions may be an artefact of the coarse resolution of the model palaeolatitudes and longitudes, which are 3.75 • × 2.5 • , and which cannot therefore distinguish correctly between land and ocean near the coast. Some sites are characterised by relatively stable modelled temperatures (and therefore small adjustment factors) over 10's of millions of years, for example Demerara Rise and Blake Nose during the late Cretaceous.
It is not possible with our current experimental design to partition the component of the adjustment factor due to solar constant from that due to palaeogeography; however, it is possible to partition the effect due to plate movements. Figure S12 shows the modelled temperature evolution over the CPE at each site in Fig. 9a, assuming either that the climate stays constant through the CPE while the site location moves ( Fig. S12a and b), or that the site location stays constant while climate varies ( Fig. S12c and d), with the constant being either that of the late Eocene ( Fig. S12a and c) or early Cretaceous (Fig. S12b and d). Comparison of Fig. 9b or Fig. S12e with Fig. S12a and Fig. S12c shows that for many sites the majority of the SST change is due to changing palaeogeography. However, in the Early Cretaceous, site movement appears to be playing a role for Falkland Plateau and Saxony Basin. Here we illustrate the adjustment process in the context of the Falkland plateau (Site 511) site. Temperature data for this site have been reconstructed using TEX 86 by Liu et al. (2009) for the latest Eocene, and by Jenkyns et al. (2012) for the middle of the early Cretaceous (Fig. 10a). For the purposes of this example, we do not also include δ 18 O data to avoid complications from comparison of different proxies. This indicates significantly warmer temperatures in the Cretaceous than in the Eocene, and could be interpreted as indicating a global cooling over this period. However, our model output over this timescale, in which the global mean temperature is almost constant, also indicates a significant cooling at this site (Fig. 10a). If this cooling related to local palaeogeography is used to adjust the proxy data, then the apparent cooling in the proxy record is greatly reduced (Fig. 10b). As such, any inferred atmospheric CO 2 changes implied by the adjusted proxy temperature record would be of lesser magnitude than that implied by the unadjusted record. It is clear that it is crucial to take into account the magnitude of this non-CO 2 component of local climate change, before proxies from single sites are interpreted in a global context. This analysis can be summarised on a global scale, indicating regions where this adjustment process is small. Figure 11a shows the total change in temperature (maximum temperature through the CPE minus minimum temperature through the CPE), for all modern locations. Small values are where the adjustment factor is constant (note that some of the locations with large values are associated with transitions between continental and marine settings). This can be interpreted as a map indicating places where future proxies could be targeted in order to reconstruct the purest CO 2 -induced temperature change, where the complicating contributions of other (palaeogeographical, solar, and plate movement) processes are minimised (low values in Fig. 11 represent the "best" regions in this context). White regions indicate that the modern crust was not present at the beginning of the CPE. Figure 11b-d are similar, but show the total change in temperature across the Eocene (d), the Eocene-Palaeocene (c), and the Eocene-Palaeocene-early Cretaceous (b). Not all of the "best" regions will have suitable sedimentary material, obviously, but, combined with all other considerations, this work could provide useful information for supporting targets for drilling localities and outcrop studies.

Discussion
It is anticipated that the results from these simulations will be of interest to the palaeoclimate data community; as such, we make the results available on our website: http://www.bridge. bris.ac.uk/resources/simulations, including variables not discussed in this paper. In addition, in the Supplement we provide the raw data which underlies Figs. 9 and 11, in netcdf (.nc) and Excel (.xlsx) format, which will allow others to develop their own adjustments, over any period in the CPE, for any site in the world.
The Eocene simulations (Ypresian, Lutetian, Bartonian, and Priabonian) described in this paper have been discussed in a previous publication (Inglis et al., 2015), as has a lower CO 2 simulation of the Priabonian simulation (Kennedy et al., 2015), and a less spun-up version of the Maastrichtian simulation (Brown, 2013). In addition, in future studies we expect to investigate many aspects of the simulations which have not been possible in the scope of this study, including the evolution of monsoon systems, ENSO, vegetation, and atmospheric circulation. Furthermore, we intend to carry out sensitivity studies, especially to CO 2 in order to investigate the evolution of climate sensitivity through geological time.
However, there are some aspects of the simulations which could be modified and improved, although we do not think that they will have a first-order effect on our results. The fol- lowing is not a complete list, but includes the main aspects that we intend to explore as we commence Phase 5 of the ensemble and beyond.
The version of the model used in this study has received little or no tuning. The internal model parameters in the atmosphere are identical to those in HadCM3 (Gordon et al., 2000), which did receive some (largely undocumented) tuning at the UK Met Office. However, compared with HadCM3, our model has a lower resolution ocean and a different land-surface scheme. In addition, the ozone correction discussed in Sect. 2.5 cools the climate somewhat. Furthermore, the subgridscale parameters derived from the Getech 0.5 • × 0.5 • resolution palaeogeographies are not necessarily consistent with those derived from higher resolution observational data sets. As such, the modern climate for this version of the model has greater biases than the HadCM3 model from which it is derived. Future work will involve tuning the model, using techniques such as those developed by Irvine et al. (2013).
In order to maintain stability in the atmosphere and ocean, some Stages received more or less smoothing of topography or bathymetry than others (see Table 1). In Phase 5, we will be more consistent, and apply the same amount of smoothing to all Stages.
There are still trends in the ocean temperatures, at all depths (see Fig. 4). Although computational constraints mean that no GCM of this complexity could currently be run to full equilibrium, and we argue that the main findings presented here will not be affected significantly by further spinup, we do aim to run Phase 5 for a further 1000 years in order to further approach equilibrium. As such, at present the adjustment factors are presented only for the sea surface temperatures, and the exact mechanisms associated with ocean overturning changes should be regarded as hypotheses at this stage.
Getech Plc provide maps of runoff basins and nodes, but these are not currently used. Instead, as discussed in Sect. 2.5, water is routed downhill according to the model resolution topography. In Phase 5 we will make use of the observational constraints on river basins and river mouths by using the Getech maps.
We have not been consistent in our definition of islands for the purposes of the barotropic circulation calculation. For example, in some Stages single gridcell islands are defined as such, and in others they are not (see Fig. S4). Furthermore, we have not defined the continent of Antarctica and Australasia as an island in the mid-Cretaceous simulations, which could affect the flow through the Tethys-Pacific seaway and Drake Passage (see discussion in Kennedy et al., 2015).
The model does not rigorously conserve water, due to the build-up of snow on polar continents, the loss of water in inland endorheic regions, and a salinity cap which affects inland basins. In the modern, a prescribed freshwater flux is applied in polar regions, in an attempt to mitigate against salinity drift in long simulations. However, in these simulations we do not apply such a correction. In Phase 5 we will diagnose a freshwater flux in order to maintain constant ocean mean salinity.
The palaeogeographies themselves have associated uncertainties; for example, the timing of the evolution of key ocean gateways, and uplift of mountain ranges, is not always well constrained. As such, it will be important to investigate the sensitivity of our results to these palaeogeographic uncertainties, in particular as it is known that they can potentially have a significant effect on surface temperatures (e.g. Foster et al., 2010;Zhang et al., 2011).
All our simulations are carried out at 4× pre-industrial CO 2 values. It has been shown that the ocean circulation is a function of CO 2 value (e.g. Lunt et al., 2010); this, coupled with likely thresholds in the system means that it is possible that the adjustment factor we present is dependent on the background CO 2 concentration. As discussed above, this is currently being explored by carrying out additional simulations at 2× pre-industrial CO 2 values.
There is undoubtedly some degree of model dependency to our results, although the extent of this is somewhat uncertain. Previous work has shown that different models can give variable results for the Eocene time period , although that was a study in which the different models were constrained by different boundary conditions. The extent of model dependency in the simulation of CPE climates is currently being explored in a consistent fashion in the frame-work of the "DeepMIP" intercomparison project (see https: //wiki.lsce.ipsl.fr/pmip3/doku.php/pmip3:wg:ppc:index).

Conclusions
1. We have carried out a set of 19 GCM simulations covering 115 million years, one for each Stage, from the earliest Cretaceous to the latest Eocene, with constant CO 2 but varying palaeogeography and solar constant ( Table 1). All simulations are within 1 W m −2 of equilibrium after more than 1400 years of simulation.
2. The global mean temperatures across the simulations are remarkably constant, with a trend of only 0.004 • C million years −1 (Fig. 5). The lack of trend results from a cancelling of effects due to changing solar constant with effects due to changing palaeogeography.
3. There is also little scatter around the trend, ∼ ±0.5 • C (Fig. 7). The scatter correlates weakly with changing land area, indicating the albedo contrast between land and ocean may play a role; continental temperatures correlate weakly with mean orography, indicating lapse rate and area above snowline also may be playing a role (Fig. 8). Energy balance analysis indicates that the solar and palaeogeographic forcing is amplified by planetary albedo and emissivity feedbacks.
4. The largest Stage-Stage transitions through the CPE are associated with changes in the mode of ocean circulation. For example, the largest transition, Berriasian-Valanginian, is associated with a reduction in deepwater formation in the North Pacific, and a reduction in the meridional positively overturning cell; the second largest transition, Aptian-Albian, is associated with a reduction in deepwater formation off the coast of Antarctica, and a reduction in the negatively overturning cell. In some cases, these ocean circulation changes can be directly related to palaeogeographic change associated with gateway opening or closing, for example the isolation of the Arctic at the Berriasian-Valanginian transition.
5. Although the global mean changes are relatively small across the CPE, local temperature changes are much larger (Fig. 11). This has implications for interpretations of proxy records. In particular, our results allow the non-CO 2 (i.e. palaeogeographic and solar constant) components of proxy records to be removed, through the application of an adjustment factor, leaving a global component associated with carbon cycle change. This adjustment factor is illustrated for seven key sites in the CPE (Fig. 9), and applied to proxy data from Falkland plateau, and data provided so that similar adjustments can be made to any site and for any time period within the CPE.
6. Regions where the adjustment factor is constant throughout the CPE could indicate places where future proxies could be targeted in order to reconstruct the purest CO 2 -induced temperature change, where the complicating contributions of other processes are minimised. Therefore, combined with other considerations, this work could provide useful information for supporting targets for drilling localities and outcrop studies.

Data availability
In the supplement to this paper we provide the following: the file values.nc is a netcdf file which contains a 96 × 73 × 19 level field. The 96 × 73 are the longitudes × latitudes of the model, which have a 3.75 • × 2.5 • resolution. The 19 represents the 19 stratigraphic stages of the CPE. The field itself represents the modelled annual mean 1.5 m near-surface air temperature of each modern grid cell, for that particular stratigraphic stage. The plots in Figs. 9b and 10 in the main paper can be generated from the data in this file. The file masks.nc is very similar, but instead of temperature it has 0 for ocean and 1 for land, to indicate whether a particular site in the modern was marine or terrestrial at a certain period in the past. The file values.xlsx is an Excel version of values.nc. Each time period is a separate sheet, and the overview sheet allows you to extract the values and make plots of the adjustment factor very quickly for any location. The default is for Falkland plateau, which reproduces the data in Figs. 9b and 10. To choose a different site, just edit cell C4. The files land_ypr.nc, orog_ypr.nc, and bath_ypr.nc contain the model-resolution land-sea mask, orography, and bathymetry respectively, in digital form, for the Ypresian stage (early Eocene). If you use these as boundary conditions for model simulations, then please cite this paper, and acknowledge Getech PlC, who produced the original palaeogeographies. In addition, those interested may explore our results in further detail here: http://www.bridge.bris.ac. uk/resources/simulations, following the link to "Simulations featured in papers".
The Supplement related to this article is available online at doi:10.5194/cp-12-1181-2016-supplement.