Large-scale features of Last Interglacial climate: Results from evaluating the lig127k simulations for CMIP6-PMIP4

The modeling of paleoclimate, using physically based tools, is increasingly seen as a strong out-of-sample 50 test of the models that are used for the projection of future climate changes. New to CMIP6 is the Tier 1 lig127k https://doi.org/10.5194/cp-2019-174 Preprint. Discussion started: 21 January 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
Quaternary interglacials can be thought of as natural experiments to study the response of the climate system to variations in forcings and feedbacks (Tzedakis et al., 2009). The current interglacial (Holocene, the last 11,600 years) and the last interglacial (LIG; ~129,000-116,000 years before present) are well represented in the geological record and provide an opportunity to study the impact of differences in orbital forcing. Two interglacial timeslices, 75 the mid-Holocene (midHolocene or MH, ~6,000 years before present) and the early part of the LIG (lig127k; 127,000 thousand years before present), are included as Tier 1 simulations in the Coupled Model Intercomparison Project (CMIP6) and Paleoclimate Modeling Intercomparison Project (PMIP4). These equilibrium simulations are designed to examine the impact of changes in the Earth's orbit on the latitudinal and seasonal distribution of incoming solar radiation (insolation) at times when atmospheric greenhouse gas levels were similar to those of the preindustrial 80 period and the continental configurations were also very similar to modern. They test our understanding of the interplay between radiative forcing and atmospheric circulation, and the connections between large-scale and regional climate changes giving rise to phenomena such high-latitude amplification in temperature changes, and responses of the monsoons, as compared to today.
were identified as important time periods to compare data reconstructions and model simulations (Joussaume et al., 1999;Braconnot et al., 2000). A novel aspect in CMIP5 is having the same models and configurations used in the 90 paleo-climate simulations as with the transient 20th century and future simulations, implying that consistency -both in the overall forcings and in how they are imposed -with those other experiments is a very strong requirement. In addition to MH and LGM experiments, CMIP5 and PMIP3 included coordinated protocols for Last Millennium (LM, 850-1850 CE) and the mid-Pliocene Warm Period (mPWP, 3.3-3.0 million years ago) experiments.

95
The LIG is recognized as an important period for testing our knowledge of climate and climate-ice sheet interactions to forcing in warm climate states. Although the LIG was discussed in the First Assessment Report of the IPCC (Folland et al., 1990), it gained more prominence in the IPCC Fourth and Fifth Assessments (AR4 and AR5) (Jansen et al., 2007;Masson-Delmotte et al., 2013). Evidence in the geologic record indicate a warm Arctic (Turney et al., 2010) and a global mean sea level highstand at least 5 m higher (but probably no more than 10 m higher) than present 100 for several thousand years (Dutton et al., 2015) during the LIG. The ensemble of LIG simulations examined in the AR5 (Masson-Delmotte et al., 2013) was not wholly consistent; the orbital forcing and GHG concentrations varied between the simulations. While it had been suggested that differences in regional temperatures between models might reflect differences in cryosphere feedback strength (Yin and Berger, 2012;Otto-Bliesner et al., 2013) or differences in the simulation of the Atlantic Meridional Overturning Circulation (AMOC) (Bakker et al., 2013), differences 105 between models could also have arisen because of differences in the experimental protocols. Furthermore, the LIG simulations were mostly made with older and/or lower-resolution versions of the models than were used for future projections, making it more difficult to use the results to assess model reliability .
For the first time the LIG has been designated as a CMIP6 simulation, setting a common experimental protocol and 110 asking modeling groups to run with the same model and at the same resolution as the DECK simulations (Otto-Bliesner et al., 2016). At the PAGES QUIGS workshop in Cambridge in 2015, the community identified the 127ka time slice for the CMIP6-PMIP4 LIG experiment for several reasons: large Northern Hemisphere seasonal insolation anomalies, no (or little) remnants of the North American and Eurasian ice sheets, and sufficient time to allow for dating uncertainties to minimize the imprint of the previous deglaciation and the Heinrich 11 (H11) meltwater event 115 (Marino et al., 2015). The Tier 1 lig127k experiment addresses the climate responses to stronger orbital forcing, relative to the midHolocene. It also provides a basis to address the linkages between ice sheets and climate change in collaboration with the Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6) (Nowicki et al., 2016).
In this paper, we provide a brief review of the experimental design of the lig127k simulation (Otto-Bliesner et al., 120 2016), an overview of the models included in the ensemble, and details of the new syntheses of surface temperature and precipitation proxies, targeted for 127 ka, developed for this project. We present the multi-model ensemble mean and model spread for surface temperature, precipitation, and sea ice responses as compared to the CMIP6 DECK piControl simulations from the same models and as compared to the data reconstructions. We also explore differences in the responses of surface temperature, monsoon precipitation, and Arctic sea to the different magnitudes 125 https://doi.org/10.5194/cp-2019-174 Preprint. Discussion started: 21 January 2020 c Author(s) 2020. CC BY 4.0 License. and seasonal character of the insolation anomalies at 127 ka versus 6 ka. We conclude with a discussion of possible reasons for the model-data differences.

