Large-scale features of Last Interglacial climate: results from evaluating the lig127k simulations for the Coupled Model Intercomparison Project (CMIP6)–Paleoclimate Modeling Intercomparison Project (PMIP4)

The modeling of paleoclimate, using physically based tools, is increasingly seen as a strong out-of-sample test of the models that are used for the projection of future climate changes. New to the Coupled Model Intercomparison Project (CMIP6) is the Tier 1 Last Interglacial experiment for 127 000 years ago (lig127k), designed to address the climate responses to stronger orbital forcing than the midHolocene experiment, using the same state-of-the-art models as for the future and following a common experimental protocol. Here we present a first analysis of a multi-model ensemble of 17 climate models, all of which have completed the CMIP6 DECK (Diagnostic, Evaluation and Characterization of Klima) experiments. The equilibrium climate sensitivity (ECS) of these models varies from 1.8 to 5.6 ^∘C. The seasonal character of the insolation anomalies results in strong summer warming over the Northern Hemisphere continents in the lig127k ensemble as compared to the CMIP6 piControl and much-reduced minimum sea ice in the Arctic. The multi-model results indicate enhanced summer monsoonal precipitation in the Northern Hemisphere and reductions in the Southern Hemisphere. These responses are greater in the lig127k than the CMIP6 midHolocene simulations as expected from the larger insolation anomalies at 127 than 6 ka.

New synthesis for surface temperature and precipitation, targeted for 127 ka, have been developed for comparison to the multi-model ensemble. The lig127k model ensemble and data reconstructions are in good agreement for summer temperature anomalies over Canada, Scandinavia, and the North Atlantic and for precipitation over the Northern Hemisphere continents. The model–data comparisons and mismatches point to further study of the sensitivity of the simulations to uncertainties in the boundary conditions and of the uncertainties and sparse coverage in current proxy reconstructions.

The CMIP6–Paleoclimate Modeling Intercomparison Project (PMIP4) lig127k simulations, in combination with the proxy record, improve our confidence in future projections of monsoons, surface temperature, and Arctic sea ice, thus providing a key target for model evaluation and optimization.

