The relative role of direct orbital forcing versus CO_2 and ice feedbacks on Quaternary climate

Domain: biology
Published: 2026-03-19
Editorial outcome: major_revision

Abstract
During the Quaternary (the last 2.58 million years), Earth’s climate has fluctuated between glacials and
interglacials, paced by external forcings and mediated by internal feedbacks. However, General
Circulation Models (GCMs), essential for addressing the mechanisms associated with these fluctuations,
require substantial computational resources, meaning they are unsuitable for exploring orbital-scale
variability on million-year timescales. Here, we use a GCM to calibrate a faster statistical model, or
emulator, and apply this to the Quaternary.
We show a good agreement between the emulated climate and proxy data over the last 800,000 years,
especially the timing of glacial-interglacial cycles. A series of sensitivity experiments allows us to
identify the dominant components driving long-term climate change. The results show that a
combination of the CO and ice sheet feedbacks provide the dominant contribution to the annual mean
2
temperature signal, with the direct orbital radiative forcing playing only a minor role.

Introduction
For the last few million years (Myr), the Earth’s climate has fluctuated between cooler (glacial) and
warmer (interglacial) intervals, as revealed by proxy records of palaeoclimate1. The timing and
magnitude of these glacial-interglacial variations are forced solely by changes in the Earth’s orbit and
axis inclination around the sun, the three relevant parameters of which are precession, obliquity and
eccentricity (which have main periods of ~23 thousand years (kyr), ~41 kyr, and ~96 and 400 kyr,
respectively2,3. The direct orbital forcing is mediated by internal feedbacks in the climate system,
including changes in the concentration of atmospheric CO and variations in the extent and thickness of
2
global ice sheets themselves. For the Quaternary (2.58 million years ago (Ma) to present day),
comprising the Pleistocene, Holocene and, as some argue, the Anthropocene4), glacial-interglacial
fluctuations are evident in a relatively large number of palaeoclimate records, constructed from various
climate proxies3. However, although proxies for the Quaternary climate provide a clear picture of
climate variability, especially over the last 800 kyr, the interplay between forcings and feedbacks
(associated with, for example, changes in the extent and thickness of ice sheets), is not well known5.
Importantly, it is still not clear what proportion of the signals of Quaternary climate change is driven
directly by the radiative effect of orbital forcing, and how much is associated with CO and/or ice sheet
2
feedbacks which mediate and amplify that orbital forcing. Although this question of external forcings
versus internal feedbacks has been investigated by previous studies, because of computational
constraints (discussed below) almost all have only looked at certain periods during the Quaternary, such
as the Last Glacial Maximum6 (LGM, ~21 kyr ago (ka), or have focused on shorter timespans, such as
the last ~400 kyr7.
To simulate changes in past climate that are physically-based and spatially resolved, climate models are
required. More specifically, fully-coupled General Circulation Models (GCMs), or Earth System
Models (ESMs) have proved useful for investigating the driving mechanisms, dynamics, feedbacks, and
sensitivity of the climate system associated with variations in orbital forcing, CO , and ice sheets8.
2
However, because of their structural complexity, along with their relatively high spatial and temporal
resolutions, these models require substantial computational resources and time to run. As a result, they
often cannot be run in fully-coupled transient mode over the long timescales (> 100 kyr) over which
these forcings and feedbacks vary. Some studies have applied GCMs to conduct transient simulations
of past climatic changes, but this has generally been in the form of a number of snapshot equilibrium
simulations (the assumption being that the system is never far from equilibrium with the orbital
forcings) in order to build a transient simulation9; however, even these simulations are typically limited
to a single glacial-interglacial cycle. Alternatively, a limited number of fully transient GCM simulations
have been conducted, such as in the framework of the Deglaciation PMIP experiment10 or the PalMod
project11, but these in general simulate at most a single glacial-interglacial cycle. As such, simulations
of multiple glacial-interglacial cycles with a full GCM or ESM are generally not practical, and the
number of sensitivity studies that can be performed is limited. As far as we are aware, to date only one
GCM simulation can be considered fully transient12, which involved running the Community Earth
System Model version 1.2 (CESM1.2) for the entire Quaternary. This itself was an extension of earlier
work13, which ran the same model for the last 2 Myr to examine climate-driven habitat suitability of five
hominin species. However, even this cannot be considered as a classic transient simulation, because the
period was split into 42 chunks (of between 32-125 kyr each) that were run in parallel, and then
concatenated afterwards to create a continuous timeseries12. In addition, the simulation also used the
acceleration technique14, in which some of the boundary conditions (such as the orbital parameters) are
accelerated by a certain factor, resulting in a compressed simulation12; in this example, the acceleration
technique resulted (using a factor of 5) in the 3 Myr simulation being compressed to 600,000 model
years12. Arguably, a different (but complementary) approach is to use models with a lower complexity,
such as Earth system Models of Intermediate Complexity (EMICs), to simulate long-term transient past
climate changes15. However, EMICs provide a more schematic and/or lower spatial representation than
GCMs, often with a relatively idealised representations of atmospheric processes that do not capture
enough detail to understand climate variability, even at these long timescales.
Therefore, to resolve this competing desire for maximum process complexity and sufficient spatial
resolution versus practical computing resource requirements and running time, statistical methods are
being developed and applied to simulate long-term changes in climate. One emerging approach is that
of emulators, which are statistical models that are calibrated on data from a more complex climate
model, such as a physics-based GCM. They give a projection of the climate resulting from a certain set
of input conditions (climate drivers), along with an estimation of the uncertainty (as shown by the
variance at every timestep) associated with the projection. It should be noted, as a point of clarification,
that while we refer to the emulator simulations in this paper, strictly data are interpolated rather than
simulated. However, whilst the climate projection is based on the physics of the GCM and is at the
same spatial resolution, the requirements of the emulator in terms of computational resources and time
are a fraction of those required for the full GCM, once the emulator has been trained on existing GCM
simulations. This makes them useful for investigating processes and climatic changes occurring over
long time scales of hundreds of kyr or longer, thus allowing for comparison to proxies for climate16 and
model validation17. Emulators have been applied in sensitivity analyses of climate to orbital,
atmospheric CO , and ice sheet configurations18,19 and to investigate parametric uncertainty in models20.
2
Here, we apply a climate emulator to simulate and understand the evolution of long-term climate change
in the Quaternary. In short, the emulator is initially trained on a number of GCM simulations designed
to capture the full range of climate conditions (including global sea level, GSL) from very cold (e.g. the
LGM) to very warm (e.g. the Pliocene, 5.33-2.58 Ma), and is driven by five components, each relatively
well understood global metrics: three orbital parameters (based on astronomical predictions of planetary
motions), a proxy-derived CO component, and a proxy-derived GSL component, used to represent ice
2
sheets. Many other feedbacks/processes are represented in the emulator because they are included in the
GCM itself (e.g. atmosphere and ocean dynamics (including the Atlantic Meridional Overturning
Circulation, AMOC) and thermodynamics, clouds, sea ice, and vegetation). However, some
feedbacks/processes are not represented in the GCM, and are therefore not included in the emulator,
including dust, methane, lakes and freshwater fluxes associated with changes in ice sheet volume.
Likewise, it should be noted that the emulator is a purely statistical model designed to replicate the
output of the GCM. The emulator does not explicitly include any representation of any physical
processes, but, as explained above, because they are represented in the GCM, the output of the emulator
does implicitly take account of these processes. In principle, any variable that can be output by the
GCM can be emulated, but here we choose to emulate near-surface air temperature (SAT) and mean
annual precipitation (hereafter PREC), because these are the variables that can most robustly be
evaluated with proxy data. A full description of the emulator, the process of creating the training data,
and then building, optimising, testing and running the emulator, is described in the Methods. After a
brief validation exercise in which the emulator output is compared to proxy data, the main purpose here
is to use the emulator to answer the above question i.e. what proportion of Quaternary climate change is
directly orbitally-forced versus CO and/or ice sheet feedbacks, something which the efficiency of this
2
approach allows given that the emulator takes several orders of magnitude less computational cost (see
Methods, in particular Efficiency and computational cost) than if a climate model alone were used.

Results
All drivers simulation: Tmodel-data comparison

Figure 1 compares the proxy and emulated temperature data over the last 2.58 Myr when all forcings
and feedbacks are included in the emulator (simulation E ; see Methods, and in particular Section 5.3
11111
for more details, but in short this is when all forcings and feedbacks in the emulator are switched on, in
comparison to e.g. E where they are all switched off), at one location in Antarctica from the Dome
00000
C ice core, providing δD-based annual-mean SAT, and four Ocean Drilling Program (ODP) sites
providing sea surface temperature (SST) from the UK ’-alkenone proxy (see Table 1 for details,
37
Methods for the selection criteria used to obtain the sites, and Figure S1 in the Supplementary Material
for the geographical location of the sites).

The results from Figure 1 are discussed further below in quantitative terms but, qualitatively, the timing
and duration of many of the interglacials and glacial maxima are similar in the model and proxy
reconstructions, and the agreement between the magnitude of variations is generally very good (Fig. 1).
These results (focusing on just the last 800 kyr) are also shown as a scatterplot in the Supplementary
Material (Fig. S2), where it can be seen that the relationship between emulated and proxy temperature
varies according to location but is strongest at Dome C (r2 = 0.73). At many of the locations there are
instances when the relative maxima and minima are not well reproduced by the emulator, such as the
interglacials at approximately 330 kyr, 240 kyr, and 125 kyr before present (BP), corresponding to
Marine Isotope Stages (MIS) 9, 7 and 5, respectively (Fig. 1). Again, this is further demonstrated by
Figure S2, where these interglacials appear warmer in the proxy data relative to the emulated
temperature, matching the biased periods identified above in Figure 1. At Dome C (Fig. 1a) there is a
relatively large amount of SAT variability in the glacial-interglacial cycles, which continues back in
time until approximately 1.5 Myr and then decreases even further back in time towards the early
Quaternary (in line with the other ODP sites, Fig. 1b-e); possible reasons for this are discussed below.
At Dome C in particular (Fig. 1a), where the temporal range of proxy data is currently limited by
existing drilling methods, the emulator gives an indication of what deeper (i.e. older) temperature
reconstructions (as pushed for by consortia such as Beyond EPICA, https://www.beyondepica.eu/en/)
throughout the Quaternary might resemble there, if the CO used to drive the emulator21, is realistic.
2

The results presented here are in line with the previously-discussed simulation12. Concerning
magnitudes of change (relative to the preindustrial, PI), they found a global-mean glacial cooling of
5.7°C during the LGM in their simulation, consistent with their combined proxy-model estimate of 6.3-
5.6°C cooling12. This is similar, albeit slightly colder, to our results; here the emulator shows a global
mean cooling of 5.1°C during the LGM. This is discussed further in the Supplementary Material
(Section S1 and Fig. S3). In Antarctica, however, the results presented here are in better agreement with
the previously-discussed simulation; an LGM cooling of 9.5°C12 compared to ~10°C here (Fig. 1a).
Concerning temporal trends, a clear decrease in global mean temperature anomalies and interannual
variability in shown in the previous simulation12 until the Mid-Pleistocene Transition (MPT, 1.25-0.7
Ma), after which the trend stabilises and variability increases12. Here, although we show the same
change in variability before and after ~1 Myr, the emulated temperatures are not showing the long-term
decline at any of the record locations (Fig. 1). Possible reasons for this are discussed below.

Figure 1 demonstrates that the proxy data show more glacial-interglacial amplitude than the emulator,
especially in SST prior to 1 Myr (e.g. Fig. 1c-e). In fact, for the emulator to reproduce the variability
shown by proxy data (at each individual location), a more extreme synthetic CO curve is needed as
2
input to the emulator. To obtain this (and again taking advantage of its efficiency), the emulator was
run multiple times with idealised (constant) CO ; the resulting temperature, in conjunction with the
2
proxy data (using the same five locations as in Figure 1), was used to generate timeseries of CO via
2
interpolation of the available data; for details see Methods (Section 5.2), and Figure S4 in the
Supplementary Material for the temperatures associated with the idealised CO . The resulting synthetic
2
CO curves are shown in Figure S5 of the Supplementary Material (for just the last 1600 kyr,
2
corresponding to the maximum extent of the available proxy data), where the synthetic CO can be seen
2
against available reconstructed and modelled CO , at each of the five locations on which the synthetic
2
CO is based. In short, a comparison between the synthetic CO and actual CO suggests that the
2 2 2
emulator would need input CO that, in some locations, is more than 1.5 times larger than its current
2
driving CO to match the variability in temperature shown by the proxy data (Fig. S5). At Dome C, the
2
synthetic CO corresponds well in terms of variability to the observations, but the magnitudes are often
2
greater than the proxy data, especially during the interglacials discussed above (Fig. S5a). Elsewhere, in
some locations (such as in the central South Pacific) the synthetic CO is much greater in both
2
magnitude and variability relative to either reconstructed or modelled CO (Fig. S5c). The observation
2
that this synthetic CO is very different in places to the driving CO might be explained by: i) the low
2 2
variability shown by the emulator is correct but the proxy temperature data are too variable; ii) the
variability of the proxy temperature data is correct but the CO records driving the emulator are causing
2
the resulting temperature to underestimate the correct amount of variability (this can only apply prior to
the ice core record, because the ice-core CO is assumed to be accurate); or iii) using other additional
2
driving components (such as those discussed above that are not included in the GCM) might improve
this match.

Concerning PREC, the composite Chinese speleothem δ18O data records (we acknowledge that this is
not the same as PREC, but is used here as a first-order comparison for precipitation) are compared with
emulated precipitation in Figure 2. As with temperature, the emulator reproduces the major variations
evident in the δ18O record reasonably well, particularly with regards to the precession-driven timing of
maxima and minima in precipitation. Unlike SAT and SST, however (where the level of variability
decreases further back in time), a similar level of variability continues when the emulator is run further
back beyond the proxy data records, and this is true at the three sites considered (Fig. 2).

To quantify the model-proxy comparison, we calculate the arcsin Mielke measure, or M score22 to
assess the goodness of fit between two variables with the same units, and therefore the skill of the
emulator. The resulting skill score has a maximum possible value of 1000, while a score of 0 or
negative values indicates no skill22. See Methods, in particular Metric for quantifying model
performance. We calculate an M score between the emulated temperature and proxy data, for the all
drivers simulation (in other words, a simulation in which all five drivers - three orbital parameters, CO
2
and GSL - varied according to proxy data, denoted E ) throughout the 2.58 Myr; see Table 1. SAT
11111
over Dome C is shown to be the best match with the proxy data, with an M score of 620 (out of a
possible 1000, which represents an exact agreement). For SST at ODP locations, M scores are 328,
200, 129 and 128 (Table 1).

It is important to put these M scores into context. A point of reference is to calculate the M score for a
constant climate signal with the same mean as the proxy record (i.e. keeping the driving CO
2
, orbit, and
ice as constant, averaged over the last 800 kyr; see Methods and in particular Section 5.1 for details).
As expected, at all sites the modelled variations perform substantially better than a constant signal; for
example, at Dome C the original emulated results gave an M score of 620 whereas the constant climate
signal based on the mean gave an M score of 0.65, the constant climate signal based on the mean LGM
gave an M score of 0.48, and a random climate signal gave an M score of -0.04. As a further
contextualisation, we calculate an M score between the model and proxy at every possible combination
of sites (see Supplementary Material, Table S1). We assess agreement as good for a particular site if the
best M score is obtained when the model results at that site are compared with proxies from the same
site (i.e the diagonals in Table S1). Although some sites (such as Dome C and ODP 1239) show good
agreement, other sites do not. However, this can be explained by the relative geographical proximity of
some of the sites; proxy data at ODP 1012, for example, are giving the best match to emulated data at
ODP 1239 (M score = 360), both of which are in the eastern Pacific. Likewise, agreement is sometimes
best when the sites are geographically distant but, for example, at a similar latitude.