Experimental design
The CMIP DECK piControl for 1850 CE (see Eyring et al. 2016 for description of this experiment) is the reference 130 simulation to which the lig127k paleo-experiment is compared. The modeling groups were asked to use the same model components and follow the same protocols for implementing external forcings as are used in the piControl.
The boundary conditions for the lig127k, and piControl experiments are given in Table 1, and more detailed information is given below.

135
Earth's orbital parameters (eccentricity, longitude of perihelion, and obliquity) are prescribed following Berger and Loutre (1991). The DECK piControl simulations use the orbital parameters appropriate for 1850 CE (Table 1, Fig. 1) (Eyring et al., 2016), when perihelion occurs close to the boreal winter solstice. The orbit at 127 ka is characterized by larger eccentricity than at 1850 CE., with perihelion occurring close to the boreal summer solstice (Table 1, Fig.   1). The tilt of the Earth's axis was maximal at 131 ka and remained higher than in 1850 CE through 125 ka; obliquity 140 at 127 ka was 24.04° (Table 1). The solar constant for the lig127k simulations is prescribed to be the same as in the DECK piControl simulation, which is fixed at the mean value for the first two solar cycles of the historical simulation (i.e. 1850-1871) (Eyring et al., 2016).
The orbital parameters affect the seasonal and latitudinal distribution and magnitude of solar energy received at the 145 top of the atmosphere, resulting in large positive insolation anomalies during boreal summer at 127 ka as compared to 1850 CE (Fig. 1). Positive insolation anomalies are present from April to September and from 60°S to 90°N. These anomalies peak at over 70 W m -2 in June at 90°N. Insolation in the Arctic (defined here as 60-90°N) is more than 10% greater at 127 ka than 1850 CE during May through early August. The higher obliquity at 127 ka contributes to a small but positive annual insolation anomaly compared to preindustrial at high latitudes in both hemispheres. The 150 global difference in annual insolation forcing between the lig127k and piControl experiments is negligible.
Ice-core records from Antarctica provide measurements of the well-mixed GHGs: CO2, CH4, and N2O. By 127 ka, the concentrations of atmospheric CO2 and CH4 had increased from their respective levels during the previous glacial period, to values comparable to, albeit somewhat lower than, preindustrial levels (Table 1).

155
Natural aerosols show large variations on glacial-interglacial time scales, with low aerosol loadings during interglacials compared to glacials, and during the peak of the interglacials compared to present day (Albani et al., 2015;deMenocal et al., 2000;Kohfeld and Harrison, 2000). Modeling groups were asked to implement changes in atmospheric dust aerosol in their lig127k simulations following the treatment used for their DECK piControl 160 simulations (see Table 2 for details). The background volcanic stratospheric aerosol used in the CMIP6 DECK https://doi.org/10.5194/cp-2019-174 Preprint. Discussion started: 21 January 2020 c Author(s) 2020. CC BY 4.0 License. piControl was to be used for the lig127k simulation. Other aerosols included in the DECK piControl simulations should similarly be included in the lig127k simulations.
There is evidence for changes in vegetation distribution during the LIG (e.g. LIGA Members, 1991;CAPE, 2006; 165 Larrasoana, 2013). However, there is insufficient data coverage for many regions to be able to produce reliable global vegetation maps. Furthermore, given the very different levels of complexity in the treatment of vegetation properties in the current generation of climate models, paleo-observations do not provide sufficient information to constrain their behavior in a comparable way. The treatment of natural vegetation in the lig127k simulations was therefore to be the same as in the DECK piControl simulation. Accordingly, depending on what was done in the DECK 170 piControl simulation, vegetation could either be prescribed to be the same as in that simulation, or prescribed but with interactive phenology, or predicted dynamically (Table 2).
Paleogeography and ice sheets were to be kept at their present-day configuration.

Models
The lig127k simulations performed with 17 models are evaluated in this paper (

Calendar adjustments
The output is corrected following Bartlein and Shafer (2019), to account for the impact that the changes in the length of months or seasons over time, related to changes in the eccentricity of Earth's orbit and precession, have on the analysis (i.e 'celestial' calendar). Not considering the "paleo-calendar effect" can prevent the correct interpretation of data and model comparisons at 127ka, with the largest problems occurring in boreal fall/austral spring (Bartlein and 185 Shafer, 2019).

Preindustrial evaluation
While it is not appropriate to provide a complete evaluation of all the models' pre-industrial simulations here, it is nevertheless important to demonstrate that they are suitable for paleoclimate purposes. Figure 2 presents the 190 climatological temperature  and precipitation (Adler et al., 2003) determined over the historical, observed period. The ensemble mean difference between each of the model's pre-industrial climatology and these observations is also shown (termed here as ensemble mean bias). The climatological (reanalyzed) surface temperatures are taken from the Early Industrial (1871-1900Compo et al., 2011) to coincide with the reference frame used for proxy reconstructions. It is clear that the PMIP4-CMIP6 models are in general cooler that than the 195 observations, most noticeably in the poles, over land and the NH oceans (Fig. 2). The models are too warm over the https://doi.org/10.5194/cp-2019-174 Preprint. Discussion started: 21 January 2020 c Author(s) 2020. CC BY 4.0 License. eastern boundary ocean upwelling regions, Atlantic and Indian Ocean sectors of the Southern Ocean, and the Labrador Current -indicating a failure to sufficiently resolve regional ocean circulation features.
The climatological zonal mean temperature profile is compared to the simulations (Fig. 3). Whilst the model ensemble compares favorably to the HadCRUT4 database (Morice et al., 2012;Ilyas et al. 2017) in the Tropics, the 200 polar latitudes are still too cold in the simulations (Fig. 3).
The piControl simulations of minimum (August-September) Arctic sea ice extent (Figs. 4a and S2), i. e. area where sea ice concentration is at or above 15%, show good agreement with the 15% contour from the HadISST data averaged over the 1870-1920 period . Only one of the models, FGOALS-g3, shows greater 2). Rainfall associated with the mid-latitude storm tracks is generally too little, and the precipitation over land can show consistent wet and dry biases (Fig. 2).