1 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 time slices, the mid-Holocene (midHolocene or MH, ∼ 6000 years before present) and the early part of the LIG (lig127k; 127 000 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 and hence the latitudinal and seasonal distribution of incoming solar radiation (insolation) at times when atmospheric greenhouse gas levels and continental configurations were similar to those of the preindustrial period. 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.

The modeling of paleoclimate, using physically based tools, has long been used to understand and explain past environmental and climate changes (Kutzbach and Otto-Bliesner, 1982; Braconnot et al., 2012; Harrison et al., 2015; Schmidt et al., 2014). In the first phase of PMIP, the MH and the Last Glacial Maximum (LGM, ∼ 21 000 years ago) 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 was applying the same models and configurations used in the paleoclimate simulations as in the transient 20th-century and future simulations, providing consistency – both in the overall forcings and in how they are imposed – between experiments. In addition to MH and LGM experiments, CMIP5 and PMIP3 included coordinated protocols for the last millennium (LM, 850–1850 CE) and the mid-Pliocene Warm Period (mPWP, 3.3–3.0 × 10⁶ years ago) experiments.

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 indicates a warm Arctic (CAPE, 2006; Turney and Jones, 2010) and a global mean sea level highstand at least 5 m higher (but probably no more than 10 m higher) than the present for several thousand years during the LIG (Dutton et al., 2015). The ensemble of LIG simulations examined in the AR5 (Masson-Delmotte et al., 2013) was not wholly consistent; the orbital forcing and greenhouse gas (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 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 (Lunt et al., 2013).

For the first time an LIG experiment is included as a CMIP6 simulation, setting a common experimental protocol and asking modeling groups to run with the same model and at the same resolution as the DECK simulations (Otto-Bliesner et al., 2017). At the PAGES QUIGS workshop in Cambridge in 2015, the community identified the 127 ka 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 (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 start with a brief overview of the experimental design of the lig127k (Otto-Bliesner et al., 2017). We briefly summarize the simulation of temperature, precipitation, and sea ice, in the subset of CMIP6 piControl simulations that have a corresponding lig127k simulation, as compared to observational datasets. We then provide an initial analysis of the multi-model ensemble mean and model spread in the lig127k surface temperature, precipitation, and Arctic sea ice responses as compared to the CMIP6 DECK piControl simulations. A new syntheses of surface temperature and precipitation proxies, targeted for 127 ka, is used for comparison to the model simulations. We also explore differences in the responses of surface temperature, monsoon precipitation, and Arctic sea ice to the different magnitudes and seasonal character of the insolation anomalies at 127 ka versus 6 ka. We then conclude with a discussion of possible reasons for the model–data differences and implications for future projections.

2 Methods

2.1 Experimental design

The CMIP DECK piControl for 1850 CE (see Eyring et al., 2016, for description of this experiment) is the preindustrial (PI) reference 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 used in the piControl. The boundary conditions for the lig127k and piControl experiments are described in Otto-Bliesner et al. (2017) and the Earth System Documentation (2019). More detailed information is given below and in Table 1.

Table 1 Protocols: forcings and boundary conditions.

^∗ The term “orbital parameters” is used to denote the variations in the Earth's eccentricity and longitude of perihelion as well as changes in its axial tilt (obliquity).

Figure 1 (a, b) Orbital configurations for the piControl and lig127k experiments. The number of days between the vernal equinox and summer solstice, summer solstice and autumnal equinox, etc., are indicated along the periphery of the ellipse. Latitude–month insolation anomalies 127 ka–1850 in (c) W m⁻² and (d) percentage change from PI.

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 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.

The orbital parameters affect the seasonal and latitudinal distribution and magnitude of solar energy received at the 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⁻² 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 1850 CE at high latitudes in both hemispheres and negative annual insolation anomaly at tropical latitudes. The 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: CO₂, CH₄, and N₂O. By 127 ka, the concentrations of atmospheric CO₂ and CH₄ had increased from their minimum levels during the previous glacial period to values comparable to, albeit somewhat lower than, preindustrial levels (Table 1).

Table 2 Summary of CMIP6–PMIP4 models in this intercomparison.

¹ ECS uses the Gregory method from a 150-year run of an instantaneously quadrupled CO₂ simulation (Meehl et al., 2020; Wyser et al., 2020). ² Unless otherwise noted, prescribed vegetation and aerosols are as in each model's piControl simulation.

Natural aerosols show large variations on glacial–interglacial timescales, with low aerosol loadings during interglacials compared to glacials and during the peak of the interglacials compared to the 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 simulations (see Table 2 for details). The background volcanic stratospheric aerosol used in the CMIP6 DECK piControl was also 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; 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 piControl simulation, vegetation could either be prescribed to be the same as in that simulation, prescribed but with interactive phenology, or predicted dynamically (see Table 2 for implementations in the models).

Paleogeography and ice sheets were to be kept at their present-day configuration.

2.2 Model evaluation

The 17 modeling groups that have completed CMIP6 lig127k simulations are presented in this paper (Table 2). All used the CMIP6 version of their model also used for their DECK experiments. The equilibrium climate sensitivity (ECS) varies from 1.8 to 5.6 ^∘C. The years analyzed for each model and DOIs for each of the simulations are given in Table S1 in the Supplement. The analysis uses data available on the CMIP6 ESGF (Earth System Grid Federation) for surface air temperature (tas), precipitation (pr), and sea ice concentration (siconc).

2.3 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 have on the analysis (Fig. 1). This correction is necessary to account for the impact of the changes in the eccentricity of the Earth's orbit and the precession when using the “celestial” calendar. Not considering the “paleo-calendar effect” can prevent the correct interpretation of data and model comparisons at 127 ka, with the largest problems occurring in boreal fall/austral spring (Joussaume and Braconnot, 1997; Bartlein and Shafer, 2019). A more detailed discussion of the application of the PaleoCalAdjust software to past time periods with strong orbital forcing can be found in Bartlein and Shafer (2019) and Brierley et al. (2020).

3 Simulation results

3.1 Preindustrial simulations

Brierley et al. (2020) provide an extensive evaluation of the CMIP6 preindustrial simulations as compared to observational datasets: reanalyzed climatological temperatures (between 1871–1900 CE; Compo et al., 2011) for the spatial patterns, zonal averages of observed temperature for the period 1850–1900 CE from the HadCRUT4 dataset (Morice et al., 2012; Ilyas et al., 2017), and climatological precipitation data for the period between 1970 and the present day (Adler et al., 2003). In summary, they find that the PMIP4–CMIP6 models are in general cooler than the observations, most noticeably at the poles, over land, and over the NH oceans. The poleward extent of the North African monsoon, in particular, is underestimated in the CMIP6 preindustrial simulations.

Figure 2 (a) Comparison of the preindustrial zonal mean temperature profile of individual climate models and MMM to the 1850–1900 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., 2017; Morice et al., 2012). (b) Changes in zonal average, mean annual surface air temperatures (lig127k minus piControl).

The CMIP6 midHolocene and lig127k have 14 models in common (see Fig. 4a in Brierley et al., 2020, and Fig. 2a in this paper). The piControl multi-model mean (MMM), zonal-average temperature is slightly cooler than observed at high (60–90^∘ N) Northern Hemisphere (NH) latitudes (Fig. 2a). There is a large spread across the models though, with eight of the models simulating colder (up to 4 ^∘C) than observed temperatures and nine of the models simulating warmer (up to 2 ^∘C) than observed temperatures. The piControl MMM, zonal-average temperature is noticeably warmer than observed at high (60–90^∘ S) Southern Hemisphere (SH) latitudes, again with a large spread across the models. Two models – MIROC-ES2L and EC-Earth3-LR – have biases in excess of 5 ^∘C. Hajima et al. (2020) attribute the MIROC-ES2L piControl warm bias over the Southern Ocean to it being mainly associated with the model's representation of cloud radiative processes. The spread of the piControl simulations is smaller at low and midlatitudes (Fig. 2a).

Figure 3 Comparison of the piControl sea ice distributions (a) in the Northern Hemisphere for August–September and (b) in the Southern Hemisphere for February–March. For each 1^∘ × 1^∘ longitude–latitude grid cell, the figure indicates the number of models that simulate at least 15 % of the area covered by sea ice. The observed 15 % concentration boundaries (black lines) are the 1870–1919 CE interval based on the Hadley Centre Sea Ice and Sea Surface Temperature (HadISST; Rayner et al., 2003) dataset. See Figs. S2 and S4 for individual model results.

We adopt the definition of sea ice area of the Sea Ice Model Intercomparison Project (SIMIP; SIMIP Community, 2020), i.e., sea ice concentration times the cell area. The multi-model ensemble of piControl simulations of minimum (August–September) Arctic sea ice distribution (Figs. 3a, S2) show good agreement with the 15 % contour from the HadISST data averaged over the 1870–1920 period (Fig. S1) (Rayner et al., 2003). Two models – FGOALS-g3 and EC-Earth3-LR – show noticeably greater minimum summer sea ice extent in the Nordic Seas as compared to the HadISST period (Fig. S2). Further, evaluation of the piControl simulations can be found in Kageyama et al. (2021). In particular, they find that in comparison to sea ice reconstruction sites, the models generally overestimated sea ice cover at sites close to the sea ice edge.

Figure 4 The simulated Arctic (a, c, e) and Antarctic (b, d, f) annual cycle of sea ice area (10⁶ km²) for the (a, b) PI, (c, d) LIG, and (e, f) LIG minus PI. The monthly mean sea ice areas from the NOAA_OI_v2 dataset for 1982–2001 (Reynolds et al., 2002) are shown in panels (a) and (b).

Figure 4 shows the seasonal cycle of Arctic sea ice area in the piControl simulations for each model and the MMM. These are compared to the NOAA OI_v2 observational dataset, with higher temporal and spatial coverage than the HadISST dataset. The NOAA_OI_v2 dataset (Reynolds et al., 2002), also used in Kageyama et al. (2021), only extends back to 1981. It should be noted that atmospheric CO₂ concentrations had already risen to 340 ppm by 1981, as compared to 284.7 ppm specified in the piControl simulations. We find a large spread across the piControl simulations. The range in March is 12.27 to 19.16 × 10⁶ km² and the MMM is 15.30 ± 1.89 × 10⁶ km². The range in September is 3.56 to 9.73 × 10⁶ km², and the MMM is 6.13 ± 1.66 × 10⁶ km². Generally, those models with less sea ice in March than the MMM also have less sea ice in September than the MMM. Observed estimates of sea ice area from the NOAA-OI_v2 dataset for 1982–2001 are 14.7 × 10⁶ km² for March and 5.1 × 10⁶ km² for September.

The MMM piControl simulations of austral summer minimum (February–March) sea ice distribution around Antarctica, however, show less consensus among the models and less agreement with the HadISST data, with many models significantly underestimating the observed austral summer minimal extent (Figs. 3b, 4b, S4). The range in February is 0.02 to 3.82 × 10⁶ km² and the MMM is 1.65 ± 1.21 × 10⁶ km². Antarctic sea ice melts back largely to the continent's edge in February–March in four models (AWI-ESM-2-1-LR, EC-Earth3-LR, MIROC-ES2L, and MPI-ESM1-2-LR) (Fig. S5). The spread of models is even greater in their simulations of piControl austral winter sea ice area around Antarctica, ranging from 3.27 × 10⁶ km² in September in MIROC-ES2L to over 19 × 10⁶ km² in IPSL-CM6-LR and FGOALS-g3. The September MMM is 17.13 ± 5.21 × 10⁶ km². Observed estimates of sea ice area from the NOAA-OI_v2 dataset for 1982–2001 are 2.7 × 10⁶ km² for February and 16.5 × 10⁶ km² for September.

3.2 Surface temperature responses

The seasonal character of the insolation anomalies results in warming and cooling over the continents in the lig127k ensemble (relative to the piControl) in June–July–August (JJA) and December–January–February (DJF), respectively, except for the African and southeast Asian monsoon regions in JJA. These patterns of seasonal, continental warming and cooling are a robust feature across the models, with more than 70 % of the models agreeing on the sign of the temperature change (Fig. 5a, c).

Figure 5 Multi-model ensemble average changes (a, c, e) and across-ensemble standard deviations (b, d, f) of surface air temperatures (^∘C) for lig127k minus piControl. Shown are June–July–August (a, b), December–January–February (c, d), and annual mean (e, f) changes. Dots indicate where less than 12 (70 %) of the 17 models agree on the sign of the change.

The warming during JJA is greater than 6 ^∘C at midlatitudes in North America and Eurasia (Fig. 5a), though with significant differences in the magnitude of the warming in the southeast US, Europe, and eastern Asia among the models (Fig. 5b). Further investigation of the effects of preindustrial vegetation, including crops, for these regions in the lig127k protocol would be useful (Otto-Bliesner et al., 2020). Subtropical land areas in the Southern Hemisphere (SH) also respond to the positive (but more muted) insolation anomalies, with JJA temperature anomalies more than 2 ^∘C warmer than PI. JJA warming over most of the oceans is a robust feature across the models. This warming is greatest in the North Atlantic and the Southern Ocean, though with large differences across the ensemble of models (Fig. 5b). Cooling over the Sahel and southern India in JJA is associated with the increased cloud cover associated with the enhanced monsoons (see Sect. 3.4).

In response to the negative insolation anomalies at all latitudes (Fig. 1), the lig127k MMM simulates cooling during DJF over the continental regions of both hemispheres and low and midlatitude oceans (Fig. 5c). The largest DJF temperature anomalies occur over southeastern Asia and northern Africa. Ocean memory has been shown to provide the feedback to maintain positive or neutral DJF temperature anomalies in the Arctic and North Atlantic (see Serreze and Barry, 2011, for further discussion). As indicated by the standard deviations of the ensemble changes, large differences in the magnitude of the DJF high-latitude, surface temperature responses and feedbacks exist among the models (Fig. 5d). Understanding these differences warrants further analyses in future studies.

These seasonal patterns of change are similar to those found in Lunt et al. (2013), though the warming is larger in the CMIP6 simulations. It should be noted that the MMM in Lunt et al. (2013) includes simulations that have varying orbital years (between 125 and 130 ka) and greenhouse gas concentrations.

Table 3 Metrics for surface air temperature change (^∘C) for CMIP6–PMIP4 lig127k simulations.

¹ 60–90^∘ N minus 0–30^∘ N. ² 60–90^∘ S minus 0–30^∘ S.

Annually, the MMM surface temperature changes between the lig127k and piControl are generally less than 1 ^∘C over most of the globe, with two exceptions (Fig. 5e): greater negative surface temperature anomalies across the North African and Indian monsoon regions and positive surface temperature anomalies in the Arctic. Although more than 70 % of the models agree on the sign of the changes in these regions, as well as in the Indian sector of the Southern Ocean (Fig. 5e), the across-ensemble standard deviations indicate differences in the magnitudes of the annual surface temperature responses (Fig. 5f). Globally, the MMM change in annual surface air temperature is close to zero (−0.2 ± 0.32 ^∘C), though with a large spread among the models (−0.48 to 0.56 ^∘C) (Table 3). Conclusions about the land versus ocean or NH versus SH annual temperatures changes are complicated by mean changes being close to zero and not consistently positive or negative (Table 3).

The large spread of mean annual surface temperature change among the models in the polar regions (60–90^∘ latitude) is further illustrated in Fig. 2b. Annual Arctic surface temperature changes in the lig127k simulations range from −0.39 to 3.88 ^∘C. The MMM is 0.82 ± 1.20 ^∘C. Notably, EC-Earth3-LR and HadGEM3-GC3.1-LL have anomalies greater than 3 ^∘C in their lig127k simulations as compared to their piControl simulations, while AWI-ESM-1-1-LR and FGOALS-f3-L are cooler in their lig127k simulations as compared to their piControl simulations. The spread (and magnitude) of mean annual temperature change for the SH polar region is less, with 7 of 17 models simulating a modest warming of 0–1 ^∘C and 3 models simulating a cooling of the mean annual surface temperature (Fig. 2b). The MMM is 0.38 ± 0.63 ^∘C. The change in the NH latitudinal gradient is positive from all models: 1.27 ^∘C in the MMM though ranging quite significantly among models for 0.30 ^∘C in FGOALS-f3-L and 3.94 ^∘C in EC-Earth3-LR (Table 3). The change in the SH latitudinal gradient is smaller (0.47 ^∘C in the MMM), reflecting the prescription of a modern Antarctic ice sheet in the lig127k experiment (Table 3). Changes in the size of the Antarctic ice sheet during the Last Interglacial would be expected to result in warming at polar latitudes in the SH and an increase in the SH latitudinal gradient (Bradley et al., 2012; Otto-Bliesner et al., 2013; Stone et al., 2016)

3.3 Sea ice responses

Boreal insolation anomalies at 127 ka enhance the seasonal cycle of Arctic sea ice (Fig. 4c). There is a ∼ 50 % reduction and shift of minimum area in the MMM from 6.1 × 10⁶ km² in August–September for PI to 3.1 × 10⁶ km² in September for lig127k, with a range of 0.22 to 7.47 × 10⁶ km² in the individual lig127k simulations. The lig127k MMM maximum winter sea ice area in the Arctic in March is 15.68 ± 2.08 × 10⁶ km² with a range of 12.27 to 20.28 × 10⁶ km². The INM-CM4-8 and AWI-ESM2-1-LR have small reductions in sea ice area in all seasons with the largest decrease in October (Fig. 4e). HadGEM3-G31-ll and EC-Earth3-LR have large reductions in minimum Arctic sea ice area. HadGEM3-GC31-ll simulates an ice-free Arctic in August–September–October, with the largest reduction in October (Fig. 4c, e). EC-Earth3-LR has the largest reduction of March sea ice area for lig127k as compared to its piControl, and AWI-ESM2-1-LR has a notable increase (Fig. 4e). As shown also in Kageyama et al. (2021), PI biases in simulation of the minimum Arctic sea ice are not always a good predictor of reductions at lig127k (Fig. 4c).

The individual model lig127k minimum (August–September) Arctic sea ice area anomalies show negative correlations (−0.65) with the Arctic (60–90^∘ N) annual surface temperature anomalies from their respective piControl simulations and negative correlation (−0.53) with the corresponding JJA temperature anomalies, both significant at the 0.05 significance level (Fig. 7). Memory in the ocean and cryosphere memory provide feedbacks to maintain positive temperature anomalies, DJF and annually, in the Arctic (Fig. 5). Analyzing the summer atmospheric heat budgets across the models, Kageyama et al. (2021) find that the different Arctic sea ice responses can be related to the sea ice albedo feedback, i.e., phasing of the downward solar insolation changes associated with the orbital forcing and reflected upward shortwave flux changes associated with the sea ice cover changes. As has been done for evaluating simulations of present sea ice distributions, it would be useful for further studies to also explore model differences in the simulated changes in high-latitude cloudiness, boundary layer, winds, and ocean processes (Kattsov and Källén, 2005; Arzel et al., 2006; Chapman and Walsh, 2007).

Previous studies suggest that the mean-ice state in the control climate can influence the magnitude and spatial distribution of warming in the Arctic in future projections (Holland and Bitz, 2003). Thinner Arctic sea ice is more susceptible to summer melting than thicker Arctic sea ice. Arctic sea ice thickness varies substantially across the 1850 CE ensemble, ranging from 1–1.5 m in CNRM-CM6-1 and NESM3 to ∼ 7.5 m in MIROC-ES2L (not shown). No robust relationship to the August–September lig127k minimum Arctic sea ice area anomaly is present. This is also true for the CMIP6–PMIP4 mid Holocene simulations (Brierley et al., 2020). One reason for a lack of any relationship may be the seasonal nature of the lig127k and midHolocene insolation forcings as compared to the annual forcing by greenhouse gas changes in future projections.

The lig127k austral summer sea ice around Antarctica has a minimum in February in the MMM of 1.84 ± 1.42 × 10⁶ km² (Fig. 4d). This is similar to the MMM of the piControl simulations (Fig. 4b). In both the lig127k and piControl, the models exhibit widely different sea ice areas (0.06 to 4.65 × 10⁶ km²) and distributions for the austral summer (Fig. S5). Those models that simulate summer sea ice in the Weddell Sea in the piControl (Fig. S4) retain this sea ice in their lig127k simulation. The maximum austral winter sea ice around Antarctica also varies widely among the models, with the MIROC-ES2L simulating the smallest area (and seasonal cycle) and IPSLCM6 simulating the highest areal extent (and seasonal cycle) (Fig. 4b, d) in the piControl and lig127k simulations. ACCESS-ESM1-5 has the greatest sensitivity to the lig127k forcings (Fig. 4f).

Figure 6 Comparison of the lig127k sea ice distributions (a) in the Northern Hemisphere for August–September and (b) in the Southern Hemisphere for February–March. For each 1^∘ × 1^∘ longitude–latitude grid cell, the figure indicates the number of models that simulate at least 15 % of the area covered by sea ice. The average 15 % concentration boundaries (black lines) are averaged for 2000–2018. See Figs. S3 and S5 for individual model results.

Figure 7 (a) lig127k August–September sea ice NH area anomaly (10⁶ km²) versus lig127k annual 60–90^∘ N surface air temperature anomaly (^∘C); (b) lig127k August–September NH sea ice area anomaly (10⁶ km²) versus lig127k JJA 60–90^∘ N surface air temperature anomaly (^∘C).

The consensus from the lig127k sea ice distributions is a reduced minimum (August–September) summer sea ice extent (defined as 15 % concentration) in the Arctic (Fig. 6) as compared to the piControl simulations (Fig. 3). It is interesting to compare the MMM simulated summer sea ice extents in the lig127k simulations to the observed sea ice extents for 2000–2018 (black lines in Fig. 6). More than half of the models simulate a retreat of the Arctic minimum (August–September) ice edge at 127 ka, similar to the average of the last 2 decades. The pattern of February–March Southern Ocean sea ice extent is broadly similar in the lig127k simulations to 2000–2018, though four models simulate a larger sea ice area.

3.4 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 southeast Asia, and northwestern Mexico/the southwestern US (Fig. 8a). In contrast, summer monsoonal precipitation decreases over South America, southern Africa, and Australia. The spread among models is large, however, as shown by the across-ensemble standard deviations (Fig. 8b, d) and percentage changes in area-averaged precipitation during the monsoon season for seven different regional monsoon domains for the individual lig127k simulations (Fig. 16a). The models generally agree on the sign of the percentage changes in the area-averaged precipitation rate during the monsoon season for the monsoon regions, except for the East Asian, South Asian, and Australian–Maritime Continent monsoons where some models simulate increased monsoonal precipitation whereas others show decreases.

Figure 8 Multi-model ensemble average changes (a, c) and across-ensemble standard deviations (b, d) of precipitation (mm d⁻¹) for lig127k minus piControl. Shown are June–July–August (a, b) and December–January–February (c, d) changes. Dots indicate where less than 12 (70 %) of the 17 models agree on the sign of the change.

Over the tropical Pacific Ocean, reduced DJF precipitation over the Intertropical Convergence Zone (ITCZ) is a robust feature across the ensemble of lig127k simulations (Fig. 8c). The models simulate a shift of the tropical Atlantic ITCZ northward in JJA and southward in DJF, though with significant differences among the models of the ensemble (Fig. 8a, b). Over the Indian Ocean, the ensemble mean indicates more precipitation in DJF over the entire basin and less in JJA, particularly in the central and eastern basin, though again with large standard deviations (Fig. 8).

Figure 9 Ensemble-averaged Last Interglacial change in monsoon-related rainfall rate (in mm d⁻¹). Red and blue contours show the boundaries of lig127k and piControl monsoon domains, respectively, using the definitions of Wang et al. (2011).

Figure 9 shows the ensemble-averaged lig127k change in monsoon-related rainfall rate and global monsoon domain. Increases in the summer rainfall rate and areal extent of the North Africa and East Asia monsoon are clear and are robust across the multi-model ensemble. The spread across the multi-model ensemble is considerable, though, for the North African (NAF) monsoon, with the percentage change in the areal extent varying from ∼ 40 %–120 % (Fig. 16b) and the percentage change in the total amount of water precipitated in each monsoon season varying from ∼ 70 %–140 % (Fig. 16c). 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. 16b) and the percentage change in the total amount of water precipitated in each monsoon season varying from ∼ 25 %–40 % (Fig. 16c). The lig127k and piControl simulations produce more muted changes for the other monsoon regions in the MMM, with regards to the regional monsoon-related rainfall rate and the monsoon domains (Fig. 9). Four models (AWI-ESM-1-1-LR, AWI-ESM-2-1-LR, MPI-ESM1-2-LR, NESM3) in the lig127k ensemble include interactive vegetation. Even then, these four models generally fall within the spread of the models with prescribed vegetation for the three metrics and seven monsoon regions (Fig. 16).

4 Data reconstructions

4.1 Marine temperatures

The lig127k climate model simulations are assessed using two complementary compilations of sea surface temperature (SST) anomalies at 127 ka (Tables S3–S5, S7), which are both individually based on stratigraphically consistent chronologies (Capron et al., 2017; Hoffman et al., 2017).

The multi-archive high-latitude compilation by Capron et al. (2014, 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 five 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 records calculated for 41 sites as the average of the summer and winter records with a model- and observation-consistent 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. (2017) age models. SSTs 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., 2014, 2017; Hoffman 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 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) (Bazin et al., 2013; Veres et al., 2013) in the compilation by Capron et al. (2014, 2017) and a chronology based on millennial-scale variations observed in independently dated Asian speleothem records (Speleo-Age) (Barker et al., 2011) in the compilation by Hoffman et al. (2017). Note that the two reference chronologies diverge by about 1 kyr at 127 ka (Capron et al., 2017).

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 in the temperature reconstruction's method, and both produce 1000 possible chronologies to 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 time series interpolated every 1 kyr. In the compilation of Capron et al. (2014, 2017), the temperature at 127 ka is taken as the median temperature averaged over the 128–126 ka period. Finally, temperature anomalies relative to the preindustrial period are calculated in both cases for marine sites using the HadISST dataset (Rayner et al., 2003), over the intervals 1870–1899 and 1870–1889 CE, in the compilations by Capron et al. (2017) and Hoffman et al. (2017), respectively. For both compilations, the provided 2σ uncertainties integrate errors linked to relative dating and surface temperature reconstruction methods.

Nevertheless, because of (1) the different reference chronologies used, (2) the different tie points and associated relative age uncertainties defined to derive the chronology of each site, and (3) the different calculation methods (Bayesian statistics versus linear interpolation between tie points) used in the Monte Carlo age model analysis of each site (despite apparently relatively similar approaches), the two compilations by Capron et al. (2014, 2017) and by Hoffman et al. (2017) are listed as such in Tables S2–S5, S7. Implications of these methodological differences in the inferred 127 ka values are best illustrated when comparing the surface temperature time series deduced from the two different approaches for a same North Atlantic (62^∘ N) site: at 127 ka, a temperature offset of ∼ 2 ^∘C is observed between the two reconstructions (see Fig. 4 of Capron et al., 2017).

4.2 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, 2017) (Tables S2 and S4). For ice cores, preindustrial conditions are estimated using borehole temperature measurements for Greenland and 1870–1899 CE water isotopic profiles for Antarctica (Capron et al., 2017). 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).