Single driver experiments
Following the model-data comparison, the relative impacts of the individual drivers on emulated SAT
were investigated. To overcome potential disadvantages of linear factorisation (see Section S2 in the
Supplementary Material text, and Figures S6 and S7, for the results), a linear-sum/shared-interaction
factorisation was also conducted23, again globally. This is shown in Figure 3, with the results being
very similar (but more robust) compared with the linear factorisation (Figure S7).

The key finding is that the direct impact of the orbital parameters (i.e. obliquity, eccentricity and
precession) on the M scores is negligible and they are thus contributing little to the overall signal (Fig.
3b-d). Of the three orbital parameters, obliquity provides the highest contribution to the overall signal,
particularly over high latitudes (e.g. the Southern Ocean) with M scores reaching ~200-300 as opposed
to ~100 or less from eccentricity and precession (Fig. 3b-d). It should be noted that eccentricity and
precession are combined as two forcings, one of ecosω and one of esinω, consistent with previous
work18 and motivated by the fact that eccentricity modulates the precession forcing. Also, as discussed
in more detail below (see Discussion), we are not suggesting that the overall orbital forcing is
negligible, but that the direct orbital radiative forcing on the annual mean temperature response is
negligible. The orbital forcing does have a very significant indirect impact, through controlling the CO
2
and ice sheet forcings (which, given that these are used to drive the emulator, therefore account for the
indirect impact of orbital forcing), and also has a significant direct impact on seasonal precipitation (e.g.
the West African monsoon). In contrast to this negligible orbital contribution, the results also show that
CO and ice are contributing the most to the overall signal, particularly over the tropical landmasses (for
2
CO , Fig. 3a) and polar regions (for ice, Fig. 3e). When averaged globally CO is contributing 54% to
2 2
the overall signal and ice is contributing 35%, whereas the orbital parameters’ contributions are all < 5%
(Fig. 3f).

Methods
An emulator is a statistical representation of a more complex model, in this case a GCM. The
underlying principle is that a relatively small number of simulations are carried out initially using the
GCM, which sample the multidimensional climate driver input space (using Latin hypercube sampling,
see below); in this case, five dimensions consisting of three direct orbital forcing dimensions, a CO
2
feedback dimension, and a GSL feedback dimension. The emulator is calibrated on (or tuned/trained to)
these GCM simulations, including the full two-dimensional latitude-longitude global coverage of the
grid GCM data, with the aim of being able to interpolate the GCM results such that a prediction is made
of the output that the GCM would produce if it were run using any particular combination of input
climate drivers17. The emulator also produces an estimation of the uncertainty associated with the
prediction. For any single combination of values of these climate forcings, the emulator simulates the
corresponding state of a particular climate variable, in this case global two-dimensional fields of the
SAT and PREC. As such, temperature and precipitation emerge as output from the emulator for any
combination of input forcings. The emulator can be validated by comparing results with additional
GCM results not included in the original training stage17. The combination of the emulator and the
GCM itself can be validated by comparing the resulting output to existing proxy reconstructions (i.e. a
model-proxy validation). The process of building, optimising, testing and running the emulator is
described below, but first a description of the emulator’s efficiency and low computational cost is
provided.

Efficiency and computational cost
To run over the full 2.58 Myr at a temporal resolution of 1 kyr, calculate comparisons with the proxy
data, and produce output files and figures, the emulator ran at several orders of magnitude faster than a
standard GCM (minutes, as opposed to months or years). As a comparison, for a GCM such as
HadCM3 to carry out the same process (i.e. produce a quasi-transient simulation of 2,580,000 years),
even at its optimal speed of approximately 120 model years/day, it would take ~59 years to run. Even
the shorter, 800,000-year simulation would take just over 18 years to run. Therefore, carrying out a
single quasi-transient simulation of the entire Quaternary using this model, let alone the 31 additional
sensitivity simulations of the late Pleistocene, would be impossible. It should be noted that prior to
running the emulator, 182 HadCM3 simulations were undertaken to provide the training data for the
emulator. The computationally expensive part of the process, therefore, is the creation of the training
data (i.e. the HadCM3 simulations), but this only needs to be done once, at the beginning; thereafter, for
any quasi-transient simulation or sensitivity tests, the emulator is several orders of magnitude faster than
a GCM. Moreover, the emulator can be run, at this speed, on a standard desktop PC, whereas HadCM3
requires resource-intensive supercomputing infrastructure, clearly demonstrating the efficacy of the
emulator and its suitability for subsequent work.

Step 1: Creation of training data
The emulator was calibrated on simulations run using the HadCM3 climate model, a fully coupled
atmosphere-ocean GCM developed by the UK Met Office. The specific model setup of the GCM is
denoted HadCM3B-M2.1aE35; hereafter referred to as HadCM3. The horizontal resolution of the
atmospheric component is 2.5° latitude by 3.75° longitude with 19 vertical levels, whilst the ocean has a
resolution of 1.25° by 1.25° and 20 vertical levels. Please see Supplementary Material (Section S4) for
more details on HadCM3.

Overall, a 182-member ensemble of HadCM3 simulations was used to train the emulator, with input
parameters varying across five dimensions: three main orbital parameters of longitude of perihelion (ϖ),
obliquity (ε) and eccentricity (e), atmospheric CO concentration (from ~180 ppmv to ~1900 ppmv) and
2
GSL (from -100 m to +30 m relative to modern). This is one more dimension than a previous study
based on these simulations, which did not include GSL17. The 182-member ensemble consisted of two
60-member sub-ensembles (modice and lowice) and a 62-member sub-ensemble (highice).

The two 60-member ensembles had ranges of orbital and CO values (see Table S3) that were
2
appropriate for a wide range of palaeo and future climate states17. Eccentricity and longitude of
perihelion were combined under the forms esinϖ and ecosϖ given that, in general at any point in the
year, insolation can be approximated as a linear combination of these two terms and obliquity36. Latin
hypercube sampling was used to optimally fill the 4-dimensional input space of orbit and CO . The two
2
60-member sub-ensembles differed solely by having two different fixed ice sheet configurations,
consisting of changes to the Greenland and West Antarctic ice sheets (GrIS and WAIS, respectively)
only. One sub-ensemble had present-day ice sheet extents, denoted the modice sub-ensemble, and the
other had reduced extents, based on the PRISM4 mid-Piacenzian (a warm interval, ~3.3-3.0 Ma, during
the Pliocene) reconstructions37, denoted the lowice sub-ensemble. The lowice sub-ensemble was
created from the modice sub-ensemble using the anomaly approach17. A third sub-ensemble consisting
of 62 GCM simulations was also run with HadCM3; a series of snapshot simulations covering the last
glacial cycle (LGC, ~115-11 kyr before present (BP)), forced by appropriate changes in orbit,
atmospheric CO concentration and ice sheet evolution, denoted the highice sub-ensemble24. All of the
2
HadCM3 simulations were initialised from preindustrial climatological conditions and run for several hundred
years, by the end of which they had all reached a quasi-equilibrium state at the surface and subsurface;
the training data used here, therefore, can be considered independent of initial conditions, as each
ensemble member represents a stable climate state.

Step 2: Building the emulators
Creating the interglacial and glacial training data
The emulator described here is an updated version of a previous emulator17, updated firstly by including
an additional driving term of ice sheet volume (using GSL as a proxy for this), and secondly by adding
additional HadCM3 simulations to the calibration that included various ice sheet configurations (sub-
ensemble highice, as described in Step 1 above). In practice two separate emulators were used, one
calibrated on the lowice and modice 60-member sub-ensembles and used for interglacial conditions
(defined as GSL >= 0 m), and one calibrated on the 60-member modice and 62-member highice sub-
ensembles and used for glacial conditions (defined as GSL < 0 m). Taken together, therefore, the
HadCM3 simulations cover the full range of GSL including LGM, glacial-interglacial cycles and
Pliocene-like conditions. Please see Supplementary Material (Section S3) for more details on creating
the training data.

The approach of having two separate emulators (shown in the Supplementary Material, Fig. S9) with
separate training data for the interglacial emulator applied for GSL >=0 m compared with the glacial
emulator for GSL < 0 m, was taken to ensure that there was no leakage of state-specific climate
variations across the glacial-interglacial transition. The concept and effect of this leakage can be
illustrated (see Figure S10 in the Supplementary Material) by considering surface air temperature for
two climate states, one with an emulated climate with GSL = +5 m (Fig. S10a and b), and one with an
emulated climate with sea level = -5 m (Fig. S10c and d), firstly as produced by the separate emulators
(Fig. S10a and c, as used in this paper) and secondly produced by a combined emulator, in which all
training data are used to train a single emulator (Fig. S10b and d). For the +5 m climate state, both the
separate and combined emulators show clear warming associated with melt of the Greenland and
Antarctic ice sheets (Fig. S10a and b). This is correct, and appears because the lowice training data
include ice sheet reduction in these regions, associated with the PRISM4 mid-Pliocene ice sheet
reconstruction. However, for the -5 m climate state, there are clear differences between the results of
the separate and combined emulator versions in these regions (Fig. S10c and d). For the combined
emulator (Fig. S10d), the algorithm is still influenced by the Pliocene lowice training data, producing
anomalous cooling in the regions that show melt in the lowice training data. However the two separate
emulators are, for this GSL, only influenced by the highice and modice training data, and so do not have
this unphysical leakage across the glacial-interglacial transition (Fig. S10c). In essence, the use of
separate emulators is required because in the real climate system, warming relative to modern results in
ice sheet loss primarily in the Greenland/Antarctica region, whereas cooling relative to modern results
in ice sheet growth primarily in the Laurentide/Fennoscandian region. The use of two separate
emulators therefore ensures that the different states are characterized separately and correctly. In
addition, checks were carried out to ensure that there were no discontinuities across the glacial-
interglacial transition between the climate projections from the two different emulators. This lack of
discontinuity is unsurprising, given that the modice ensemble was common to both emulators, and sits at
the surface GSL=0, on the boundary of the two emulators. We therefore consider it appropriate to use
the two separate emulators, as combining them into one leads to higher uncertainties at the boundary
between glacials and interglacials, due to the leakage effect discussed above.

The reconstructed ice sheet extents in the highice ensemble are based on the ICE-5G model32. These
reconstructions were primarily selected because the associated HadCM3 model simulations already
existed24. The reconstructions are based on palaeo data (GSL, local sea level, and ice-sheet extent) and
include the range of associated data that was required for this work. In addition, for glaciations prior to
the last glacial cycle, there is very little or no palaeo data available that enables the global 3-D
reconstructions which are required here. There are, however, other reconstructions (see Section S3 in
the Supplementary Material); although we acknowledge the existence of these, the lowice, modice and
highice ensembles used here, when combined, do capture changes in climate and ice sheet extent
ranging from interglacial states to full glacial conditions

Input parameter details and optimisation
In this study, a principal component analysis (PCA) Gaussian process (GP) emulator is used38, with the
subsequent Bayesian treatment39,40 and a PCA approach41. Climate data for the entire global grid, from
the aforementioned HadCM3 simulations, are used to calibrate the emulators, as opposed to calibrating
separate emulators based on data for individual proxy locations. This approach is taken because, for
past climate, the global response overall is of interest, rather than just the response at specific locations
individually. It also means that the results are consistent across all locations.

The driving parameters for the interglacial emulator can be expressed as a 120 × 5 matrix (n × p),
representing n = 120 HadCM3 simulations and p = 5 input forcing parameters, as discussed above.
These 120 × 5 input driving parameters are single values; the emulator is trained on the corresponding
120 spatial patterns of temperature and precipitation from the 120 simulations coming out of HadCM3,
expressed as PCAs (see below).

The matrix containing the training data from HadCM3 had dimensions of longitude × latitude × n = 96
× 73 × N, for the interglacial (N=120) and glacial (N=122) emulator, respectively. A PC analysis was
performed on each of the n GCM simulations, the results of which were then used to calibrate the two
emulators. Five correlation length hyperparameters (δ) were used, one for each of the five input factors,
which describe the smoothness of the climate response in the HadCM3 data to the input conditions17. A
nugget term (υ) was also used, which accounts for any non-linearity in the output response to the inputs,
non-explicitly specified inactive inputs and the effects of lower-order PCs that are excluded from the
emulator. The optimal values for these hyperparameters and the number of PCs retained were
calculated during calibration and evaluation of the emulator.

Four emulators were created in total, two of which were trained on SAT data and two of which were
trained on PREC data, for an interglacial (trained on modice and lowice ensembles) and glacial (trained
on modice and highice) version, respectively. Before being used, the emulator configurations were
optimised, which involved generating an ensemble of emulators with different numbers of PCs retained
(the majority of the PCs were discarded as they only account for a very small amount of total variation
in the HadCM3 data) and different values for a number of hyperparameters used in the emulator,
including the correlation length hyperparameters (δ) and the nugget parameter (ν). The performances of
the different emulators were compared in order to identify the optimal number of PCs to retain and the
hyperparameter values to use. It has been demonstrated that the optimisation can be carried out on
SAT, and then the optimised configuration applied to other climate variables (e.g. precipitation), with no
significant loss of performance17. Further details on the optimisation process are discussed in Section
S5 of the Supplementary Material.

Step 3: Testing the optimized emulator
Leave-one-out approach
The performance of the optimised emulators was evaluated, firstly by using a leave-one-out cross-
validation approach. Here, a series of emulators was constructed and used to predict one left out
HadCM3 simulation each time. For example, for the interglacial emulator (normally based on 120
HadCM3 simulations), the emulator was built on 119 simulations and used to reproduce the 120th, with
the results then compared to the actual 120th HadCM3 simulation. This was then repeated for every
possible combination (N=120) of left-out simulations. The number of grid boxes for each experiment
calculated to lie within different SD bands, and the root mean squared error (RMSE) averaged across all
the emulators were used as a performance indicator. The results suggest that that the emulators are able
to reproduce the left out ensemble simulations reasonably well (see Figure S11 in the Supplementary
Material for individual measures of this e.g. the percentage of grid boxes for which the emulator
predicts the SAT of the left out experiment to within 1, 2, 3 and >3 SD), with no obvious systematic
errors in their predictions.

Last Glacial Maximum approach
One of the input training simulations in the highice ensemble is an LGM simulation. Because extensive
proxy data exist for the LGM, this allows us to evaluate the HadCM3 climate model itself for this time
period. At the same time, we can evaluate how well the emulator reproduces the HadCM3 LGM
simulation (a less stringent test than the leave-one-out test described above, as the target simulation is
part of the training data). This evaluation is discussed in Section S1 of the Supplementary Material.

Step 4: Running the emulators
Model-data comparison
In addition to evaluating the emulator/HadCM3 output with the aforementioned LGM reconstructions,
another way of testing the emulator is to compare the output with timeseries of existing palaeo proxy
data.