220
The seasonal character of the insolation anomalies results in strong warming and cooling over the Northern Hemisphere (NH) continents in the lig127k ensemble (relative to the piControl) in JJA and DJF, respectively. The warming during JJA is greater than 6°C at mid-latitudes in North America and Eurasia (Fig. 5a), though with significant differences among the models (Fig. 5b). Subtropical land areas in the Southern Hemisphere (SH) respond to the positive (but more muted) insolation anomalies, with JJA temperature anomalies more than 2°C. The North 225 Atlantic, North Pacific and Southern Ocean also exhibit warming, with the latter associated with the simulation of winter sea ice around Antarctica, though with large differences across the ensemble of models. Cooling over the Sahel and southern India in JJA is associated with the increased cloud cover associated with the enhanced monsoons (see Section 3.4).

230
The lig127k multi-model ensemble simulates cooling during DJF relative to the piControl, in response to the negative insolation anomalies at all latitudes (Fig. 5c). The largest DJF temperature anomalies occur over southeastern Asia and northern Africa. Cryosphere (and ocean) memory, particularly lower sea ice area allowing enhanced air-sea https://doi.org/10.5194/cp-2019-174 Preprint. Discussion started: 21 January 2020 c Author(s) 2020. CC BY 4.0 License. interaction, provide the feedback to maintain positive or neutral temperature anomalies annually in the Arctic, northern Europe, and over the Southern Ocean. As indicated by the standard deviations of the ensemble changes, 235 large differences in the magnitude of the DJF surface temperature responses, and thus cryosphere feedbacks, exist across the multi-mean ensemble (Fig. 5d).
Annually, the multi-model ensemble changes in surface temperature between the lig127k and piControl are generally less than 1°C, except for negative surface temperature anomalies across the North African and Indian monsoon 240 regions, and positive surface temperature anomalies in the Arctic (Fig. 5e). Large differences exist among the models in the magnitude of the seasonal and annual surface temperature responses in the polar regions.
The large spread of mean annual surface temperature change among the models in the Arctic (60°-90°N) is further illustrated in Figure 3. The FGOAL-g3 and HadGEM3 simulate 3°C and over 4°C warmer annual surface 245 temperatures in their lig127k simulations than piControl, respectively. Two of the models simulate no change or cooling annually for the Arctic: AWI and ACCESS, respectively. The spread (and magnitude) of mean annual temperature change for the SH polar region (60-90N) is less, with 7/10 models simulating a modest warming of 0-1°C, and three models simulating a cooling of the mean annual surface temperature. In the NH monsoon zonal band from 1-30°N, 7/10 models simulate colder mean annual surface temperature, with the other 3 simulating little or no 250 change.

