Articles | Volume 17, issue 5
Research article
21 Oct 2021
Research article |  | 21 Oct 2021

Data-constrained assessment of ocean circulation changes since the middle Miocene in an Earth system model

Katherine A. Crichton, Andy Ridgwell, Daniel J. Lunt, Alex Farnsworth, and Paul N. Pearson

Since the middle Miocene (15 Ma, million years ago), the Earth's climate has undergone a long-term cooling trend, characterised by a reduction in ocean temperatures of up to 7–8 C. The causes of this cooling are primarily thought to be due to tectonic plate movements driving changes in large-scale ocean circulation patterns, and hence heat redistribution, in conjunction with a drop in atmospheric greenhouse gas forcing (and attendant ice-sheet growth and feedback). In this study, we assess the potential to constrain the evolving patterns of global ocean circulation and cooling over the last 15 Ma by assimilating a variety of marine sediment proxy data in an Earth system model. We do this by first compiling surface and benthic ocean temperature and benthic carbon-13 (δ13C) data in a series of seven time slices spaced at approximately 2.5 Myr intervals. We then pair this with a corresponding series of tectonic and climate boundary condition reconstructions in the cGENIE (“muffin” release) Earth system model, including alternative possibilities for an open vs. closed Central American Seaway (CAS) from 10 Ma onwards. In the cGENIE model, we explore uncertainty in greenhouse gas forcing and the magnitude of North Pacific to North Atlantic salinity flux adjustment required in the model to create an Atlantic Meridional Overturning Circulation (AMOC) of a specific strength, via a series of 12 (one for each tectonic reconstruction) 2D parameter ensembles. Each ensemble member is then tested against the observed global temperature and benthic δ13C patterns. We identify that a relatively high CO2 equivalent forcing of 1120 ppm is required at 15 Ma in cGENIE to reproduce proxy temperature estimates in the model, noting that this CO2 forcing is dependent on the cGENIE model's climate sensitivity and that it incorporates the effects of all greenhouse gases. We find that reproducing the observed long-term cooling trend requires a progressively declining greenhouse gas forcing in the model. In parallel to this, the strength of the AMOC increases with time despite a reduction in the salinity of the surface North Atlantic over the cooling period, attributable to falling intensity of the hydrological cycle and to lowering polar temperatures, both caused by CO2-driven global cooling. We also find that a closed CAS from 10 Ma to present shows better agreement between benthic δ13C patterns and our particular series of model configurations and data. A final outcome of our analysis is a pronounced ca. 1.5 ‰ decline occurring in atmospheric (and ca. 1 ‰ ocean surface) δ13C that could be used to inform future δ13C-based proxy reconstructions.

1 Introduction and background

Since the middle Miocene (∼15 Ma), the Earth has experienced a pronounced and quasi-continuous period of global cooling, characterised by an expansion of ice sheets over Antarctica and the subsequent establishment of the Greenland ice sheet (Zachos et al., 2008; Cramer et al., 2011). The Pleistocene (2.6 Ma to present) also saw the intensification of glacial–interglacial cycles and the episodic establishment of the North American ice sheet. Terrestrial temperature proxies indicate that the Miocene was significantly warmer than the present day (Pound et al., 2012, and references therein). Marine data also indicate a significantly warmer-than-present Miocene climate, with global surface ocean temperatures 6 C warmer than present (Stewart et al., 2004; Herbert et al., 2016) and global deep-ocean temperature 4 to 6 C warmer than present (Cramer et al., 2011). However, early proxy estimates of atmospheric CO2 in the mid-Miocene show that levels may have been similar to or even lower than present (Pagani et al., 2005, Pearson and Palmer 2000). Climate modelling studies have tended to apply a forcing of around 400 ppm for the general Miocene period (but not necessarily to represent peak mid-Miocene warmth) (Burls et al., 2021; Londoño et al., 2018, and references therein). More recent work places mid-Miocene CO2 higher still, from ∼500 ppm to ∼1000 ppm, although uncertainty still exists (Londoño et al., 2018; Sosdian et al., 2018; Stoll et al., 2019). Regardless, falling atmospheric CO2 is thought to have been a driver of global cooling since the middle Miocene (Rae et al., 2021).