4.3 Terrestrial temperatures

Calibrated, well-dated reconstructions of Last Interglacial temperatures over the continents are quite limited. We have assembled two distinct compilations of continental air temperature reconstructions: a dataset of air temperatures over Europe at 127 ka based on Brewer et al. (2008) and a compilation of peak Last Interglacial summer temperatures reconstructed at Arctic sites from pollen and insect assemblages (Table S6). For both we report anomalies comparing 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 insight into the LIG vegetational development across Europe. A number of pollen-based climate reconstructions based 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.

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 transferring the marine chronology to land sequences for certain pollen stratigraphical events on the basis of joint pollen and paleoceanographic 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 time slice, 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 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-year 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 the 15 terrestrial sites used by Brewer et al. (2008), 8 were excluded due to uncertainties over their chronostratigraphical or chronological assignments or because they did not extend to 127 ka.

The Arctic dataset compiles the most stratigraphically complete, best time-constrained, and calibrated summer temperature reconstructions published from above 65^∘ N latitude. We report the mean of the two 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. This differs from error reporting for the European dataset above. Importantly, the Arctic compilation also differs from the 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 ¹⁴C (non-finite ¹⁴C ages) and in some cases luminescence or tephrochronology. In contrast to the North American Arctic sites, in northern Finland (Sokli) and northeast Russia (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., 2012; Salonen et al., 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.

Despite an abundance of LIG pollen records from Eurasia and various attempts at pollen-based climate reconstructions (e.g., compilations Velichko et al., 2008; Turney and 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 time slice had not been correctly identified. An alternative approach would have been to select peak temperatures from a wider interval (e.g., 127 ± 2 ka) and assume that these are quasi-synchronous. In addition, the Arctic reconstruction may be skewed towards warmer temperatures than the models, given that by definition this reconstruction reports the warmest period from each Arctic site rather than the 127 ka time slice.

4.4 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. (2021). The sea ice records are derived from dinoflagellate cysts, subpolar foraminers, and ostracods.

4.5 Precipitation

Compilations 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 comparison with results from the lig127k simulations. The main patterns that emerge, about precipitation change between the LIG and the preindustrial/recent past, are near-ubiquitous higher LIG annual precipitation over the NH (Fig. 13). Exceptions to this are individual sites in western North Africa, the Levant, northern South America, Borneo, the northwest of the 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 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.

5 Model–data comparisons

5.1 Temperature

Figures 10 to 12 compare the 127 ka temperature reconstructions discussed in Sect. 4 to the MMM and individual models. Details can be found in Tables S2–S7.

Figure 10 High-latitude surface temperature anomaly between 127 ka and the preindustrial from models (ensemble average in colors) and proxies (filled markers): circles for the compilation by Hoffman et al. (2017), squares and diamonds for marine sites and ice cores, respectively, of the compilation by Capron et al. (2014, 2017), pluses for the compilation of Brewer et al. (2008), and triangles for the Arctic compilation. (a, b) June–July–August; (c, d) December–January–February; and (e, f) annual. The preindustrial reference is 1850 CE for model anomalies, and for the data it is 1870–1899 for Capron, Brewer, and Arctic; 1870–1889 for Hoffman.

Figure 11 Same as Fig. 10 but for low-latitude (40^∘ S–40^∘ N) surface temperature.

Figure 12 Comparison of proxy estimates of surface temperature anomalies (± 1 standard deviation) with modeled temperature anomalies at the locations of the proxy data. Annual anomalies for (a) 40–80^∘ N, North Atlantic and Pacific SST and Greenland; (b) 40^∘ S–40^∘ N SST; and (c) 40–80^∘ S SH ocean SST and Antarctic surface air temperature. Seasonal anomalies for (d) 40–70^∘ N, JJA NH oceans; (e) 60–80^∘ N, JJA NH terrestrial; and (f) 40–60^∘ S, DJF SH oceans. All units are degrees Celsius. Data and model values supporting this figure can be found in Tables S2–S7.

NH high-latitude terrestrial temperature proxies for the boreal summer (JJA) match the large warming in the lig127k MMM for most sites (Fig. 10a), except for Lake CF8 on Baffin Island and Wax Lips Lake in northeastern Greenland (Fig. 12e, Table S6). These estimates are from subfossil midges and use published climatic and biogeographic calibration for calculating the mean temperature of the warmest month, rather than JJA, and represent the peak LIG temperatures and not necessarily 127 ka. The only model that simulates the warming reconstructed for these two sites is EC-Earth3-LR, though its lig127k simulation overestimates the warming farther south. Over Europe, the temperature proxies show generally positive anomalies for JJA, but these are often smaller than those of the lig127k MM (Fig. 10a). The lig127k MMM DJF temperatures over North America and Eurasia are significantly colder with respect to PI, except over western Europe (Fig. 10c). The proxies for the latter show a mixed signal. The MMM indicates much warmer surface temperatures in DJF over the Arctic Ocean, Baffin Bay, and Labrador and Greenland Seas, which cannot be evaluated given the available reconstructions (Fig. 10c). Annually, the MMM shows notable warming for Greenland and the ocean surrounding it (Fig. 10e). The range of warming is significant for sites poleward of 50^∘ N (Fig. 11a, Table S2). For the marine sites south of Greenland and near Iceland, the warming simulated by the individual models bracket the proxy estimate. For Greenland, all models are within the 2σ uncertainty for the NEEM ice core.

Over the North Atlantic, the MMM and proxy JJA temperature anomalies are generally in good agreement (Fig. 10a). The exceptions are in the northwestern North Atlantic and the Nordic Seas, where the Capron data suggest significant cooling. This mismatch could be associated with meltwater from potentially remnant ice sheets over Canada and Scandinavia, ice sheets that are not incorporated by the lig127k simulations. EC-Earth3-LR, HadGEM3-GC31-LL, and ACCESS-ESM1-5 simulate the greatest warming at the three northernmost sites (poleward of 68^∘ N) in the Norwegian Sea, with EC-Earth3-LR warming outside the 2σ uncertainty range of the proxy JJA temperature anomalies (Fig. 12d, Table S5).

The marine reconstruction of Capron et al. (2017) provides evidence of significant LIG warm temperature anomalies for the austral summer (DJF) near New Zealand, which is neither exhibited by the lig127k MMM (Fig. 10d) nor the individual models which all cluster around little or no change in DJF temperature change (Fig. 12f, Table S7). This discrepancy suggests regional circulation changes not resolved by the models. The multi-model ensemble indicates austral winter (JJA) warming over the Southern Ocean and Antarctica, but the lack of proxies does not allow an assessment (Fig. 10b). The simulated annual temperature anomalies for the Antarctic ice cores are cooler than the reconstructed values but generally fall within the 2σ uncertainties (Figs. 10f, 12c, Table S4).

At lower latitudes (40^∘ S–40^∘ N), marine proxy data from the Hoffman reconstruction are available (Fig. 11). They generally correspond with the MMM changes. The SST proxies from the tropical Atlantic match the colder MMM lig127k SSTs in DJF (Fig. 11b). The reconstructed cooling there in JJA is not captured in the MMM, leading to a failure to also capture the annual mean signal (Fig. 11a, c). Proxy indications of much warmer SSTs in the upwelling regions off the west coasts of southern Africa, North America, and South America are not simulated by the models (Figs. 11, 12b, Table S3). The resolution of CMIP6 models is generally not adequate to properly simulate these narrow coastal upwelling regions.

5.2 Precipitation

As shown in a comparison with a smaller ensemble of 127 ka simulations (Scussolini et al., 2019), precipitation proxies from the global compilation largely match the annual precipitation from the models included in the lig127k MMM (Fig. 13a), with the overall hit rate comparing matches between the sign of the anomaly in the models and in the proxies of 65 % (Fig. 13b). The agreement between the MMM and NH proxies is even higher over North Africa–the Middle East (hit rate of 76 %), North America–Greenland (hit rate of 78 %), and South Asia (hit rate of 73 %). It should be noted that the range across the individual model is quite large for North America–Greenland (hit rates of 45 % to 90 %) and South Asia (hit rates of 40 % to 87 %). Proxies and MMM weakly disagree over much of Europe, central Asia, and the region between them, where proxies indicate wetter LIG conditions or no change, and the MMM indicates somewhat drier conditions or no change (Fig. 13a). The overall MMM hit rate for Europe (68 %) is much improved as compared to the smaller ensemble analyzed by Scussolini et al. (2019), but the range across the models is quite large (36 % to 77 %). Other instances of more regional disagreements 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 MMM. 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. The agreement between the MMM and SH proxies is noticeably less than for the NH, with hit rates of less than 50 % except for South America with a hit rate of 89 % (Fig. 13b). 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 the MMM indicate wetter and drier LIG, respectively (Fig. 13a). The hit rates for individual models show that some models perform significantly better over Australia (Fig. 13b).

Figure 13 (a) Annual precipitation anomaly between 127 ka and PI, from model simulations (ensemble average in contoured colors) and from proxies (filled markers). Green colors indicate higher precipitation at 127 ka from models or proxies, and the opposite is true for brown colors. Proxy anomalies are on a semi-quantitative scale: dark green (much wetter LIG), light green (wetter), white (no noticeable anomaly), light brown (drier), and dark brown (much drier). Different markers represent different types of proxy records as specified in the legend. Proxy reconstruction from Scussolini et al. (2019). (b) Annual precipitation anomaly between LIG and PI; comparison of models and proxies. The hit rate is the percentage of matches between the sign of the anomaly in the models and in the proxies. N is the number of model–proxy comparisons per region.

6 Comparison of the model sensitivities to the insolation anomalies at 127 and 6 ka

The large-scale features and evaluation of the CMIP6–PMIP4 midHolocene simulations in comparison to data 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 ice to the different magnitudes and seasonal character of the insolation anomalies at 127 ka versus 6 ka.

6.1 Orbital forcing

The orbit at 6 ka was characterized by a smaller eccentricity than at 127 ka, similar to 1850 CE (Fig. 14). Perihelion occurred near the boreal autumn equinox as compared to close to the boreal summer solstice at 127 ka and near aphelion at 1850 CE. NH summer insolation anomalies at 6 ka, ∼ 5 %–10 % greater than at 1850 CE, are considerably less than at 127 ka (Figs. 1 and 14). In addition, the positive insolation anomalies of greater than 10 % in the Arctic occur in July–August at 6 ka 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 127 ka than at 6 ka.

Figure 14 (a) Orbital configuration for the midHolocene (6 ka) experiment. (b) Latitude–month insolation anomalies 6 ka–1850 as percentage change from PI.

6.2 Surface temperature responses

Figure 15 compares the MMM changes and standard deviations of ensemble changes of surface air temperatures for lig127k and midHolocene simulations. In the tropics and the Southern Hemisphere, the JJA zonal-average temperature anomaly is positive ($\sim +$0.5 ^∘C) for the lig127k ensemble but negative ($\sim -$0.5 ^∘C) in the midHolocene ensemble. The maximum JJA surface temperature anomalies occur at ∼ 40–65^∘ N for both time periods but are significantly greater at 127 ka (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 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 MMM and up to 3 ^∘C higher in the lig127k MMM than the piControl.

Figure 15 Multi-model ensemble mean and 2 standard deviation, zonal surface air temperature anomalies (^∘C) for midHolocene and lig127k simulations for JJA and DJF (see Brierley et al., 2020, for more details on midHolocene simulations). Note that 14 models completed both the midHolocene and lig127k experiments. Three models (ACCESS-ESM1-5, AWI-ESM-2-1-LR, CNRM-CM6-1) completed only the lig127k experiment, while three models (MRI-ESM2-0, UofT-CCSM-4, BCC-CSM1-2) completed only the mid-Holocene experiment.

6.3 Precipitation responses

The signs of the percentage changes in the areal extents of the regional monsoon domain (Fig. 16b) and the percentage changes in the total amount of water precipitated in each monsoon season (Fig. 16c) are similar for the lig127k and midHolocene simulations as compared to piControl simulations, but the responses are generally enhanced in the lig127k simulations as compared to the midHolocene simulations. Both time periods show greater areal extent and the total amount of water precipitated for the NAF and EAS monsoons, with the lig127k MMMs outside the midHolocene quartile range. Similarly, the Australian–Maritime Continent, southern Africa, and South America monsoons show greater reductions in areal extent and total water precipitated in the lig127k simulations than in the midHolocene simulations as compared to the piControl simulations. Both time periods show a mixed simulated response of the North American monsoon (NAMS).

Figure 16 Relative changes in MMM and individual lig127k monsoons. Three different monsoon diagnostics as computed for each of seven different regional domains for the individual CMIP6 lig127k simulations. The comparable results from the midHolocene simulations are shown with boxes and whiskers (for details, see Brierley et al., 2020). (a) The percentage changes in area-averaged precipitation rate during the monsoon season; (b) the percentage change in the areal extent of the regional monsoon domain; (c) the percentage change in the total amount of water precipitated in each monsoon season (computed as the precipitation rate multiplied by the areal extent). The abbreviations used to identify each regional domain are as follows: North America Monsoon System (NAMS), North Africa (NAF), southern Asia (SAS), and East Asia (EAS) in the Northern Hemisphere and South America Monsoon System (SAMS), southern Africa (SAF), and Australian–Maritime Continent (AUSMC) in the Southern Hemisphere. Note that 14 models completed both the midHolocene and lig127k experiments. Three models (ACCESS-ESM1-5, AWI-ESM-2-1-LR, CNRM-CM6-1) completed only the lig127k experiment, while three models (MRI-ESM2-0, UofT-CCSM-4, BCC-CSM1-2) completed only the mid-Holocene experiment.

6.4 Arctic sea ice responses

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. Similar to the analyses of the ensemble of lig127k simulations (see Sect. 3.3), Brierley et al. (2020) also found that in the ensemble of CMIP6 midHolocene simulations, the summer sea ice reduction in the Arctic is correlated to the magnitude of annual warming over the Arctic but has little Arctic-wide relationship with the simulated PI sea ice extents.

Figure 17 The Arctic annual cycle of area of sea ice greater than 15 % (10⁶ km²) for the (a) PI, (b) MH for a subset of the models. Note that 14 models completed both the midHolocene and lig127k experiments. Three models: ACCESS-ESM1-5, AWI-ESM-2-1-LR, CNRM-CM6-1 completed only the lig127k experiment.

7 Discussion

The Tier 1 lig127k experiment was designed to address the climate responses to stronger orbital forcing (relative to 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 ones. 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. The 17 CMIP6 models that have completed the lig127k experiment are presented.

The CMIP6–PMIP4 lig127k simulations show warming and cooling over the continents during JJA and DJF, respectively, in response to the seasonal character of the insolation anomalies. The JJA MMM warming is greater than 6 ^∘C at midlatitudes in North America and Eurasia, though with across-ensemble standard deviations in excess of 2 ^∘C over the eastern US and central Europe. The simulations exhibit a ∼ 50 % reduction and shift of Arctic minimum summer sea ice area to 3.1 × 10⁶ km² in September for lig127k, though with a large range of 0.22 to 7.47 × 10⁶ km². Positive temperature anomalies are present in the lig127k simulations annually in the Arctic and over the Southern Ocean, though with across-ensemble standard deviations in excess of 2 ^∘C. The large spread across the models in their simulations of Arctic sea ice, even now with all models adopting a common experimental protocol, points to the need to better diagnose the atmosphere and ocean feedbacks that differ across the lig127k ensemble (Kageyama et al., 2021). As 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/the southwestern US. 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, with the percentage change in its areal extent ranging from less than 50 % to more than 150 % and total amount of water precipitated during the monsoon season ranging from ∼ 65 % to more than 200 %. The four models with interactive vegetation fall within the spread of the models with prescribed vegetation for the three metrics and seven monsoon regions. The lig127k individual monsoon changes are of a similar sign, but a greater magnitude, to those seen in the midHolocene simulations (Brierley et al., 2020).

New syntheses for surface temperature and precipitation, targeted for 127 ka, have been developed for comparison to the CMIP6–PMIP4 lig127k simulations. The surface temperature reconstructions include two complementary 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 based on pollen records, and peak LIG summer temperatures in the Arctic inferred from pollen and insect assemblages. Anomalies were consistently computed comparing the reconstructed temperatures with observationally based preindustrial climate estimates from the end of the 19th century. A new precipitation reconstruction has expanded from previous regional compilations to now near-global coverage.

Over Canada, Scandinavia, parts of midlatitude Europe, and the North Atlantic, the proxy and lig127k positive JJA temperature anomalies are in good agreement. The exceptions are in the northwestern North Atlantic and the Nordic Seas, where the Capron reconstruction (Capron et al., 2017) suggest significant cooling. The 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 regional ocean circulation changes not resolved by the models as well as meltwater from potential remnant ice sheets over Canada and Scandinavia as well as memory in the ocean of the H11 event (Govin et al., 2012; Marino et al., 2015), which the lig127k simulations do not incorporate. The latter would lead to cooling in the North Atlantic and warming in the Southern Ocean (Stone et al., 2016; Holloway et al., 2018).

The simulated annual temperature anomalies for the Greenland and Antarctic ice cores are cooler than the reconstructed values but generally fall within the 2σ uncertainties. The lig127k 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) demonstrated that the distinctive peak in δ¹⁸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, 2018). The CMIP6–PMIP4 Tier 2 LIG 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., 2017).

Other reasons for mismatches between the models and the reconstructions for temperature and precipitation will also be explored with CMIP6–PMIP4 Tier 2 LIG experiments. 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 evapotranspiration on the surface 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). Recent results show large changes in hydrology, with, e.g., the possible existence of an extensive river network across the Sahel and Sahara, therefore also pointing to the need to prescribe or model vegetation changes to capture regional feedbacks (Scussolini et al., 2020). 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).