For SAT and PREC separately, the emulator was therefore run for the last 800 kyr (i.e. the late
Pleistocene and Holocene). This was achieved by running the emulator 801 times, once for each 1 kyr
timeslice between the present (i.e. 0 kyr) and 800 ka, each with appropriate values for the 5 driving
parameters. For timeslices that had a GSL value of equal to or higher than 0 m (equivalent to present-
day interglacial conditions), the interglacial emulator was used, whilst for timeslices that had a GSL
value lower than 0 m the glacial emulator was used. Hereafter, and throughout the main manuscript,
these 4 emulators (two SAT, two PREC; two interglacial, two glacial) are known collectively as the
emulator. For each timeslice, the input to the emulator is the five input forcing parameters, either
calculated from astronomical solutions2 (for the 3 orbital parameters) or from proxy data (for CO and
2
GSL), as described below. The resulting emulated latitude-longitude SAT and PREC were then
compared with the proxy data at the location of the various sites, as described in the Results section.
The last 800 kyr were initially selected because a range of high-resolution palaeo records exist, of
multiple variables and in multiple locations, providing both forcing data for the emulator and proxy
climate reconstructions which can be compared to the emulator results. However, some of the sites
provide reconstructions going further back in time and so, once the above model-data comparison had
been undertaken, an additional emulator simulation of 2.58 Myr (i.e. the entire Quaternary) was
conducted.

Forcing parameters used when running the emulator
The palaeo data used to force the emulator are illustrated in Figure S8 of the Supplementary Material.
For the first 800 kya, orbital data (Fig. S8a) were calculated at 1 kyr resolution2. For CO , a composite
2
record of observed atmospheric CO was used (Fig. S8b), which was measured from the Dome C ice
2
cores from Antarctica42,43, interpolated to 1 kyr resolution. The sea level stack reconstructed from δ18O
data from ocean sediment cores44 was used in order to provide the GSL index used as the fifth driving
parameter to the emulator (Fig. S8c), which represents changes in global ice sheet volume, again
interpolated to 1 kyr resolution. We equate a GSL index value of 0.0 with the modern ice sheet volume,
and a GSL value of -125 m with the ice volume of the LGM. Prior to 800 ka, the emulator was run with
orbital parameters calculated in the same way as in the 800 kyr simulations2. However, CO and GSL
2
were from different sources: i) CO was taken from a record of model-derived CO 21 that is consistent
2 2
with global sea level reconstructions45; and ii) GSL was taken from a 1-D ice sheet model to reconstruct
sea level over the last 20 Myr45. Although there are more recent proxy-based estimates of CO 46, this
2
cannot be used in our study as we require an orbitally-resolved record through the entire 2.58 Myr that
captures that relationship between the CO feedback and the direct orbital forcing. Moreover, this other
2
CO record is not orbitally resolved prior to the ice core record46. Both CO and GSL were interpolated
2 2
to 1 kyr resolution here, and these two records were spliced with the existing 800 kyr records (see
above) during the overlapping 5 kyr (i.e. 795-800 ka).

Proxy data compared with output from the emulator
The palaeoclimate data used to compare with the emulated SAT and PREC were taken from several
different locations (Table 1, Figure S1), and were selected due to their high temporal resolution and
coverage of the Quaternary. For temperature, a quality control exercise was undertaken which resulted
in the elimination of any records with: i) a temporal resolution > 5 kyr, which may not capture orbital-
scale variability; ii) a starting point > 3 kyr BP, as the emulated anomalies were calculated using data
between 0-3 kyr; and iii) < 500 kyr worth of data during the Quaternary, to allow a reasonable model-
data comparison. This exercise yielded one record of reconstructed SAT, four records of reconstructed
SST and three records of reconstructed precipitation:
• Reconstructed Antarctic SAT, derived from ẟD data from the Dome C ice core42,43.
• Reconstructed SST data are based on a number of proxy methodologies, including the UK ’
37
index, Mg/Ca ratios from foraminiferal calcite, and the modern analogue technique (MAT),
listed in Table 1.
• δ18O data measured from cave speleothems, which is taken to be a first-order approximation for
intensification of East Asian monsoon rainfall. Records from three caves were used to provide
full coverage of the last 640 kyr, namely the Hulu Cave in eastern China47, the Sanbao Cave in
central China48, and the Dongge cave in southern China49.
Palaeo SST data is used for comparison with emulated SAT because multiple records exist from varying
global locations, which cover hundreds of thousands of years at a sufficient resolution to capture orbital
cycles. Outside regions of sea ice cover, the offset between SST and SAT at a single location is small
compared with the local glacial-interglacial cyclicity. Each proxy data set was compared to emulated
mean annual SAT or PREC for the appropriate grid box for each site. Temperature and precipitation
data are generally shown as anomalies compared to PI conditions at the different locations, taken as the
difference with the average of the most recent 3 kyr in the case of SST, and the PI control simulation
(i.e. 0 kyr) in the case of SAT and PREC.

Step 5: Sensitivity simulations
Once the above model-data validation process was complete, the final stage was to undertake a number
of sensitivity simulations. Firstly to investigate the sensitivity to the various forcings of the emulated
data at the proxy sites, through constant-driver simulations (Section 5.1), secondly to identify a
synthetic CO forcing that results in the best-fit of emulated climate to proxy data (Section 5.2), and
2
thirdly to separate out the dominant driving mechanisms from among the five input parameters, on a
global scale (Section 5.3).

(5.1) Sensitivity of emulated results at proxy locations
For the first investigation, a number of idealised sensitivity simulations, covering just the most recent
800 kyr, were carried out in order to investigate the sensitivity of the sites. To do this, instead of using
actual realistic values for all five driving parameters (used as input when running the emulator), these
were set as either (i) a constant average, (ii) constant PI, (iii) constant first year (because of missing
data, the first year is not always the PI), (iv) constant LGM and (v) randomly generated values.

(5.2) Sensitivity of emulated results to CO
2
For the second investigation, a synthetic timeseries of CO was generated, by initially running the
2
emulator multiple times using constant CO as input, to provide a range of temperature values that sit
2
above and below all possible proxy reconstructions. The other four input parameters were kept as
naturally/realistically varying. Specifically, the emulator was run seven times for the full 2.58 Myr,
with the CO held constant in each simulation as multiples of the PI (280.4 ppm): 0.125x, 0.25x, 0.5x,
2
1x, 2x, 4x and 8x (corresponding to 35.05, 70.1, 140.2, 280.4, 560.8, 1121.6 and 2243.2 ppm,
respectively). Therefore, in these simulations, any variability in the resulting emulated temperature is
associated with ice and/or orbital parameters only. Figure S4 in the Supplementary Material shows the
resulting temperature from these simulations at one example location (in the central South Pacific),
alongside the proxy data. These constant-CO simulations (Fig. S4a) were then used to estimate the
2
CO that would give the best fit of emulated temperature to proxy data. At every location and at every
2
timeslice, the proxy temperature was identified, as was the emulated temperature immediately above
and below and the constant CO values used to create these temperatures. These values were then used
2
to interpolate the exact CO (initially in log space, and then converted back to actual CO ) 
2 2
corresponding to each temperature. As a test, the emulator was then run again with these synthetic CO
2
curves, and the corresponding emulated temperatures were, as expected, identical to the proxy
temperatures on which they were calculated (see Fig. S4b). Finally, for each proxy record, the synthetic
CO was compared with CO reconstructions, including the observed42 and modelled21 CO used in the
2 2 2
model-data comparison, as well as two other CO records further back into the Quaternary: CO at ODP
2 2
999, spanning 1088-1244 kyr50 and CO at ODP 668, spanning 1384-1535 kyr51. All CO records were
2 2
firstly interpolated to 1 kyr resolution to enable comparison with the emulated temperature and
optimised CO (see Fig. S5). Note that although other CO records exist even further back in time, such
2 2
as for the Pliocene52, these were not included here because they are not orbitally resolved over multiple
successive glacial-interglacial cycles, and additionally no SAT proxy data (which passed the above
quality-control tests) were available for this time period.

(5.3) Sensitivity of emulated results to dominant driving input
For the third investigation, the emulator was run with a single forcing or subset of the forcings being
incrementally varied, following the multi-variate factorisation method23, whilst the other forcings were
held constant at PI values. For these, the nomenclature uses a binary system of 5 digits representing (in
order) CO , obliquity, eccentricity, precession and ice, where 0 = constant PI values and 1 = varying
2
values. Hence E (the no drivers simulation) has all drivers set to constant PI values, E has
00000 10000
varying CO but everything else is set to constant PI values, and so on, up to E where all forcings
2 11111
are varying (the all drivers simulation). This gave a total of 32 simulations, capturing all possible
combinations of forcings (see Table S2 in the Supplementary Material). To reiterate, the 32 simulations
refer to the sensitivity simulations using the emulator, not the HadCM3 simulations used to train the
emulator (of which there were 120 and 122 simulations for the interglacial and glacial emulators,
respectively, as discussed above). For the factorisation, only the most recent 800 kyr were considered,
as the input drivers (i.e. CO and GSL) are most reliable during this time period.
2