Sea ice responses
Boreal insolation anomalies at 127 ka enhance the seasonal cycle of Arctic sea ice ( There appears to be a relationship between the Equilibrium Climate Sensitivity (ECS) ( Table 2)  Previous studies suggest that the mean-ice state in the control climate can influence the magnitude and spatial 290 distribution of warming at high-latitudes in future projections (Holland and Bitz, 2003). Thinner Arctic sea ice is more susceptible to summer melting than thicker Arctic sea ice. Figure 8 investigates if this relationship holds for the models in the multi-model ensemble. Arctic sea ice thickness, averaged for model grid cells with at least 15% concentration, varies substantially across the 1850 CE ensemble, ranging from 1-1.5m in CMRM-CM6-1 and NESM3 to ~7.5m in MIROC-ES2L. No robust relationship to the August-September lig127k minimum Arctic sea ice 295 area anomaly is present. One reason for a lack of any relationship may be the seasonal nature of the lig127k insolation forcing as compared to the annual forcing by greenhouse gas changes in future projections.

Precipitation responses
The seasonal character of the insolation anomalies results in enhanced summer monsoonal precipitation in the lig127k ensemble (relative to the piControl ensemble) over northern Africa, extending into Saudi Arabia, India and 300 southeast Asia, and northwestern Mexico/southwestern U.S (Figs. 9, 10). In contrast, summer monsoonal precipitation decreases occur over South America, southern Africa, and Australia. The spread among models is large, however, as shown by the standard deviations of ensemble changes of the precipitation changes ( Fig. 9 and percentage changes in area-averaged precipitation during the monsoon season for seven different regional monsoon domains for the individual lig127k simulations (Fig. 16). The models generally agree on the sign of the percentage 305 changes in the area-averaged precipitation rate during the monsoon season for the monsoon regions, except for the https://doi.org/10.5194/cp-2019-174 Preprint. Discussion started: 21 January 2020 c Author(s) 2020. CC BY 4.0 License. varying from ~70-140% (Fig. 18c). The models are in closer agreement for the East Asian monsoon (EAS), with the percentage change in the areal extent varying from ~10-35% (Fig. 18b) and the percentage change in the total amount of water precipitated in each monsoon season varying from ~25-40% (Fig. 18c). The lig127k and piControl simulations produce more muted changes for the other monsoon regions in the multi-model ensemble, with regards to the regional monsoon-related rainfall rate and the monsoon domains (Fig. 10). 325 4 Data reconstructions

Marine temperatures
The lig127k climate model simulations are assessed using two complementary compilations of sea surface temperature (SST) anomalies at 127 ka (Tables S1 and S2), which are both individually based on stratigraphically consistent chronologies Hoffman et al., 2017).

330
The multi-archive high-latitude compilation by Capron et al. (2014Capron et al. ( , 2017 includes 42 sea surface annual/summer temperature records with a minimum temporal resolution of 2 kyr for latitudes above 40°N and 40°S, along with 5 ice core surface air temperature records. In contrast, the global marine compilation by Hoffman et al. (2017) includes 186 annual, summer and winter SST records from the Atlantic, Indian and Pacific oceans, with a minimum temporal resolution of 4 kyr on their published age models. Note that, in addition to the annual microfossil assemblage SST 335 records calculated for 41 sites as the average of the summer and winter records with a model-and observationconsistent correction for annual offsets (Hoffman et al., 2017), we also provide here for these specific sites the updated seasonal (summer and winter) SST estimates on the Hoffman et al. age models. SST from marine cores are reconstructed in both compilations from foraminiferal Mg/Ca ratios, alkenone unsaturation ratios or microfossil faunal assemblage transfer functions (Capron et al., 2014Hoffman et al., 2017). To derive the LIG marine chronologies, both compilations make use of the climate model-supported hypothesis that surface-water temperature changes in the sub-Antarctic zone of the Southern Ocean (respectively in the North Atlantic) occurred simultaneously with air temperature variations above Antarctica (respectively Greenland) (Capron et al., 2014;Hoffman et al., 2017). The compilation by Hoffman et al. (2017) then uses basin-synchronous LIG changes in the oxygen isotopic composition of benthic foraminifera, as observed in previous studies of benthic 345 foraminiferal isotope changes across glacial terminations (Lisiecki and Raymo, 2009) within the same ocean basins, to align intra-basin chronologies. However, a major difference is the underlying reference chronology used in both compilations: the Antarctic Ice Core Chronology 2012 (AICC2012) Veres et al., 2013) in the compilation by Capron et al. (2014Capron et al. ( , 2017, and a chronology based on millennial-scale variations observed in independently-dated Asian speleothem records (Speleo-Age) (Barker et al., 2011) in the 350 compilation by Hoffman et al. (2017). Note that the two reference chronologies diverge by about 1 ka at 127 ka .
The two compilations then follow quite similar Monte Carlo approaches to propagate temperature and chronological uncertainties. Indeed, both compilations generate 1000 realizations of the site-specific surface temperature records to integrate the uncertainty on the temperature reconstruction's method, and both produce 1000 possible chronologies to 355 propagate the relative age uncertainty related to alignment of records. For a given site, the temperature at 127 ka is the temperature value directly taken at 127 ka in the compilation by Hoffman et al. (2017), using dated temperature timeseries interpolated every 1 ka. In the compilation of (Capron et al., 2014, the temperature at 127 ka is taken as the median temperature averaged over the 128-126 ka period. Finally, temperatures anomalies relative to preindustrial are calculated in both cases for marine sites using the HadISST dataset

Ice core temperatures
Surface air temperature records for one site (NEEM) on the Greenland ice sheet and four sites on the Antarctic ice sheet are deduced from ice core water isotopic profiles (Capron et al., 2014 (Table S3). For ice cores, preindustrial conditions are estimated using borehole temperature measurements for Greenland, and 1870-1899 CE 375 https://doi.org/10.5194/cp-2019-174 Preprint. Discussion started: 21 January 2020 c Author(s) 2020. CC BY 4.0 License. water isotopic profiles for Antarctic . Temperatures are again the median for the 126-128 ka period, and are considered to represent annual averages. Uncertainty is estimated using the same Monte Carlo procedure as was used for the marine cores in the compilation of Capron et al. (2017). Because it uses the same reference timescale (AICC2012) the ice core dataset can be considered coherent with the marine SST dataset of Capron et al (2017).  (Table S4). For both we report anomalies comparing 385 reconstructed temperatures with preindustrial climate estimated from 1871-1900.
In Europe, favorable geological conditions have led to the accumulation of numerous LIG sediment sequences from a variety of depositional environments (Tzedakis, 2007). These include former kettle lakes overlying late Saalian (MIS 6) till, depressions left by the penultimate alpine glaciation or local ice caps, and volcanic crater lakes or tectonic grabens mainly in the unglaciated south. Over several decades, a substantial body of pollen evidence has provided an 390 insight into the LIG vegetational development across Europe. A number of pollen-based climate reconstructions on reference sequences have been attempted, using a variety of methods. However, differences between underlying assumptions and data employed (e.g. taxon presence-absence versus abundance) mean that results have been difficult to compare.

395
Here, we include data from one study that has applied a multi-method approach to assess combined uncertainties of reconstruction and age models on a set of reference pollen records (Brewer et al., 2008). The reconstruction methods used are (i) partial least squares, (ii) weighted average partial least squares, (ii) generalized additive models, (iv) artificial neural network, (v) unweighted modern analogue technique, (vi) weighted modern analogue technique, and (vii) revised analogue method using response surfaces. Timescales for the pollen sequences were developed by 400 transferring the marine chronology to land sequences for certain pollen stratigraphical evens on the basis of joint pollen and palaeoceanographic analyses in deep-sea sequences on the Portuguese Margin and Bay of Biscay (Shackleton et al., 2003;Sánchez Goñi et al., 2008). With particular reference to constraining the 127 ka timeslice, the pollen stratigraphical events used were the onset of the Quercus (128.8±1 ka) and Carpinus (124.77±1 ka) expansions (Brewer et al., 2008). For each site, chronological uncertainties were estimated at each sample by 405 randomly sampling an age from the range around each control point, fitting a linearly interpolated age model and repeating this 1000 times (Brewer et al., 2008). Reconstructions were made at 500 yr intervals by randomly sampling within the chronological uncertainties and reconstruction errors for each method, resulting in 1400 estimates for each time t (Brewer et al., 2008). Here we present the mean value and standard deviation for mean annual temperature, mean temperature of the coldest month and mean temperature of the warmest month across all sites for 127 ka. Of 410 https://doi.org/10.5194/cp-2019-174 Preprint. Discussion started: 21 January 2020 c Author(s) 2020. CC BY 4.0 License. the seventeen sites used by Brewer et al. (2008), four were excluded because they either did not extend to 127 ka or were from marine pollen records.
The Arctic dataset compiles the most stratigraphically complete, well-constrained calibrated summer temperature reconstructions published from each subregion of the Arctic above 65°N latitude. We report the mean of the two 415 warmest consecutive reconstructions at each site, utilizing the original published models and reconstructions. For sites where both insect-and pollen-based temperature reconstructions have been published, or where multiple models have been applied to the same proxy, we report here the average of those reconstructions. We report the original published model uncertainties (e.g., root mean square error of prediction for weighted averaging models), including the most conservative (largest) model uncertainties for sites where multiple proxies/models are applied. Importantly, 420 the Arctic compilation differs from all other paleotemperature datasets used here, in that it reports the warmest LIG conditions registered at each site rather than temperatures at 127 ka. This approach was necessitated by the coarse temporal resolutions and chronologies of the North American Arctic reconstructions, which come from stratigraphically discontinuous deposits dated by 14 C (non-finite 14 C ages) and in some cases luminescence or tephrochronology. In contrast to the North American Arctic sites, in northern Finland (Sokli) and northeast Russia 425 (El'gygytgyn) correlative dating provides continuous chronologies. The reported peak warmth at those sites occurred at ~125 and 127-125 ka, respectively (Melles et al. Science 2012;Salonen et al. Nature Comms 2018). Reconstructed temperature at Sokli at 127 ka was ~1°C lower than the peak temperature reported here from that site. The Greenland ice core-derived temperature reconstruction from NEEM complements the Arctic terrestrial dataset, but it reflects annual rather than summer-specific climate.

Arctic sea ice
A summary of LIG sea ice data obtained from marine cores in the Arctic, Nordic Seas and northern North Atlantic, their interpretation, and comparison to the lig127k simulations can be found in Kageyama et al. (2019). The sea ice records are derived from dinoflagellate cysts, subpolar foraminers, and ostracods.

435
Compilation of the existing proxy evidence for LIG precipitation have been presented for the Northern Asian and circum-Arctic region (CAPE, 2006;Kim et al., 2010;Velichko et al., 2008). Recently, a compilation with near-global coverage was presented in Scussolini et al. (2019), including 138 proxy sites based on different types of proxies and archives, mostly from pollen, lacustrine sediment composition, speleothem, and multi-proxy reconstructions. This, in contrast to previous work, aimed to select proxy signals approximately corresponding to 127 ka, in order to facilitate 440 comparison with results from the lig127k simulations of PMIP4. The main patterns that emerge, about precipitation change between the LIG and the pre-industrial/recent past, are near-ubiquitous higher LIG annual precipitation over the NH (Fig. 14). Exception to this are individual sites in western North Africa, the Levant, northern South America, Borneo, the northwest of modern United States and Alaska, northern Scandinavia, and northern Siberia. Over the SH, the proxy signal is more irregular: Australia and the west coast of South America have proxies predominantly 445 https://doi.org/10.5194/cp-2019-174 Preprint. Discussion started: 21 January 2020 c Author(s) 2020. CC BY 4.0 License.
indicating higher precipitation in the LIG, sites in the rest of South America indicate lower precipitation or no change and over southern Africa changes are geographically more heterogeneous.

Temperature
NH high-latitude terrestrial temperature proxies for the boreal summer (JJA) match the large warming from the 450 lig127k ensemble for the Canadian, Greenland and Scandinavian sites ( Figure 11). On the other hand, although the models also simulate large warming over Alaska and Siberia in JJA, these sites have temperature anomalies near zero. Over Europe the temperature proxies can indicate substantial warmings right next to substantial coolings making a credible comparison difficult, though a much weaker lapse rate at 127 ka (Brewer et al., 2008)  Labrador and Greenland Seas, which cannot be evaluate given the available reconstructions. Simulated LIG DJF 465 temperatures over North America and Eurasia are significantly colder with respect to preindustrial, except over western Europe. The proxies for the latter again show a mixed signal.
The most consistent picture from the temperature proxies representing annual conditions is warmer LIG temperatures over Greenland and Antarctica. The ensemble mean matches the pattern but not the magnitude of LIG temperature 470 anomalies over the polar ice sheets reconstructed from the ice cores. However, the model spread in the temperature response averaged over the 60-90 o latitude range is large as shown in Fig. 2, and thus, some models may compare more favorably with the temperature proxies than others, which pulls down the mean response. It should also be noted that internal variability is large in high latitudes, which can also hinder data-model comparisons. Elsewhere, at high latitudes in both hemispheres, the temperature proxy reconstructions indicate a mixed response at LIG.

475
At low latitudes (40°S -40°N), marine proxy data from the Hoffman reconstruction are available (Fig. 12). They generally correspond with the ensemble mean changes. The SST proxies from the tropical Atlantic match the colder CMIP models is generally not adequate to properly simulate these narrow coastal upwelling regions, and their warm bias in the preindustrial simulations (Fig. 2) leaves little potential for strong SST increases.

485
As shown in a comparison with a similar ensemble (Scussolini et al., 2019), precipitation proxies from the global compilation largely match the higher precipitation from the models included in the lig127k ensemble over much of the NH continents: over most of northern Africa, the Middle-East, the Mediterranean, South and East Asia, northeast Asia, and North America (Fig 10). Proxies and models weakly disagree over much of Europe, Central Asia and the 490 region between them, where proxies indicate wetter LIG conditions or no change, and the ensemble average indicates somewhat drier conditions or no change. Other instances of disagreement in the NH are over the southern side of Northern Africa, with drier proxies and wetter models, and over the Mississippi basin, with a wetter proxy site and somewhat drier model average. However, the coastal proxy sites near the Bay of Bengal, which show strongly drier conditions, are near the region of strongly drier conditions over the Atlantic suggesting a northward shift in the ITCZ.

495
In the SH, proxies and models mostly agree over South America, while they disagree over Australia and in several locations over southern Africa, where many proxies and models indicate wetter and drier LIG, respectively. More details of a comparison with a similar model ensemble can be found in Scussolini et al. (2019).

Comparison of the model sensitivities to the insolation anomalies at 127 ka and ka
The large-scale features and evaluation of the CMIP6/PMIP4 midHolocene simulations in comparison to data 500 reconstructions and in the CMIP5/PMIP3 endeavor can be found in Brierley et al., 2020). In this section, we briefly explore differences in the responses of surface temperature, monsoon precipitation, and Arctic sea to the different magnitudes and seasonal character of the insolation anomalies at 127 ka versus 6 ka.

Orbital forcing
The orbit at 6 ka was characterized by a smaller eccentricity that at 127 ka, similar to 1850 CE. Perihelion occurred 505 near the boreal autumn equinox as compared to close to the boreal summer solstice at 127ka and near aphelion at 1850CE. NH summer insolation anomalies at 6ka, ~5-10% greater than at 1850 CE, are considerably less that at 127ka ( Fig. 1 and Fig. 14). In addition, the positive insolation anomalies of greater than 10% in the Arctic occur in July-August at 6ka as compared to May-August at 127 ka. At SH mid-and high latitudes, the anomalous insolation anomalies are shifted to boreal fall/austral spring. As such, the orbital forcing on climate is expected to be stronger at 510 127 ka than at 6 ka. but are significantly greater at 127ka (over 3°C at 127 ka as compared to ~1°C at 6 ka). The DJF zonal-average surface temperature anomalies are near zero or slightly negative south of 65°N for both time periods. Cryosphere and ocean feedbacks provide the memory for positive surface temperature anomalies in DJF, even with negative insolation anomalies, with DJF Arctic surface temperatures averaging about 0.5°C higher in the midHolocene 520 ensemble mean and up to 3°C higher in the lig127k ensemble mean than the piControl.

Precipitation responses
The sign of the percentage change in the areal extents of the regional monsoon domain (Aav) and the percentage change in the total amount of water precipitated in each monsoon season (Totwater) are similar for the lig127k and  Bartlein et al, (2017) in the mid-Holocene appears to also occur at 127 ka, where simulations indicate a weak drying or no change while proxies suggest it was wetter (Fig. 13).

535
The boreal insolation anomalies at 6 ka enhance the seasonal cycle of Arctic sea ice, though much less so than in the lig127k simulations (Fig. 17). None of the models currently in this analysis have the Arctic becoming ice-free in their midHolocene simulations.

Conclusions
The Tier 1 lig127k experiment was designed to address the climate responses to stronger orbital forcing (relative to 540 the midHolocene experiment) using the same state-of-the-art models and following a common experimental protocol.
At 127 ka, atmospheric greenhouse gas levels were similar to those of the preindustrial period, land ice likely only remained over Greenland, and Antarctica, and the continental configurations were almost identical to modern. In addition, within uncertainties in chronology and dating, this time period allows data reconstructions for comparison to the model simulations allowing an assessment of responses to the large insolation changes, yet does not include the 545 effects of the Heinrich 11 meltwater event. The 17 models that have completed the lig127k experiment are presented.
The CMIP6-PMIP4 lig127k simulations show strong warming and cooling over the NH continents during JJA and DJF, respectively, in response to the seasonal character of the insolation anomalies. The JJA ensemble average warming is greater than 6°C at mid-latitudes in North America and Eurasia, though with significant differences 550 https://doi.org/10.5194/cp-2019-174 Preprint. Discussion started: 21 January 2020 c Author(s) 2020. CC BY 4.0 License. among the models. The simulations exhibit reduced minimum (August-September) summer sea ice extent (defined as 15% concentration or greater) in the Arctic, but with substantial differences among the ensemble and one model having the Arctic becoming ice-free in late summer. Cryosphere and ocean memory and increased air-sea interaction provide the feedback to maintain positive temperature anomalies in the lig127k simulations annually in the Arctic, northern Europe, and over the Southern Ocean, though again with a large spread across the model ensemble. As 555 expected from the larger insolation anomalies in the lig127k than midHolocene simulations, the boreal summer responses in NH surface temperature and Arctic sea ice are amplified.
The CMIP6-PMIP4 lig127k simulations produce enhanced summer monsoonal precipitation and areal extent over northern Africa, which extends into Saudi Arabia, India and southeast Asia, and northwestern Mexico/southwestern 560 U.S. In contrast, summer monsoonal precipitation decreases over South America, southern Africa, and Australia.
The spread across the multi-model ensemble is particularly large for the North African monsoon, for the percentage change in its areal extent and total amount of water precipitated during the monsoon season. The lig127k individual monsoon changes are of a similar sign, but a greater magnitude, to those seen in the midHolocene simulations (Brierley et al., 2019).

565
New syntheses for surface temperature and precipitation, targeted for 127ka, have been developed for comparison to the CMIP6-PMIP4 lig127k simulations. The surface temperature reconstructions include two complimentary compilations of SST based on stratigraphically consistent chronologies, surface air temperatures from the Greenland and Antarctic ice sheets deduced from the ice core water isotopic profiles, continental air temperatures for Europe 570 based on pollen records, and peak LIG summer temperatures in the Arctic. Anomalies were consistently computed comparing the reconstructed temperatures with preindustrial climate estimate from the end of the 19 th century. A new precipitation reconstruction has expanded from previous regional compilations to now near-global coverage.
Over Canada, Scandinavia, and the North Atlantic, the proxy and lig127k positive JJA temperature anomalies are in 575 good agreement. The exceptions are in the northwestern North Atlantic and the Nordic Sea, where the Capron reconstruction suggest significant cooling. Capron reconstruction also provides evidence of significant positive DJF temperature anomalies over the Southern Ocean, which is not exhibited by the ensemble mean. These mismatches could be associated with meltwater from potentially remnant ice sheets over Canada and Scandinavia as well as memory in the ocean of the H11 event (Marino et al., 2015), which the lig127k simulations do not incorporate, 580 leading to cooling in the North Atlantic and warming in the Southern Ocean (Stone et al., 2016).
The most consistent picture from the temperature proxies representing annual conditions is warmer LIG temperatures over Greenland and Antarctica. The ensemble mean matches the pattern but not the magnitude of LIG temperature anomalies over the polar ice sheets reconstructed from the ice cores. However, the model spread is large. The lig127k 585 Tier 1 experiment protocol prescribed modern Greenland and Antarctic ice sheets rather than allowing them to evolve to smaller and lower ice sheets of the lig127k climate. A modeling study with the HadCM3 (a CMIP3 model) https://doi.org/10.5194/cp-2019-174 Preprint. Discussion started: 21 January 2020 c Author(s) 2020. CC BY 4.0 License. demonstrated that the distinctive peak in d 18 O observed in Antarctic ice cores at 128 ka was likely due to the loss of winter sea ice in the Atlantic, Indian, and Pacific sectors of the Southern Ocean. To achieve this winter sea ice extent required forcing by the H11 meltwater event (Holloway et al., 2017(Holloway et al., , 2018. The CMIP6-PMIP4 Tier 2 LIG 590 experiments (lig127k-H11, lig127-gris, lig127k-ais) will allow modeling groups to explore the effects of the H11 meltwater event and the Antarctic and Greenland ice sheets at their minimum LIG extent and lower elevations (Otto-Bliesner et al., 2016).
Other reasons for mismatches between the models and the reconstructions for temperature and precipitation will also 595 be explored with CMIP6-PMIP4 Tier 2 LIG experiments. Most of the models contributing to the lig127k still do not include interactive vegetation. An intercomparison of the responses of the lig127k simulations with interactive vegetation would be useful.The CMIP6-PMIP4 Tier 2 lig127k-veg experiments will consider the sensitivity of the responses to prescribed boreal forests in the Arctic and shrub/savanna over the Sahara, separately and together.
Incorporating these vegetation changes has been shown to impact the albedo and evapotransporation on the surface 600 energy and water budgets, reducing model and data mismatches at high latitudes (Swann et al., 2010) and for the North African monsoon (Pausata et al., 2016). Additionally, the CMIP6 models do not currently simulate changes to soil texture or color for different climate states. A previous modeling study suggests that soil feedback can drive the African monsoon northward during interglacials (Levis et al., 2004).

605
Despite an abundance of LIG pollen records from Eurasia and various attempts at pollen-based climate reconstructions (e.g. compilations Velichko et al., 2007;Turney & Jones 2010), chronological and methodological uncertainties continue to complicate comparisons with climate model outputs. The lack of spatial coherence in the European temperature reconstructions may reflect depth-age model issues at individual sites, which implies that the 127 ka timeslice had not been correctly identified. An alternative approach would have been to select peak 610 temperatures from a wider interval (e.g. 127±2 ka) and assume that these are quasi-synchronous. However, closer inspection of the Brewer et al. (2008) data reveals that the timing of peaks in annual and JJA temperatures are often not synchronous. In addition, the Arctic reconstruction may be skewed a bit warmer than the models, given that by definition this reconstruction is reporting the warmest period from each Arctic site, rather than the 127 ka timeslice. Climate Sensitivity (ECS) varying from 2.1 to 5.3°C. There appears to be a clear relationship between the ECS of each model and its simulation of August-September lig127k minimum Arctic sea ice extent. With very limited Arctic sea ice proxies for 127 ka, and with evolving interpretation of the relationships of these proxies with sea ice coverage (Stein et al., 2017;Kageyama et al., 2020), it is difficult to rule out the high or low values of ECS from the proxy 630 data. Warmer summer temperatures over Greenland, warmer oceans year-round surrounding Greenland, and reduced Arctic summer sea ice all have the potential to force a retreat of the ice sheet in the future. The lig127k results can be used to force Greenland ice sheet models, both one-way as included in the ISMIP6 protocols (Nowicki et al., 2016) and fully coupled to a climate model as is now being done by several modeling groups. With the availability of LIG ice and marine core records, LIG simulations with an evolving Greenland ice sheet will allow an assessment of the 635 corresponding future projection simulations.

660
This work is a contribution to the PAGES QUIGS Working Group and the PMIP project. We acknowledge the World Climate Research Programme's Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups for producing and making available their model output.  observations. The area-averaged, annual mean surface air temperature for 30° latitude bands in the CMIP6 models and a spatially complete compilation of instrumental observations over 1850-1900(black, Ilyas et al., 2017Morice et al., 2012). (bottom) Changes in latitudinal temperature gradients (lig127k minus piControl). The simulated annual mean temperature change in averaged over broad zonal bands is presented for each of the individual CMIP6 models.
https://doi.org/10.5194/cp-2019-174 Preprint. Discussion started: 21 January 2020 c Author(s) 2020. CC BY 4.0 License.            (Table 1) on occurs close to the boslightly varies depending the number of days used the fact that the length of recession and eccentricity the vernal equinox must e simulations (piControl, bit at 6 ka was character-, similar to 1850 CE (Ta-5 ) and perihelion at 6 ka uinox. The orbit at 127 ka tricity than at 1850 CE, he boreal summer solstice s was maximal at 131 ka E through 125 ka; obliq-The different orbital cond lig127k result in differtion of top-of-atmosphere piControl (Fig. 3). Both nsolation anomalies durnomalies between 40 and t 127 ka and 25 W m 2 at and 6 ka contributes to a n anomaly compared to n both hemispheres and a opics in the annual mean. forcing between the interstrial is negligible. or the midHolocene and in the DECK piControl ean value for the first two ulation (i.e., 1850-1871) (1360.7 W m 2 ) is lower nt used by some models leads to a global reducmpared to the PMIP3 exerences in orbital paramce periods to be used for eference used for PMIP3 ing with a slight decrease oreal autumn. The comto an overall reduction: real spring and is about N. and Greenland provide HGs: CO 2 , CH 4 , and N 2 O iven as mole fractions in illion (ppm) or parts per Figure 2. Orbital configurations for the piControl, midHolocene, and lig127k experiments. Note that the aspect ratio between the two axes of the ellipse has been magnified to better highlight the differences between the periods. However, the change in ratio between the different periods is proportional to the real values. In these graphs VE stands for vernal equinox, SS for summer solstice, AE for autumnal equinox, and WS for winter solstice. The numbers along the ellipse are the number of days between solstices and equinoxes.
billion (ppb), respectively. For simplicity, we use the term "concentration" for these GHG levels. By 6 and 127 ka, the concentrations of atmospheric CO 2 and CH 4 had increased from their respective levels during the previous glacial periods, the Last Glacial Maximum and the penultimate glaciation, to values comparable to preindustrial levels. midHolocene. In PMIP4-CMIP6, we use a revised version of an earlier trace gas reconstruction (Joos and Spahni, 2008). The CO 2 concentration for the mid-Holocene is derived from ice-core measurements from Dome C (Monnin et al., 2001(Monnin et al., , 2004) and dated using the AICC2012 age scale . A smoothing spline (Bruno and Joos, 1997;Enting, 1987) with a nominal cut-off period of 3000 years was used to produce a continuous CO 2 record. This yields a CO 2 concentration of 264.4 ppm at