Associated with this interval of global cooling, fundamental changes have also occurred in large-scale ocean circulation patterns (Butzin et al., 2011). Specifically, the Atlantic Meridional Overturning Circulation (AMOC), which today redistributes heat to the Northern Hemisphere at the expense of the South, is thought to have established in its current form sometime after the middle Miocene (Sepulchre et al., 2014; Bell et al., 2015) (although this was not necessarily its first appearance; see Abelson and Eres 2017). This “switching on” of a prominent AMOC is suspected to have been linked to the closing of the seaway between the Atlantic and the Pacific with the creation of the Isthmus of Panama (Lunt et al., 2008; Montes et al., 2015; O'Dea et al., 2016; Jaramillo et al., 2017) that finally closed the Central American Seaway (CAS). Other seaways have also been tectonically transformed since the mid-Miocene, particularly the disappearance of the Tethys Sea due to the northward movement of Africa (Rogl et al., 1999; Hamon et al., 2013; Bialik et al., 2019), the restriction of the Indonesian seaway with the northward movement of Australia (Srinivasan et al., 1998), and the widening of the Drake Passage driven by the northward movement of South America relative to Antarctica (Lagabrielle et al., 2009). Associated with these plate movements, the mid-Miocene to Holocene was characterised by significant mountain building which may have played a direct role in the drawdown of atmospheric CO2 via weathering and, hence, a key role in the progressive cooling (Filippelli et al., 1997; Raymo et al., 1988), although this is far from certain (Steinthorsdottir et al., 2020).

Putting all of these changes together, Miocene climate modelling efforts have found that reconciling the combined constraints of ocean temperature, CO2 indicators, and Antarctic ice-sheet dynamics is a non-trivial task (Micheels et al., 2009; Henrot et al., 2010, Bradshaw et al., 2012; Sijp et al., 2014). Vegetation feedbacks have also likely been integral in creating the Miocene climatic conditions (Henrot et al., 2010; Knorr et al., 2011; Micheels et al., 2011; Bradshaw et al., 2015) in addition to altered bathymetry, topography, and atmospheric CO2 as drivers of a warmer climate.

In this paper, we aim to step back from the spatial and dynamical complexities of the land surface and terrestrial biosphere and their interaction with the atmosphere, and we instead focus on exploring the extent to which proxy data can constrain changing ocean circulation and temperature patterns in conjunction with a simplified forcing of surface climate, using an Earth system model (of “intermediate complexity”). We do not provide an exhaustive sensitivity study of Miocene seaways or details of the marine carbon cycle here. For an overview of Miocene circulation and the effects of seaways on that circulation, the reader is referred to studies such as Butzin et al. (2011) and Sijp et al. (2014). A full review of δ13C and δ18C as ocean circulation tracers used to infer seaway and global ice volume and climate changes in the Cenozoic is available in MacKensen and Schmiedl (2019), and a review of Miocene climate and data indicators is given in Steinthorsdottir et al. (2020). Instead, we aim to combine constraints of temperature and ocean circulation to create plausible and self-consistent palaeo-realisations of climate and ocean circulation for seven respective time slices spanning the middle Miocene to the late Holocene. Hence, our focus in this paper is on the model–data methodology as well as the self-consistency and plausibility of the outcome.

2 Model–data methodology

Our overall methodology is based on creating a suite of data sets that can help constrain the large-scale patterns of ocean circulation and temperature (and, hence, climate in a broad sense). We compile and employ foraminifera proxy data for surface temperature, benthic δ18O (to calculate benthic temperature), and benthic δ13C. The surface and benthic (seafloor) temperature data allow us to evaluate the model skill in reproducing ocean heat distribution and, hence, the model's simulated global-scale pattern of ocean circulation. It also provides a means of determining the atmospheric forcing (in terms of CO2) that produces surface and deep-ocean temperatures that best match the data. We further compare modelled benthic global mean ocean temperature to a published estimate from Cramer et al. (2011) that used the Mg/Ca thermometer. The isotope of carbon (carbon-13, or δ13C) has been used as a tracer for palaeo-ocean circulation for many years (Lynch-Stieglitz, 2003), and we also employ it here as a further circulation constraint. Different water masses have characteristic δ13C signatures, which depend on biological activity (carbon pumps in the ocean) and, critically for this study, also on ocean circulation patterns. Rather than attempt a continuous model–data analysis for the entire 15 Ma long interval, we discretise the mid-Miocene to the present (late Holocene) (Table 1) in a series of seven time slices.

Table 1Data slices for the approximate carbon cycling locations (in time), and the HadCM3L set-ups used to create cGENIE bathymetry. The 7.5 and 2.5 Ma time slices are just outside the ages to which they are assigned; however, in order to not have two points in the other ages, they are placed in the Messinian and Piacenzian respectively.

Download Print Version | Download XLSX

For convenience, uncertainty in model climate forcing is assumed to be represented by atmospheric CO2 – this should be understood as an “equivalent” CO2 forcing that encompasses changes in all atmospheric greenhouse gases (especially methane). Ocean circulation in the model is a function of surface ocean boundary conditions (especially wind stress) and climate (including atmospheric CO2 forcing), ocean bathymetry, and the existence and nature of seaways and gateways – all of which we adjust at each time slice (along with the solar constant). We further introduce a variable (parameter uncertainty) to the model that alters the salinity transfer between the North Pacific and North Atlantic – a classic “flux adjustment”. This represents the effect of atmospheric moisture transport and the relative precipitation that east- vs. west-draining watersheds receive on the North American continental land mass, none of which can be resolved well in the simple 2D energy–moisture balance model-based approach employed here (Edwards and Marsh, 2005; Marsh et al., 2011). This “flux correction” results in an Atlantic Meridional Overturning Circulation (AMOC) being induced (and/or strengthened) in cGENIE and, hence, represents a parameter “knob” in trying to find a fit to ocean proxy data. (Indeed, it is required in the standard present-day configuration (Ridgwell et al., 2007; Marsh et al., 2011) in order to rectify the aforementioned simplification in modelled atmospheric moisture transport and dynamics.)

2.1 Data sets

We selected published surface ocean temperature data based on either alkenones or TEX86 for all seven slices (with the exception of two data points at 15 Ma), noting that our selected proxies, like any proxy, are still subject to uncertainties and limitations, such as seasonality effects, surface/subsurface water representation, or non-thermal influences (discussed in Richey and Tierney, 2016, for U37K and TEX86). Figure 1 shows surface temperature data locations, and Supplement A contains the data set. For slices 2.5 to 15 Ma, published benthic δ13C and δ18O data were selected only for the Cibicidoides and Planulina foraminifera species. These species were selected so that temperature could be calculated from shell δ18O using the linear model from Marchitto et al. (2014) in order to reduce uncertainty driven by species-specific effects on shell isotopic signatures. We ensured that at least 15 data points were available for each time slice and that coverage included the Atlantic and Pacific basins (Fig. 1 shows the benthic data locations, Fig. 2 shows the δ13C data plotted by palaeolatitude, and Supplement B provides the data set). Final benthic temperatures calculated from δ18O consider the effect of benthic water salinity on δ18Osw which, in turn, is influenced by ocean circulation (the calculated temperatures in Table S2 in the Supplement are uncorrected for salinity), and we use the modelled benthic salinity field to create the correction for the local water δ18Osw. The method used in this calculation is described in detail in Appendix A. Palaeo-locations for each data point were found using the reconstructions from (last access: 11 October 2021) (Urban and Hardisty 2013) that provide Zanclean (4.466 Ma), Tortonian (9.427 Ma), and Burdigalian (18.2 Ma) palaeo-locations using PLATES reconstructions (University of Texas, 2015). Thereafter, the palaeo-locations for each data point at our time slices (Table 1) were interpolated using these three points and the present-day location in a cubic model regression.

Figure 1Present-day location of surface ocean temperature data points (a) and benthic δ18O and δ13C data points (b) collated for this study (the Holocene δ13C data set from Peterson et al., 2014, is not plotted).

Figure 2Benthic δ13C data and palaeo-locations overlaid on palaeo-coastlines (left column). Benthic δ13C data plotted by latitude for each time slice, showing the linear regression and R2 for the Atlantic and Pacific basins (right column). The colour scale for the map is also the y axis scale for the per-latitude data, and all axes have a range of 1.8 ‰.


To account for age model uncertainties in the isotope data, a window of ±1 Myr was inspected for each δ18O and δ13C data point to ensure that the data value was not unrepresentative of the general value at that site around the target age (rather than taking a mean value over a certain window of time). This approach was applied to try to ameliorate the effect of uncertainties in age models, particularly for periods in which isotope values change significantly over the 2 Myr interval. In data sets with high resolution (i.e. showing possible orbital-type climate cycling) and generally better age models, an intermediate value between the high and low nearest to the target age was selected (thereby targeting mid-climate states). All of the data points and their palaeo-locations as well as the ± million-year time series plotted for each data location are available in Supplement C, Fig. S1.

We created data sets for the late-Holocene time slice slightly differently. For this, we took the Holocene benthic δ13C from Peterson et al. (2014). For ocean temperature, we used data from the World Ocean Atlas (WOA) 2009 (Levitus et al., 2010), assuming that deep-ocean temperature has changed little between the late Holocene and the present. Data points for the deep ocean (at 4000 m) were extracted from WOA 2009 to create the benthic temperature data set for model–data fit statistics (only data locations where the cGENIE model's seafloor depth closely matched the WOA 2009 seafloor depth were used).

2.2 The cGENIE.muffin Earth system model

We run the open-source intermediate-complexity Earth system model cGENIE (“muffin” release) (see the “Code availability” section for details). As implemented in our study, this comprises (1) a 3D ocean circulation model configured on a 36×36 equal area grid with 16 non-equally spaced vertical levels in the ocean, (2) a 2D energy–moisture balance atmospheric model (EMBM) component, and (3) a 2D dynamic–thermodynamic sea-ice model component. These three individual components and their coupling are described in Marsh et al. (2011, and references therein). The basic physics parameter calibration of the climate model component is as per Cao et al. (2009) unless described otherwise. For the modern Atlantic Basin configuration (and, as per the findings of this paper, also over the past 15 Ma), the lack of both atmospheric dynamics and explicit topography on land in our configuration model requires that a freshwater flux adjustment is applied. This is implemented in the modern configuration of the model, principally by transferring salinity to the North Atlantic region from the North Pacific as described in Edwards and Marsh (2005), although we deviate from this flux pattern as described in the following.

Our implementation of the cGENIE model includes a relatively complete description of the cycling of carbon and oxygen in the ocean as well as exchange with the atmosphere, as described in Ridgwell et al. (2007) and Cao et al. (2009). In addition, the carbon isotopic (δ13C) composition of all of the carbon pools as well as associated fractionations between them are explicitly account for, as described in Ridgwell et al. (2007) and with additional description and evaluation in Kirtland Turner and Ridgwell (2016). Given the paucity of constraints on the evolving patterns of aeolian fluxes to the ocean surface, we omit a marine iron cycle and include only a single nutrient (phosphate) as potentially limiting to biological productivity in the ocean in our configuration. We also employ an idealised organic matter export scheme and, hence, biological carbon pump (Cao et al., 2009) and do not attempt to simulate ecosystem composition and dynamics (as per, for example, Wilson et al., 2018).

In the original modern configurations of the cGENIE model, ocean bathymetry and land–sea mask grids as well as files supporting the calculation of ocean circulation files (defining “islands” and circulation paths around islands) are derived from filtered global topographic observations (ETOPO5) (Edwards and Marsh, 2005). Furthermore, because the atmospheric EMBM component lacks clouds and dynamics (e.g. winds), additional spatial fields must be provided as fixed (annual average) boundary conditions. Firstly, a zonally averaged planetary albedo profile is applied (originally as a simple cosine function of latitude as per Edwards and Marsh, 2005). Secondly, fields for (a) vectors of wind stress on the ocean surface, (b) vectors of surface wind velocity in the atmosphere, and (c) short-term wind speed are derived from modern observations and applied to (a) driving surface ocean circulation, (b) transporting heat and moisture in the atmosphere and for calculating heat and moisture exchange between the ocean surface and atmosphere, and (c) in calculating air–sea gas exchange respectively.

For this study, we create a series of new continental and surface boundary condition configurations of the cGENIE model for the seven time slices spanning the late Holocene through to mid-Miocene. As per previous (deeper time) palaeo-applications of the cGENIE.muffin model (e.g. Ridgwell and Schmidt, 2010), rather than observations, we derive the required boundary conditions from a representative fully coupled GCM (general circulation model) experiment. To create the Miocene slices' boundary conditions, we re-grid the GCM HadCM3LM2.1aE (hereafter referred to as HadCM3L) model configuration to the cGENIE grid using the method described in Appendix B, including bathymetry, wind stress and velocity, and planetary albedo. Due to the re-gridding process, this results in an open Central American Seaway (CAS) from 15 Ma up to and including the present (late Holocene). To test for the role of closing the CAS, we create a parallel series of an additional five time slices spanning 10 Ma to present, in which we force a fully closed CAS in the cGENIE model grid. (Although the CAS is not open today, we include an open CAS for the late-Holocene time slice for completeness and to test for the ability of benthic δ13C observations to distinguish between the two alternatives.) The resultant cGENIE bathymetry for each time slice is shown in Fig. 3 for all 12 (7 open-CAS and 5 closed-CAS) continental configurations, with a close-up section around the CAS and a cross section of the shallowest CAS depths for each open CAS configuration. We compare our PLATES-derived palaeo-locations to those derived from the HadCM3L in Fig. B1.

Figure 3Model bathymetry and coastlines for each time slice, and a close-up of the region covering the Central American Seaway (CAS), including a transect showing the most depth-restricted pathway for the open-CAS configurations and the ocean depths for both open- and closed-CAS configurations.

The impact of changing grid resolution and swapping the HadCM3L physics for that of the cGENIE model on the ocean circulation is included in Fig. B2 in the form of a comparison of modern (late-Holocene) simulations. HadCM3L generally shows much stronger currents, such as the Antarctic Circumpolar Current (ACC), than cGENIE, with the finer-resolution HadCM3L enabling features like the western boundary current in the Atlantic to be reasonably resolved, whereas it appears weak and diffuse at much lower resolution in cGENIE. In addition to the lower-resolution model grid employed in this study, the frictional geostrophic approximation physics employed in cGENIE (Edwards et al., 1998) also tends to dampen the wind-driven circulation. However, the fidelity of the large-scale patterns of ocean circulation in cGENIE is sufficient to support the simulation of anthropogenic carbon uptake and deep-ocean radiocarbon on par with other 3D ocean-circulation-based Earth system models (Cao et al., 2009) as well as generally plausible distributions of dissolved nutrients and oxygen (e.g. Crichton et al., 2021).

2.3 Model experimental design

In a “perfect” palaeo-model (which does not exist) and under the correct greenhouse gas boundary conditions (which are poorly constrained), climate and ocean circulation would be a correct and emergent property of the model and would provide an exact (within proxy measurements and calibration error) match to the data. Rather, the premise of our experimental design and study is that we can find a specific combination of climate state and pattern of global ocean circulation that reasonably accounts for the observed distribution of ocean temperature proxies and benthic δ13C. In cGENIE, we consider CO2 as a primary uncertainty in the model – not only in terms of uncertainty in the real past value but also with respect to the radiative forcing required in cGENIE to generate appropriate ocean temperatures. Additionally, because there is no prior expectation that the emergent state of ocean circulation is at all correct for a “correct” surface climate in cGENIE, we consider the salinity flux adjustment (hereafter “FwF”) as an additional “unknown” and also as a means to an end in adjusting basin-scale (to global-scale) patterns and the strength of circulation (specifically Atlantic overturning circulation). By varying both controlling factors (parameters), we seek to find the circulation pattern and climate (temperature) state that best reproduces the data. It should be noted that we ignore the strength of the biological pump as an additional and independent control on benthic δ13C for this study, and in our particular model configuration, remineralisation rates of organic matter in the ocean interior are not dependent on local temperature. Similarly, we assume the modern ocean nutrient (phosphate) inventory throughout all time slices.

For each of the 12 (7 open-CAS and 5 closed-CAS) time slices and model configurations, we carry out a 2D parameter space (CO2 vs. FwF) sweep via an ensemble of model experiments. In each ensemble, we test (radiative forcing equivalent) atmospheric CO2 values of 280, 400, 560, 800, 1120, and 1600 ppm, and vary the salinity flux correction applied to the North Atlantic (FwF) between 0.0 and 0.7 Sv, in increments of 0.1 Sv, for a total of 48 members in each ensemble. The CO2 values are chosen as either simple multiples of 280 ppm (pre-industrial/late Holocene), or commonly assumed GCM values (e.g. the 400 ppm in Farnsworth et al., 2019) and multiples thereof. We chose the range of FwF values to encompass the equivalent tuned modern cGENIE model value of 0.32 Sv (Edwards and Marsh, 2005). Furthermore, because the land–sea mask progressively changes across the seven time slices, we create a common mask for the purpose of salt transfer: first, we identify all of the grid points that are ocean (“wet” points) across all configurations and that lie northwards of the northernmost extent of the CAS, and we add salinity to these; second, we identify all of the grid points outside of the North Atlantic that are ocean in all configurations, and we remove salinity from these; finally, the total FwF value is divided evenly across all (equal area) grid points in the North Atlantic (positive salinity input) and “elsewhere” (negative salinity input) such that global ocean salinity is always conserved. The mask is shown in Fig. B3.

For each ensemble member and each time slice (a total of 48×12=576 model simulations), we spin-up the cGENIE model for 10 000 years. In the absence of appreciable inter-annual variability and unlike in fully coupled GCM experiments, multi-decadal averaging is not necessary in the cGENIE model; thus, we take the last annual average (year 10 000) of the simulation in order to carry out the model–data comparison.

The one caveat and complication regarding how we conduct the model–data comparison is that, in running the parameter ensembles, we identified self-sustained oscillations in global ocean circulation in a subset of the simulations that affected both mean benthic temperature (up to several C) and δ13C (several tenths of a ‰). These oscillations were of varying period and magnitude and occurred only in simulations with sea-ice cover present (and, hence, in the lower range of eq.CO2 forcing), and generally for the low to middle range of FwF values (see Appendix C and Fig. C1 for more details). Whilst this raises extremely interesting questions about past ocean circulation dynamics, it is not the focus of this study (although it will be followed up in a subsequent study). Here, we need to identify a “representative” ocean state of the model in order to carry out the model–data comparison consistently. Therefore, to create a representative state for each CO2–FwF combination, we parsed the ensemble cGENIE model output, identifying ensemble members characterised by self-sustained oscillations. For these ensemble members, we identified the period of the oscillation and averaged the model output over one full period, starting the average from the end of the 10 000-year spin-up and working back (towards the start of the model experiment). Both unmodified run-end and reconstructed mean annual averages were then treated exactly the same in terms of carrying out model–data comparison.

2.4 Model–data comparison

The collated data sets were used as an observational constraint to determine which combination of CO2 and North Atlantic salinity flux correction (FwF) produces the best fit ocean state to the data. This is performed quantitatively by statistically comparing local model output with the data points for all three data sets and combining these to produce a final “best fit” model setting for each time slice. The statistical methods applied are as follows: (1) the difference between the mean of the data set and the mean of the model output at the data locations, producing an offset value or mean bias, and (2) an overall measure of goodness of fit of the model to the data set, known as the “M-score” (Watterson, 1996). Where several data points are located within one model grid square, the mean of the data values is used for comparison with the model value. For the benthic data, each data point is assumed to be on the ocean floor, so the model's bathymetry determines the data depth, and the data value is compared with the model value in the deepest water in that location.

3 Results

3.1 Global mean and spatial data constraints

We start by assessing the model ensemble members against global mean observational constraints. In Fig. 4, we show the global mean benthic temperature data estimate from Cramer et al. (2011) plotted on top of the modelled global mean benthic (seafloor arithmetic mean) ocean temperature for our ensemble. Both CO2 and FwF have a strong effect on benthic temperature. Bathymetry, albedo, and wind fields (which differ in our time slices; see Appendix B) also have a direct effect on benthic ocean temperature, with 15 Ma having a tendency toward a warmer deep ocean compared with other time slices, even for the same CO2 and FwF in the open-CAS cases. A range of combinations of CO2 and FwF would satisfy the global mean benthic temperature constraint (shown as the dashed white line in Fig. 4, with the range shown using thin white lines). The FwF affects North Atlantic surface salinity and, therefore, has an impact on thermohaline circulation; higher FwF values create more saline surface North Atlantic waters of lower buoyancy, transporting more of the surface heat to the deep ocean. This results in combinations of higher CO2–low FwF or lower CO2–high FwF being equally possible to satisfy the benthic global mean temperature constraint, although with a general trend from 15 Ma towards the Holocene of decreasing CO2 and FwF (as the dashed line moves towards the top left in Fig. 4) for both the open-CAS and closed-CAS cases. For almost all slices, an elbow-type shape is present in the modelled benthic temperatures. As CO2 decreases for a given FwF, the global mean benthic temperature is either stable or even increases at a certain point (Fig. 4). At CO2 levels lower than this “elbow”, the cooler surface ocean favours the sinking of waters at high latitudes, supporting the thermohaline circulation. At 12.5 Ma, the elbow turning-point CO2 value is around 800 ppm, and this gradually decreases to about 400 ppm in the Holocene (for mid-range FwF). At CO2 higher than this elbow, there is only a weak (or shallow) AMOC. In the closed-CAS cases, this elbow is less evident than for the open-CAS cases, which may suggest that an open CAS results in a climate more (or differently) sensitive to CO2 changes than a closed CAS, via AMOC.

Figure 4Modelled benthic temperature with the global mean temperature estimate (and minimum and maximum range) (Cramer et al., 2011) overlaid (dashed line and solid white lines respectively). The global mean temperature shown for the Holocene time slice denotes 0.1 Ma, representing a mean value for the Pleistocene glacial–interglacial climate.


The results of the M-score spatial statistical comparison (the closer the value is to 1, the better the model can simulate the data) for the three individual data sets and all model ensemble members are shown in Fig. 5, where the FwF from 0.0 to 0.7 is plotted sequentially for each eq.CO2 forcing.

Figure 5The M-score for the fit of the model to data for surface temperature, benthic temperature, and benthic δ13C. The model settings are plotted on the y axis (totalling 48 simulations per chart).


The surface temperature constraint is rather insensitive to FwF, with eq.CO2 being the main controller on the fit to data. Whether the CAS is open or closed also has very little effect on the ability of the model to reproduce the data. The highest M-scores indicate that a general fall in eq.CO2 forcing results in the best fit to the surface temperature data, with generally higher M-scores achieved in the more recent time slices. The 12.5 and 15 Ma time slices have low scores (less than 0.4) compared with other time slices for surface ocean temperature, indicating that, even at 1600 ppm, the model cannot reproduce the data-indicated surface ocean temperatures very well.

For the benthic temperature data set, as for the global mean ocean temperature (Fig. 4), both eq.CO2 and FwF control deep-ocean heat distribution, with the highest M-scores again generally showing decreasing CO2 and decreasing FwF towards the present. The open-CAS case results in higher M-scores (i.e. better fit to data) than the closed-CAS case, even for the Holocene. As we know that the CAS is closed in the Holocene, these higher M-scores may be the result of a model bias.

The δ13C M-scores for the Holocene show higher values for a closed CAS than for an open CAS, indicating that δ13C may be a better indicator than benthic temperature for whether CAS is open or closed in our cGENIE model configurations. For the 10 Ma slice, a closed CAS better fits benthic δ13C data for a higher CO2 value, and if gradually decreasing CO2 is assumed (from the surface temperature data) toward the present, a closed CAS generally results in higher M-scores for all other time slices as well. A shift occurring in the patterns of benthic δ13C can also be seen in the raw data. The latitudinal trend for δ13C data (dotted lines) in Fig. 2 shows that the Atlantic appears to switch trend between 12.5 and 10 Ma. Prior to 12.5 Ma, the North Atlantic data indicate a generally more negative δ13C than in the South Atlantic – a situation that reverses from 10 Ma onwards when more positive δ13C values tend to occur in the North Atlantic (although R2 is low). In the Pacific Ocean, the δ13C gradient from south to north tends to intensify towards the present, with the largest gradient seen at 2.5 Ma. Overall, the range of δ13C values increases from 15 Ma towards the present. (Note that the colour bar values change, but the range – the difference between the minimum and maximum value – is 1.8 ‰ for all time slices in Fig. 2.)

3.2 Combining the constraints

All proxy data constraints combined are shown in Fig. 6, the surface and the deep-ocean temperature constraints (showing both M-scores and regions where the bias in deep-ocean temperature is less than 2 C), and the global mean benthic temperature estimates are mapped onto the statistical fit surface (M-score) of model vs. observed δ13C. We mark the best fit selection with a star and include a range of estimates delineated by whiskers on the best fit CO2 and FwF values. Note that the best fit is only identified from those specific values for which we ran the model ensemble members (i.e. we do not interpolate between settings).

Figure 6All proxy data constraints combined. The M-score for the modelled benthic δ13C fit to the data set is shown using filled contours (grey). The modelled benthic global mean temperature nearest to that identified in Cramer et al. (2011) is shown using a dashed white line, with the maximum and minimum estimates shown as solid thin white lines. The benthic temperature constraint is shown using blue; the blue filled region is where modelled benthic temperature has a bias of less than <2C with respect to the δ18O-derived temperature data (or the WOA 2009 for the Holocene), and the blue contour indicates the M-score, with the value given in blue. The M-score contour for the modelled surface ocean temperature fit to data is shown using red shading, with the value marked on the edge of the contour. The overall best fit settings are shown using a white star, with the upper and lower range estimates shown using whiskers.


It is clear from this that the primary constraint on atmospheric eq.CO2 in the model in this study is the proxy-reconstructed surface ocean temperature. This is much less sensitive to the FwF value than the other benthic ocean constraints. This is because surface ocean circulation patterns are largely dictated by the wind stress forcing, which we do not vary within any single time slice. Given a surface ocean circulation pattern that is relatively immune to changes in the applied FwF, the ensemble member eq.CO2 value and the surface climate state control the mean and pole-to-Equator gradient of sea surface temperature (SST) and, hence, the model–data surface temperature fit. Thus, for all of the following model–data fit analyses, we start by considering SSTs and the choice of eq.CO2 for each respective time slice and then bring in additional constraints and consideration of the value of FwF.

A priority for selecting good combinations of eq.CO2 and FwF values is to ensure that the global ocean heat distribution shows reasonable agreement with our proxy temperature data. The benthic global mean temperature estimate of Cramer et al. (2011) is then used as an ancillary guide for our benthic temperature set derived from δ18O and the salinity correction; in most cases, these show quite good agreement with one another (with 15 Ma and both of the Holocene cases showing the lowest agreement; note that the Cramer et al. (2011) estimate applied for the Holocene actually represents a mid-climate state rather than the interglacial, which explains the difference there). It should be noted that the global ice-volume-linked global δ18Osw value that we use in our temperature calculation is derived from the Mg/Ca temperature data set from Cramer et al. (2011) (Table S2 in the Supplement), so these are not fully independent data sets. The M-score for the δ13C is used to determine the direction of any adjustment needed when SST and benthic temperature constraints have low agreement for the best eq.CO2–FwF combination, especially for particularly high δ13C M-scores (e.g. for the 10 Ma open-CAS case, where benthic temperature suggests an FwF value of 0.1 Sv and an eq.CO2 value of 1120 ppm, but the δ13C M-score strongly increases in the direction of lower CO2 and higher FwF values). More details on the selected combinations are available in Appendix D. All of the selected settings for eq.CO2 and FwF and their ranges are summarised in Fig. 7 (and marked in Fig. 5), along with some recent data estimates of atmospheric CO2. Also plotted in Fig. 7 are global mean ocean surface and benthic temperatures as well as sea-ice extent and thickness for each time slice, showing the modelled evolution from 15 Ma to the Holocene.

Figure 7The eq.CO2 and FwF settings (and estimated ranges) that agree with temperature and benthic δ13C data constraints according to this study. Plotted with eq.CO2 are estimates for atmospheric CO2 levels from two recent studies on tropical alkenones (Stoll et al., 2019) and boron (Sosdian et al., 2018). Also shown are modelled global mean surface ocean temperature, global mean benthic ocean temperature (and the Mg/Ca thermometer estimate from Cramer et al., 2011), modelled global sea-ice extent, and modelled global mean sea-ice thickness for the best fit settings.


Modelled surface ocean temperature drops with decreasing eq.CO2; deep-ocean (benthic) temperatures also fall, but they show some dependence on FwF at the same time (Fig. 7). In the highest FwF case (4.5 Ma, open CAS), despite a halving of eq.CO2 compared with 7.5 Ma, deep-ocean temperatures are only slightly lower than 7.5 Ma. This is due to more saline North Atlantic waters and cooler polar temperatures (demonstrated by the presence of more sea ice) promoting the sinking of these surface waters to the deep, thereby delivering relatively more surface heat to the deep ocean at 4.5 than at 7.5 Ma. Conversely, in the 2.5 Ma open-CAS case, FwF is lower than at 4.5 Ma, but eq.CO2 is the same. Here, less of the surface heat is delivered to the deep ocean, resulting in a cooler deep ocean but a warmer surface ocean and lower sea ice than at 4.5 Ma (for the best fit setting). Sea ice is present from eq.CO2 levels of 800 ppm from 10 Ma, generally increasing with respect to both extent and thickness with decreasing eq.CO2 (with the exception of the previous example), irrespective of whether the CAS is open or closed.

It should be noted that the imposed eq.CO2 does not change due to changes in ocean carbon distribution or carbon exchange with the atmosphere; we apply a flux correction such that atmospheric CO2 is always restored to the eq.CO2 value that we initially impose; therefore, we do not attempt to attribute the source of CO2 decline since 15 Ma to any particular mechanism.

3.3 Atmospheric and global mean benthic carbon-13 ratios

Using the CO2 FwF settings identified from the combined constraints, we can derive model-predicted estimates of the trends in atmospheric δ13CO2 and of global mean benthic δ13C over the past 15 Myr. In the ensembles, the model atmosphere is forced with a δ13CO2 value of −6.5 ‰ for all simulations; thus, by determining the mean bias between the modelled benthic δ13C points and the data benthic δ13C points, we can derive the δ13CO2 value that would produce the lowest overall model–data error (thereby essentially inverse modelling the δ13CO2 value). This is shown in Fig. 8. This diagnosed history of atmospheric δ13CO2 can be used as a means of identifying changes in the global carbon cycle (e.g. Hilting et al., 2008) or as initial condition values for future model-based (and model–data-based) studies.

Figure 8Inferred global benthic mean δ13C, global mean surface ocean δ13C, and atmospheric δ13CO2 calculated using the mean bias of the best fit model to the data-indicated benthic ocean δ13C in each time slice. The data points from the δ13C benthic stack are from Westerhold et al. (2020). The shaded regions for predicted δ13CO2 show the range of δ13CO2 using the best ranges identified in Fig. 6.


Similarly, using this model–data benthic δ13C bias we can derive the global surface and benthic mean δ13C trend using the local bias of our model to the benthic δ13C data. This global mean deep-ocean δ13C is shown in Fig. 8 and is plotted along with the benthic δ13C stack from Westerhold et al. (2020). The data points from our benthic δ13C data set are also shown in Fig. 8 and tend to be higher (especially in younger time slices) than the benthic stack data; unlike the benthic stack data, our benthic data set is not limited to low latitudes and mid-latitudes (see Fig. 2). The North Atlantic data points generally have more positive δ13C values than other regions, and a large proportion of our data points are from this region. We use local model–data benthic δ13C differences to constrain the global benthic δ13C mean, but we still find this mean to be on the higher (heavier) end of the benthic δ13C stack from Westerhold et al. (2020). The modelled surface and deep-ocean mean δ13C falls by around 0.8 ‰ between 15 Ma and the Holocene, whereas atmospheric δ13CO2 falls by up to 1.5 ‰ at the same time.

3.4 Best fit model–data maps

In general, in the older time slices, cGENIE underestimates temperatures in the North Atlantic (Fig. 9), and this results in a higher M-score for the highest CO2 settings for the 15, 12.5, and 10 Ma slices (see Fig. 5). In these older time slices, the model high-latitude temperatures tend to be too low compared with data: this is an established characteristic of warm climates, where climate models tend to struggle to reproduce the flatter latitudinal temperature gradients seen in data (Goldner et al., 2014). Whether the CAS is open or closed makes little difference to the model–data fit, as our data points are not in regions that are very strongly affected by any resulting change in surface heat distribution. However, the temperature difference map shows that a closed CAS results in a colder modelled South Atlantic surface ocean than an open CAS as well as a slightly warmer Pacific surface ocean (note that although the FwF values are not the same between the cases, it is the CAS configuration that dominates the surface ocean temperature differences). These temperature differences of 1 or 2 C at high latitudes would have implications for sea-ice and ice-shelf growth in the cooler time slices.

Figure 9Modelled and data-indicated surface ocean temperatures for the best fit model settings, with data overlaid as shaded circles (temperature maps), and cross plots of data points (x axis) and model points (y axis) for both CAS cases. The surface temperature difference is shown on the right. The model–data fit statistic (M-score) is shown in Fig. 5.

Figure 10Benthic temperature for the best fit model settings, with data overlaid as shaded circles, and cross plots of data points (x axis) and model points (y axis), with temperature calculated with no salinity correction (semi-transparent) and with salinity correction (opaque). Noted above each map is the FwF setting (in Sv). The model–data fit statistic (M-score) is shown in Fig. 5.

Taking account of the deep-ocean salinity in the temperature calculation with δ18O data makes a significant difference to deep-ocean temperature patterns (Fig. 10). When flux corrections are higher (mainly from 10 Ma), the saltier North Atlantic waters are less buoyant and sink more readily. This higher salinity results in higher calculated temperatures (as δ18Osw is corrected for local salinity), with data points in shallower waters most strongly affected (see Appendix A). Accounting for water salinity increases some data points' temperature by more than 3 C (Fig. 10) in the North Atlantic. The inverse is true for data locations in the deep Pacific, where some temperatures are corrected lower. In both cases, this affects the model–data fit, likely resulting in higher CO2 and/or higher FwF combinations than would otherwise be the case if we had conducted the model–data comparison using only uncorrected δ18O-derived temperatures. The open-CAS case for the Holocene shows a higher M-score than the closed-CAS case for the model fit to benthic temperature (which we know should not be the case, as the CAS is actually closed in the Holocene) (Fig. 5). The Indian Ocean appears to be a location where modelled benthic temperature is too high compared with our δ18O-derived data, whereas the Pacific is generally too low (and the open-CAS case shows a better fit to data). The 2.5 Ma modelled benthic temperature is slightly too high, and this is probably due to the 400 ppm forcing, where 280 ppm would also have been a good fit (see Appendix D). The two oldest time slices show modelled benthic temperature with a much smaller range than the data indicates, which may be a result of assuming all benthic data points are located on the ocean floor in the cGENIE model grid (where ocean bathymetry is smoothed out compared with the real-world case). In general, a closed CAS results in a relatively warmer deep Atlantic and cooler Pacific than an open CAS for the best fit settings.

Figure 11Modelled δ13C of dissolved inorganic carbon (DIC) for the best fit model settings, with data overlaid as shaded circles, and cross plots of mean-adjusted data points (x axis) and mean-adjusted model points (y axis). Note that shading is with respect to an atmospheric δ13CO2 forcing of −6.5 ‰. This shading allows the comparison of the benthic ocean δ13C patterns between the time slices. The data-indicated δ13C value is shown to the left of each colour bar. The model–data fit statistic (M-score) is shown in Fig. 5.

The range of modelled δ13C values are in general agreement with the data-indicated range (Fig. 11). At 15 Ma, the modelled higher δ13C values in the high southern latitudes are not clearly shown in the data, which particularly affects the Indian Ocean. The 12.5 Ma slice shows a very low correlation for the distribution of δ13C. This may reflect the large climate and/or circulation changes occurring at the time (the mid-Miocene climate transition; Mackensen and Schmiedl, 2019) and the fidelity of the complied data set, especially with respect to age or, again, perhaps partly due to the smoothed model bathymetry. From 10 Ma to the Holocene, all of the closed-CAS cases show better fit to data than open-CAS cases (see also Fig. 5; the linear regression for each is shown in the cross plots (Fig. 11), and the closer this line is to diagonal, or the 1:1 line, the better the fit to data). Importantly, the model–data benthic δ13C correlation for 0 Ma is notably better for the closed-CAS cases than the open-CAS cases. This can also be seen in Fig. 5 where, regardless of the parameter combination, a better M-score is always obtained for a closed CAS (0 Ma) (with the exception of a single ensemble member for 0 Sv FwF and 280 ppm CO2).

Every ensemble member, regardless of the specific atmospheric CO2 assumption, was driven with a δ13CO2 value of −6.5 ‰. The model output in Fig. 11 is colour-coded according to this value (coloured shading) so that each model time slice can be compared with any other; the coloured shading indicates deep-ocean δ13C relative to a universal atmospheric δ13CO2 value (the data-indicated actual benthic δ13C value is shown to the left of each colour bar, and the scale for δ13CO2 at −6.5 ‰ is shown to the right of each colour bar). As a general rule, during the warmest earlier periods, the surface ocean southern high-latitude signal dominates the deep ocean from the south, and as time progresses and the climate cools, the surface North Atlantic becomes dominant over the deep Atlantic. This is accentuated for the closed-CAS cases.

3.5 Best fit ocean cross sections

The modelled Pacific and Atlantic ocean δ13C cross sections are presented in Fig. 12. (Note that for benthic δ13C data points, we always assume that they are at the modelled seafloor, so they are not shown here.) The cross sections are for the transects shown on the inset map in Fig. 12; the mean value of model output from a swath of 30 of latitude about the centre of the transect line is used to create the model ocean cross section data. Again, the coloured shading is with respect to a δ13CO2 of −6.5 ‰ in order to be able to compare time slices (as in Fig. 11). In the older slices, there is a larger offset between the whole-ocean δ13C and atmosphere than in the younger slices (the ocean is relatively heavier in δ13C in older time slices). The Pacific Ocean becomes progressively more negative (lighter) in δ13C through time, accentuated by a closed CAS. The mid-depth mid-latitude Pacific appears sensitive to the closing or opening of CAS, with modelled δ13C indicating some 0.2‰–0.5‰ differences between the open- or closed-CAS cases.

Figure 12Pacific and Atlantic ocean cross sections of δ13C, where the shading and contour labelling are with respect to an atmospheric δ13CO2 forcing of −6.5 ‰. This shading allows the comparison of the inner-ocean δ13C between the time slices, for example, 15 Ma shows values further from the atmospheric value than the Holocene. The data-indicated δ13C value is shown on the left of each colour bar. The inset map shows the transect locations used to make the cross sections.

Figure 13Modelled barotropic stream function (left) and global stream function (right) for each time slice for the best fit model settings. Contours for the barotropic stream function are at 20 Sv intervals, and contours for the global stream function are at 10 Sv intervals. Note the difference in colour scale between the barotropic and global stream functions.

The modelled barotropic and global stream functions for the best fit settings are shown in Fig. 13. Surface wind fields derived from HadCM3L exert some control on the barotropic streamfunction, as do thermohaline circulation and bathymetry. When the CAS is open, a circum-South America circulation is in operation, allowing the mixing of Pacific and Atlantic water masses. The Antarctic Circumpolar Current (ACC) is already well established by 15 Ma, and it strengthens towards 2.5 Ma with cooling and with the northern movement of South America (and associated widening of the Drake Passage). In the Holocene, the ACC is weaker compared with 2.5 Ma; this may be a similar effect to the weakening of the ACC in cold glacials relative to interglacials (Toyos et al., 2020) and may possibly be related to increased sea-ice extent in the lower-eq.CO2 Holocene (280 ppm) compared with the 2.5 Ma (400 ppm forcing), among other things (such as bathymetry and wind field).

The global streamfunction shows the increasing dominance of Northern Hemisphere deep water through the cooling from the Miocene to the Holocene, with a gradually deepening and strengthening AMOC. This is the case for both open-CAS and closed-CAS sets, although the open CAS-set requires a higher FwF in the North Atlantic to induce it than the closed-CAS set. The strongest AMOC is achieved in the coldest Holocene time slice, with eq.CO2 of 280 ppm, cool poles, and the largest sea-ice extent and thickness.

4 Discussion

Overall, we find that global surface and deep-ocean cooling since the mid-Miocene can be accounted for in cGENIE by the combined effects of changing bathymetry, eq.CO2 forcing, and salinity of the North Atlantic Ocean. The middle-Miocene eq.CO2 required in the model is relatively high (1120 ppm), and a small (0.1 Sv) salinity forcing of the North Atlantic can reproduce surface to deep-ocean heat distribution fairly well. However, latitudinal temperature gradients are not well reproduced, with high latitudes being too cool compared with data in the warmer time slices. The eq.CO2 forcing decreases steadily from the mid-Miocene to the present to satisfy surface temperature data, while the required salinity forcing of the North Atlantic increases to a maximum in the late Miocene/early Pliocene and then decreases again towards the present. Our modelled closed-CAS cases from 10 Ma to present all show better agreement with δ13C data distributions than the open-CAS cases. (Note that, due to data-indicated δ13C Atlantic Ocean latitudinal gradient trends in Fig. 2, we did not simulate a closed CAS at 12.5 or 15 Ma.)

In the next sections, we compare findings from this study against published evidence for the plausibility of the climate forcings and conditions that we have identified.

4.1 Heat distribution and eq.CO2

Previous climate modelling studies of the Miocene have generally found that proxy data for an atmospheric CO2 concentration of around 400pm are insufficient to explain the observed warmth (Greenop et al., 2014; Goldner et al., 2014; Bradshaw et al., 2012; Krapp and Jungclaus, 2011). Estimates of Miocene CO2 in more recent studies have trended higher compared with the earlier estimates (e.g. ∼220 ppm at 15 Ma; Pagani et al., 1999). Our required eq.CO2 forcing in the oldest, 15 Ma time slice is 1120 ppm. This contrasts with recent estimates of 450–550 ppm using a multi-method multi-taxon pCO2 reconstruction (Steinthorsdottir et al., 2021), 470 to 630 ppm from a boron CO2 proxy (Sosdian et al., 2018), 528 or 912 ppm from Neotropical fossil leaf assemblages (Londoño et al., 2018), or in the region of 1000 ppm in a recalculation of the alkenone εp CO2 proxy (Stoll et al., 2019). However, despite our eq.CO2 forcing being higher than these recent estimates, we still obtain a low M-score for the surface ocean temperature distribution.

A portion of the CO2 forcing model–data discrepancy may be explained by the contribution of other greenhouse gases to climate forcing during the Miocene. In the configuration of the cGENIE model used here, only CO2 is imposed as a greenhouse gas, so the imposed atmospheric eq.CO2 concentration should be understood as an “equivalent” climate forcing encompassing other greenhouse gases. Wetlands are a significant source of atmospheric methane in the present day as well as representing large terrestrial carbon stores. Extensive wetlands existed during the middle Miocene (Eronen and Rossner, 2007; Hoorn et al., 2010; Morley and Morley, 2013), with the generally warmer, wetter climate and lower-elevation topography being favourable to their initiation and persistence. These conditions would have been conducive to elevated methane production in the past (Dean et al., 2018). A large Antarctic ice sheet has been identified for the period soon after the Miocene climatic optimum (Hautvogel et al., 2012; Badger et al., 2013; Pierce et al., 2017), which would have caused a fall in sea level. Gas hydrates were found to have likely destabilised during the sea level lowering of the Miocene in the (present-day) Apennines in Italy (Argentino et al., 2019) and in the Black Sea basin (Kitchka et al., 2016), which may have had a transient impact on atmospheric methane levels and climate fluctuations. A larger contribution of methane to eq.CO2 forcing in the Miocene relative to the Holocene is certainly plausible (atmospheric chemistry and baseline methane levels also affect the lifetime of atmospheric methane; Schmidt and Shindell, 2003).

During the warmest intervals – 12.5 and 15 Ma – we find our largest differences in model–data surface temperature in the North Atlantic (Fig. 8), with modelled North Atlantic temperature being much lower than sea surface temperature (SST) indicators. The model overestimates the latitudinal temperature gradient in the warmer climates, with heat transport to the poles likely being too low. Increased ocean heat transport was found to reduce the latitudinal temperature gradient and the location of Hadley and Ferrel cells in a modelling study by Knietzsch et al. (2015). Changes in land surface cover and increased northern transient eddy atmospheric heat transport (linked to storm tracks at mid-latitudes) were found to increase heat transport to high latitudes in a modelled late Miocene in a fully coupled atmosphere–ocean general circulation model (AOGCM) by Micheels et al. (2011). These indicate the importance of atmospheric circulation, which is highly idealised in the 2D energy–moisture balance atmosphere configuration used here.

Equilibrium climate sensitivity is almost the same in all of our modelled climates, with a value of 2.9 (±0.05 at 1σ) and a range of 0.174. However, in more complex models, the interaction between vegetation and palaeogeography has been shown to give rise to a higher climate sensitivity in Miocene climate modelling (Bradshaw et al., 2015; Micheels et al., 2010), with warming occurring independent of CO2 increase due to vegetation and lower-elevation topography (Henrot et al., 2010) and their atmosphere feedbacks (Knorr et al., 2011). Indeed, chemistry–climate feedbacks linked to vegetation were also found to be as strong as or more important than CO2 forcing in the Pliocene (Unger and Yue, 2014). The significant changes occurring in global vegetation distribution since the middle Miocene (Pound et al., 2012) may then be critical to fully reproducing observed SST patterns. The absence of any representation of vegetation and associated feedbacks in the version of cGENIE that we employ may then help account for the low M-scores in model–data fit at a high mid-Miocene CO2 forcing of 1120 ppm as well as some of the discrepancies between model eq.CO2 and data-indicated CO2 in the older time slices (Fig. 7). Therefore, a higher climate sensitivity in the Miocene, vegetation–atmospheric feedbacks, and increased heat transport to the poles may be able to improve the model–data fit.

Miocene climatic optimum (MCO) proxy data suggest global surface temperatures that are 7 to 8 C warmer than modern temperatures (Steindthorsdottir et al., 2020). Our modelled mid-Miocene is around 6 C warmer than the Holocene, but this is in response to a relatively high estimate of eq.CO2 forcing (1120 ppm). The Mg/Ca thermometer estimates deep-ocean temperatures that are around 6 to 7 C warmer than present (Cramer et al., 2011, Lear et al., 2015), whereas we model temperatures that are around 3 to 4 C warmer (as a deep-ocean global mean). Indian Ocean benthic temperatures were found to be 9 C warmer than present at the Miocene climate optimum and 6 C warmer than present after the Miocene climate transition at 14.7 to 13 Ma by Modestou et al. (2020), which are far warmer than we model (Fig. 10). If high-latitude surface waters were warmer in the model (as data suggest), this would not only raise the global mean (and possibly allow for a lower eq.CO2 forcing) but would also drive a globally warmer deep ocean, as meridional overturning circulation transports high-latitude warm waters to the deep; thus, this may then explain some of the disagreement in surface and deep-ocean T constraints at 15 Ma (Fig. 6). However, it would probably not raise the deep Indian Ocean temperatures enough to be in line with the recent estimates from Modestou et al. (2020), suggesting a possible model circulation bias for the Indian Ocean in the mid-Miocene.

The feedbacks affecting climate sensitivity and the latitudinal temperature gradient are likely climate-sensitive, affected by topography, and are perhaps also directly temperature-sensitive. Topography, climate, and temperature approach those seen in the present over each successive time slice, suggesting that these climate feedbacks that we do not capture may become less important with each time slice as we approach the Holocene.

4.2 North Atlantic salinity and mountain uplift

The salinity flux correction that we apply to the North Atlantic can be considered to primarily represent an incompletely resolved runoff pattern across North America, which arises from (poorly resolved) atmospheric moisture transport and (unresolved) topography. In a warm, higher-CO2 climate like the Miocene, the atmosphere is able to hold more water vapour, resulting in a globally wetter climate. However, the North American plains saw the rise of grasslands that favour a drier climate (Janis et al., 2002) through the Miocene, and this is linked to the uplift of mountain ranges in the west that created a rain shadow on the central plains. Over the last 12 Myr in the Sierra Nevada, a rain shadow similar to present was identified by Mulch et al. (2008), indicating that mountain uplift in this area occurred prior to 12 Ma. Further north in the Cascades, high uplift and erosion rates were dated to the late Miocene (12 to 6 Ma) and the Coast Mountains and British Colombia uplift from 10 Ma (Reiners et al., 2002). The flux correction implicitly includes the impact of mountain building from around 12.5 to 6 Ma in addition to changing atmospheric moisture content (and after 6 Ma, a global cooling reduces the energy in the hydrological cycle). The modelled CAS being open or closed has an effect on the required FwF, with a closed CAS reducing the required salinity correction by isolating the Atlantic from the Pacific water masses and affecting North Atlantic water buoyancy. Our best fit FwF increases from 15 Ma, with the highest values seen at 4.5 Ma for an open CAS and 7.5 Ma for a closed CAS (Fig. 7), which is generally in agreement with evidence on the timing of North American mountain building and changing atmospheric moisture capacity, and with open-CAS cases requiring a higher FwF than the closed-CAS cases.

4.3 Ocean gateways and circulation

The Antarctic Circumpolar Current (ACC) strengthens through time up to 2.5 Ma in our best fit model time slices, with an ACC already clearly established at 15 Ma, when the Drake Passage was already open and fairly deep (Figs. 3, 13). Some evidence suggests that a volcanic arc may have inhibited ACC flow until after the mid-Miocene (Dalziel et al., 2013) or even later (Pérez et al., 2021), although other studies find an earlier onset for the ACC (e.g. Pfuhl et al., 2005). However, even if a volcanic arc was present, the upscaling to the cGENIE grid would smooth out this feature. A full consideration of changing bathymetry and gateways and how this impacts the ACC, Antarctic bottom water production and properties, and any knock-on effect on global climate would need to be the subject of a separate study.

Evidence shows a global cooling and an increase in ice growth at the mid-Miocene climate transition (MMCT, occurring from 14.5 to 12.5 Ma), which we do not identify in our 12.5 Ma time slice. The MMCT has been linked to the closing of the eastern gateway of the Tethys Sea, which is already closed at 15 Ma in our model set-up. An open eastern Tethys gateway allowed Tethyan Indian Saline Water (TISW) flow into the Indian Ocean that may have inhibited East Antarctic sea-ice growth until the gateway's closure at 14 Ma (MacKensen and Schmiedl, 2018). Surface North Atlantic temperature data lend some support to an enhanced AMOC or strengthened North Atlantic Current during the MMCT (and that heat distribution patterns are not consistent with solely CO2-driven cooling in this location) (Super et al., 2020). Furthermore, for the 12.5 Ma time slice, our δ13C data set provides only a very weak constraint on the forcings. It is certainly possible that with a better (less noisy) δ13C data set, the 12.5 Ma best fit could target a lower-CO2 value (more in line with CO2 proxy data, such as the data shown in Fig. 7) and a higher-FwF value, with a more saline North Atlantic very possibly linked to a salty-water contribution from the Tethys Sea that is now closed off from the Indian Ocean. However, this lower CO2 value would then result in lower modelled surface ocean temperatures at low latitudes, which are already slightly too low at 1120 ppm in cGENIE according to data (Fig. 9).

We identify a possible early CAS closure/restriction from 10 Ma according to the combined temperature and δ13C constraint. A “washhouse” (warm and wet) climate was identified in Europe in certain periods of the Miocene, from 10.2 to 9.8 and from 9.0 to 8.5 Ma, which appeared to be linked to temperatures in the deep North Atlantic Ocean (Böhme et al., 2008). This was subsequently linked to a possible temporary restriction of the CAS and greater northward heat transport. A middle-Miocene (even if temporary; Jaramillo et al., 2017) closure of the CAS was identified by Montes et al. (2015), further supporting our findings. The CAS is thought to have been definitely closed by 2.5 Ma (O'Dea et al., 2016), and evidence shows that some flow between the Pacific and Atlantic was still occurring up until that point (O'Dea et al., 2017, Jaramillo et al., 2017). A closed CAS from 10 Ma shows better agreement with δ13C data in our study, but we do not consider a restricted CAS – only open or closed (Fig. 3 shows CAS depths of at least 1000 m for our open-CAS cases); periodic CAS closing or restriction from 10 up until 2.5 Ma is consistent with our findings. A shallow CAS still allowed the formation of deep water in the Miocene North Atlantic in a GCM model study by Nisancioglu et al. (2003), and even an open CAS in this study (at >1000 m; Fig. 3) does not impede the establishment of the AMOC if the surface North Atlantic is sufficiently saline. Our modelled global overturning circulation shows a gradual increase in the dominance of Northern Hemisphere deep water from the mid-Miocene as well as the development of the AMOC (Fig. 13). The onset of the AMOC was thought to be a key factor in the expansion of ice at the poles. However, Bell et al. (2015) found that an early Pliocene (4.7 to 4.2 Ma) shoaling of CAS had no profound impact on climate evolution, as North Atlantic deep-water formation was found to already be vigorous by 4.7 Ma. Vigorous North Atlantic deep-water formation appears to have probably started by 4.5 Ma in our study (Fig. 13), when CO2 drops to 400 ppm. A definitive closure of the CAS dates to the late Pliocene (∼3 Ma; O'Dea et al., 2017). This final Pliocene closure of CAS had been thought to be linked to conditions conducive to Northern Hemisphere glaciation, although Lunt et al. (2008) found that an open or closed CAS made little difference to ice-sheet size in their model study of the Pliocene. The closed CAS has a lesser effect on sea ice compared with CO2 in our study; sea-ice growth intensifies when CO2 drops to 280 ppm in the Holocene time slice (Fig. 7) for both open or closed CAS, indicating that it is CO2 that would probably drive Northern Hemisphere glaciation (in cGENIE), rather than a CAS closure.

4.4 Benthic temperatures and the AMOC

The pattern of temperature in the deep ocean arises as a combination of the pattern of ocean surface temperature (and surface climate in general) and the large-scale circulation of the ocean. In our model, increasing the net salinity transport into the North Atlantic region induces and strengthens a meridional overturning circulation that, in turn, drives a warming of the deep ocean, independently of atmospheric CO2 forcing (Fig. 4). Local water salinity has a strong control on δ18O. For example, an increase in measured benthic δ18O in benthic foraminifera from the North Atlantic during the late Miocene could be interpreted as evidence of a cooling or, equally, it could be attributable to the increased salinity of the sea water, where salinity (rather than temperature) dominates the δ18O signal. With the strengthening of Atlantic overturning circulation during the Miocene, the increased salinity of deep North Atlantic waters exerts a relatively stronger control on temperature calculated from shell δ18O (than for a weaker AMOC). When accounting for local salinity in the temperature calculation, increases of up to 3C in some locations are seen when compared with the temperature that is uncorrected for salinity (Fig. 10). Without this correction to the δ18O temperature calculation, the North Atlantic would appear to be cooler (and the Pacific warmer), and this would have a large impact on what the best fit settings would be for our study (see Appendix A).

4.5 Caveats for δ13C as an ocean circulation tracer

Changes in deep-ocean circulation patterns and specifically the ageing of water masses and progressive accumulation of isotopically light respired carbon should, in theory, be identifiable in global patterns of δ13C data from benthic foraminifera. However, the processes that control δ13C are complex, involving both circulation and ocean carbon pumps, which are also not independent; hence, large-scale changes in circulation that affect nutrient return to the surface will also modulate the strength of the biological pump in the ocean. Changing ocean interior temperature patterns may also influence where carbon is respired via a temperature control on the rate of carbon respiration (John et al., 2014), further affecting δ13C.

In this study, we have applied two of the circulation-linked controllers of δ13C by applying palaeobathymetry and altering thermohaline circulation by forcing eq.CO2 and North Atlantic salinity (FwF). There are uncertainties in the applied palaeobathymetries as well as in the re-gridding method that up-scales to the cGENIE grid (as demonstrated in the re-gridding resulting in an open and fairly deep CAS up to and including the Holocene, which is an issue equally applicable to GCMs given the relatively small width of the Isthmus of Panama). Re-gridding also affects the depth at which benthic data are assumed to be located, as it smoothes out peaks and troughs in the ocean floor; thus, some data locations may have been forced deeper or shallower in the ocean than they would have actually been.

Orbital variations will also have an effect on recorded δ13C via changes in climate and ocean circulation. Some of the benthic data were high-resolution data sets where it was possible to select a mid-climate value (i.e. between the highest and lowest values in the orbit oscillation). However many of the data were not of this type; therefore, the data set may not represent the same point in time well (in an orbital cycle). This is also generally the case for uncertainties in age models of the benthic data, and it would particularly affect periods in which climate is likely to have been changing over the 2 Myr window (that we used to select data) – such as the 12.5 Ma slice. This may explain why δ13C distributions show very little model–data agreement for any eq.CO2–FwF combination at 12.5 Ma. These combined effects create noise in the δ13C data set (on top of uncertainties in the data itself), with some time slices showing higher M-scores (in general) than others (Fig. 5).

Despite these limitations, the δ13C model–data comparison combined with temperature data provided higher M-scores for more plausible forcings at each time slice and also showed agreement with published work on Miocene to Pleistocene climate and circulation patterns. We identify a likely restricted CAS at 10 and 4.5 Ma as well as a closed CAS at 2.5 Ma and in the Holocene. At 7.5 Ma, we identify a possible intermediate state (where data may be a mix of both open and closed cases), with low M-scores from δ13C in both open-CAS and closed-CAS “best fit” cases. At 12.5 Ma, our combined δ13C–temperature approach showed very poor agreement with other indicators of climate changes, for example, in reconciling surface ocean temperature reconstructions with evidence of extensive cooling and ice growth at this time. Importantly, benthic δ13C patterns are able to distinguish between open and closed CAS for the present-day in the model, which we found is not the case for benthic temperature proxies.

5 Conclusions and summary

In this study, we used proxy data estimates for both surface and benthic temperature (δ18O) as well as benthic δ13C in order to constrain the evolution of atmospheric eq.CO2 and large-scale ocean circulation in the “cGENIE.muffin” Earth system model and, hence, identify plausible climatic states for seven respective time slices spanning the mid-Miocene to Holocene cooling. Constrained by changes in the absolute magnitude and pattern of benthic δ13C, we also diagnosed a plausible history of atmospheric δ13CO2 over this time interval for use as boundary conditions in future modelling studies or for use as a data target for assimilation in geochemical box models.

In the cGENIE model, we diagnose a progressive reduction in atmospheric greenhouse gas forcing since the mid-Miocene, driving global cooling. Simultaneously, from the middle Miocene, we diagnose a gradual strengthening of overturning circulation in the Atlantic that transports heat to the deep Atlantic Ocean over the cooling period and leads to a stronger cooling in surface waters (at ∼6C) than in deep waters (at ∼3C). This onset and strength of the AMOC in cGENIE is controlled by the combined effects of the progressive restriction of the Central American Seaway together with our two free parameters: a salinity adjustment (representing mountain building in North America and an increasing Atlantic–Pacific salinity gradient) and declining atmospheric CO2. Declining CO2 drives cooling which helps to promote the sinking of salty waters in the North Atlantic. The net result in the model is a strong and deep AMOC in the Holocene when the CAS is closed and atmospheric CO2 is low.

Appendix A: Correcting benthic temperature calculated from δ18O for local salinity

The Cramer et al. (2011) δ18Osw estimate for sea water is a global mean value. Therefore, applying the Marchitto et al. (2014) palaeotemperature calculation using one global mean δ18Osw results in an uncertainty in temperature depending on location (Fig. A1) due to local differences in δ18Osw values (Fig. A2).

Figure A1Example palaeotemperature calculation of Marchitto et al. (2014) for a δ18OCibicidoides, Planulina of 3.3 ‰ with a variable δ18Osw. The calculated temperature using the global mean δ18Osw from Cramer et al. (2011) for the Holocene is marked, and selected ocean region mean δ18Osw values have been labelled on the figure (values evident in Fig. A2 for ocean basin means). When applying a local correction to δ18Osw, the North Atlantic calculated temperature is ∼0.5C warmer, and the North Pacific calculated temperature is ∼0.7C cooler, than the when applying the global mean δ18Osw.


The seawater δ18Osw is determined by global ice volume (which we get from Cramer et al., 2011), local temperature, and local salinity (Rohling, 2013). In the present day, with an active AMOC, the North Atlantic benthos have a more positive δ18Osw than, for example, the North Pacific. This is due to both temperature and salinity, with the salty North Atlantic waters transported to the deep by the AMOC. In our model ensembles, the benthic salinity is affected by both CO2 and flux correction, as these affect ocean circulation. Therefore, for each simulation, we apply a δ18Osw-driven correction to the palaeotemperatures due to their locally modelled salinity. To do this, we use present-day deep-water (2500 m) δ18Osw (LeGrande and Schmidt, 2006) and salinity (WOA 2013; Zweng et al., 2013) to create a general linear model (Eq. A1):

(A1) δ 18 O sw = 0.8 S - 27.7 ,

where S is salinity.

The North Atlantic is the region that is most affected by changes in salinity; therefore, the greatest temperature offset will be in this location. We adjusted the salinity model to get best fit values for the North Atlantic as well as a good fit for the difference in δ18Osw between the North Atlantic and the Pacific. All ocean locations with δ18Osw between −0.3 ‰ and 0.3 ‰ are shown in a cross plot of data and salinity-model-derived δ18Osw in Fig. A3. The grouping of points clearly offset from the 1:1 line are the high southern latitudes (also clearly visible in Fig. A2). The accuracy of the salinity model in finding δ18Osw is ±0.03 ‰ at 1 standard deviation and ±0.06 ‰ at the 95 % confidence level excluding latitudes higher than 70 (where we have no foraminifera shell δ18O data points in any case).

The palaeotemperature equation that we apply is the linear model from Marchitto et al. (2014, Eq. 8 therein), using the global ice volume estimate from Cramer et al. (2011). The δ18Osw model from salinity is also linear, so we apply a simple linear correction to the calculated temperature (Table S2). To obtain the δ18Osw offset, we first find the modelled global mean salinity (not including latitudes higher than 70) and then subtract it from the modelled benthic salinity, giving a ΔS field (an offset from the mean). We apply Eq. (A2) to find ΔT, where 0.8 is the gradient of the linear δ18Osw model (Eq. A1), and 0.224 is the gradient of the linear palaeotemperature model (Marchitto et al., 2014), and we use this to correct T for salinity. The uncertainty for the temperature correction is ±0.13C at 1 standard deviation and 0.23 C at the 95 % confidence level (based on the δ18O model uncertainties). The range of temperature correction corresponding to the benthic δ18Osw spread (of ∼0.4 ‰) in the present day is 1.8 C.

(A2) Δ T = 0.8 Δ S 0.224

Figure A2δ18Osw data from LeGrande and Schmidt (2006) at 2500 m depth (at depths below 2500 m, the data set shows little change in δ18Osw compared with the values shown) (a) and salinity-derived δ18Osw (b), where salinity data from WOA 2013 (Zweng et al., 2013) are used in the model in Eq. (A1) (all in VSMOW).

Figure A3Cross plot of data δ18Osw and salinity-model-derived δ18Osw as shown in Fig. A2. The offset region (where data shows low δ18Osw compared with the derived salinity model) is the Southern Ocean, where the polar front more strongly affects ocean water δ18Osw than salinity/thermohaline circulation.


Appendix B: Model boundary conditions

In this paper, as per previous (deeper time) palaeo-applications of the cGENIE.muffin model (e.g. Ridgwell and Schmidt, 2010), we derive the required boundary conditions from a representative fully coupled GCM experiment using the “muffingen” open-source software version v0.9.20 (which is assigned the following DOI:, rather than using observations. This software takes a specific GCM experiment as input and carries out the following:

  1. It creates output (here, 36×36 with 16 levels in the ocean) grids.

  2. Using (1), it derives a land–sea mask and ocean bathymetry on the output grid. The land–sea mask is filtered to prevent the occurrence of isolated inland seas and single-cell-width coastal embayment, while the ocean bathymetry is filtered to avoid single-cell “holes” occurring in the ocean floor.

  3. The model then generates drainage basins determining where precipitation is directed towards the ocean (and, hence, runoff direction). The specific scheme used here is known as a “roofing scheme” and operates to create a watershed approximately equidistant from the coast.

  4. It next derives island and ocean path files required by the ocean circulation model.

  5. The model then re-grids GCM wind stress and (10 m) wind velocity to the output grid (for wind stress, this means re-gridding to both u- and v-edge grids). Wind speed is calculated from the mean annual wind velocity components.

  6. Finally, it re-grids the GCM planetary albedo and converts to a zonally averaged profile.

The GCM simulations underlying our cGENIE model configurations were carried out using HadCM3LM2.1aE, which is described in detail in Valdes et al. (2017). The models are constrained with palaeogeographies and a solar constant appropriate for each geological stage in the Miocene as well as a CO2 mixing ratio of 400 ppmv. The experimental design is described in detail in Farnsworth et al. (2019). The specific GCM experiments that we used from Farnsworth et al. (2019) are those for the Holocene (0 Ma), Piacenzian (2.58–3.6 Ma), Zanclean (3.6–5.333 Ma), Messinian (5.333–7.246 Ma), Tortonian (7.246–11.63 Ma), Serravallian (11.63–13.82 Ma), and Langhian (13.82–15.97 Ma); we simplified and adopted these experiments to approximately evenly spaced slices at 0.0, 2.5, 4.5, 7.5, 10.0, 12.5, and 15.0 Ma (Table 1 main text).

For all time-slice configurations we generally follow the default re-gridding algorithm of the muffingen software in order to avoid imposing prior assumptions regarding the importance of specific ocean features, meaning that a Mediterranean Sea is not present for most of the reconstructed continental configurations in cGENIE. Only in the 12.5 and 15.0 Ma time slices is the remnant Tethys Ocean sufficiently expansive to re-grid as an ocean basin at the selected 36×36 (16 levels) model resolution. We do, however, make the following manual interventions in the generation of the land–sea mask (but not in ocean bathymetry):

  • for 0 Ma (late Holocene), the Panama Isthmus as well as the tip of South America is made continuous;

  • for 10.0, 12.5, and 15.0 Ma, the Arctic is opened up (turning land cells to ocean) in order to approximately preserve the global land fraction of 0.29 in the underlying GCM across all seven reconstructions;

  • in a second set of simulations, the Panama Isthmus is made continuous for the 10 to 2.5 Ma time slices and the late-Holocene Panama Isthmus is left open.

The resulting cGENIE bathymetry for each time slice is shown in Fig. 3. The Drake Passage is already open at 15 Ma, but it widens by 12.5 Ma and gradually deepens. The Atlantic Ocean widens and deepens throughout, and the Pacific generally deepens. The Australian land mass moves north, and the Indonesian seaway gradually reduces. Africa moves north, reducing the Tethys Sea area, which is already closed off from the Indian Ocean at 15 Ma. The Tethys disappears by 10 Ma, and the Mediterranean Sea is not included in the simulations. Greenland is isolated from North America until 7.5 Ma, and the Bering Strait is closed until 7.5 Ma. Only at the Holocene is the Bering Strait open and is Greenland once again isolated from North America, allowing mixing between the Arctic Ocean, the North Pacific, and North Atlantic at the same time.

We chose and calculate a zonally average (GCM-derived) planetary albedo profile (rather than a 2D re-gridded one) in order to retain closer back-compatibility with the original GENIE configuration in which an idealised zonal profile is applied (e.g. Edwards and Marsh, 2005). Different GCMs average and save wind speed differently (or only as velocity vectors), meaning that the final re-gridded wind speed product can differ substantially between GCMs and with modern observations. Because of this, we re-scale air–sea gas exchange in the late-Holocene configuration in order to give a mean global and annual average modern air–sea coefficient value for CO2 of approximately 0.058 molm-2yr-1µatm-1. This same air–sea gas exchange scaling is then applied to all older time slices. As compared to Cao et al. (2009), we also forego the high southern latitude zone of reduced atmospheric diffusivity, previously used for present-day model configurations (described further in Marsh et al., 2011). Initial mean salinity was reduced by 1 to 33.9 PSU in all time-slice configurations for simplicity and consistency (although we recognise that mean ocean salinity should realistically progressively decrease from a modern value of 34.9 back in time with progressively decreasing global land ice volume). The orbital configuration was kept at its modern settings throughout all time slices (0.0167 eccentricity, 0.397789 for sine of obliquity, and 102.92 for the longitude of perihelion; longitude is in degrees). However, we did vary the following model parameter values and initial conditions as a function of (geological) time (time-slice age):

  • The solar constant is assumed to change with time and to follow Gough (1981) (see also Feulner, 2012), resulting in a small increase between 15.0 and 0 Ma, from 1366.09 W m−2 (a reduction of 0.14 % compared to modern) at 15.0 Ma to 1368.0 W m−2 by the late Holocene.

  • The mean Mg/Ca ratio of the ocean is also assumed to change with time, following Tyrrell and Zeebe (2004). The corresponding range is then from 13.15 mmolkg-1Ca2+ and 41.21 mmolkg-1Mg2+ at 15.0 Ma to 10.28 mmolkg-1Ca2+ and 52.81 mmolkg-1Mg2+ in the present-day ocean. (Changing ocean Mg/Ca with time influences the calculation of carbonate saturation as well as dissociated constants.)

All boundary configurations and relevant parameter settings (as well as the cGENIE.muffin model code itself) are open-source. Refer to our “Code availability” statement for information regarding obtaining and running any or all of these model configurations.

Figure B1Comparison of PLATES-derived rotation (applied to the full data set) to the HadCM3 rotation (used as a basis for the cGENIE model). The mean difference between PLATES and HadCM3 palaeo-locations for latitude is 0.216 (standard deviation of 1.14), whereas for longitude, the mean difference is −0.498 (standard deviation of 1.80).


Figure B2Comparison of cGENIE to HadCM3 near-surface and deep-ocean currents and to data-indicated near-surface currents from OSCAR (Ocean Surface Current Analysis in Real-time) (ESR, 2009). Note the scale difference for the surface and deep-ocean colour bars and vectors.

Figure B3Location of salinity addition and salinity removal common to all configurations for the salinity flux adjustment controller “FwF”.


Appendix C: Sustained climate oscillations in the ensemble

The settings for which sustained climate oscillations were identified are shown in Fig. C1. These oscillations were present in a subset of ensemble members in all time-slice ensembles except the Holocene, and they were more prevalent for the open-CAS configurations. They ranged in amplitude from 0.3 to 2.5 C in benthic temperature with periods of 1.7 to 4.5 kyr.

Figure C1Ensemble settings that displayed sustained oscillations in benthic ocean temperature (top). Model output for atmospheric temperature, sea ice, and benthic temperature for the longest- and shortest-period oscillations (bottom panels).


Appendix D: Details on selecting the “best fit” settings

D1 Holocene, 0 Ma

For the Holocene slice, we have a priori knowledge that CO2 during this warm interglacial period was around 280 ppm (Indermühle et al., 1999) – a value which in cGENIE agrees with the surface ocean temperature data. This gives us some confidence in the cGENIE Earth system model and the methodology; however, although we use a new and different GCM-derived (Farnsworth et al., 2019) modern continental configuration here, cGENIE has already been calibrated against present-day observed ocean temperatures at an atmospheric CO2 value of ca. 280 ppm (e.g. Price et al., 2009; Ridgwell et al., 2007), so an acceptable fit to SSTs for an eq.CO2 value of 280 ppm is not necessarily unexpected. The global benthic temperature estimate from Cramer et al. (2011) represents a mean climate for both glacial and interglacials, which is around 2 C cooler than the warmer Holocene. For a global benthic temperature that is 2 C warmer than that from Cramer et al. (2011), a relatively low flux correction of 0.1 to 0.2 Sv is required to fit the data (in comparison, the present-day calibrated value in cGENIE is 0.32 Sv). However, the benthic δ13C data constraints (Fig. 6) tend towards a higher M-score for a higher flux correction, so we chose a best fit value of FwF 0.2 Sv for the closed-CAS case. In an open-CAS case, the δ13C constraint suggests a higher FwF, so the open-CAS best fit is set at 0.3 Sv.

D2 Piacenzian, 2.5 Ma

For the 2.5 Ma slice, surface and benthic temperatures both suggest a higher (and larger range for a fit to) eq.CO2 than for the Holocene, so 400 ppm is selected as the best fit eq.CO2 value (although 280 ppm would also work well). Similar to the Holocene, a higher flux correction (FwF) at 2.5 Ma tends to show a better model–data benthic δ13C agreement for the open-CAS case, although this constraint is less strong than the benthic temperature constraint (maximum M-score for benthic temperature is ∼0.4, whereas it is ∼0.2 for δ13C at 0.5 Sv). Hence, for the open CAS, we select 0.3 Sv as the FwF value, which is at the higher end of the benthic temperature (δ18O-derived) data set but is in agreement with global mean benthic temperature from Cramer et al. (2011). The closed CAS shows a generally cooler benthic ocean temperature and a lower FwF, so an eq.CO2 value of 400 ppm is combined with a FwF of 0.1 Sv for the closed CAS.

D3 Zanclean, 4.5 Ma

At 4.5 Ma, the surface temperature data support an even higher eq.CO2 value in cGENIE, but here the δ13C constraint is stronger (compared with that at 2.5 Ma) for both CAS cases (and instead supports a lower CO2 and higher FwF) as the M-score increases in this direction of parameter space (Fig. 6). Hence, we place the best fit in this direction for the open CAS, with an eq.CO2 value of 400 ppm and an FwF value of 0.5 Sv. For the closed CAS, the benthic δ13C constraint is stronger still, and, together with a benthic temperature constraint indicating a lower FwF than for the open CAS, an eq.CO2 of 400 ppm is combined with a FwF of 0.2 Sv for the closed CAS.

D4 Messinian, 7.5 Ma

At 7.5 Ma, surface temperature data require a significantly increased eq.CO2 compared with 4.5 Ma, with 800 ppm clearly being the best fit CO2 value. A similar tendency (to 4.5 Ma) toward lower CO2 but higher FwF is also apparent in the δ13C constraint for this time slice, although at 800 ppm, this is less strong (with a low M-score for the δ13C constraint at 800 ppm and high FwF). Again, with a closed CAS, the benthic temperature constraints indicate lower FwF than for the open CAS. The flux correction is set at 0.4 Sv, which is near the centre of the benthic temperature maximum M-score (shown as a dashed blue contour in Fig. 6) for the open-CAS case. For the closed-CAS case, 800 ppm is combined with an FwF of 0.3 Sv, which is towards the higher end of the benthic temperature constraint, but in the direction of increasing δ13C M-score.

D5 Tortonian, 10 Ma

For 10 Ma and older time slices, the model–data M-score for surface temperature declines, although with higher scores for higher CO2 in all three cases. At 10 Ma, for open CAS, the benthic (δ18O) temperature data set suggests a CO2 of 800 ppm and an FwF of 0.3 to 0.4 Sv. The δ13C fit surface tends towards higher M-score values for higher flux, so we set FwF at 0.4 Sv for open CAS. For closed CAS, both benthic constraints suggest a lower FwF, so 800 ppm combined with an FwF of 0.3 Sv is selected for closed CAS.

D6 Serravallian, 12.5 Ma

The 12.5 Ma slice seems to show a transition state for the trends in δ13C, with an overall low M-score for all combinations of CO2 and FwF. As δ13C provides a weak constraint, the CO2 and FwF values are selected as 1120 ppm and 0.2 Sv respectively as a compromise between the high-CO2 requirement for surface temperature and the lower-CO2–higher-FwF for the benthic temperature.

D7 Langhian, 15 Ma

At 15 Ma, the ensemble surface of the M-score for δ13C is inversed compared with other time slices, showing a higher score for (generally) higher eq.CO2 combined with a lower FwF. Although we have fewer surface temperature data points, they suggest a high eq.CO2, which is somewhat in disagreement with the benthic temperature data set that favours a mid-range eq.CO2 and a higher FwF. As the δ13C constraint tends towards a lower FwF (compared with time slices younger than 10 Ma), we select an eq.CO2 value of 1120 ppm and a low flux correction of 0.1 Sv.

Code availability

The code for the version of the “muffin” cGENIE Earth system model used in this paper is tagged as v0.9.22 and is assigned the following DOI: (Ridgwell et al., 2021a). Configuration files for the specific experiments presented in the paper can be found in the following directory: genie-userconfigs/MS/crichtonetal.CP.2021. Details of the experiments, as well as the command line needed to run each one, are given in the readme.txt file in that directory. All other configuration files and boundary conditions are provided as part of the code release.

A manual detailing code installation, basic model configuration, tutorials covering various aspects of model configuration and experimental design, and the output and processing of results is assigned the following DOI: (Ridgwell et al., 2021b).

Data availability

All data sets used in this study and derived from previously published data are available as tables in the Supplement. The code required to run the simulations is freely available (see “Code availability” above).


The Supplement to this paper provides data tables of surface ocean temperature data, benthic ocean δ18O data and calculated non-salinity-corrected benthic temperature, and benthic ocean δ13C data. Time-slice maps showing the location of the benthic data data points and the ±1 Myr time series for each are also given. The supplement related to this article is available online at:

Author contributions

KAC collated the temperature, δ18O, and δ13C data; devised the salinity correction for temperature; and carried out the model–data comparison and analyses. AR produced the cGENIE time-slice configurations, ran the ensembles, and extracted mean-climate data. DJL and AF provided the HadCM3L configurations. PNP was principal investigator of the project under which this research was carried out. All authors wrote the paper.

Competing interests

The authors declare that they have no conflict of interest.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Katherine A. Crichton (NERC; grant no NE/N001621/1 to Paul N. Pearson) and Daniel J. Lunt and Andy Ridgwell (grant no. NE/K014757/1) were supported by the Natural Environment Research Council. Andy Ridgwell acknowledges further support from the National Science Foundation (grant nos. 1702913 and 1736771).

Financial support

This research has been supported by the Natural Environment Research Council (grant nos. NE/N001621/1 and NE/K014757/1) and the National Science Foundation (grant nos. 1702913 and 1736771).

Review statement

This paper was edited by Yannick Donnadieu and reviewed by Anta-Clarisse Sarr and two anonymous referees.


Abelson, M. and Erez, J.: The onset of modern-like Atlantic meridional overturning circulation at the Eocene-Oligocene transition: Evidence, causes, and possible implications for global cooling, Geochem. Geophy. Geosy., 18, 2177–2199, 2017. 

Argentino, C., Conti, S., Fioroni, C., and Fontana, D.: Evidences for Paleo-Gas Hydrate Occurrence: What We Can Infer for the Miocene of the Northern Apennines (Italy), Geosciences, 9, 134,, 2019. 

Badger, M. P., Lear, C. H., Pancost, R. D., Foster, G. L., Bailey, T. R., Leng, M. J., and Abels, H. A.: CO2 drawdown following the middle Miocene expansion of the Antarctic Ice Sheet, Paleoceanography, 28, 42–53,, 2013. 

Bell, D. B., Jung, S. J., Kroon, D., Hodell, D. A., Lourens, L. J., and Raymo, M. E.: Atlantic deep-water response to the early Pliocene shoaling of the Central American Seaway, Sci. Rep., UK, 5, 12252,, 2015. 

Bialik, O. M., Frank, M., Betzler, C., Zammit, R., and Waldmann, N. D.: Two-step closure of the Miocene Indian Ocean Gateway to the Mediterranean, Sci. Rep., 9, 8842,, 2019. 

Böhme, M., Ilg, A., and Winklhofer, M.: Late Miocene “washhouse” climate in Europe, Earth Planet Sc. Lett., 275, 393–401,, 2008. 

Bradshaw, C. D., Lunt, D. J., Flecker, R., Salzmann, U., Pound, M. J., Haywood, A. M., and Eronen, J. T.: The relative roles of CO2 and palaeogeography in determining late Miocene climate: results from a terrestrial model–data comparison, Clim. Past, 8, 1257–1285,, 2012. 

Bradshaw, C. D., Lunt, D. J., Flecker, R., and Davies-Barnard, T.: Disentangling the roles of late Miocene palaeogeography and vegetation–Implications for climate sensitivity, Palaeogeogr. Palaeocl., 417, 17–34,, 2015. 

Burls, N. J., Bradshaw, C. D., De Boer, A. M., Herold, N., Huber, M., Pound, M., Donnadieu, Y., Farnsworth, A., Frigola, A., Gasson, E., and von der Heydt, A. S.: Simulating Miocene warmth: insights from an opportunistic Multi-Model ensemble (MioMIP1), Palaeogeogr. Palaeocl., 36, e2020PA004054,, 2021. 

Butzin, M., Lohmann, G., and Bickert, T.: Miocene ocean circulation inferred from marine carbon cycle modeling combined with benthic isotope records, Paleoceanography, 26, PA1203,, 2011. 

Cao, L., Eby, M., Ridgwell, A., Caldeira, K., Archer, D., Ishida, A., Joos, F., Matsumoto, K., Mikolajewicz, U., Mouchet, A., Orr, J. C., Plattner, G.-K., Schlitzer, R., Tokos, K., Totterdell, I., Tschumi, T., Yamanaka, Y., and Yool, A.: The role of ocean transport in the uptake of anthropogenic CO2, Biogeosciences, 6, 375–390,, 2009. 

Cramer, B. S., Miller, K. G., Barrett, P. J., and Wright, J. D.: Late Cretaceous–Neogene trends in deep ocean temperature and continental ice volume: Reconciling records of benthic foraminiferal geochemistry (δ18O and Mg/Ca) with sea level history, J. Geophys. Res-Oceans, 116, C12023,, 2011. 

Crichton, K. A., Wilson, J. D., Ridgwell, A., and Pearson, P. N.: Calibration of temperature-dependent ocean microbial processes in the cGENIE.muffin (v0.9.13) Earth system model, Geosci. Model Dev., 14, 125–149,, 2021. 

Dalziel, I. W. D., Lawver, L. A., Pearce, J. A., Barker, P. F., Hastie, A. R., Barfod, D. N., Schenke, H. W., and Davis, M. B.: A potential barrier to deep Antarctic circumpolar flow until the late Miocene?, Geology, 41, 947–950, 2013. 

Dean, J. F., Middelburg, J. J., Röckmann, T., Aerts, R., Blauw, L. G., Egger, M., Jetten, M. S., de Jong, A. E., Meisel, O. H., Rasigraf, O., and Slomp, C. P.: Methane feedbacks to the global climate system in a warmer world, Rev. Geophys., 56, 207–250,, 2018. 

Edwards, N. R. and Marsh, R.: Uncertainties due to transport-parameter sensitivity in an efficient 3-D ocean-climate model, Clim. Dynam., 24, 415–433,, 2005. 

Edwards, N. R., Willmott, A. J., and Killworth, P. D.: On the role of topography and wind stress on the stability of the thermohaline circulation, J. Phys. Oceanogr., 28, 756–778,<0756:OTROTA>2.0.CO;2.1998. 

Eronen, J. T. and Rössner, G. E.: Wetland paradise lost: Miocene community dynamics in large herbivorous mammals from the German Molasse Basin, Evol. Ecol. Res., 9, 471–494, 2007. 

ESR: OSCAR third deg., Ver. 1, CA, USA, Dataset, PO.DAAC, (last access: 27 April 2021), 2009. 

Farnsworth, A., Lunt, D. J., Robinson, S. A., Valdes, P. J., Roberts, W. H., Clift, P. D., Markwick, P., Su, T., Wrobel, N., Bragg, F., and Kelland, S. J.: Past East Asian monsoon evolution controlled by paleogeography, not CO2, Science Advances, 5, eaax1697,, 2019. 

Feulner, G.: The faint young Sun problem, Rev. Geophys., 50, RG2006,, 2012. 

Filippelli, G. M.: Intensification of the Asian monsoon and a chemical weathering event in the late Miocene–early Pliocene: implications for late Neogene climate change, Geology, 25, 27–30,<0027:IOTAMA>2.3.CO;2, 1997. 

Goldner, A., Herold, N., and Huber, M.: The challenge of simulating the warmth of the mid-Miocene climatic optimum in CESM1, Clim. Past, 10, 523–536,, 2014. 

Greenop, R., Foster, G. L., Wilson, P. A., and Lear, C. H.: Middle Miocene climate instability associated with high-amplitude CO2 variability, Paleoceanography, 29, 845–853,, 2014. 

Gough, D. O.: Solar interior structure and luminosity variations, Sol. Phys., 74, 21–34,, 1981. 

Hamon, N., Sepulchre, P., Lefebvre, V., and Ramstein, G.: The role of eastern Tethys seaway closure in the Middle Miocene Climatic Transition (ca. 14 Ma), Clim. Past, 9, 2687–2702,, 2013. 

Hauptvogel, D. W. and Passchier, S.: Early–Middle Miocene (17–14 Ma) Antarctic ice dynamics reconstructed from the heavy mineral provenance in the AND-2 A drill core, Ross Sea, Antarctica, Global Planet. Change, 82, 38–50,, 2012. 

Henrot, A.-J., François, L., Favre, E., Butzin, M., Ouberdous, M., and Munhoven, G.: Effects of CO2, continental distribution, topography and vegetation changes on the climate at the Middle Miocene: a model study, Clim. Past, 6, 675–694,, 2010. 

Herbert, T. D., Lawrence, K. T., Tzanova, A., Peterson, L. C., Caballero-Gill, R., and Kelly, C. S.: Late Miocene global cooling and the rise of modern ecosystems, Nat. Geosci., 9, 843,, 2016. 

Hilting, A. K., Kump, L. R., and Bralower, T. J.: Variations in the oceanic vertical carbon isotope gradient and their implications for the Paleocene-Eocene biological pump, Paleoceanography, 23, PA3222,, 2008. 

Hoorn, C., Wesselingh, F. P., Hovikoski, J., and Guerrero, J.: The development of the Amazonian mega-wetland (Miocene; Brazil, Colombia, Peru, Bolivia), in: Amazonia: Landscape and Species Evolution: a look into the past, 123, 142,, 2010. 

Indermühle, A., Stocker, T. F., Joos, F., Fischer, H., Smith, H. J., Wahlen, M., Deck, B., Mastroianni, D., Tschumi, J., Blunier, T., and Meyer, R.: Holocene carbon-cycle dynamics based on CO2 trapped in ice at Taylor Dome, Antarctica, Nature, 398, 121,, 1999. 

Janis, C. M., Damuth, J., and Theodor, J. M.: The origins and evolution of the North American grassland biome: the story from the hoofed mammals, Palaeogeogr. Palaeocl., 177, 183–198,, 2002. 

Jaramillo, C., Montes, C., Cardona, A., Silvestro, D., Antonelli, A., and Bacon, C. D.: Comment (1) on “Formation of the Isthmus of Panama” by O'Dea et al., Science Advances, 3, e1602321,, 2017. 

John, E. H., Wilson, J. D., Pearson, P. N., and Ridgwell, A.: Temperature-dependent remineralization and carbon cycling in the warm Eocene oceans, Palaeogeogr. Palaeocl., 413, 158–166,, 2014. 

Kirtland Turner, S. and Ridgwell, A.: Development of a novel empirical framework for interpreting geological carbon isotope excursions, with implications for the rate of carbon injection across the PETM, Earth Planet. Sc. Lett., 435, 1–13,, 2016. 

Kitchka, A. A., Tyshchenko, A. P., and Lysenko, V. I.: May. Mid-late Miocene sea level falls, gas hydrates decay, submarine sliding, and tsunamites in the Black Sea Basin, in: 78th EAGE Conference and Exhibition 2016, Vienna, Austria, 30 May–2 June 2016,, 2016. 

Knietzsch, M.-A., Schröder, A., Lucarini, V., and Lunkeit, F.: The impact of oceanic heat transport on the atmospheric circulation, Earth Syst. Dynam., 6, 591–615,, 2015. 

Knorr, G., Butzin, M., Micheels, A., and Lohmann, G.: A warm Miocene climate at low atmospheric CO2 levels, Geophys. Res. Lett., 38, L20701,, 2011. 

Krapp, M. and Jungclaus, J. H.: The Middle Miocene climate as modelled in an atmosphere-ocean-biosphere model, Clim. Past, 7, 1169–1188,, 2011. 

Lagabrielle, Y., Goddéris, Y., Donnadieu, Y., Malavieille, J., and Suarez, M.: The tectonic history of Drake Passage and its possible impacts on global climate, Earth Planet. Sc. Lett., 279, 197–211,, 2009. 

Lear, C. H., Coxall, H. K., Foster, G. L., Lunt, D. J., Mawbey, E. M., Rosenthal, Y., Sosdian, S. M., Thomas, E., and Wilson, P. A.: Neogene ice volume and ocean temperatures: Insights from infaunal foraminiferal Mg/Ca paleothermometry, Paleoceanography, 30, 437–1454,, 2015. 

LeGrande A. N. and Schmidt G. A.: Global gridded data set of the oxygen isotopic composition in seawater, Geophys. Res. Lett., 33, L12604,, 2006. 

Levitus, S., Locarnini, R. A., Boyer, T. P., Mishonov, A. V., Antonov, J. I., Garcia, H. E., Baranova, O. K., Zweng, M. M., Johnson, D. R., and Seidov, D.: World ocean atlas 2009, National Oceanographic Data Center (US), Ocean Climate Laboratory, United States, National Environmental Satellite, Data, and Information Service, Silver Spring, MD, USA, 2010. 

Londoño, L., Royer, D. L., Jaramillo, C., Escobar, J., Foster, D. A., Cárdenas-Rozo, A. L., and Wood, A: Early Miocene CO2 estimates from a Neotropical fossil leaf assemblage exceed 400 ppm, Am. J. Bot., 105, 1929–1937, 2018. 

Lunt, D. J., Valdes, P. J., Haywood, A., and Rutt, I. C.: Closure of the Panama Seaway during the Pliocene: implications for climate and Northern Hemisphere glaciation, Clim. Dynam., 30, 1–18,, 2008. 

Lynch-Stieglitz, J.: Tracers of past ocean circulation, Treatise on Geochemistry, 6, 433–451,, 2003. 

Mackensen, A. and Schmiedl, G.: Stable carbon isotopes in paleoceanography: atmosphere, oceans, and sediments, Earth-Sci. Rev., 197, 102893, 2019. 

Marchitto, T. M., Curry, W. B., Lynch-Stieglitz, J., Bryan, S. P., Cobb, K. M., and Lund, D. C.: Improved oxygen isotope temperature calibrations for cosmopolitan benthic foraminifera, Geochim. Cosmochim. Ac., 130, 1–11,, 2014. 

Marsh, R., Müller, S. A., Yool, A., and Edwards, N. R.: Incorporation of the C-GOLDSTEIN efficient climate model into the GENIE framework: ”eb_go_gs” configurations of GENIE, Geosci. Model Dev., 4, 957–992,, 2011. 

Micheels, A., Bruch, A. A. and Mosbrugger, V., 2009. Miocene climate modelling sensitivity experiments for different CO2 concentrations, Palaeontol. Electron., 12, 1–19, 2009. 

Micheels, A., Bruch, A. A., Eronen, J., Fortelius, M., Harzhauser, M., Utescher, T., and Mosbrugger, V.: Analysis of heat transport mechanisms from a Late Miocene model experiment with a fully-coupled atmosphere–ocean general circulation model, Palaeogeogr. Palaeocl., 304, 337–350,, 2011. 

Modestou, S. E., Leutert, T. J., Fernandez, A., Lear, C. H., and Meckler, A. N.: Warm middle Miocene Indian Ocean bottom water temperatures: Comparison of clumped isotope and Mg/Ca-based estimates, Palaeogeogr. Palaeocl., 35, e2020PA003927,, 2020. 

Montes, C., Cardona, A., Jaramillo, C., Pardo, A., Silva, J. C., Valencia, V., Ayala, C., Pérez-Angel, L. C., Rodriguez-Parra, L. A., Ramirez, V., and Niño, H.: Middle Miocene closure of the Central American seaway, Science, 348, 226–229,, 2015. 

Morley, R. J. and Morley, H. P.: Mid Cenozoic freshwater wetlands of the Sunda region, J. Limnol., 72, 18–35,, 2013. 

Mulch, A., Sarna-Wojcicki, A. M., Perkins, M. E., and Chamberlain, C. P.: A Miocene to Pleistocene climate and elevation record of the Sierra Nevada (California), P. Natl. Acad. Sci. USA, 105, 6819–6824,, 2008. 

Nisancioglu, K. H., Raymo, M. E., and Stone, P. H.: Reorganization of Miocene deep water circulation in response to the shoaling of the Central American Seaway, Paleoceanography, 18, 1006,, 2003. 

O'Dea, A., Lessios, H. A., Coates, A. G., Eytan, R. I., Restrepo-Moreno, S. A., Cione, A. L., Collins, L. S., De Queiroz, A., Farris, D. W., Norris, R. D., and Stallard, R. F.: Formation of the Isthmus of Panama, Science Advances, 2, e1600883,, 2016. 

Pagani, M., Arthur, M. A., and Freeman, K. H.: Miocene evolution of atmospheric carbon dioxide, Paleoceanography, 14, 273–292,, 1999. 

Pagani, M., Zachos, J. C., Freeman, K. H., Tipple, B., and Bohaty, S.: Marked decline in atmospheric carbon dioxide concentrations during the Paleogene, Science, 309, 600–603, 2005. 

Pearson, P. N. and Palmer, M. R.: Atmospheric carbon dioxide concentrations over the past 60 million years, Nature, 406, 695,, 2000. 

Pérez, L. F., Martos, Y. M., García, M., Weber, M. E., Raymo, M. E., Williams, T., Bohoyo, F., Armbrecht, L., Bailey, I., Brachfeld, S., Glüder, A., Guitard, M., Gutjahr, M., Hemming, S., Hernández-Almeida, I., Hoem, F.S., Kato, Y., O'Connell, S., Peck, V. L., Reilly, B., Ronge, T. A., Tauxe, L., Warnock, J., and Zheng, X.: Miocene to present oceanographic variability in the Scotia Sea and Antarctic ice sheets dynamics: Insight from revised seismic-stratigraphy following IODP Expedition 382, Earth Planet. Sc. Lett., 553, 116657,, 2021. 

Peterson, C. D., Lisiecki, L. E., and Stern, J. V.: Deglacial whole-ocean δ13C change estimated from 480 benthic foraminiferal records, Paleoceanography, 29, 549–563,, 2014. 

Pfuhl, H. A. and McCave, I. N., Evidence for late Oligocene establishment of the Antarctic Circumpolar Current, Earth Planet. Sc. Lett., 235, 715–728, 2005. 

Pierce, E. L., van de Flierdt, T., Williams, T., Hemming, S. R., Cook, C. P., and Passchier, S.: Evidence for a dynamic East Antarctic ice sheet during the mid-Miocene climate transition, Earth Planet. Sc. Lett., 478, 1–13,, 2017. 

Pound, M. J., Haywood, A. M., Salzmann, U., and Riding, J. B.: Global vegetation dynamics and latitudinal temperature gradients during the Mid to Late Miocene (15.97–5.33 Ma), Earth-Sci. Rev., 112, 1–22,, 2012. 

Price, A. R., Myerscough, R. J., Voutchkov, I. I., Marsh, R., and Cox, S. J.: Multi-objective optimization of GENIE Earth system models, Philos. T. R. Soc. A, 367, 2623–2633, 2009. 

Rae, J. W., Zhang, Y. G., Liu, X., Foster, G. L., Stoll, H. M., and Whiteford, R. D.: Atmospheric CO2 over the Past 66 Million Years from Marine Archives, Annu. Rev. Earth Pl. Sc., 49, 609–641,, 2021. 

Raymo, M. E., Ruddiman, W. F., and Froelich, P. N.: Influence of late Cenozoic mountain building on ocean geochemical cycles, Geology, 16, 649–653,<0649:IOLCMB>2.3.CO;2, 1988. 

Reiners, P. W., Ehlers, T. A., Garver, J. I., Mitchell, S. G., Montgomery, D. R., Vance, J. A., and Nicolescu, S.: Late Miocene exhumation and uplift of the Washington Cascade Range, Geology, 30, 767–770,<0767:LMEAUO>2.0.CO;2, 2002. 

Richey, J. N. and Tierney, J. E.: GDGT and alkenone flux in the northern Gulf of Mexico: Implications for the TEX86 and UK'37 paleothermometers, Paleoceanography, 31, 1547–1561, 2016. 

Ridgwell, A. and Schmidt, D. N.: Past constraints on the vulnerability of marine calcifiers to massive carbon dioxide release, Nat. Geosci., 3, 196,, 2010. 

Ridgwell, A., Hargreaves, J. C., Edwards, N. R., Annan, J. D., Lenton, T. M., Marsh, R., Yool, A., and Watson, A.: Marine geochemical data assimilation in an efficient Earth System Model of global biogeochemical cycling, Biogeosciences, 4, 87–104,, 2007. 

Ridgwell, A., Reinhard, C., van de Velde, S., Adloff, M., Monteiro, F., DomHu​​​​​​​, Wilson, J., Ward, B., Vervoort, P., kirtlandster​​​​​​​, and Li, M.,: derpycode/cgenie.muffin: Crichton et al. [revised for CP] (v0.9.22), Zenodo [code],, 2021a. 

Ridgwell, A., DomHu​​​​​​​, Peterson, C., Ward, B., sjszas​​​​​​​, evansmn, and Jones, R.: derpycode/muffindoc: (v0.9.20), Zenodo [code],, 2021b. 

Rögl, F.: Mediterranean and Paratethys. Facts and hypotheses of an Oligocene to Miocene paleogeography (short overview), Geol. Carpath., 50, 339–349, 1999. 

Rohling E. J.: Oxygen isotope composition of seawater, in: The Encyclopedia of Quaternary Science, edited by: Elias S. A., Elsevier, Amsterdam, 915–922, 2013. 

Schmidt, G. A. and Shindell, D. T.: Atmospheric composition, radiative forcing, and climate change as a consequence of a massive methane release from gas hydrates, Paleoceanography, 18, 1004,, 2003. 

Sepulchre, P., Arsouze, T., Donnadieu, Y., Dutay, J. C., Jaramillo, C., Le Bras, J., Martin, E., Montes, C., and Waite, A. J.: Consequences of shoaling of the Central American Seaway determined from modeling Nd isotopes, Paleoceanography, 29, 176–189,, 2014. 

Sijp, W. P., Anna, S., Dijkstra, H. A., Flögel, S., Douglas, P. M., and Bijl, P. K.: The role of ocean gateways on cooling climate on long time scales, Global Planet. Change, 119, 1–22,, 2014. 

Sosdian, S. M., Greenop, R., Hain, M. P., Foster, G. L., Pearson, P. N., and Lear, C. H.: Constraining the evolution of Neogene ocean carbonate chemistry using the boron isotope pH proxy, Earth Planet. Sc. Lett., 498, 362–376,, 2018. 

Srinivasan, M. S. and Sinha, D. K.: Early Pliocene closing of the Indonesian Seaway: evidence from north-east Indian Ocean and Tropical Pacific deep sea cores, J. Asian Earth Sci., 16, 29–44,, 1998. 

Steinthorsdottir, M., Coxall, H. K., de Boer, A. M., Huber, M., Barbolini, N., Bradshaw, C. D., Burls, N. J., Feakins, S. J., Gasson, E., Henderiks, J., and Holbourn, A.: The Miocene: the Future of the Past, Palaeogeogr. Palaeocl., e2020PA004037,, 2020. 

Steinthorsdottir, M., Jardine, P. E., and Rember, W. C.: Near-Future pCO2 During the Hot Miocene Climatic Optimum, Palaeogeogr. Palaeocl., 36, e2020PA003900,, 2021. 

Stewart, D. R., Pearson, P. N., Ditchfield, P. W., and Singano, J. M.: Miocene tropical Indian Ocean temperatures: evidence from three exceptionally preserved foraminiferal assemblages from Tanzania, J. Afr. Earth Sci., 40, 173–189,, 2004. 

Stoll, H. M., Guitian, J., Hernandez-Almeida, I., Mejia, L. M., Phelps, S., Polissar, P., Rosenthal, Y., Zhang, H., and Ziveri, P.: Upregulation of phytoplankton carbon concentrating mechanisms during low CO2 glacial periods and implications for the phytoplankton pCO2 proxy, Quaternary Sci. Rev., 208, 1–20, 2019. 

Super, J. R., Thomas, E., Pagani, M., Huber, M., O'Brien, C. L., and Hull, P. M.: Miocene Evolution of North Atlantic Sea Surface Temperature, Palaeogeogr. Palaeocl., 35, e2019PA003748,, 2020. 

Toyos, M. H., Lamy, F., Lange, C. B., Lembke-Jene, L., Saavedra-Pellitero, M., Esper, O., and Arz, H. W.: Antarctic circumpolar current dynamics at the Pacific entrance to the Drake Passage over the past 1.3 million years, Palaeogeogr. Palaeocl., 35, e2019PA003773,, 2020. 

Tyrrell, T. and Zeebe, R. E.: History of carbonate ion concentration over the last 100 million years, Geochim. Cosmochim. Ac., 68, 3521–3530,, 2004. 

Unger, N. and Yue, X.: Strong chemistry-climate feedbacks in the Pliocene, Geophys. Res. Lett., 41, 527–533, 2014. 

University of Texas (Institute for Geophysics): PLATES, (last access: 14 December 2019), 2015. 

Urban, T. and Hardisty, F.: Developing PAMS – A Paleolocation Web Service, in: International Cartographic Conference, Dresden, Germany, 30 May–2 June 2016 128, 2013. 

Valdes, P. J., Armstrong, E., Badger, M. P. S., Bradshaw, C. D., Bragg, F., Crucifix, M., Davies-Barnard, T., Day, J. J., Farnsworth, A., Gordon, C., Hopcroft, P. O., Kennedy, A. T., Lord, N. S., Lunt, D. J., Marzocchi, A., Parry, L. M., Pope, V., Roberts, W. H. G., Stone, E. J., Tourte, G. J. L., and Williams, J. H. T.: The BRIDGE HadCM3 family of climate models: HadCM3@Bristol v1.0, Geosci. Model Dev., 10, 3715–3743,, 2017. 

Watterson, I. G.: Non-dimensional measures of climate model performance. Int. J. Climatol., 16, 379–391,<379::AID-JOC18>3.0.CO;2-U, 1996. 

Westerhold, T., Marwan, N., Drury, A. J., Liebrand, D., Agnini, C., Anagnostou, E., Barnet, J. S., Bohaty, S. M., De Vleeschouwer, D., Florindo, F., and Frederichs, T.: An astronomically dated record of Earth's climate and its predictability over the last 66 million years, Science, 369, 1383–1387, 2020.  

Wilson, J. D., Monteiro, F. M., Schmidt, D. N., Ward, B. A., and Ridgwell, A.: Linking marine plankton ecosystems and climate: A new modeling approach to the warm early Eocene climate, Paleoceanography and Paleoclimatology, 33, 1439–1452,, 2018. 

Zachos, J. C., Dickens, G. R., and Zeebe, R. E.: An early Cenozoic perspective on greenhouse warming and carbon-cycle dynamics, Nature, 451, 279,, 2008. 

Zweng, M. M., Reagan, J. R., Antonov, J. I., Locarnini, R. A., Mishonov, A. V., Boyer, T. P., Garcia, H. E., Baranova, O. K., Johnson, D. R., Seidov, D., and Biddle, M. M.: World ocean atlas 2013, 2, Salinity, edited by: Levitus, S. and Mishonov, A., NOAA Atlas NESDIS 74, 2013. 

Short summary
The middle Miocene (15 Ma) was a period of global warmth up to 8 °C warmer than present. We investigate changes in ocean circulation and heat distribution since the middle Miocene and the cooling to the present using the cGENIE Earth system model. We create seven time slices at ~2.5 Myr intervals, constrained with paleo-proxy data, showing a progressive reduction in atmospheric CO2 and a strengthening of the Atlantic Meridional Overturning Circulation.