Metric for quantifying model performance
For the model-data comparison, idealised sensitivity simulations and factorisation, we generate a skill
score for temperature (SAT and SST) of the emulator compared with the proxy data, using the non-
dimensional arcsin Mielke measure, M22. This is the mean square error (MSE), non-dimensionalised by
the variance22. For emulated values x and proxy data y: 
M = (2 / π)arcsin{1 - mse / [𝑽 + 𝑽 + (𝑮 - 𝑮 )2]} * 1000
𝒙 𝒚 𝒙 𝒚
Eq. (1)
where V is the variance and G is the mean. For the model-data comparison, M is calculated between the
emulated (using the E simulation) and proxy data at the same site (such as ODP 1012 versus proxy
11111
ODP 1012) and at every possible combination of sites (such as emulated ODP 1012 versus proxy ODP
1123). For the idealised sensitivity simulations, M is calculated between the emulated (using the E
11111
simulation) and the modified (to be constant) proxy data at the same site. For the linear and linear-
sum/shared-interaction factorisation, M is calculated between the E simulation and all other
11111
simulations (i.e. treating the E simulation as observations). Whereas the model-data comparison
11111
figures presented here show anomalies (relative to the PI), the M scores were calculated using the
emulated and proxy absolute values.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Williams et al 2024 Nat Comm. Review\nThis study uses two statistical climate (temperature and precipitation) emulators a stand-in for a forward climate model in order to\ndo a simulation of the last 2.6 Million years while circumventing long integration times.\nOverall, I do like the concept of model emulators to facilitate researchers parsing through large ensemble sets, filling in the\nblanks between known forward model integrations, mining through large parameter spaces, etc. It is highly likely that this\ntechnique will become a vital supplemental technique to forward models in the future.\nHowever, although the application of this concept towards doing long paleoclimate simulations is somewhat novel (barring the\nheavily self-cited Lord et al 2017), the implementation here seems incomplete (could be done better).\nAlso—I do not understand why the authors (some of whom I SAW on the call on the PMIP WINGS seminar in fall 2023 actively\nlistening to Axel Timmermann present) did not cite the Yun et al 2023 paper out of that group using the CESM1.2 model to do\ntransient simulations of the exact same 3 Million years to present time period. https://cp.copernicus.org/articles/19/1951/2023/cp-\n19-1951-2023.pdf Its a pretty big glaring miss. And this is THE simulation to compare to now. (And maybe in another decade\nwe’ll have a few other models who have managed such a transient, but for now…)\nFirst, the methods presented in the main paper are totally insufficient to understand what is going on. You must dig into the\nsupplemental, or happen to catch in the discussion, to realize they are using HadCM3 – a CMIP3 class model – to do the\ntraining of their statistical emulator. (They should not use ‘GCM’ in the text—they should define HadCM3 is the GCM that they\nuse, as opposed to taking ALL the various flavors of PMIP GCMs and training on those GCMs, the general term whereas they\nuse HadCM3, a specific model.)\nThe emulator training should be described concisely in this paper so that it can stand alone. I found myself unable to read the\ntext without going back and reading Lord 2017. After reading Lord 2017, I am still a little confused at a few things. Why are there\nonly 4 parameters? What does it mean to do that with two distinct emulators? (Why isn’t ice sheet extent and amoc response a\n5th / n-th parameter?) Vegetation extent parameter? methane separate from co2 as some sort of proxy for veg extent perhaps?\nLake extent parameter? Why did they only train for a global composite number? Why not regional? If you’re going to train for a\nglobal composite number, why train on (?deep ocean?) temperature instead of training on what the archives actually are d18Oc\nor MgCa or UK37 or whatever. In the Lord 2017 paper, there is a thoughtful review exchange between Ganopolski and Crucifix.\n“In our case, as in most applications we have seen so far, the most important judge- ment is that the GCM response is smooth,\nbut it does not need to be linear. Another important judgement is that the GCM internal variability is Gaussian. ” But is this true.\nDo ice sheets respond smoothly? Does AMOC? When you look at the emulators being done by folks at MIT, this part of the\nworld they don’t do so well. (https://www.youtube.com/watch?v=PbcFWN5dtJc 35:59) This project of course is much more\nexpansive, but -BUT- I think it illustrates how tricky it IS to get AMOC and ice sheet responses right.\nNext, is the parameter space of the Lord 2017 paper – which was made to look at high CO2 environments applicable here? (see\nnext comment)\nHadCM3 is a fine model to use for making an emulator—however, this model is also at this point relatively inexpensive to use.\nWhen they name drop on page 7, they admit to do the actual simulation using the real forward model would only take 280 days.\nFor a paleoclimate simulation, 9 months of integration really doesn’t seem that bad. That seems pretty par for the course. For this\nstudy, they probably *should* go back and run many more simulations of HadCM3. The code repository says only 32 runs were\nused. That seems pretty tiny. And, the parameter space of the Lord 2017 paper, we later learn is not good enough, so they\nactually train another emulator for more ice climates. So this paper is really an amalgamation of two emulators. Which is a weird\n(not 100% sure defensible, especially as buried and not stress tested enough as this info is currently in the paper) place to be in\nif the times in between deglaciated and glaciated state is what they are trying to use the emulator to fill in.\nPlus, there have been deglacial HadCM3 runs ‘out there’ for probably 20 years. Why not use them more. And there IS Ice6G from the\nEemian onwards. Why NOT use lots of time slices there.\nFor the results pieces, the most interesting part – the Mid-Pleistocene-Transition doesn’t really show up. Its disappointing that\nthere are two end-state emulators and nothing in between. Isn’t the MPT all about oscillating in a stable way in the in-between\npart? I don’t quite follow why they are showing individual ODP cores when they are training to global numbers. Wouldn’t some\nsort of global stack be better (actually, I think they should at the minimum follow the Yun23 example and use regional stacks)?\nOr, alternatively, train to all the various ODP and speleothem and ice records that are available. And then do tons of\ntraining/testing period exercises to quantify how putting in/out the various archives influences the answer.\nIt feels supremely hand-wavy to say that variability seems a combination of “CO2 forcing and ice sheet feedbacks” because they\nhave emulated the former well while the latter is left in this weird zone between the two emulators.\nIn summary—I think this could be a great paper if the idea was fleshed out better, but I do not think my above concerns can be\ncorrected with minor revisions. Its probably 6-12 months’ work.\nBut the author should indeed persevere because this type of approach will become increasingly important, and it is a challenge\nto implement novel applications.", "reviewer_2": "Review of “The relative role of orbital forcing, CO2 and ice sheet feedbacks on Quaternary climate” by Williams et al.\nThe manuscript of Williams et al. shows a long simulation of glacial and interglacial cycles over the Quaternary. For that\npurpose, they use an emulator calibrated from a coupled model. They show that CO2 and ice sheet feedbacks are the main\ndriver of the temperature variations over those cycles and that the orbital forcing effect is minimal.\nOverall, the manuscript of Williams et al. is truly fascinating, well-written and it was a pleasure to review it. However, I do think\nsome of the most interesting features are the least highlighted. I also have a few comments to improve the clarity of the\nmanuscript which in some places is not so easy to understand, in particular for someone who is not into the field.\nMajor comments:\nL93-94-95 (and in the whole manuscript in general): I think there is a big lack of explanations regarding what an emulator is. I\nhave never heard of them before, and without looking at the methods, I still had absolutely no idea what it is and how it works by\nthe end of the paper. I understand that Nature’s papers have a structure such that methods will be after the core of the\nmanuscript, but I think a few lines need to be added in the main text to explain more in details what is an emulator. Does it have\nany kind of physics? Is it some kind of extrapolation method? Calling it just “a statistical model which estimate uncertainty” is just\nvague and means everything and nothing at the same time. It is very difficult for a reader to constantly jump back and forth\nbetween the text and the methods while keeping motivated to continue reading. For instance, some of the sentences around\nL387 were quite helpful for me to understand.\nL166-167-168 (and in the whole manuscript in general): This whole paragraph is crazy to believe and I found it to be the most\ninteresting part of the paper. While I understand the importance of simulating glacial/interglacial cycles, I could not say I was\nsurprised at all CO2 and ice sheet feedbacks are the main drivers. The main topic that amazed me in this paper is the use of a\nmethod, which seems like a massive step forward in climate modelling to simulate those cycles. After reading this, I was just\nasking myself “Why do we bother using GCMs then?”. I feel that the efficacy of this emulator and its potential use in any field in\npaleoclimate modelling is as important, if not more, than the other results of this paper. Yet, I think the authors did not discuss\nthis aspect so much. I feel the paper’s highlight is really this role of CO2 and ice sheets, while the authors basically showed that\nan emulator could replace a GCM in any kind of sensitivity tests and many complex situations, all of that while saving an\nenormous amount of resources. So on one hand I have the role of CO2 feedback in glacial cycles, and on the other hand a\ncritical tool that could be a revolution in paleoclimate modelling…\nL262: Since orbital forcing seems to contribute so little to temperature signal, does that imply constraining orbital forcing in any\nPMIP runs is pointless? PlioMIP3 is exploring sensitivity runs with different orbital values. Does that make any sense? Also, is\nthis result strongly model dependent?\nL322: Are there any kind of vegetation feedbacks in the emulator which could also influence the temperature signal?\nL333-334: I wish the authors would talk more on this point, because from what I understood from the rest of the paper, this aspect\ndoes not seem very difficult to test considering the performance of an emulator. I think one needs some kind of large ensembles,\nbut those large ensembles already exist. Is the issue that you need large ensembles in both glacial and interglacial conditions?\nMinor comments:\nL83: Too many commas? “A different, but complementary, approach is…”\nFig.1: I am not sure about Nature’s requirements regarding figure but if the figure was in PDF it would really improve it… The\ncaptions talk about “Light blue error bands” but with this poor figure quality it is impossible to see much.\nL152: Maybe a very brief explanation of what an M score is in this caption is required (I know it is in the method but it is the same\nissue as the lack of explanation regarding an emulator in the core of the manuscript)\nL206: Typo “underestimate the show the correct…”"}, "author_response": "Responses to reviewers following submission of manuscript\n“The relative role of orbital forcing, CO and ice sheet feedbacks on\nREVIEWER COMMENTS\nWe extend our sincere appreciation to both Reviewers for a thorough examination of our\nmanuscript. Here, we address their suggestions, comment-by-comment. In the below,\nreviewers’ comments are blue and in a smaller font, and our corresponding response follows\nin a standard font. Line numbers in the following refer to the manuscript.\nReviewer #1 (Remarks to the Author):\nWilliams et al 2024 Nat Comm. Review\nThis study uses two statistical climate (temperature and precipitation) emulators a stand-in for a forward climate\nmodel in order to do a simulation of the last 2.6 Million years while circumventing long integration times.\nOverall, I do like the concept of model emulators to facilitate researchers parsing through large ensemble sets,\nfilling in the blanks between known forward model integrations, mining through large parameter spaces, etc. It\nis highly likely that this technique will become a vital supplemental technique to forward models in the future.\nWe thank the reviewer for their comments.\nHowever, although the application of this concept towards doing long paleoclimate simulations is somewhat\nnovel (barring the heavily self-cited Lord et al 2017), the implementation here seems incomplete (could be done\nbetter).\nThe methods section, and associated description of the implementation, has been rewritten\n(see Methods, lines 453-814) and comments below.\nAlso—I do not understand why the authors (some of whom I SAW on the call on the PMIP WINGS seminar in\nfall 2023 actively listening to Axel Timmermann present) did not cite the Yun et al 2023 paper out of that group\nusing the CESM1.2 model to do transient simulations of the exact same 3 Million years to present time\nperiod. https://cp.copernicus.org/articles/19/1951/2023/cp-19-1951-2023.pdf Its a pretty big glaring miss. And\nthis is THE simulation to compare to now. (And maybe in another decade we’ll have a few other models who\nhave managed such a transient, but for now…)\nThis has been done, discussing the Yun et al. 2023 paper extensively in the Introduction (lines\n83-95), Results (lines 206-219) and Discussion (lines 313-330).\nFirst, the methods presented in the main paper are totally insufficient to understand what is going on. You must\ndig into the supplemental, or happen to catch in the discussion, to realize they are using HadCM3 – a CMIP3\nclass model – to do the training of their statistical emulator. (They should not use ‘GCM’ in the text—they\nshould define HadCM3 is the GCM that they use, as opposed to taking ALL the various flavors of PMIP GCMs\nand training on those GCMs, the general term whereas they use HadCM3, a specific model.)\nThe methods section has been rewritten (see Methods, lines 453-814), such that it is now\nunderstandable as a stand-alone document. Terminology and definitions have been refined\nthroughout the Methods section.\nThe emulator training should be described concisely in this paper so that it can stand alone.\nSee above response.\nstill a little confused at a few things. Why are there only 4 parameters? What does it mean to do that with two\ndistinct emulators? (Why isn’t ice sheet extent and amoc response a 5th / n-th parameter?) Vegetation extent\nparameter? methane separate from co2 as some sort of proxy for veg extent perhaps? Lake extent parameter?\nWhy did they only train for a global composite number? Why not regional? If you’re going to train for a\nglobal composite number, why train on (?deep ocean?) temperature instead of training on what the archives actually\nare d18Oc or MgCa or UK37 or whatever.\nFirstly, see above response about the methodology being rewritten, so that it is now a stand-\nalone document and understandable without reading Lord et al. (2017). This now explains\nthe 5 input parameters (see Methods, lines 477-483), and why there are two distinct emulators\n(lines 485-546). Secondly, as discussed in the Methods (lines 569-574), the parameters\nmentioned in the comment above (e.g. vegetation extent) are not included in the emulator but\nare accounted for by the underlying GCM simulations used to train the emulator. Thirdly, the\nemulator is not trained on a global composite number or regional data, but rather it is trained\non the spatial patterns from the GCM simulations i.e. global gridded data (lines 563-567).\nIn the Lord 2017 paper, there is a thoughtful review exchange between Ganopolski and Crucifix. “In our case, as\nin most applications we have seen so far, the most important judge- ment is that the GCM response is smooth,\nbut it does not need to be linear. Another important judgement is that the GCM internal variability is Gaussian. ”\nBut is this true. Do ice sheets respond smoothly? Does AMOC? When you look at the emulators being done by\nfolks at MIT, this part of the world they don’t do so well.\n(https://www.youtube.com/watch?v=PbcFWN5dtJc 35:59) This project of course is much more\nexpansive, but -BUT- I think it illustrates how tricky it IS to get AMOC and ice sheet responses right.\nThis question is addressed by the (rewritten) evaluation of the emulator, discussed in the\nMethods (lines 615-628). The training GCM data includes non-linearities such as changes in\nthe AMOC. The leave-one-out the analysis demonstrates that the emulator can capture these\nnon-linearities.\nNext, is the parameter space of the Lord 2017 paper – which was made to look at high CO2 environments\napplicable here? (see next comment)\nThe input parameter space has now been explained more clearly (see Methods, lines 477-\n483).\nHadCM3 is a fine model to use for making an emulator—however, this model is also at this point relatively\ninexpensive to use. When they name drop on page 7, they admit to do the actual simulation using the real\nforward model would only take 280 days. For a paleoclimate simulation, 9 months of integration really doesn’t\nseem that bad. That seems pretty par for the course.\nWe thank the reviewer for highlighting this error. This is now been corrected (see Results,\nlines 140-163), showing that actually it would take many decades to run HadCM3 for the\nequivalent amount of years.\nFor this study, they probably *should* go back and run many more simulations of HadCM3. The code\nrepository says only 32 runs were used. That seems pretty tiny. And, the parameter space of the Lord 2017\npaper, we later learn is not good enough, so they actually train another emulator for more ice climates. So this\npaper is really an amalgamation of two emulators. Which is a weird (not 100% sure defensible, especially as\nburied and not stress tested enough as this info is currently in the paper) place to be in if the times in between\ndeglaciated and glaciated state is what they are trying to use the emulator to fill in.\nRather than 32 runs being used (in the comment above), there were actually 182 GCM\nsimulations (of 500 years each, giving a total of 91 kyr) that were used to train the emulator.\nfactorisation, not the number of GCM simulations used to train/calibrate the emulator (see\nMethods, lines 781-794).\nPlus, there have been deglacial HadCM3 runs ‘out there’ for probably 20 years. Why not use them more. And\nthere IS Ice6G from the Eemian onwards. Why NOT use lots of time slices there.\nSee above comment. We acknowledge that there are more model data available but, as\ndiscussed in the Methods (lines 548-551), the GCM simulations used here, taken together,\ncapture the full range from interglacial states to full glacial conditions, and the simulations of\nSingarayer and Valdes (2010) that are used as part of the training data do include the\ndeglaciation.\nFor the results pieces, the most interesting part – the Mid-Pleistocene-Transition doesn’t really show up. Its\ndisappointing that there are two end-state emulators and nothing in between. Isn’t the MPT all about oscillating\nin a stable way in the in-between part? I don’t quite follow why they are showing individual ODP cores when\nthey are training to global numbers. Wouldn’t some sort of global stack be better (actually, I think they should at\nthe minimum follow the Yun23 example and use regional stacks)? Or, alternatively, train to all the various\nODP and speleothem and ice records that are available. And then do tons of training/testing period exercises to\nquantify how putting in/out the various archives influences the answer.\nFirstly, a discussion about why the MPT doesn’t obviously show up has now been included\n(see Discussion, lines 332-340). Secondly, we are not training the emulator on global\nnumbers but rather on globally gridded data i.e. capturing the spatial patterns from the GCM;\ntherefore a comparison of the emulator results with individual ODP cores is appropriate. This\nis discussed more clearly in the Methods (lines 563-567).\nIt feels supremely hand-wavy to say that variability seems a combination of “CO2 forcing and ice sheet\nfeedbacks” because they have emulated the former well while the latter is left in this weird zone between the\ntwo emulators.\nNeither CO nor ice sheets (as represented by global sea level) are produced by the emulator,\n2\nbut rather both are used as inputs to the emulator and are based on proxy data. Therefore,\ngiven that the final part of the Results section discusses the output of the emulator (i.e. based\non temperature only), and identifying which of the 5 input parameters are providing the\nlargest signal (via the multi-variate factorisation), we consider this observation to be\nappropriate.\nIn summary—I think this could be a great paper if the idea was fleshed out better, but I do not think my above\nconcerns can be corrected with minor revisions. Its probably 6-12 months’ work.\nBut the author should indeed persevere because this type of approach will become increasingly important, and it\nis a challenge to implement novel applications.\nWe thank the reviewer for their comments, and apologise for any confusion caused by a\ndisjointed and unclear Methods section. We very much hope that the methodology is now\nmuch clearer to understand, and without referencing multiple other sources.\nReviewer #1 (Remarks on code availability):\nThey definitely have code there. I didn't look through it closely, but it has comments and is 'clean' enough.\n(I'm a little surprised that R has the memory for this task.)\nWhen You open up one of the files, it has dependencies.These need to be packaged up better because I couldn't\nget the links to work. (Tarball all of these dependent, non public source R files)\n-->\nsource('C:/Users/cw18831/OneDrive - University of\nBristol/Documents/Research/SKB/SKB_Alan_Code_DJL/PosivaSKB/PosivaSKB/Emulator/2015_Bristol_5D_\nv001/R/col_parula.R')\nsource('C:/Users/cw18831/OneDrive - University of\nBristol/Documents/Research/SKB/SKB_Alan_Code_DJL/PosivaSKB/PosivaSKB/Emulator/2015_Bristol_5D_\nv001/R/col_bwr.R')\nsource('C:/Users/cw18831/OneDrive - University of\nBristol/Documents/Research/SKB/SKB_Alan_Code_DJL/PosivaSKB/PosivaSKB/Emulator/2015_Bristol_5D_\nv001/R/col_wr.R')\nThe dependencies mentioned above are not standard R routines, publicly available from any\nR repository (e.g. https://www.r-project.org/). All R routines that are personal to this work\nhave been uploaded.\nReviewer #2 (Remarks to the Author):\nReview of “The relative role of orbital forcing, CO2 and ice sheet feedbacks on Quaternary climate” by\nWilliams et al.\nThe manuscript of Williams et al. shows a long simulation of glacial and interglacial cycles over the Quaternary.\nFor that purpose, they use an emulator calibrated from a coupled model. They show that CO2 and ice sheet\nfeedbacks are the main driver of the temperature variations over those cycles and that the orbital forcing effect is\nminimal.\nOverall, the manuscript of Williams et al. is truly fascinating, well-written and it was a pleasure to review it.\nHowever, I do think some of the most interesting features are the least highlighted. I also have a few comments\nto improve the clarity of the manuscript which in some places is not so easy to understand, in particular for\nsomeone who is not into the field.\nWe thank the reviewer for their comments.\nMajor comments:\nL93-94-95 (and in the whole manuscript in general): I think there is a big lack of explanations regarding what an\nemulator is. I have never heard of them before, and without looking at the methods, I still had absolutely no idea\nwhat it is and how it works by the end of the paper. I understand that Nature’s papers have a structure such that\nmethods will be after the core of the manuscript, but I think a few lines need to be added in the main text to\nexplain more in details what is an emulator. Does it have any kind of physics? Is it some kind of extrapolation\nmethod? Calling it just “a statistical model which estimate uncertainty” is just vague and means everything and\nnothing at the same time. It is very difficult for a reader to constantly jump back and forth between the text and\nthe methods while keeping motivated to continue reading. For instance, some of the sentences around L387\nwere quite helpful for me to understand.\nThe methods section has been rewritten (see Methods, lines 453-814). In particular, to\nanswer the above questions about what exactly the emulator is and whether it contains any\nphysics, please see Methods (lines 569-574), in combination with the rest of that section.\nL166-167-168 (and in the whole manuscript in general): This whole paragraph is crazy to believe and I found it\nto be the most interesting part of the paper. While I understand the importance of simulating glacial/interglacial\ncycles, I could not say I was surprised at all CO2 and ice sheet feedbacks are the main drivers. The main topic\nthat amazed me in this paper is the use of a method, which seems like a massive step forward in climate\nmodelling to simulate those cycles. After reading this, I was just asking myself “Why do we bother using GCMs\nthen?”. I feel that the efficacy of this emulator and its potential use in any field in paleoclimate modelling is as\nimportant, if not more, than the other results of this paper. Yet, I think the authors did not discuss this aspect so\n4\nmuch. I feel the paper’s highlight is really this role of CO2 and ice sheets, while the authors basically showed\nthat an emulator could replace a GCM in any kind of sensitivity tests and many complex situations, all of that\nwhile saving an enormous amount of resources. So on one hand I have the role of CO2 feedback in glacial\ncycles, and on the other hand a critical tool that could be a revolution in paleoclimate modelling…\nSee above comment about the Methods rewrite, and in particular the Results (lines 140-163)\nwhere we discuss in more detail the efficiency of the emulator compared to a GCM.\nL262: Since orbital forcing seems to contribute so little to temperature signal, does that imply constraining\norbital forcing in any PMIP runs is pointless? PlioMIP3 is exploring sensitivity runs with different orbital\nvalues. Does that make any sense? Also, is this result strongly model dependent?\nSome text on the above question has now been included in the Discussion (lines 382-390).\nBut in short, no, constraining orbital forcing in PMIP runs is not pointless, because here we\nare looking at annual means only, for which orbital changes only have a small impact. If we\nwere to look at seasonal changes, a separate emulator would be needed and then the orbital\ncomponent would have more of an impact.\nL322: Are there any kind of vegetation feedbacks in the emulator which could also influence the temperature\nsignal?\nThis is now more clearly discussed in the Methods (lines 569-574). But in short , no,\nparameters such as vegetation are not included in the emulator but are accounted for by the\nunderlying GCM simulations used to train the emulator.\nL333-334: I wish the authors would talk more on this point, because from what I understood from the rest of the\npaper, this aspect does not seem very difficult to test considering the performance of an emulator. I think one\nneeds some kind of large ensembles, but those large ensembles already exist. Is the issue that you need large\nensembles in both glacial and interglacial conditions?\nThis has now been done as part of the methodology rewrite (see Methods, lines 453-814),\nin particular lines 485-551 where we discuss the different ensembles going into the glacial and\ninterglacial versions of the emulator.\nMinor comments:\nL83: Too many commas? “A different, but complementary, approach is…”\nThis has been addressed.\nFig.1: I am not sure about Nature’s requirements regarding figure but if the figure was in PDF it would really\nimprove it… The captions talk about “Light blue error bands” but with this poor figure quality it is impossible to\nsee much.\nThe caption here has been modified and the light blue error bands have been removed, as\nthey are not relevant to the study.\nL152: Maybe a very brief explanation of what an M score is in this caption is required (I know it is in the\nmethod but it is the same issue as the lack of explanation regarding an emulator in the core of the manuscript)\nThis has been addressed.\nL206: Typo “underestimate the show the correct…”\n6"}, {"round": 2, "reviewer_comments": {"reviewer_1": "Williams et al 2024 Nature\nSecond phase review\nThe manuscript is much improved from what was submitted previously.\nHowever, I reiterate here as I have to the editor that I am an emulator adjacent scientist. This manuscript absolutely must be\nevaluated by a researcher who actively uses emulators and is familiar with the quickly changing field of AI/ML in climate\nscience. Although I work with colleagues who are experts in emulators and feel very familiar with their use and applications, I do\nnot use them in my own research and thus do not feel qualified to comment on the particulars of the technique applied here.\nAs in my previous comments, there are many important feedbacks within the earth’s climate in modulating the climate response.\nThe PMIP/CMIP literature acknowledges this fact with its ‘spur’ experiments: lgm lgm+dust midholocene midholocene+veg (plus\nmethane lakes etc) … ismip experiment … volmip experiment … nahosmip … In the first paragraph of the intro, the things your\nstudy can address are listed, but the things you did not consider are not discussed. Be clear to set up the limitations and\nexpectations here and realistic about what the implication of these exclusions. These caveats belong up top, in the first\nparagraph or two, and you can even set up a nice resolution of the outstanding question- whether the exclusion meant you\nmissed out on important features in the transient. It is especially disappointing to have 3 of the forcings being the orbital\nparameters—which are then determined to be unimportant using the M score (which isn’t quite believable). And it isn’t clear to\nme that using the parameter coefficients of orbital forcing versus some derived quantity for insolation is correct. More details on\norbital forcing are required for such a ringer conclusion (contrary to pretty much all the paleoclimate literature out there) about\norbital forcing not being the driver and most important forcing. And as a set, since CO2 is how orbital forcing gets parsed through\nthe carbon cycle and ice sheets are how that gets parsed into the cryosphere, perhaps it IS better to use other derived quantities\nlike vegetation or dust or 1st derivative of sea level (ice melt) or similar as the forcings instead of the raw orbital parameters.\nThe methods section is much better now. But the main text readability isn’t that great. Too much in the weeds about certain\nthings without enough detail to follow the argument for others (details below). Stylistically, I think the main text needs a strong\nediting hand.\nI also did ask that the authors make plenty of comparisons to the Yun et al 2023 simulation which is the most similar to their\nefforts. They have done that—but that is not to say that this should be to the exclusion of other comparisons and efforts.\nThere is an implicit question through the text that needs explicit address:: hysteresis. That is to say, does the climate response\ndepend on the initial state. This question comes up in the discussion about the relative merits of the speeding-up technique used\nin the yun 2023 manuscript versus the emulator. Does the emulator include a phase delayed temperature (co2 et al) forcing field\nor anything else to include initial state? Or is it judging each climate state independently.\nIt STILL appears that there are two separate emulators (glacial and interglacial) that are stitched together. Is it any wonder that\nthe main finding of ice sheet is most important is arrived upon when its baked in to the process? I remain concerned, even after\nreading the responses and rewritten methods section, that there are two distinct emulators. Without this continuum, what chance\nis there to capture changes at the mid-pleistocene transition? Authors aren’t using a linear regression that might suffer from\nassigning warming of similar magnitude to cooling from G to IG states. That’s the point of this method. I am not convinced about\nthe use of two separate emulators. What happens if they use just one. L529 ‘steps were taken’ -> this is the methods section to\ndescribe those steps.\nI still remain concerned about emulating AMOC as per the first review. I do not see this addressed in the main text, not relegated\nto a mention deep-deep in the methods L542. AMOC remains one of the confounding factors for folks working on climate model\nemulators and it is odd to not explicitly address this piece.\nMajor revisions are still needed, although somewhat less major than last round. This emulator exercise is important for the\npaleoclimate community, although I am not convinced the present work is complete enough as presently written (and possibly\ndone for the forcing comments above).\nL86: This is where you should explicitly state that this yun23 type of simulation assumes that the model simulation is relatively\ninsensitive to initial conditions – you need to also discuss your own initial conditions. Was the training HadCM3 mindful of initial\nconditions? (Did you start your Pliocene simulations from an existing Miocene simulation or was it initialized from a default\nconfiguration.) If it was not, then you also have this same potential <lack of> initial condition bias built in.\nHysteresis (aka tipmip and the like) is a major point of concern for future comparisons. iow: I would frame this comparison a bit\ndifferently. It does nothing to degrade your work by having another group achieve a similar result—rather it gives an excellent\nbasis for comparison with strengths and weaknesses of your (their) approach, (in)consistencies, etc.\nL103-5: “simulate” => “interpolate”   this is what you are doing, you are not simulating. Emulators interpolate points determined\nby data or physics-based models, and totally dependent on the training data “points”. The word “emulator” denotes the relatively\nnew application of new AI/ML techniques instead of traditional to the statistical models reproduction—please introduce this\nconcept in this paragraph.\nL119: “proxy” => “past climate” using “proxy” as a noun becomes ambiguous as you cross fields.\nL120: “physics-based”, more complex GCM…\nL121: the number of simulations used in the training data is related to the number of parameters explored – use the specific\nlanguage. What did you use to properly sample the phase-space? Latin Hypercube? Something else? Introduce here. How\nmany simulations are required per parameter?\nL126-128: In reading this section, the word parameter feels muddied with the word parameter as it is used in perturbed physics\nensembles—those take on the PPE acronym. Perhaps there is a use-specific acronym? Boundary-Condition-Parameter (BCP)?\nThis technique has bing applied to sampling the parameter space used within the model physics (eg\ndoi:10.22541/essoar.172745119.96698579/v1), to interpolate between future earth scenarios (eg\nhttps://arxiv.org/html/2404.13227v2)\nL141: This stat of ‘5 minutes to run XX years’ will not age well and there is not enough info about the emulator in-line to know\nhow complex it is (Also- is the emulator full atmosphere 3d, ocean 3d?, land surface 3d?). Framing in comparison to the amount\nof time required to run the full model is sufficient.\nFurther, the emulator is solving each state as if it is totally independent (my take after reading the methods but this should be\nspelled out briefly in the text) — which is equivalent to doing time slices. Authors use 1kyr resolution, which is equivalent to\ndoing 2580 time slice simulations. The emulator was trained on 182, 500 year time slices. Since they provide 120 years / day as\nthe length scale of these simulations, these training simulations took <5 days. If the runs were done not overlapping, that is 29\nyears. Running the ensemble at a more typical few dozen at a time, that would be a run time on order of a year or so. In any\ncase, there is too much in the weeds on this run time justification in the main text and not enough details on the emulator in the\nmain text.\nL151: Does HadCM3 really output subdaily diagnostics and then average those into monthly averages? That is an interesting\nchoice. I would say most groups do not do that I/O step because I/O tends to slowdown the simulation A LOT because you are\nstopping and writing to the disk and depending on the filesystem that can be rather slow. Many groups keep up with tracers in\nmemory at the normal timestep (20-30 min), and then only output at the temporal resolution required, vastly decreasing I/O. Are\nthe authors positive that HadCM3 does not work that way? The next few sentences are too much in the weeds for main text and\nshould be moved to the methods or supplemental.\nL189: The discussion of figure 1 is too qualitative. Did you plot and calculate T_emulator versus T_proxy? The figure 1 plot here\ntakes up too much space and yet the curves are really too tiny and delicate to ‘tell’ what the authors are talking about.\nThere also needs to be a context of traditional ‘model’ built in. If you simply regress your forcing fields on (say) the 3d SST field,\nhow does this statistic compare with the emulator’s? T_regress versus T_proxy? For the d18Oice record, how do your results\nstack up against the EPICA Grand Challenge results (where simple and fast box models were applied to link CO2 and\ntemperature for Antarctic ice core).\nWith these other methods, do you have the same bias periods that you lay out on L193?\nOn L211: 4degC is FAR too little cooling for last glacial maximums. This magnitude should be set by data assimilation methods\neg https://doi.org/10.1038/s41586-021-03984-4 and in context of PMIP3/PMIP4 LGM models.\nL 204/701: You talk about the interpolated CO2 up to 800k in this study and then punt to the methods. Probably better to put\nconcisely enough in the text to follow what you’ve done (we interpolated 10k resolution to 1k resolution etc). You use an\noutdated source for CO2 prior to 800k. I recommend https://doi.org/10.1126/science.adi5177 which is the recent Dec2023\nCenCO2PIP consensus paper.\nL217: Are you saying that you have the same as the Yun paper or the same variability. It isn’t clear.\nL221: Again here, this would show up better in a scatter plot to augment figure 1 as the highest and lowest values having the\ngreatest error.\nL232: Can you say the interpretation something like “G-IG CO2 would need to be 1.5 times larger to match the proxy G-IG\ntemperatures (Figure S6)” instead of this hand waving that requires digging into the supplemental. After looking at Figure S6, I\nam not entirely sure what you’ve done. Do you not have actual CO2 in 800k to 1.1k? Are you comparing optimized CO2 to actual\nCO2? Can you put numbers on this?\nL239: OR your mismatch could be that you have not included enough forcing fields. Like vegetation or etc written above.\nL243: MAP is defined above, but it is an uncommon shorthand and makes the paper hard to read. SAT, in contrast, is a common\nacronym. Why not use PREC to improve readability?\nAlso, Chinese speleothem d18O are not MAP. They are changes in d18O that memorializes monsoon strength.\nL262: ‘all drivers’ simulation (E11111) Each subscript means forcing on or off? Which is which? Is this spelled out somewhere?\nAlthough this shorthand likely works great in the lab, it is hard to read.\nL265: as above, some sort of basic regression model comparison/context should be added.\nL269: I don’t follow what this means. Constant climate signal of what? Do you mean “equilibrium” climate experiment like PMIP?\nIf that’s what you mean, what were their scores?\nL289: You cannot punt folks to the supplemental text and then descend into unfamiliar jargon like “where each driving\ncomponent is averaged over all 120 pathways” I have no idea what you are talking about. The supplemental should AUGMENT\nthe main text, not render it unreadable. When you read in the supplemental the full suite of runs—decode them here otherwise\nreaders have to go back and dig into the methods (as per comment above about E11111), and then dig back into the original\npoint in the paper. This isn’t readable. Please add a supplemental table what E10000, E01000… E00001 mean.\nL293: If you get a result that the impact of the orbital forcing on M scores is negligible—perhaps this is an indicator that the orbital\n‘forcing’ the way you’ve implemented it isn’t quite right. It is the impact that these have on insolation seasonality that is important.\nAnd it seems weird to have eccentricity separate from precession as two separate forcings.\nWhat is the M score of the orbital parameters to CO2 and ice sheet volume? These need to be plugged in to a function to yield\ninsolation. Does the M score technique capture that?\nPerhaps you should use a variation of the original specmap methods (insolation at 65N) Perhaps max-min insolation at mid\nlatitudes of northern and southern hemisphere (or something to indicate seasonality magnitude) plus a baseline total insolation\nat (say) 65N?\nGiven that the emulator is so inexpensive to run, it would seem that you could experiment with other forcing parameters.\nOr perhaps this is an indication that the orbital forcing is already built in to the CO2 and ice volume signal so well that they are\nredundant? And maybe these 3 parameters should be swapped out with vegetation. Vegetation could be distilled down to the\nmidHolocene pmip variations (fraction of ‘green sahara’ and shrubs at high latitudes). Or maybe other things like aerosol loading\nmight be used as a forcing (as CO2, easier done 800k to present than before). Or perhaps use the first derivative of the sea level\n(something like amount of ice melt in NH and SH).\nIn short, do not think this paragraph L297 with the orbital parameters causing <5% of the signal is defensible as written and goes\nagainst 35 years of prevailing thought.\nL346: Why are there only modern and lgm ice? I know that there are HadCM3 deglacial experiments. I realize that I made this\ncomment before, and there is a response to reviewer, and yet I am not satisfied given that the lack of ice variability is one of the\nissues discussed in the discussion.\nL642: You should use BOTH the updated Hargreaves/Annan estimates (Renoult) generally on the low-end of accepted across\nthe community and the Osman/Tierney estimates generally on the higher end of accepted across the community from data\nassimilation methods here instead of this older work.", "reviewer_2": "No more comments (all comments have been answered in the new version of the paper)", "reviewer_3": "The idea of a paleoclimate emulator is a good one. These simulations are so long and expensive that we already (often) use\npared down climate models, and this is taking it another step to pare things down even further, which is a worthy goal. The\nemulator seems to be well constructed and carefully validated. In general, I really like this manuscript. However, I have a\nquestion about whether it fits the scope of the journal.\nWhere I struggle with this paper is that it seems to be focused on describing and validating the emulator’s performance, which is\nwhat I might expect for a more methods-focused journal. What I think is lacking are insights that the emulator can provide.\nExamples could include figuring out what processes are important to represent in models for understanding paleoclimate\nvariability or revealing drivers of past changes. The authors do point out that CO2 concentration and ice are more important than\nthe inputs corresponding to orbital parameters, but given the way the results are presented and caveated, I can’t tell whether\nthose results are new or surprising, or even whether they are simply a product of what was done.\nMy hope is that these issues can be taken care of via wording and do not require much new analysis. The paper is quite well\ndone, and the tool the authors have produced is quite valuable.\nThe authors also did an excellent job of responding to the previous review comments, some of which were quite challenging."}, "author_response": "Responses to reviewers following submission of manuscript\n“The relative role of direct orbital forcing, and CO and ice sheet feedbacks\n2\non Quaternary climate” by Williams et al.\nREVIEWER COMMENTS\nWe extend our sincere appreciation to all Reviewers for their examination of our manuscript.\nHere, we address their suggestions, comment-by-comment. In the below, reviewers’\ncomments are blue and in a smaller font, and our corresponding response follows in a\nstandard font. Line numbers in the following reviewer comments refer to the previously-\nsubmitted version of the manuscript, whereas line numbers in our responses refer to the new\nversion of the manuscript. We have also have numbered the comments below and provided\nreferences to these as comments next to the text (main manuscript and supplementary\nmaterial) to aid readability.\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nWilliams et al 2024 Nature\nSecond phase review\nThe manuscript is much improved from what was submitted previously.\n0. We thank the reviewer for their comments.\nHowever, I reiterate here as I have to the editor that I am an emulator adjacent scientist. This manuscript\nabsolutely must be evaluated by a researcher who actively uses emulators and is familiar with the quickly\nchanging field of AI/ML in climate science. Although I work with colleagues who are experts in emulators and\nfeel very familiar with their use and applications, I do not use them in my own research and thus do not feel\nqualified to comment on the particulars of the technique applied here.\n0. A third reviewer was appointed, and their comments addressed (see below). The second\nreviewer was satisfied with our previous responses (see below).\nAs in my previous comments, there are many important feedbacks within the earth’s climate in modulating the\nclimate response. The PMIP/CMIP literature acknowledges this fact with its ‘spur’ experiments: lgm lgm+dust\nmidholocene midholocene+veg (plus methane lakes etc) … ismip experiment … volmip experiment …\nnahosmip … In the first paragraph of the intro, the things your study can address are listed, but the things you\ndid not consider are not discussed. Be clear to set up the limitations and expectations here and realistic about\nwhat the implication of these exclusions. These caveats belong up top, in the first paragraph or two, and you can\neven set up a nice resolution of the outstanding question- whether the exclusion meant you missed out on\nimportant features in the transient. It is especially disappointing to have 3 of the forcings being the orbital\nparameters—which are then determined to be unimportant using the M score (which isn’t quite believable). And\nit isn’t clear to me that using the parameter coefficients of orbital forcing versus some derived quantity for\ninsolation is correct. More details on orbital forcing are required for such a ringer conclusion (contrary to pretty\nmuch all the paleoclimate literature out there) about orbital forcing not being the driver and most important\nforcing. And as a set, since CO2 is how orbital forcing gets parsed through the carbon cycle and ice sheets are\nhow that gets parsed into the cryosphere, perhaps it IS better to use other derived quantities like vegetation or\ndust or 1st derivative of sea level (ice melt) or similar as the forcings instead of the raw orbital parameters.\n1. Firstly, a sentence has been added towards the end of the introduction (lines 135-141),\ndetailing which feedbacks/processes are included in the HadCM3 model (and hence included\nin the emulator), and which are not included.\n2. Secondly, we have clarified what we mean when we say that the CO and ice feedbacks\n2\nare the dominant drivers (rather than the orbital forcings), both in the Abstract (lines 41-45)\nand the Discussion (lines 403-418). To clarify, we are not saying that the net overall effect of\nthree orbital forcings are unimportant – indeed, because CO , ice and other feedbacks are\n2\nultimately driven by these orbital forcings, the net overall orbital forcing is, of course, very\nimportant. However here we are looking at the direct effect of the orbital forcing, which,\ndespite being very large seasonally, is relatively small on the annual mean, resulting in a\nrelatively small impact on annual mean temperatures.\nThe methods section is much better now. But the main text readability isn’t that great. Too much in the weeds\nabout certain things without enough detail to follow the argument for others (details below). Stylistically, I think\nthe main text needs a strong editing hand.\nWe have made a general edit of all the main text, including throughout the Introduction,\nResults, Discussion and Conclusions.\nI also did ask that the authors make plenty of comparisons to the Yun et al 2023 simulation which is the most\nsimilar to their efforts. They have done that—but that is not to say that this should be to the exclusion of other\ncomparisons and efforts.\n3. We have added in more studies undertaking similar comparisons and efforts, in the\nIntroduction (lines 67-71 and 93-95).\nThere is an implicit question through the text that needs explicit address:: hysteresis. That is to say, does the\nclimate response depend on the initial state. This question comes up in the discussion about the relative merits of\nthe speeding-up technique used in the yun 2023 manuscript versus the emulator. Does the emulator include a\nphase delayed temperature (co2 et al) forcing field or anything else to include initial state? Or is it judging each\nclimate state independently.\n4. Here, we judge each climate state independently, therefore the issue of hysteresis is not\nconsidered. We do, however, acknowledge that this is an interesting aspect, and could be\nexplored in future work. A sentence has been added in the Discussion to highlight this\nlimitation (lines 435-439).\nIt STILL appears that there are two separate emulators (glacial and interglacial) that are stitched together. Is it\nany wonder that the main finding of ice sheet is most important is arrived upon when its baked in to the process?\nI remain concerned, even after reading the responses and rewritten methods section, that there are two distinct\nemulators. Without this continuum, what chance is there to capture changes at the mid-pleistocene transition?\nAuthors aren’t using a linear regression that might suffer from assigning warming of similar magnitude to\ncooling from G to IG states. That’s the point of this method. I am not convinced about the use of two separate\nemulators. What happens if they use just one. L529 ‘steps were taken’ -> this is the methods section to describe\nthose steps.\n5. As we explain in the Methods (lines 576-593), we were concerned about using a single\nemulator due to possible unphysical ‘leakage’ effects from the glacial to the interglacial,\nresulting in unphysical warming in the region of the Laurentide, Fennoscandian, and\nAntarctic ice sheets during interglacials. In order to test this, we created a combined\nemulator, i.e. the glacial and interglacial conditions were considered as a single emulator.\nThis combined emulator is illustrated in the following figure (now also included in the\nSupplementary Material, Figure S1e), which shows each of the HadCM3 training simulations\naccording to the relative sea level, log CO and insolation at 65°N (combining the three\n2\norbital parameters into a single parameter, just for the purposes of the figure). Each training\nsimulation, represented by the dots, has been categorised according to whether relative sea\nlevel is at (red), above (yellow), or below (blue) 0 m.\nWe repeated the ‘all drivers’ simulation (i.e. where all five drivers varied according to proxy\ndata) but using this combined emulator. This can be seen in the figure below.\nWhen compared to Figure 1 in the manuscript, there is little difference at the SST locations\n(panel b-e) but substantially more uncertainty in SAT over Antarctica (panel a) when a\ncombined emulator is used.\nWe therefore consider it appropriate to use the two separate emulators, as combining them\ninto a single emulator does indeed lead to higher uncertainties at the boundary between\nglacials and interglacials, likely due to the ‘leakage’ effect discussed in the text. All of this is\nnow discussed in the Methods (lines 576-593).\nI still remain concerned about emulating AMOC as per the first review. I do not see this addressed in the main\ntext, not relegated to a mention deep-deep in the methods L542. AMOC remains one of the confounding factors\nfor folks working on climate model emulators and it is odd to not explicitly address this piece.\n6. The emulator is a purely statistical model designed to replicate the output of HadCM3, and\nas such does not explicitly include any representation of any physical processes, including the\nAMOC (or indeed any other process, e.g. atmospheric circulation, clouds, se up a ice etc.).\nHowever, because physical processes are represented in HadCM3, the output of the emulator\ndoes implicitly account for these processes, including the AMOC. This has been made more\nclear in the Introduction (lines 141-145).\nMajor revisions are still needed, although somewhat less major than last round. This emulator exercise is\nimportant for the paleoclimate community, although I am not convinced the present work is complete enough as\npresently written (and possibly done for the forcing comments above).\n7. Revisions have now been made throughout the main text, both according to comments\nhere (and Reviewer 3) and comprising more general edits.\nL86: This is where you should explicitly state that this yun23 type of simulation assumes that the model\nsimulation is relatively insensitive to initial conditions – you need to also discuss your own initial conditions.\nWas the training HadCM3 mindful of initial conditions? (Did you start your Pliocene simulations from an\nexisting Miocene simulation or was it initialized from a default configuration.) If it was not, then you also have\nthis same potential <lack of> initial condition bias built in.\n8. All of the HadCM3 training simulations were initialised from preindustrial conditions and\nthen spun-up to quasi-equilibrium; this means that the training data can considered\nindependent of the initial conditions of the underlying GCM simulations. This is now\ndiscussed in the Methods (lines 556-560). The glacial-interglacial timeseries presented in the\npaper are also independent of the initial emulated simulation in the timeseries, because each\nemulated time period within the timeseries is an independent snapshot, stitched together to\nmake a timeseries.\nHysteresis (aka tipmip and the like) is a major point of concern for future comparisons. iow: I would frame this\ncomparison a bit differently. It does nothing to degrade your work by having another group achieve a similar\nresult—rather it gives an excellent basis for comparison with strengths and weaknesses of your (their) approach,\n(in)consistencies, etc.\n9. This has already been addressed by the above comment (point 4) and in the Discussion\n(lines 435-439), where we discuss in the main text that hysteresis is not considered here.\nL103-5: “simulate” => “interpolate”  this is what you are doing, you are not simulating. Emulators interpolate\npoints determined by data or physics-based models, and totally dependent on the training data “points”. The\nword “emulator” denotes the relatively new application of new AI/ML techniques instead of traditional to the\nstatistical models reproduction—please introduce this concept in this paragraph.\n10. A sentence has been added in the Introduction to clarify the above point (lines 116-118).\nWe have also been through the manuscript and ensured that when we are talking about the\nresults from the emulator, we are clear on this e.g. by using ‘emulated’ instead of ‘simulated’.\nL119: “proxy” => “past climate” using “proxy” as a noun becomes ambiguous as you cross fields.\n11. This has been removed as part of other edits (line 129).\nL120: “physics-based”, more complex GCM…\n12. This has been added (line 113).\nL121: the number of simulations used in the training data is related to the number of parameters explored – use\nthe specific language. What did you use to properly sample the phase-space? Latin Hypercube? Something else?\nIntroduce here. How many simulations are required per parameter?\n13. This paragraph has now been moved into the Methods (lines 486-506), in line with\nReviewer 3’s comments that the main text should focus more on the results. But we have\nalso added to this sentence, to clarify that Latin hypercube sampling was indeed used and is\ndiscussed in more detail below (line 489-490).\nL126-128: In reading this section, the word parameter feels muddied with the word parameter as it is used in\nperturbed physics ensembles—those take on the PPE acronym. Perhaps there is a use-specific acronym?\nBoundary-Condition-Parameter (BCP)? This technique has bing applied to sampling the parameter space used\nwithin the model physics (eg doi:10.22541/essoar.172745119.96698579/v1), to interpolate between future\nearth scenarios (eg https://arxiv.org/html/2404.13227v2)\n14. This sentence has now been moved into the Methods (lines 492, 496), in line with\nReviewer 3’s comments that the main text should focus more on the results. But we have\nalso removed the word ‘parameter’ from this line, to avoid confusion (line 496). We have\nalso been through the manuscript and ensured we are consistent with our terminology, such as\nthe use of ‘driver’ versus ‘forcing’.\nL141: This stat of ‘5 minutes to run XX years’ will not age well and there is not enough info about the emulator\nin-line to know how complex it is (Also- is the emulator full atmosphere 3d, ocean 3d?, land surface 3d?).\nFraming in comparison to the amount of time required to run the full model is sufficient.\nFurther, the emulator is solving each state as if it is totally independent (my take after reading the methods but\nthis should be spelled out briefly in the text) — which is equivalent to doing time slices. Authors use 1kyr\nresolution, which is equivalent to doing 2580 time slice simulations. The emulator was trained on 182, 500 year\ntime slices. Since they provide 120 years / day as the length scale of these simulations, these training\nsimulations took <5 days. If the runs were done not overlapping, that is 29 years. Running the ensemble at a\nmore typical few dozen at a time, that would be a run time on order of a year or so. In any case, there is too\nmuch in the weeds on this run time justification in the main text and not enough details on the emulator in the\nmain text.\n15. This paragraph has now also been moved into the Methods (lines 508-524), in line with\nReviewer 3’s comments that the main text should focus more on the results. But we have\nalso rephrased this sentence and taken out any reference to an exact number of minutes (lines\n508-524).\nL151: Does HadCM3 really output subdaily diagnostics and then average those into monthly averages? That is\nan interesting choice. I would say most groups do not do that I/O step because I/O tends to slowdown the\nsimulation A LOT because you are stopping and writing to the disk and depending on the filesystem that can be\nrather slow. Many groups keep up with tracers in memory at the normal timestep (20-30 min), and then only\noutput at the temporal resolution required, vastly decreasing I/O. Are the authors positive that HadCM3 does not\nwork that way? The next few sentences are too much in the weeds for main text and should be moved to the\nmethods or supplemental.\n16. This sentence been removed, and the whole paragraph has been moved into the Methods\n(lines 508-524). HadCM3, as with all GCMs, makes calculations and updates the system\nstate every model timestep (15 minutes in the atmosphere). However (also like all GCMs), it\ndoes not output all these 15-minute states to disk. Instead it stores these in memory and then\ncalculates the monthly-mean at the end of the month, and outputs this to disk.\nL189: The discussion of figure 1 is too qualitative. Did you plot and calculate T_emulator versus T_proxy? The\nfigure 1 plot here takes up too much space and yet the curves are really too tiny and delicate to ‘tell’ what the\nauthors are talking about.\n17a. Firstly, we have increased the size of Figure 1, and we now feel the curves are easier to\nsee. Secondly, we have added a sentence (line 184) in the Results pointing the reader to the\nquantitative discussion of this figure, further down in the results when the M scores are\nintroduced. Finally, we have added a figure into the Supplementary Material (Figure S0),\nshowing a scatterplot of emulated temperature versus proxy temperature for the same five\nlocations as Figure 1, and discussed this in the Results (lines 187-190).\nThere also needs to be a context of traditional ‘model’ built in. If you simply regress your forcing fields on (say)\nthe 3d SST field, how does this statistic compare with the emulator’s? T_regress versus T_proxy? For the\nd18Oice record, how do your results stack up against the EPICA Grand Challenge results (where simple and fast\nbox models were applied to link CO2 and temperature for Antarctic ice core).\n17b. This sounds interesting, but correlating just one of the forcings (e.g. CO ) with the site-\n2\nspecific proxy temperature is too generous to the emulator, because the emulator ‘sees’ all 5\nforcings (CO , sea-level and the 3 orbital parameters). Doing something more complex, in\n2\nwhich for example we build a multi(5)-variate regression between the forcings and the site-\nspecific temperatures (as was done in some of the responses to the EPICA Challenge; Wolff\net al. 2011), is outside the scope of this paper and, furthermore, such a regression-based\nmodel is contrary to the ‘philosophy’ of our study, in which the underlying GCM is entirely\nmechanistic, and proxy temperature observations play no part in the model development or\ntuning; the model-data comparison in our study is a purely emergent property.\nWith these other methods, do you have the same bias periods that you lay out on L193?\n17c. The scatterplot now shown in the Supplementary material (see above comment) also\nidentifies the three MISs discussed in the Results, where the same biases are evident i.e.\nduring these periods, the proxy data are generally warmer than the emulated temperatures.\nThis is now discussed in the Results (lines 193-195).\nOn L211: 4degC is FAR too little cooling for last glacial maximums. This magnitude should be set by data\nassimilation methods eg https://doi.org/10.1038/s41586-021-03984-4 and in context of PMIP3/PMIP4 LGM\nmodels.\n18. This comparison been removed, as it is merely demonstrating the observation above (line\n209). Global means from the various LGM datasets are discussed in the Supplementary\nMaterial (lines 86-90) and are as follows (all values show LGM cooling relative to the PI):\nemulator (this paper) = 5.1°C, HadCM3 (this paper) = 5.3°C, Annan and Hargreaves (2013) =\n4°C and Osman et al. (2021) = 6-7°C.\nL 204/701: You talk about the interpolated CO2 up to 800k in this study and then punt to the methods. Probably\nbetter to put concisely enough in the text to follow what you’ve done (we interpolated 10k resolution to 1k\nresolution etc). You use an outdated source for CO2 prior to 800k. I\nrecommend https://doi.org/10.1126/science.adi5177 which is the recent Dec2023 CenCO2PIP consensus paper.\n19. A sentence has been added in the Methods (lines 711-717) to include the above reference,\nto justify why we did not use it and instead used our existing CO source (Van de Wal et al.\n2\n2011), and to clarify that we have interpolated the data to 1 kyr.\nL217: Are you saying that you have the same as the Yun paper or the same variability. It isn’t clear.\n20. This has been corrected (line 216).\nL221: Again here, this would show up better in a scatter plot to augment figure 1 as the highest and lowest\nvalues having the greatest error.\n21. This has been corrected by the comment above.\nL232: Can you say the interpretation something like “G-IG CO2 would need to be 1.5 times larger to match the\nproxy G-IG temperatures (Figure S6)” instead of this hand waving that requires digging into the supplemental.\nAfter looking at Figure S6, I am not entirely sure what you’ve done. Do you not have actual CO2 in 800k to\n1.1k? Are you comparing optimized CO2 to actual CO2? Can you put numbers on this?\n22. This sentence has now been modified to incorporate the above (lines 231-234).\nL239: OR your mismatch could be that you have not included enough forcing fields. Like vegetation or etc\nwritten above.\n23. This sentence has now been modified to incorporate the above (lines 244-245).\nL243: MAP is defined above, but it is an uncommon shorthand and makes the paper hard to read. SAT, in\ncontrast, is a common acronym. Why not use PREC to improve readability?\nAlso, Chinese speleothem d18O are not MAP. They are changes in d18O that memorializes monsoon strength.\n24. This has now been done (line 247 and throughout the rest of the manuscript), and a\nsentence has been added to clarify that Chinese speleothem d18O are not mean annual\nprecipitation but that we are using it as a first-order comparison here (line 247-249).\nL262: ‘all drivers’ simulation (E11111) Each subscript means forcing on or off? Which is which? Is this spelled\nout somewhere? Although this shorthand likely works great in the lab, it is hard to read.\n25. This has now been clarified when this simulation is first mentioned at the beginning of\nthe Results (lines 162-163).\nL265: as above, some sort of basic regression model comparison/context should be added.\n26. This has been corrected by the comment above.\nL269: I don’t follow what this means. Constant climate signal of what? Do you mean “equilibrium” climate\nexperiment like PMIP? If that’s what you mean, what were their scores?\n27. This has now been clarified, directly in the Results (lines 275-276) and also linking to the\nrelevant part of the Methods where this is discussed (lines 757-763).\nL289: You cannot punt folks to the supplemental text and then descend into unfamiliar jargon like “where each\ndriving component is averaged over all 120 pathways” I have no idea what you are talking about. The\nsupplemental should AUGMENT the main text, not render it unreadable. When you read in the supplemental the\nfull suite of runs—decode them here otherwise readers have to go back and dig into the methods (as per\ncomment above about E11111), and then dig back into the original point in the paper. This isn’t readable. Please\nadd a supplemental table what E10000, E01000… E00001 mean.\n28. This sentence has now been removed from the Results (lines 295-300) and a table has\nbeen added to the Supplementary Material (Table S3) and pointed to in the Methods (line\n803).\nL293: If you get a result that the impact of the orbital forcing on M scores is negligible—perhaps this is an\nindicator that the orbital ‘forcing’ the way you’ve implemented it isn’t quite right. It is the impact that these have\non insolation seasonality that is important. And it seems weird to have eccentricity separate from precession as\ntwo separate forcings.\n29. We have added a sentence in the Results (lines 304-307), firstly to clarify that\neccentricity and precession are not two separate forcings, and secondly to reiterate (as already\ndiscussed in the Discussion, see above comment) that we are not suggesting that the overall\nnet orbital forcing is negligible, but that the direct orbital radiative forcing on the annual\nmean temperature response is very small. This result is in line with previous work (see\nDiscussion, lines 403-406).\nWhat is the M score of the orbital parameters to CO2 and ice sheet volume? These need to be plugged in to a\nfunction to yield insolation. Does the M score technique capture that?\n30. The M score measures the goodness-of-fit of emulated temperatures (in °C) relative to\nobservations of temperature (in °C). It is not possible to construct an M score of two forcings\n(with different units), relative to each other. We have added a sentence to clarify this in the\nMethods (Metric for quantifying model performance, lines 813-814).\nPerhaps you should use a variation of the original specmap methods (insolation at 65N) Perhaps max-min\ninsolation at mid latitudes of northern and southern hemisphere (or something to indicate seasonality magnitude)\nplus a baseline total insolation at (say) 65N?\nGiven that the emulator is so inexpensive to run, it would seem that you could experiment with other forcing\nparameters.\nOr perhaps this is an indication that the orbital forcing is already built in to the CO2 and ice volume signal so\nwell that they are redundant? And maybe these 3 parameters should be swapped out with vegetation. Vegetation\ncould be distilled down to the midHolocene pmip variations (fraction of ‘green sahara’ and shrubs at high\nlatitudes). Or maybe other things like aerosol loading might be used as a forcing (as CO2, easier done 800k to\npresent than before). Or perhaps use the first derivative of the sea level (something like amount of ice melt in\nNH and SH).\n31. Addressing all 3 of the above comments together, the orbital forcing is calculated from\nLaskar et al. (2004), implemented into the full HadCM3 training-data simulations. This\nincludes the complete seasonal and latitudinal evolution of insolation, giving us the direct\norbital forcing. As the reviewer states, and as we now state in the Results (lines 309-313), the\nindirect impact of orbital forcing is implicitly included in the CO\n2\nand ice sheet forcings.\nIn short, do not think this paragraph L297 with the orbital parameters causing <5% of the signal is defensible as\nwritten and goes against 35 years of prevailing thought.\n32. In our view, this statement does not go against 35 years of thought. It is very well\nestablished, and consistent with our results, that direct orbital radiative forcing leads to\nfeedbacks associated with CO and ice sheets, which in turn result in glacial-interglacial\n2\ncycling. However, the direct radiative impact of orbital forcing itself (which is very close to\nzero on an annual-global mean), and in the absence of CO and/or ice sheet feedbacks, has\n2\nrelatively little effect on temperature (in particular, on annual-mean temperature). This has\npreviously been shown in simulations of specific time periods, such as during the LGM (e.g.\nCao et al. 2019, Shi et al. 2023), but here we show it for the entire 2.58 Myr. As per the\nreviewer’s first comment, this is now discussed in the Discussion (lines 401-418).\nL346: Why are there only modern and lgm ice? I know that there are HadCM3 deglacial experiments. I realize\nthat I made this comment before, and there is a response to reviewer, and yet I am not satisfied given that the\nlack of ice variability is one of the issues discussed in the discussion.\n33. There is not only LGM and modern ice in the HadCM3 training simulations. In addition\nto modern ice (the 60-member modice ensemble), the training data includes the glacial-\ninterglacial cycle simulations of Singarayer and Valdes (2010) (the 122-member “highice”\nensemble), and also Pliocene-like simulations (the 60-member “lowice” ensemble). A\nsentence to clarify this has been added in the Methods (Creating the interglacial and glacial\ntraining data, lines 571-574) – see also the figure of the combined emulator above which\nshows this clearly.\nL642: You should use BOTH the updated Hargreaves/Annan estimates (Renoult) generally on the low-end of\naccepted across the community and the Osman/Tierney estimates generally on the higher end of accepted across\nthe community from data assimilation methods here instead of this older work.\n34. This figure has now been modified, so that it includes the more recent Osman et al.\n(2021) data assimilation method, in addition to the other 3 panels (emulator, HadCM3 and\nAnnan and Hargreaves (2013)). We have removed this discussion from the Methods (lines\n663-669) into the Supplementary Material (specifically Section S4 (SM, lines 68-96),\nwhere we discuss this addition alongside the other LGM panels.\nReviewer #1 (Remarks on code availability):\nNo I didn't look through the code this time.\n35. We believe we have satisfied all of the comments concerning code, from both the\nreviewer and the editor.\nReviewer #2 (Remarks to the Author):\nNo more comments (all comments have been answered in the new version of the paper)\n1. We thank the reviewer for their previous round of comments.\nReviewer #3 (Remarks to the Author):\nThe idea of a paleoclimate emulator is a good one. These simulations are so long and expensive that we already\n(often) use pared down climate models, and this is taking it another step to pare things down even further, which\nis a worthy goal. The emulator seems to be well constructed and carefully validated. In general, I really like this\nmanuscript. However, I have a question about whether it fits the scope of the journal.\n2. We thank the reviewer for their comments.\nWhere I struggle with this paper is that it seems to be focused on describing and validating the emulator’s\nperformance, which is what I might expect for a more methods-focused journal. What I think is lacking are\ninsights that the emulator can provide.\n3. We have modified a number of areas to address this comment. Firstly, we have rewritten\nthe abstract (lines 33-45) so that it focuses more on the results, and less on the methodology.\nSecondly we have likewise modified the Introduction (lines 64-71 and 129-155), and the\nConclusions (lines 446-447 and 469). We have also removed some of the methodology\npreviously in the Introduction, summarising this here (lines 130-155) and moving the full\ndescription to the Methods (lines 486-506). Finally, we have moved the description of the\nemulator’s efficiency out of the Results and into the Methods (lines 508-524).\nMy hope is that these issues can be taken care of via wording and do not require much new analysis. The paper\nis quite well done, and the tool the authors have produced is quite valuable.\nWe hope that all of the above changes mean that the main part of the manuscript now focuses\nless on the methodology and emulator validation, and more on the primary research question\nof the driving mechanisms of Quaternary climate change.\nThe authors also did an excellent job of responding to the previous review comments, some of which were quite\nchallenging.\nWe thank the reviewer for their comments.\nReviewer #3 (Remarks on code availability):\nThe other reviewers focused on the code review, so I didn't do this part.\nWe believe we have satisfied all of the comments concerning code, from both the previous\nreviewers and the editor."}, {"round": 3, "reviewer_comments": {"reviewer_1": "Review of Williams et al 2025\nThis paper is greatly improved from previous versions. I do not think I will need to see this again before the editor makes a decision.\nI’m still not satisfied with two stitched together emulators being used instead of one. It makes the very interesting and good work\nseem less so because the answer gets baked in to the training.\n“The approach of having two separate emulators was taken to ensure that there was no ‘leakage’\nof state-specific climate variations across different climate states (see Supplementary Material,\nFigure S1e). For example, under glacial conditions development of the Laurentide ice sheet is\naccompanied by strong cooling over northern North America. However, this does not mean that\nduring interglacial conditions there should be a strong warming in the same region, since the\ndata are presented as an anomaly from interglacial present-day conditions (where no ice sheet is\npresent). ”\nActually, the change in height alone does mean there is rightfully a huge amount of warming in this area. In anomaly space this\ndoes look odd—it just does. If you want to get rid of that, try using vertical lapse rate to apply a correction.\nIf the authors used a single emulator – show it. It makes no sense the explanation that there is leakage over North America and\nthen differences between the two emerge on the opposite side of the world “When compared to proxy data at the five sites\npresented in Figure 1, the results\nshowed no difference in SST but higher variability for SAT over Antarctica (not shown). ”. Show both emulators.\nThe “M score” is not described adequately for readability in the text. In the main text, just a few short words like “M-score (Arcsin-\nMielke measure M) quantifies the goodness of fit of the emulator to the proxy climate archives such that where no skill is 0 and\nmaximum skill is 1000 (see methods)” will allow the reader to keep reading without “flipping” to the back page.\nThe authors insist on using the 2013 Annan and Hargreaves extremely low LGM number which no one believed at the time, and\nnot even they believe today, as memorialized here\nAnnan, J. D., Hargreaves, J. C., and Mauritsen, T.: A new global surface temperature reconstruction for the Last Glacial\nMaximum, Clim. Past, 18, 1883–1896, https://doi.org/10.5194/cp-18-1883-2022, 2022.\n“Abstract\nWe present a new reconstruction of surface air temperature and sea surface temperature for the Last Glacial Maximum. The\nmethod blends model fields and sparse proxy-based point estimates through a data assimilation approach. Our reconstruction\nupdates that of Annan and Hargreaves (2013), using the full range of general circulation model (GCM) simulations which\ncontributed to three generations of the PMIP database, three major compilations of gridded sea surface temperature (SST) and\nsurface air temperature (SAT) estimates from proxy data, and an improved methodology based on an ensemble Kalman filter.\nOur reconstruction has a global annual mean surface air temperature anomaly of ∘C relative to the pre-industrial climate. This is\nslightly colder than the previous estimate of Annan and Hargreaves (2013), with an upwards revision on the uncertainty due to\ndifferent methodological assumptions. It is, however, substantially less cold than the recent reconstruction of Tierney et al.\n(2020). We show that the main reason for this discrepancy is in the choice of prior. We recommend the use of the multi-model\nensemble of opportunity as potentially offering a credible prior, but it is important that the range of models included in the PMIP\nensembles represent the main sources of uncertainty as realistically and comprehensively as practicable if they are to be used\nin this way.”\nJust take Annan and Hargreaves 2013 out. Replace it with something folks DO believe. (And many people still believe that\n-4.5degC is on the low end.)\nThe LGM reconstruction of Ice5G is not a very good one. When it came out, no one believed it either – particularly people\nstudying the Laurentide Ice Sheet. They did not believe it because it pulled all of the downsizing ice from the Fenno-Scandian\nand Antarctic and piled it high on the Laurentide. Unbelievably so. Ice 4G and 6G are perhaps defensible, but 5G is not. The\nclimate regimes by such a large ice mountain are just different from those without it.\nThe paper reads much better. And I’m glad the authors took the suggestion to cite more works beyond Yun 2023. But now the\ndiscussion slogs through Erb 2015 and it starts to read like a review instead of a nature article. There are many more models\nand simulations beyond Erb and Yun that have explored the ice sheet height and sensitivity of the climate.\nI thank the authors for their clear explanation\n“CO2, obliquity, eccentricity, precession\nand ice, where 0 = constant PI values and 1 = varying values. Hence E00000 (the ‘no drivers’\nsimulation) has all drivers set to constant PI values, E10000 has varying CO2 but everything else is\nset to constant PI values, and so on, up to E11111 where all forcings are varying (the ‘all drivers’\nsimulation). ”\nIn addition to the (correct) explanation that eccentricity and precession are convolved, why have 5 subscripts instead of 3 (let the\nlast three be ‘orbital’)? The authors state plainly that eccentricity is not tested without precession. Is obliquity tested separately?\nIf obliquity IS tested separately, it would be very interesting to know if obliquity is more important that precession/eccentricity as\nobliquity has been proposed as the cause of glacial cycles (Huybers etc). It would also be interesting to know if this was not the\nresult.\nThe presentation of figure 1 is improved—but I still find the lines too delicate – are the barely visible fuzzy nature the\nuncertainty? It is tricky to have a multipanel plot with the much larger Antarctic temperature shifts and the b-e panels then scaled\nway too zoomed out. There is still much white space and what feels like default subplots. This can still be done better – literally\nmost paleoclimate folks do this with having no extra x-axes, moving the curvy lines physically closer together, alternating y1 and\ny2 axes labels, and expanding the smaller magnitude to be visible (just make sure to note clearly that this is what you’ve done).\nFigure 3. Figure S2. Figure S4. Figure S8. Use a projection like Robinson (etc) that doesn’t blow up the poles so much. I really\nhave no idea why all the scripting languages let awful projections be the default.", "reviewer_3": "The authors have done a nice job, and I have no additional comments. Reviewer #1 did raise some important questions that\nseem to be satisfactorily answered, but I'll happily defer to them."}, "author_response": "Responses to reviewers following submission of manuscript\n“The relative role of direct orbital forcing, and CO and ice sheet feedbacks\n2\non Quaternary climate” by Williams et al.\nREVIEWER COMMENTS\nWe extend our sincere appreciation to both Reviewers for their examination of our\nmanuscript. Here, we address their suggestions, comment-by-comment. In the below,\nreviewers’ comments are blue and in a smaller font, and our corresponding response follows\nin a standard font. Line numbers in the following reviewer comments refer to the previously-\nsubmitted version of the manuscript, whereas line numbers in our responses refer to the new\nversion of the manuscript.\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nReview of Williams et al 2025\nThis paper is greatly improved from previous versions. I do not think I will need to see this again before the\neditor makes a decision.\nWe thank the reviewer for their comments.\n1) I’m still not satisfied with two stitched together emulators being used instead of one. It makes the very\ninteresting and good work seem less so because the answer gets baked in to the training. “The approach of\nhaving two separate emulators was taken to ensure that there was no ‘leakage’ of state-specific climate\nvariations across different climate states (see Supplementary Material, Figure S1e). For example, under glacial\nconditions development of the Laurentide ice sheet is accompanied by strong cooling over northern North\nAmerica. However, this does not mean that during interglacial conditions there should be a strong warming in\nthe same region, since the data are presented as an anomaly from interglacial present-day conditions (where no\nice sheet is present).”\nActually, the change in height alone does mean there is rightfully a huge amount of warming in this area. In\nanomaly space this does look odd—it just does. If you want to get rid of that, try using vertical lapse rate to\napply a correction.\nIn the previous iteration of the manuscript we qualitatively explained this leakage effect in\nterms of the impact of the Laurentide ice sheet (lines 576-593 of the previous iteration of the\nmanuscript), but the effect is even clearer over Greenland and Antarctica. We now quantify\nand illustrate (new Figure S10) this leakage effect in the revised Methods section (‘Step 2’,\nlines 585-615), showing clearly and definitively why two separate emulators are required, as\nfollows:\n“The approach of having two separate emulators (shown in the Supplementary Material, Fig.\nS9) with separate training data for the ‘interglacial’ emulator applied for GSL >=0 m\ncompared with the ‘glacial’ emulator for GSL < 0 m, was taken to ensure that there was no\n‘leakage’ of state-specific climate variations across the glacial-interglacial transition. The\nconcept and effect of this leakage can be illustrated (see Figure S10 in the Supplementary\nMaterial) by considering surface air temperature for two climate states, one with an\nemulated climate with GSL = +5 m (Fig. S10a and b), and one with an emulated climate with\nsea level = -5 m (Fig. S10c and d), firstly as produced by the separate emulators (Fig. S10a\nand c, as used in this paper) and secondly produced by a combined emulator, in which all\ntraining data are used to train a single emulator (Fig. S10b and d). For the +5 m climate\nstate, both the separate and combined emulators show clear warming associated with melt of\nthe Greenland and Antarctic ice sheets (Fig. S10a and b). This is correct, and appears\nbecause the “lowice” training data include ice sheet reduction in these regions, associated\nwith the PRISM4 mid-Pliocene ice sheet reconstruction. However, for the -5 m climate state,\nthere are clear differences between the results of the separate and combined emulator\nversions in these regions (Fig. S10c and d). For the combined emulator (Fig. S10d), the\nalgorithm is still influenced by the Pliocene “lowice” training data, producing anomalous\ncooling in the regions that show melt in the “lowice” training data. However the two\nseparate emulators are, for this GSL, only influenced by the “highice” and “modice”\ntraining data, and so do not have this unphysical “leakage” across the glacial-interglacial\ntransition (Fig. S10c). In essence, the use of separate emulators is required because in the\n“real” climate system, warming relative to modern results in ice sheet loss primarily in the\nGreenland/Antarctica region, whereas cooling relative to modern results in ice sheet growth\nprimarily in the Laurentide/Fennoscandian region. The use of two separate emulators\ntherefore ensures that the different states are characterized separately and correctly. In\naddition, checks were carried out to ensure that there were no discontinuities across the\nglacial-interglacial transition between the climate projections from the two different\nemulators. This lack of discontinuity is unsurprising, given that the “modice” ensemble was\ncommon to both emulators, and sits at the surface GSL=0, on the boundary of the two\nemulators. We therefore consider it appropriate to use the two separate emulators, as\ncombining them into one leads to higher uncertainties at the boundary between glacials and\ninterglacials.”\nFig. S10. Surface air temperature anomalies (compared to the preindustrial control (i.e. 0 kyr), °C) for two\nclimate states, and using two approaches to running the emulator: a) climate state with GSL = +5 m, separate\nemulator approach; b) climate state with GSL = +5 m, combined emulator approach; c) climate state with GSL =\n-5 m, separate emulator approach; d) climate state with GSL = -5 m, combined emulator approach.\n2) If the authors used a single emulator – show it. It makes no sense the explanation that there is leakage over\nNorth America and then differences between the two emerge on the opposite side of the world “When compared\nto proxy data at the five sites presented in Figure 1, the results showed no difference in SST but higher\nvariability for SAT over Antarctica (not shown).”. Show both emulators.\nGiven the above explanation, we have removed the sentence referring to this additional figure\n(lines 612-615) because, as discussed, using a combined emulator produces incorrect results.\n3) The “M score” is not described adequately for readability in the text. In the main text, just a few short words\nlike “M-score (Arcsin-Mielke measure M) quantifies the goodness of fit of the emulator to the proxy climate\narchives such that where no skill is 0 and maximum skill is 1000 (see methods)” will allow the reader to keep\nreading without “flipping” to the back page.\nWe have added some text into the Results section of the main manuscript (lines 264-269)\nsimilar to the above suggestion, and removed similar text from the methodology to avoid\nrepetition (lines 835-843).\n4) The authors insist on using the 2013 Annan and Hargreaves extremely low LGM number which no one\nbelieved at the time, and not even they believe today, as memorialized here\nAnnan, J. D., Hargreaves, J. C., and Mauritsen, T.: A new global surface temperature reconstruction for the Last\nGlacial Maximum, Clim. Past, 18, 1883–1896, https://doi.org/10.5194/cp-18-1883-2022, 2022.\n“Abstract: We present a new reconstruction of surface air temperature and sea surface temperature for the Last\nGlacial Maximum. The method blends model fields and sparse proxy-based point estimates through a data\nassimilation approach. Our reconstruction updates that of Annan and Hargreaves (2013), using the full range of\ngeneral circulation model (GCM) simulations which contributed to three generations of the PMIP database,\nthree major compilations of gridded sea surface temperature (SST) and surface air temperature (SAT) estimates\nfrom proxy data, and an improved methodology based on an ensemble Kalman filter. Our reconstruction has a\nglobal annual mean surface air temperature anomaly of ∘C relative to the pre-industrial climate. This is slightly\ncolder than the previous estimate of Annan and Hargreaves (2013), with an upwards revision on the uncertainty\ndue to different methodological assumptions. It is, however, substantially less cold than the recent\nreconstruction of Tierney et al.(2020). We show that the main reason for this discrepancy is in the choice of\nprior. We recommend the use of the multi-model ensemble of opportunity as potentially offering a credible prior,\nbut it is important that the range of models included in the PMIP ensembles represent the main sources of\nuncertainty as realistically and comprehensively as practicable if they are to be used in this way.”\nJust take Annan and Hargreaves 2013 out. Replace it with something folks DO believe. (And many people still\nbelieve that -4.5degC is on the low end.)\nThis has been done, removing Annan and Hargreaves (2013) and replacing with the data\nassimilation reconstruction of Tierney et al. (2020), for Figure S3 and in the text where this is\ndiscussed (Supplementary Material, Section S1, lines 11-33).\n5) The LGM reconstruction of Ice5G is not a very good one. When it came out, no one believed it either –\nparticularly people studying the Laurentide Ice Sheet. They did not believe it because it pulled all of the\ndownsizing ice from the Fenno-Scandian and Antarctic and piled it high on the Laurentide. Unbelievably so. Ice\n4G and 6G are perhaps defensible, but 5G is not. The climate regimes by such a large ice mountain are just\ndifferent from those without it.\nWe appreciate that there is some uncertainty in the configuration of the ice sheets, including\nfor the LGM. However, it is outside the scope of the manuscript to fully explore this\nuncertainty, as this would require hundreds of additional GCM simulations. However, we do\ndiscuss this uncertainty in the Methods section of the main manuscript (‘Step 2’, at lines 617-\n626) and the Supplementary Material (Section S3, lines 73-84), as follows.\nMain manuscript: “The reconstructed ice sheet extents in the “highice” ensemble are based\non the ICE-5G model of Peltier (2004). These reconstructions were primarily selected\nbecause the associated HadCM3 model simulations already existed (Singarayer and Valdes\n2010). The reconstructions are based on palaeo data (GSL, local sea level, and ice-sheet\nextent) and include the range of associated data that was required for this work. In addition,\nfor glaciations prior to the last glacial cycle, there is very little or no palaeo data available\nthat enables the global 3-D reconstructions which are required here. There are, however,\nother reconstructions (see Section S3 in the Supplementary Material); although we\nacknowledge the existence of these, the “lowice”, “modice” and “highice” ensembles used\nhere, when combined, do capture changes in climate and ice sheet extent ranging from\ninterglacial states to full glacial conditions.”\nSupplementary Material: “Other ice sheet extents reconstructions are available (see, for\nexample, Schmidt et al. 2014), as well as results from recent ice sheet models (see, for\nexample, Tarasov et al. 2025), which could result in changes to the simulated climate. For\nexample, some studies have suggested a lower maximum elevation for the Laurentide ice\nsheet during glacial conditions (e.g. Abe-Ouchi et al. 2015). Modelling studies have\nsuggested that the topography of the ice sheet can have a significant impact on mean sea\nlevel pressure, which affects the location of the wind-driven gyre circulation in the subpolar\nNorth Atlantic and on the strength of the AMOC, resulting in warming in the North Atlantic\n(Colleoni et al. 2016). Storm tracks in the North Atlantic have also been shown to be\naffected, with associated impacts on precipitation and snowfall in northern Europe (e.g.\nColleoni et al. 2016). However, we do not expect these uncertainties to impact our main\nconclusions in terms of the relative importance of the various drivers.”\n6) The paper reads much better. And I’m glad the authors took the suggestion to cite more works beyond Yun\n2023. But now the discussion slogs through Erb 2015 and it starts to read like a review instead of a nature\narticle. There are many more models and simulations beyond Erb and Yun that have explored the ice sheet\nheight and sensitivity of the climate.\nWe have made the discussion of Erb et al. (2015) more concise while keeping all the salient\npoints, and have included some other related modelling studies (lines 383-411).\n7) I thank the authors for their clear explanation “CO2, obliquity, eccentricity, precession\nand ice, where 0 = constant PI values and 1 = varying values. Hence E00000 (the ‘no drivers’\nsimulation) has all drivers set to constant PI values, E10000 has varying CO2 but everything else is\nset to constant PI values, and so on, up to E11111 where all forcings are varying (the ‘all drivers’\nsimulation). ” In addition to the (correct) explanation that eccentricity and precession are convolved, why have 5\nsubscripts instead of 3 (let the last three be ‘orbital’)? The authors state plainly that eccentricity is not tested\nwithout precession. Is obliquity tested separately? If obliquity IS tested separately, it would be very interesting\nto know if obliquity is more important that precession/eccentricity as obliquity has been proposed as the cause\nof glacial cycles (Huybers etc). It would also be interesting to know if this was not the result.\nWe need 5 subscripts because in the sensitivity studies we are individually testing the\nsensitivity to the 3 orbital components, plus ice sheets, plus CO ; therefore, 5 subscripts are\n2\nneeded to fully describe the relative contribution of each orbital individual driver and the ice\nand CO drivers, which necessitates the ability to ‘turn on or off’ any or all of the 5 drivers\n2\nindependently. We have clarified this in the Results section of the main manuscript (lines\n307-309). The obliquity is therefore tested separately, and the results of this are shown in\nFigure 3; here, compared to eccentricity and precession, obliquity is indeed providing a\nhigher contribution, particularly over high latitudes e.g. the Southern Ocean. This is\ndiscussed in the main manuscript (lines 305-307). However, we do not discuss this in the\ncontext of Huybers & Wunsch (2005) because here we are just looking at direct orbital\nforcing, rather than the impact of orbital forcing on CO and ice sheets, which is important\n2\nfor the Huybers hypothesis.\n8) The presentation of figure 1 is improved—but I still find the lines too delicate – are the barely visible fuzzy\nnature the uncertainty? It is tricky to have a multipanel plot with the much larger Antarctic temperature shifts\nand the b-e panels then scaled way too zoomed out. There is still much white space and what feels like default\nsubplots. This can still be done better – literally most paleoclimate folks do this with having no extra x-axes,\nmoving the curvy lines physically closer together, alternating y1 and y2 axes labels, and expanding the smaller\nmagnitude to be visible (just mak’s e sure to note clearly that this is what you’ve done).\nThis has been done: removing the whitespace, changing the y-axis scales to make each\nappropriate for each location, alternating the y-axes and removing the multiple x-axes.\n9) Figure 3. Figure S2. Figure S4. Figure S8. Use a projection like Robinson (etc) that doesn’t blow up the poles\nso much. I really have no idea why all the scripting languages let awful projections be the default.\nThis has been done, with all of the above figures reproduced using the Robinson projection.\nReviewer #3 (Remarks to the Author):\nThe authors have done a nice job, and I have no additional comments. Reviewer #1 did raise some important\nquestions that seem to be satisfactorily answered, but I'll happily defer to them.\nWe thank the reviewer for their previous round of comments."}]

Decision letter