Temperature reconstructions are not available for many regions where the lig127k multi-model ensemble shows interesting responses to the lig127k forcing. These include the polar regions during their respective winter seasons: Arctic and North Atlantic oceans in DJF and the Southern Ocean and Antarctica in JJA. Development of terrestrial reconstructions for most continents and marine reconstructions for the Indian and Pacific oceans would be useful for assessing the model responses.

The CMIP6–PMIP4 lig127k experiment has potential implications for confidence in future projections. More than half of the models simulate a retreat of the Arctic minimum (August–September) ice edge in their lig127k simulations that is similar to the average of the last 2 decades (Fig. 6). ECS (Table 2) and simulation of August–September lig127k minimum Arctic sea ice area across the models show a significant (at the 0.5 level) correlation of −0.62 (Fig. 18). INM-CM4-8 with the smallest ECS of 1.8 ^∘C simulates the largest August–September lig127k Arctic sea ice area. CESM2 has a high ECS of 5.2 ^∘C (Gettelman et al., 2019); HadGEM3 similarly has a high ECS of 5.6 ^∘C (Guarino et al., 2020). Both predict an almost ice-free or ice-free Arctic in their lig127k experiments. Their predicted years of disappearance of September sea ice in the SSP8-8.5 scenario is 2038 and 2035, respectively (Guarino et al., 2020). Across CMIP6 models, Kageyama et al. (2021) noted a nearly linear relationship between the simulations of Arctic summer sea ice in their 1pctCO2 simulations at the time of doubling and their lig127k simulations. 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., 2021), it is currently difficult to rule out the high or low values of ECS from the proxy data.

Figure 18 lig127k August–September NH sea ice area (10⁶ km²) versus equilibrium climate sensitivity (ECS, K).

Radiative perturbations on the Arctic system, even though related to summer insolation during the LIG and MH rather than greenhouse gas radiative forcing, might provide useful insights on the state of the future Arctic system (Schmidt et al., 2014). Using CMIP5 MH and RCP4.5 simulations from 10 climate models, Yoshimori and Suzuki (2019) examined the relevance of Arctic warming in the MH to that in the future. The radiative forcing in the RCP4.5 experiment is dominated by the elevated atmospheric CO₂ concentrations and is relatively uniform globally and seasonally. The radiative forcing in the MH associated with orbital forcing is seasonal, peaking in July–August. Yet for MH and RCP4.5, the largest Arctic warming and sea ice reduction occurs in late summer and early autumn. The surface energy balance analysis identifies local Arctic feedbacks associated with positive albedo feedback in summer and a consequent increase in heat release from the ocean to the atmosphere in autumn to be important contributors for both climate states.

Large differences exist among the models in the magnitude of the seasonal and annual surface temperature responses in the polar regions reflecting differences in the feedback processes represented by each model. These should be investigated. 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 corresponding future projection simulations.