A multi-model CMIP6-PMIP4 study of Arctic sea ice at 127 ka: sea ice data compilation and model differences

The Last Interglacial period (LIG) is a period with increased summer insolation at high northern latitudes, which results in strong changes in the terrestrial and marine cryosphere. Understanding the mechanisms for this response via climate modelling and comparing the models' representation of climate reconstructions is one of the objectives set up by the Paleoclimate Modelling Intercomparison Project for its contribution to the sixth phase of the Coupled Model Intercomparison Project. Here we analyse the results from 16 climate models in terms of Arctic sea ice. The multi-model mean reduction in minimum sea ice area from the pre industrial period (PI) to the LIG reaches 50 % (multi-model mean LIG area is 3.20×10⁶ km², compared to 6.46×10⁶ km² for the PI). On the other hand, there is little change for the maximum sea ice area (which is 15–16×10⁶ km² for both the PI and the LIG. To evaluate the model results we synthesise LIG sea ice data from marine cores collected in the Arctic Ocean, Nordic Seas and northern North Atlantic. The reconstructions for the northern North Atlantic show year-round ice-free conditions, and most models yield results in agreement with these reconstructions. Model–data disagreement appear for the sites in the Nordic Seas close to Greenland and at the edge of the Arctic Ocean. The northernmost site with good chronology, for which a sea ice concentration larger than 75 % is reconstructed even in summer, discriminates those models which simulate too little sea ice. However, the remaining models appear to simulate too much sea ice over the two sites south of the northernmost one, for which the reconstructed sea ice cover is seasonal. Hence models either underestimate or overestimate sea ice cover for the LIG, and their bias does not appear to be related to their bias for the pre-industrial period. Drivers for the inter-model differences are different phasing of the up and down short-wave anomalies over the Arctic Ocean, which are associated with differences in model albedo; possible cloud property differences, in terms of optical depth; and LIG ocean circulation changes which occur for some, but not all, LIG simulations. Finally, we note that inter-comparisons between the LIG simulations and simulations for future climate with moderate (1 % yr⁻¹) CO₂ increase show a relationship between LIG sea ice and sea ice simulated under CO₂ increase around the years of doubling CO₂. The LIG may therefore yield insight into likely 21st century Arctic sea ice changes using these LIG simulations.

1 Introduction

The Last Interglacial period (LIG) was the last time global temperature was substantially higher than the pre-industrial period (PI) at high northern latitudes. It is important in helping us understand warm-climate sea ice and climate dynamics (Otto-Bliesner et al., 2013, 2017; Capron et al., 2017; Fischer et al., 2018). Stronger LIG spring and summertime insolation contributed to this warmth, as well as feedbacks amplifying the initial insolation signal, in particular feedbacks related to the marine and land cryosphere. Previous climate model simulations of the LIG, forced by appropriate greenhouse gas (GHG) and orbital changes, have failed to capture the observed high temperatures at higher latitudes (Malmierca-Vallet et al., 2018; Masson-Delmotte et al., 2011; Otto-Bliesner et al., 2013; Lunt et al., 2013). Models used during the previous Coupled Model Intercomparison Project 5 (CMIP5) disagree on the magnitude of Arctic sea ice retreat during the LIG: the diversity of sea ice behaviour across models was linked to the spread in simulated surface temperatures and in the magnitude of the polar amplification (Otto-Bliesner et al., 2013; Lunt et al., 2013; IPCC, 2013). However it was difficult to compare some of the LIG simulations because they were not all run using identical protocol. These studies thus highlighted the need of a systematic approach to study the role of Arctic sea ice changes during the LIG.

Coupled Model Intercomparison Projects (CMIPs) coordinate and design climate model protocols for past, present and future climates and have become an indispensable tool to facilitate our understanding of climate change (IPCC, 2007, 2013; Eyring et al., 2016). The Paleoclimate Model Intercomparison Project 4 (PMIP4) is one of the individual Model Intercomparison Projects that is taking part in CMIP6 (Kageyama et al., 2018). Within this framework, a common experimental protocol for LIG climate simulation was developed by Otto-Bliesner et al. (2017). CMIP models differ among each other in their physical formulation, numerical discretisation and code implementation. However, this CMIP6-PMIP4 LIG standard protocol facilitates model inter-comparison work.

Alongside a previous lack of a common experimental protocol, our ability to evaluate CMIP models has previously been hindered by difficulties in determining LIG sea ice extent from marine core evidence (e.g. Otto-Bliesner et al., 2013; Sime et al., 2013; Malmierca-Vallet et al., 2018; Stein et al., 2017). Planktonic foraminifer assemblages that include a subpolar component suggest reduced sea ice in the Arctic Ocean (Nørgaard-Pedersen et al., 2007; Adler et al., 2009). Microfauna found in LIG marine sediments recovered from the Beaufort Sea Shelf, an area characterised by ice-free conditions during summers today, also support ice-free conditions during those times; this indicates that more saline Atlantic water was present on the Beaufort Shelf, suggesting reduced perennial Arctic sea ice during some part of the LIG (Brigham-Grette and Hopkins, 1995). On the other hand, a reconstruction of LIG Arctic sea ice changes based on sea ice biomarker proxies (see below for details) suggests that the central part of the LIG Arctic Ocean remained covered by ice throughout the year, while a significant reduction of LIG sea ice occurred across the Barents Sea continental margin (Stein et al., 2017). On the modelling side, no previous coupled climate model has simulated an ice-free Arctic during the LIG (Otto-Bliesner et al., 2006; Lunt et al., 2013; Otto-Bliesner et al., 2013; Stein et al., 2017).

Here we address the question of LIG Arctic sea ice by providing a new marine core synthesis. Additionally, the CMIP6-PMIP4 LIG experimental protocol developed by Otto-Bliesner et al. (2017) provides the systematic framework to enable us to examine the question of the simulation of LIG Arctic sea ice using a multi-model approach. This is important given the current level of interest in the ability of climate models to accurately represent key Arctic climate processes during warm periods, including sea ice formation and melting. We compare the LIG Arctic sea ice simulated by each model against our new data synthesis and investigate why different models show different Arctic sea ice behaviour.

2 Materials and methods

2.1 Current Arctic sea ice

Our main objective is to investigate LIG sea ice. However, a quick assessment of the sea ice simulated in the reference state, i.e. the pre-industrial control experiment (referred to as piControl in the CMIP6 terminology, and PI in this paper) was necessary. In the absence of extensive sea ice data for the PI, we used data for a recent period before the current sea ice cover significant decrease. We use the NOAA Optimum Interpolation version 2 data (Reynolds et al., 2002) for the period 1982 to 2001. The sea ice data in this dataset are obtained from different satellite and in situ observations. We have used the monthly time series at a resolution of 1^∘. This dataset is termed “NOAA_OI_v2” in the rest of this paper.

2.2 Marine records of LIG Arctic sea ice

We focus here on records of sea ice from marine cores. Table 1 provides a summary of LIG sea ice information and data obtained from marine sediment cores collected in the Arctic Ocean, Nordic Seas and northern North Atlantic. South of 78^∘ N, the records show ice-free conditions. Most of these sea ice records are derived from quantitative estimates of sea surface parameters based on dinoflagellate cysts (dinocysts). North of 78^∘ N the sea-ice-related records are rare and different types of indicators were used. In addition to dinocysts, the records are based on biomarkers linked to phototrophic productivity in sea ice and on foraminifers and ostracods that both provide indication on water properties and indirectly on sea ice (de Vernal et al., 2013b). Between 78 and 87^∘ N, the faunal data have been interpreted as indicating seasonal sea ice cover conditions during the LIG.

Cronin et al. (2010)Stein et al. (2017)Cronin et al. (2010)Stein et al. (2017)Nørgaard-Pedersen et al. (2007)Kremer et al.2018bMatthiessen et al.2001Matthiessen and Knies (2001)Stein et al. (2017)Stein et al. (2017)Adler et al. (2009)Cronin et al. (2010)Kremer et al.2018aCronin et al. (2010)20112013Zhuravleva et al.2017201120082008Eynaud (1999)Eynaud et al.2004Hillaire-Marcel et al.200120082011Hodell et al. (2009)Hillaire-Marcel et al. (2011)Penaud et al. (2008)Eynaud et al.2000

Table 1 Marine core records of Arctic sea ice from MIS5e. The references indicated for the dinocyst reconstructions are those for the initial core, the reconstruction itself follows de Vernal et al. (2013a, b, 2020) (see the main text for details). SIC stands for sea ice concentration. Question marks indicate that seasonal duration or annual mean SIC are uncertain and not quantitatively estimated. More general qualitative statements are still possible and are given in the “Qualitative sea ice state” column.

Among sea ice cover indicators, dinocyst assemblages have been used as quantitative proxy based on the application of the modern analogue technique applied to a standardised reference modern data base developed from surface sediment samples collected at middle to high latitudes of the Northern Hemisphere (de Vernal et al., 2005b, 2013b, 2020). The sea ice estimates from dinocysts used here are from different studies (see references in Table 1) and reconstructions based the new database, including 71 taxa and 1968 stations (de Vernal et al., 2020). The reference sea ice data used for calibration are the monthly 1955–2012 average of the National Snow and Ice Data Center (NSIDC) (Walsh et al., 2016). The results are expressed in term of annual mean of sea ice cover concentration or as the number of months with > 50 % of sea ice. The error of prediction for sea ice concentration is ±12 % and that of sea ice cover duration through the year is ±1.5 months yr⁻¹. Such values are very close to the interannual variability in areas occupied by seasonal sea ice cover (see de Vernal et al., 2013b).

Our biomarker approach for sea ice reconstruction is based on the determination of a highly branched isoprenoid (HBI) with 25 carbons (C25 HBI monoene = IP25) (Belt et al., 2007). This biomarker is only biosynthesised by specific diatoms living in the Arctic sea ice (Brown et al., 2014), meaning the presence of IP25 in the sediments is a direct proof for the presence of past Arctic sea ice. Meanwhile, this biomarker approach has been used successfully in numerous studies dealing with the reconstruction of past Arctic sea ice conditions during the late Miocene to Holocene (for a review, see Belt, 2018). By combining the sea ice proxy IP25 with (biomarker) proxies for open-water (phytoplankton productivity such as brassicasterol, dinosterol or a specific tri-unsaturated HBI, HBI-III), the so-called PIP25 index has been developed (Müller et al., 2011; Belt et al., 2015; Smik et al., 2016). Based on a comparison (“calibration”) PIP25 data obtained from surface sediments with modern satellite-derived (spring) sea ice concentration maps (Müller et al., 2011; Xiao et al., 2015; Smik et al., 2016), the PIP25 approach may allow a more semi-quantitative reconstruction of present and past Arctic Ocean sea ice conditions from marine sediments, i.e. estimates of spring sea ice concentration (or in the Central Arctic probably more the summer situation due to light limitations for algae growth in the other seasons). Based on these data, one may separate “permanent to extended sea ice cover” (> 0.75) and “seasonal sea ice cover”; (0.75–0.1), perhaps including the sub-groups “ice-edge” (0.75–0.5) and “less/reduced sea ice” (0.5–0.1), and “ice-free” (< 0.1). For pros and cons of this approach, we refer to a recent review by Belt (2018).

Figure 1 Map showing the location of the LIG Arctic sediment cores listed in Table 1. Open symbols correspond to records with uncertain chronology, and filled symbols correspond to records with good chronology. The map background has been created using http://visibleearth.nasa.gov (last access: 1 January 2020).

Based on several IP25/PIP25 records obtained from central Arctic Ocean sediment cores (see Fig. 1 for core locations and Table 1 for data), perennial sea ice cover probably existed during the LIG in the Central Arctic, whereas along the Barents Sea continental margin, influenced by the inflow of warm Atlantic Water, sea ice was significantly reduced (Stein et al., 2017). However, Stein et al. (2017) emphasise that the PIP25 records obtained from the central Arctic Ocean cores indicating a perennial sea ice cover have to be interpreted cautiously as the biomarker concentrations are very low to absent (see Belt, 2018 for further discussion). The productivity of algal material (ice and open water) must have been quite low, so that (almost) nothing reached the seafloor or is preserved in the sediments, and there must have been periods during the LIG when some open-water conditions occurred, since subpolar foraminifers and coccoliths were found in core PS51/038 and PS2200 (Stein et al., 2017). It is however unclear whether these periods equate to more than 1 month yr⁻¹ of open water (or seasonal ice conditions). This explains why some sites show both seasonal and perennial interpretations at the same site. The reader is referred to the original publications (Table 1) for more information on these data. Furthermore and importantly, a new revised ²³⁰Th chronology of late Quaternary sequences from the central Arctic Ocean (Hillaire-Marcel et al., 2017) questions the age model of some of the data listed in Table 1. Thus, further verification of age control is still needed and the data from the central Arctic Ocean should be interpreted with caution. We have therefore marked the chronological control as “uncertain” for these cores, while the chronological control is good for cores outside the Central Arctic.

The information given by the different types of sea ice indicators shows that care should be taken when comparing them with model results. We have used the qualitative information given in Table 1, taking into account the threshold given in this table. Indeed, for instance, “perennial sea ice cover” does not automatically mean 100 % sea ice cover, or a sea ice concentration (SIC) of 1.0. It means rather that there is sea ice but not necessarily at a concentration of 100 % over the core site throughout the year (i.e. the summer season is not totally ice-free). Most qualitative reconstructions cite a threshold of 75 %, which we have therefore used in our model-data comparison. We have also used the quantitative mean annual sea ice reconstructions. Finally, similar to studies for future climate, we have considered the Arctic to be ice-free when, in any given month, the total area of sea ice is less than 1×10⁶ km². This means that some marine core sites could remain ice covered for the summer, but the Arctic would nevertheless remain technically ice-free.

2.3 CMIP6-PMIP4 models

The last Coupled Model Intercomparison Project Phase 5 (CMIP5) collected climate simulations performed with 60 different numerical models by 26 research institutes around the world (IPCC, 2013). The follow-on CMIP6 archive, to be completed in 2020, is expected to gather model outputs from over 30 research institutes. Of these, currently 15 models have run the CMIP6-PMIP4 LIG simulation (Table 2). We present results here from all these models.

Ziehn et al.Sidorenko et al. (2015)Sidorenko et al. (2019, 2015)Danabasoglu et al. (2020)Voldoire et al.Hazeleger et al. (2012)Li et al. (2020)Zheng et al.Kelley et al.Williams et al.Volodin et al.Boucher et al.Goosse et al.Hajima et al.20182020Cao et al.Guo et al.Seland et al.

Table 2 Overview of models that have run the CMIP6-PMIP4 LIG simulation. For each model, denomination, physical core components, horizontal and vertical grid specifications, details on prescribed vs. interactive boundary conditions, reference publication, and LIG simulation length are shown.

Table 2 provides an overview of the models used in this study. They are state-of-the-art coupled general circulation models (GCMs) and Earth System Models (ESMs) simulating the atmosphere, ocean, sea ice and land surface processes dynamics with varying degrees of complexity. These 15 CMIP6-PMIP4 models have been developed for several years by individual institutes across the world and, in the context of CMIP6, are used in the same configuration to seamlessly simulate past, present and future climate. We have added the results from the LOVECLIM Earth System Model of Intermediate Complexity, which can be used for longer simulations.

Table 2 shows the following qualities for each model: model denomination, physical core components, horizontal and vertical grid specifications, details on prescribed vs. interactive boundary conditions, relative publication for an in-depth model description, and LIG simulation length (spin-up and production runs).

2.4 PMIP4 LIG (lig127k) simulation protocol

Results shown here are from the main Tier 1 LIG simulation, from the standard CMIP6-PMIP4 LIG experimental protocol (Otto-Bliesner et al., 2017). The prescribed LIG (lig127k) protocol differs from the CMIP6 pre-industrial (PI) simulation protocol in astronomical parameters and the atmospheric trace greenhouse gas concentrations (GHG). LIG astronomical parameters are prescribed according to Berger and Loutre (1991), and atmospheric trace GHG concentrations are based on ice core measurements. Table 3 from Otto-Bliesner et al. (2017) summarises the protocol. All models followed this protocol, except CNRM-CM6-1 for which the most important forcings for the LIG, i.e. the astronomical parameters, have been imposed at the recommended values, but the GHG have been kept at their pre-industrial values of 284.3170 ppm for CO₂, 808.2490 ppb for CH₄ and 273.0211 ppb for N₂O. All other boundary conditions, including solar activity, ice sheets, aerosol emissions, etc., are identical to PI protocol. Both the Greenland and Antarctica ice sheets are known to have shrunk during the interglacial with different timings, and therefore taking PI characteristics for the lig127k protocol is an approximation, particularly for the Antarctic ice sheet, which was possibly smaller than in the PI at that time (Otto-Bliesner et al., 2017). The Greenland ice sheet likely reached a minimum at around 120 ka and was probably still close to its PI size at 127 ka. Given the dating uncertainties and the difficulty for models to include the largest changes in ice sheets for 127 ka, i.e. changes in West Antarctica, the choice of the PMIIP4 working group on interglacials was to use the PI ice sheets as boundary conditions for the Tier 1 PMIP4-CMIP6 experiments presented here and to foster sensitivity experiments to ice sheet characteristics at a later stage. In terms of the Greenland ice sheet, the approximation is considered quite good and ideal for starting transient experiments through the whole interglacial period.

Table 3 Astronomical parameters and atmospheric trace gas concentrations used to force LIG and PI simulations.

LIG simulations were initialised either from a previous LIG run, or from the standard CMIP6 protocol pre-industrial simulations, using constant 1850 GHGs, ozone, solar, tropospheric aerosol, stratospheric volcanic aerosol and land use forcing.

Although PI and LIG spin-ups vary between the models, most model groups aimed to allow the land and oceanic masses to attain approximate steady state i.e. to reach atmospheric equilibrium and to achieve an upper-oceanic equilibrium. LIG production runs are all between 100–200 years long, which is generally within the appropriate length for Arctic sea ice analysis (Guarino et al., 2020).

The LIG orbital parameters result in modifications of the definitions of the months and seasons (in terms of start and end dates within a year). Since daily data was not available for all models to re-compute LIG-specific monthly averages, we have corrected these averages using the method of Bartlein and Shafer (2019). Unless otherwise specified, we use these results adjusted for the LIG calendar throughout this paper.

2.5 The CMIP6 1pctCO2 protocol

We compare the response to the lig127k forcings to idealised forcings for future climate. We have chosen to use the 1pctCO2 simulation from the CMIP6 DECK (Diagnostic, Evaluation and Characterization of Klima Eyring et al., 2016). These simulations start from the PI (piControl) experiment and the atmospheric CO₂ concentration is gradually increased by 1 % yr⁻¹ for at least 150 years, i.e. 10 years after atmospheric CO₂ quadrupling.

3 Results: simulated Arctic sea ice

Since all LIG production runs are at least 100 years in length, all model results are averaged over at least 100 years. We refer to the multi-model mean throughout as the MMM. We consider both the sea ice area (SIA), defined as the sum, over all Northern Hemisphere ocean cells, of the sea ice concentration × the cell area and the sea ice extent (SIE), defined as the sum of the areas of ocean cells where the sea ice concentration is larger than 0.15. Both quantities are used in sea ice studies, SIE has been used widely in IPCC AR5 (Vaughan et al., 2013), while SIA tends to be used more for CMIP6 analyses (e.g. SIMIP Community, 2020).

3.1 PI sea ice

For the present-day we have satellite and in situ observations with which to evaluate the models. The use of present-day sea ice data implies that we might expect the simulated PI sea ice to be generally somewhat larger than the observed mean. Indeed the atmospheric CO₂ levels for the years for which we chose the observation dataset (1982 to 2001) were between 340 and 370 ppm, compared to the PI level of 280 ppm. Figure 2 shows the mean seasonal cycle of the Arctic sea ice extent simulated for the PI and LIG alongside the observed Arctic sea ice extent.

Figure 2 Mean seasonal cycle of the Arctic sea ice area (SIA, left-hand side) and sea ice extent (SIE, right-hand side), in 10⁶ km², simulated for the PI and LIG periods by the PMIP4 models. The top row shows the results for PI. The grey shading shows the monthly minimum and maximum in the SIA and SIE observed over the years 1982–2001, as given by the NOAA_OI_v2 dataset. The second and fourth row show the LIG results, with no calendar adjustment and with calendar adjustment, respectively. The third and bottom row show the corresponding LIG–PI anomalies, with no calendar adjustment and with calendar adjustment, respectively.

Table 4 Sea ice area and extent (in 10⁶ km²) for the PI and LIG simulations (calendar-adjusted values). MMM stands for the multi-model mean, SD for the multi-model standard deviation. NA stands for not available.

The summer minimum monthly MMM SIA for the PI is 6.46 ± 1.41×10⁶ km², compared to the observed 1981 to 2002 mean of 5.65×10⁶ km². In terms of SIE, the summer minimum for PI is 8.89 ± 1.41×10⁶ km², to be compared to the observed 7.73×10⁶ km². Interestingly this MMM PI area and extent is a little larger than the 1981–2002 area. The majority of the simulations show a realistic representation of the geographical extent for the summer minimum (Fig. 3, Table 4), with 9 out of 16 models showing a slightly smaller area compared to the present-day observations and 7 showing an overestimated area. LOVECLIM, EC-Earth, FGOALS-g3, GISS-E2-1-G and INM-CM4-8 clearly simulate too much ice (Table 4). The other models generally exhibit realistic PI summer minimum ice conditions. The detail of the geographical distribution of sea ice for the models, the MMM and the NOAA_OI_v2 datasets (Fig. 3) confirms the results in terms of Arctic sea ice extent. Overestimations appear to be due to too much sea ice being simulated in the Barents–Kara area (LOVECLIM, FGOALS-g3, GISS-E2-1-G), in the Nordic Seas (EC-Earth, FGOALS-g3) and in Baffin Bay (LOVECLIM, INM-CM4-8, EX-Earth). MIROC-ES2L performs rather poorly for the PI, with insufficient ice close to the continents. The other models generally match the 0.15 isoline from the NOAA_OI_v2 dataset in a realistic manner. The winter maximum monthly MMM areas show little difference between the present-day and PI simulated areas. The MMM PI area is 15.16 ± 1.90×10⁶ km², compared to the observed 1981 to 2002 mean of 14.44×10⁶ km². For both the summer and winter, the simulations and observations mostly agree on the month that the minimum and maximum are attained: August–September for the minimum and February–March for the maximum for every model.

Figure 3 PI sea ice concentration for the month of minimum SIA as computed for Fig. 2. The magenta contour shows the 0.15 isocontour of the NOAA_OI_v2 observations (Reynolds et al., 2002, see the Data availability section) averaged over the years 1982–2001.

Before we carry out the comparison between model results and sea ice cover reconstructions for the LIG period, we compare the results of the models for PI to the observations at the reconstruction sites (Fig. 4 for the comparison of annual mean sea ice concentrations and Fig. 5a and b for winter and summer). Models generally overestimate sea ice cover at the three northernmost sites in summer and in annual mean and over the seven northernmost sites for the winter season. Those sites are actually very close to the sea ice edge and the overestimation could correspond to the fact that the observations are for 1981 to 2002 period, which was already warmer than the pre-industrial one.

Figure 4 Sea ice annual concentration simulated for PI, for the multi-model mean (MMM) and for each model. The colour filling of the symbols on the maps correspond to the observed values at each site, which are classified into three categories according to the NOAA_OI_v2 dataset: perennial cover (9 to 12 months), seasonal cover (3 to 9 months) and ice-free state (0 to 3 months). On the MMM panel, for each data site the colour of the symbol outline corresponds to the number of models simulating the observed ice cover. On the panels for individual models, the shape of the symbol depends on the observed result being below the simulated one (triangle down), above the simulated one (triangle up) or in the same category as the simulated one (circle). The number of sites for which reconstructions are equal to and above the number of months simulated by models are written at the bottom-right corner of each panel.

Figure 5 Model–data comparison as a function of latitude and record site for PI (a, b) and LIG (c, d). For each LIG data site, the NOAA_OI_v2 observations (PI) or reconstructions (LIG) are shown in the first column of each plot, and the model results are in the columns to the right. Both the model results and the NOAA_OI_v2 observations are shown in terms of sea ice fraction averaged over the month of minimum (a, c) or maximum (b, d) Northern Hemisphere sea ice area and the previous and following months. For the LIG, the qualitative assessments (eighth column of Table 1) have been used for records with good chronological control. The letter next to the name of the site stands for indicator used for the reconstruction: dinocysts (“d”) or IP25/PIP25 (“i”). For summer conditions, dark blue shading is used for the “no sea ice” category, light blue shading is used for “summer probably ice-free” conditions, and white shading is used for “summer sea ice concentration > 75 %” and “perennial sea ice”. For winter, dark blue shading is used for “ice-free all year round” conditions and white shading is used for “seasonal sea ice conditions” and “perennial sea ice”. The model results are averaged as they are for PI and shown following the colour scale on the right-hand side of the plots.

3.2 LIG sea ice

The models show a minimum monthly MMM SIA for the LIG of 3.20 ± 1.50×10⁶ km², and a maximum MMM SIA of 15.95 ± 2.61×10⁶ km². In terms of SIE, the minimum MMM extent is 5.39 ± 2.13×10⁶ km², while the maximum MMM extent is equal to 18.38 ± 3.12×10⁶ km². Thus, compared to the PI results, there is a reduction of ca. 50 % in the MMM minimum (summer) monthly SIA in the LIG results, and of nearly 40 % in terms of SIE, but a slight increase in the winter monthly MMM SIA and SIE. Every model shows an often substantial reduction in summer sea ice between the PI and LIG.

There is a large amount of inter-model variability for the LIG SIA and SIE during the summer (Fig. 6 and Table 4). Out of the 16 models, 1 model, HadGEM3, shows a LIG Arctic Ocean free of sea ice in summer, i.e. with an SIE lower than 1×10⁶ km². CESM2 and NESM3 show low SIA values (slightly above 2×10⁶ km²) in summer for the LIG simulation, but their minimum SIE values are below 3×10⁶ km². Both HadGEM3 and CESM2 realistically capture the PI Arctic sea ice seasonal cycle. On the other hand, NESM3 overestimates winter ice and the amplitude of the seasonal cycle in SIA and SIE, while simulating realistic PI values for both SIA and SIE (Cao et al., 2018). This seasonal cycle is amplified in the LIG simulation, with an increase in SIA and SIE in winter and a decrease in summer, following the insolation forcing. Hence, the difference in the response of these models to LIG forcing in terms of sea ice does not appear to only depend on differences in PI sea ice representation.

Figure 6 LIG sea ice concentration for the month of minimum SIA (computed with calendar adjustment) as computed for Fig. 2. The magenta contour shows the 0.15 isocontour of the corresponding PI simulation.

For the winter, only one model (EC-Earth) simulates a decrease in SIA and SIE of around 2 × 10⁶ km², two other models (ACCESS and INM-CM4-8) simulate a slight decrease in SIA and SIE, all other models simulate an increase in both SIA and SIE. All models therefore show a larger sea ice area amplitude for LIG than for PI, and the range of model results is larger for LIG than for PI. The summer season and the seasons of sea ice growth and decay are therefore key to understanding the behaviour of LIG sea ice and the inter-model differences, as will be confirmed in Sect. 4.

3.3 LIG model–data comparison

We limit our comparison to the sites for which the chronology is good. These cores mostly show ice-free conditions in summer, except for the northernmost site (core PS92/039-2), which is at least 75 % covered by ice in summer (Fig. 5c). Two other sites at high latitude (PS213861 and PS93/006-1, for which sea ice has been reconstructed based on dinocysts and IP25/PIP25), show summer conditions which are “probably ice-free”. Only four models simulate more than 75 % sea ice concentration over the northernmost site, but they also simulate more than 75 % sea ice concentration over the two following sites (in descending order of latitudes), and FGOALS-g3 simulates more than 75 % sea ice concentration for another four sites for which the reconstructions show no sea ice. On the other hand, 10 models simulate no sea ice concentration at all over the reconstruction sites in summer, and therefore probably overestimate the LIG summer sea ice reduction. From these reconstructions, we cannot distinguish the performance of the models simulating a strong reduction of sea ice from the model simulating a nearly total disappearance of summer sea ice in the Arctic. Apart from FGOALS-g3, which simulates extensive sea ice cover for both periods, there does not appear to be a strong relationship between the PI and LIG model results over the data sites: models which simulate sea ice cover over the three northernmost sites at the LIG do not necessarily simulate large sea ice concentrations over these sites for PI (e.g. LOVECLIM, AWIESM1 and AWIESM2).

For the winter season, the reconstructions show the four northernmost sites to be ice covered. The reconstructions for most other sites are qualitatively given as “nearly ice-free all year round” or “ice-free all year round”. Model results are generally in agreement with the reconstructions for the three to four northernmost sites (Fig. 5d). Most models simulate sea ice over some of the sites characterised by “nearly ice-free all year round” conditions, and only one model (IPSLCM6) simulates sea ice cover over a site for which the reconstructions show ice-free conditions. The model–data agreement is therefore quite good for the winter season. In this case, the model results for LIG appear to be related to their results for PI, with models simulating more sea ice for LIG being those simulating more sea ice for PI.

Figure 7 shows a quantitative model–data comparison in terms of annual mean sea ice concentration, which is the variable for which we have the highest number of reconstructions (Table 1). From this, we see that it is more difficult for the models to realistically capture sea ice change over the core sites near Greenland close to the sea ice edge. If we cross-compare the observation–model match for each model for both the PI (Fig. 4) and the LIG (Fig. 7) then FGOALS-g3 and NESM3 have difficulties in accurately capturing sea ice cover at the core site locations in the Nordic Seas, whilst AWIESM1 and NORESM1-F best display sea ice cover close to the sea ice edge near Greenland and in the Nordic seas for both time periods. It is these Nordic Seas sea ice edge differences (over the core sites listed in Table 1) that make the difference between the simulation–data matches for each model.

Figure 7 Sea ice annual concentration during the LIG (computed with calendar adjustment), for the multi-model mean (MMM) and for each model. The colour-filling of the symbols on the maps correspond to the reconstructed values, classified into three categories: perennial cover (9 to 12 months), seasonal cover (3 to 9 months) and ice-free state (0 to 3 months). On the MMM panel, for each data site, the colour of the symbol outline corresponds to the number of models simulating the reconstructed ice cover. On the panels for individual models, the shape of the symbol depends on the model result being below the reconstructed one (triangle down), above the reconstructed one (triangle up) or in the same category as the reconstructed one (circle). The number of data points which are above, equal and below the number of months simulated by models are written in the bottom-right corner of each panel.

4 Discussion of model differences

Whilst we cannot yet definitely establish the most likely Arctic sea ice conditions during the LIG, we can investigate sea ice differences across models when we have sufficient model data. We have first performed this analysis for the three models for which we had sufficient data: CESM2, HadGEM3, and IPSLM6. These models each represent a distinct sea ice response to the LIG forcing, i.e. summer sea ice concentration less than 0.15 everywhere (HadGEM3), significant summer sea ice retreat with concentration less than 0.8 in the Central Arctic (CESM2) and modest summer sea ice retreat with a small area with sea ice concentration close to 1 in the Central Arctic (IPSLCM6).

Sea ice formation and melting can be affected by a large number of factors inherent to the atmosphere and the ocean dynamics, alongside the representation of sea ice itself within the model (i.e. the type of sea ice scheme used). In coupled models it can be extremely difficult to identify the causes of essentially coupled model behaviour. Nevertheless, we discuss the short-wave (SW) surface energy balance, ocean, and atmosphere circulations and comment on cloudiness and albedo changes.

4.1 Atmospheric energy budget differences

The atmospheric energy budget LIG–PI anomaly (Fig. 8) is negative in winter and strongly positive in summer, following the imposed insolation anomaly. These anomalies in total heat budget are dominated by the SW budget contribution from May to August. We split the SW budget into the downward (SWdn) and upward (SWup) contributions. Both fluxes are defined to be positive when they are downward and negative when they are upward. Hence, the total SW budget (in black) is the sum of the SWdn contribution (in red) and the SWup contribution (in blue). In this figure, a positive SWup anomaly means that the SWup is less intense at LIG than at PI, hence contributing to an increase in the net SW flux.

Figure 8 The main components of the atmospheric energy budget at the surface averaged over the Arctic (70–90^∘ N) for HadGEM3, CESM2 and IPSLCM6. The LIG–PI anomalies as shown as a function of the month for the total energy budget (Ftot, black), the SW budget (SW, violet), and for the downward (SWdn, red) and upward SW (SWup, blue) fluxes. The sign convention for all fluxes is that fluxes pointing downward are positive and fluxes pointing upward are negative. Panels (b), (d) and (f) show results for which the LIG calendar has been taken into account (for the LIG simulations), while panels (a), (c) and (e) show the results averaged on the PI calendar both for PI and LIG.

For all the models, the total heat budget anomaly is due to (i) an increased downward short-wave flux in spring resulting from the insolation forcing and (ii) a decreased upward short-wave flux in summer, related to the decrease of the albedo due to the smaller sea ice cover. During summer, this decrease in upward short-wave flux more than compensates the decrease in SWdn, which is maximum in August.

Figure 9 Albedo over the Arctic for PI (a, b), LIG (c, d) and LIG–PI (e, f) for HadGEM3, IPSL-CM6 and CESM2. The albedo has been recomputed from the SWup and SWdn fluxes. Panels (a), (c) and (e) show the results for the whole Arctic, while panels (b), (d) and (f) show the results for areas where the sea ice fraction is larger than 0.9. All LIG values have been calendar adjusted.

The summer anomaly reaches 80 W m⁻² in June for HadGEM3, 60 W m⁻² for IPSLCM6 and 50 W m⁻² for CESM2. The differences between the model results are due to a different phasing of the SWdn and SWup anomalies for HadGEM3, compared to the other two models: for HadGEM3, the two fluxes peak in June, while for CESM2 and IPSLCM6, the SWdn flux peaks in May and the SWup signal peaks in July, and thus the anomaly in these fluxes partly compensate for each other. HadGEM3 shows a larger net SW increase despite a SWdn anomaly which is smaller than for the other two models. On the other hand, HadGEM3's SWup component is stronger and always positive, which is different to the other two models, which show a negative SWup contribution in April–May. These differences are associated with differences in albedo for the three models (Fig. 9). HadGEM3's sea ice and Arctic Ocean albedos are always smaller than those simulated by IPSLCM6 and CESM2 and the difference is larger for LIG than for PI. The albedo simulated by HadGEM3 in May and June is particularly low compared to the two other models, which explains why the SWup component peaks earlier. The albedo LIG–PI anomalies over the whole Arctic show that the sea ice albedo feedback is most effective in HadGEM3.

In terms of cloudiness, IPSLCM6 shows differences in the properties of clouds, in terms of optical depth, between PI and LIG, but this could not be investigated due to a lack of data (thus far) for the other models. Thus we cannot tell if LIG–PI anomalies in SWdn fluxes, i.e. differences between HadGEM3's and CESM2 flux, also have a contribution due to cloud changes.

Figure 10 The main components of the atmospheric energy budget averaged over the Arctic (70–90^∘ N) for HadGEM3, CESM2, CNRM-CM6-1, IPSLCM6, INM-CM4-8, MIROC-ES2L, NorESM1 and NorESM2. The LIG–PI anomalies are shown as a function of the month for the total energy budget (black), the SW budget (violet), and for the downward (red) and upward SW (blue) fluxes. The sign convention for all fluxes is the same as for Fig. 8.

The comparison to other model results (Fig. 10) confirm that the behaviour of HadGEM3 is unusual in terms of energy budget. It is the only model in which the anomalies in SWup and SWdn are exactly in phase and produce a much larger anomaly in total heat budget, while in other models those anomalies are not in phase and partly compensate each other.

4.2 Ocean and atmosphere circulation differences

Changes in Arctic sea ice related to ocean heat transport have been found for the CESM large ensemble (Auclair and Tremblay, 2018). The differences can then be amplified by the sea ice albedo feedbacks. We check this in our models by calculating long-term means of the maximum meridional stream function at 26^∘ N for the PI and LIG simulations. These are 19.5 and 18.7 for CESM2, 15.6 and 15.8 Sv for HadGEM3, and 12.9 and 10.4 for IPSLCM6. Thus, the CESM2 and HadGEM3 models exhibit an Atlantic Meridional Overturning Circulation (AMOC) that is almost unchanged between PI and LIG, while in the IPSLCM6 model the AMOC weakens. This implies that a reduced northward oceanic heat transport could prevent sea ice loss in the Central Arctic in some but not all models (see also Stein et al., 2017).

Some differences in the response of sea ice to LIG forcing therefore appear to be due either to differences in atmospheric response (HadGEM3 vs. IPSL-CM6 and CESM2), similar to mechanisms found for current sea ice decline (e.g. He et al., 2019; Olonscheck et al., 2019) or to changes in ocean heat transport (CESM2 vs. IPSLCM6). But while AMOC changes partially explain the differences found between IPSL (more sea ice in Central Arctic) and CESM2 and HadGEM3 (less sea ice in Central Arctic), they do not explain differences between ice-free and ice-covered conditions in HadGEM3 and CESM2.

Differences in atmospheric circulation changes could also explain difference in sea ice response to LIG forcings. We therefore investigate LIG–PI anomalies in sea level pressure (Fig. 11). Most models simulate a decrease in summer mean sea level pressure largely encompassing the Arctic Ocean and adjacent continents. This decrease is not as strong over the Nordic Seas as it is over the Arctic, and this local heterogeneity over the Nordic Seas is model dependent. However, the anomaly in atmospheric circulation is more zonal over the Nordic Seas and northern North Atlantic in HadGEM3 than in CESM2 or IPSLCM6, and therefore differences in atmospheric circulation are probably not causing more warm air to enter the Arctic for HadGEM3 and are thus not the cause of HadGEM3 being so warm over the Arctic. The mean sea level pressure winter anomaly is characterised by a deepening of the Icelandic low for all models except NESM3.

Figure 11 Anomalies (LIG–PI, hPa, shading) in mean sea level pressure for DJF (top plots) and JJA (bottom plots). Contours indicating PI values are superimposed: values every 5 hPa, 1005 hPa isobar in white, black contours for lower values and grey contours for higher values.

Other factors that remain to be investigated include clouds and ocean heat uptake in the Arctic in the different models, e.g. as a function of stratification.

4.3 Transient CO₂ forced responses: LIG vs. transient 1pctCO2

The LIG has higher insolation than PI at high northern latitudes during spring and summer and less significant changes in winter insolation. This is distinct from the increased GHG, which is the dominant forcing for future climates. However, since sea ice minimum occurs in summer, it is of interest to consider possible relationships between CMIP6 model responses for the LIG and those for the transient 1pctCO2 experiments. A total of 12 models have the LIG, PI and 1pctCO2 simulations available. These include models with large, small and intermediate responses in sea ice for the LIG.

Figure 12 suggests that there is indeed such a relationship between the summer sea ice concentration decreases for LIG and the averages from years 50 to 70 of the transient 1pctCO2 simulations: the models that respond strongly at the LIG also respond strongly for the 1pctCO2 forcing, and the model with the smallest response for the LIG (INMCM4-8) has the smallest response to the 1pctCO2 forcing. The relationship shown in Fig. 12 does not last for later periods in the 1pctCO2 runs, when the winter sea ice is also affected by the increased greenhouse gas forcing. This implies inter-comparisons between the LIG simulation and simulations with moderate CO₂ increase (during the transition to high CO₂ levels) should be investigated.

Figure 12 LIG vs. 1pctCO2 July–August–September sea ice areas (for sea ice concentrations larger than 0.15). The results for the 1pctCO2 simulations have been averaged for years 50 to 70.

5 Conclusions

The Last Interglacial period (LIG) was the last time global temperature was substantially higher than the pre-industrial period at high northern latitudes (Otto-Bliesner et al., 2013; Capron et al., 2017; Otto-Bliesner et al., 2017; Fischer et al., 2018; Otto-Bliesner et al., 2021). To help understand the role of Arctic sea ice in these changes, we present a new synthesis of LIG sea ice information using marine core data collected in the Arctic Ocean, Nordic Seas and northern North Atlantic and compare this to PMIP4-LIG simulations.

Our synthesis shows that south of 79^∘ N in the Atlantic and Nordic seas the LIG was definitely seasonally ice-free. These southern sea ice records provide quantitative estimates of sea surface parameters based on dinoflagellate cysts (dinocysts). North of 79^∘ N the sea-ice-related records are more difficult to obtain and interpret. However, the core at 81.5^∘ N brings evidence of summer being probably seasonally ice-free during the LIG from two indicators: dinocysts and IP25/PIP25. The northernmost core with good chronology is located at 81.9^∘ N and shows evidence of substantial (> 75 %) sea ice concentration all year round. Other cores, with debated chronologies, have not been used for model–data comparisons in the present study.

Model results from 16 models show a multi model mean (MMM) summer SIA LIG of 3.20 ± 1.29×10⁶ km², and a winter monthly MMM area of 15.95 ± 1.21×10⁶ km². This is a reduction in SIA of 50 % for the minimum summer month between the PI and LIG but almost no change for the winter month MMM. Every model shows an often substantial reduction in summer sea ice between the PI and LIG. For the winter, only one of the 16 models shows a (small) winter reduction in sea ice between the PI and LIG. This reinforces that the key seasons for understanding LIG warming are the spring, summer and autumn.

We investigate reasons for inter-model differences in LIG Arctic sea ice simulations: we find that the LIG total heat budget anomaly in the Arctic is due to (i) an increased downward short-wave flux in spring, resulting from the insolation forcing, and (ii) a decreased upward short-wave flux in summer, related to the decrease of the albedo due to the smaller sea ice cover. During summer, this decrease in upward short-wave flux more than compensates the decrease in the SWdn, which is at a maximum in August. Differences between the model results are due to a different phasing of the up and down short-wave anomalies in the different models and are associated with the differences in model albedo.

Analysis of IPSLCM6 results shows differences in the properties of clouds, in terms of optical depth, between PI and LIG. Further work is required to identify if this is also important for other models. Changes in Arctic sea ice may also be related to ocean heat transport. Here, we have shown that ocean circulation changes occur for some (but not all) LIG simulations. Other factors that remain to be investigated include clouds and ocean heat uptake in the Arctic in the different models.

Most models agree with the reconstructed year-round ice-free northern North Atlantic. Model–data disagreement for the LIG occur over the Nordic Seas, close to Greenland and at the boundary with the Arctic Ocean, where many models overestimate annual mean sea ice concentration. This is not fully related to the model performance for summer. Indeed, 12 of 16 models simulate little sea ice cover over the northernmost site and 10 of the models simulate less than 25 % sea ice concentration over the site at 81.5^∘ N. It is not possible, from the available data, to decide on the best models, in particular in terms of summer sea ice. The northernmost site appears to discriminate those models that simulate very little sea ice at this site. However, models which do simulate > 75 % summer sea ice concentration at this site also simulate > 75 % summer sea ice concentration for the two sites at 81.5 and 79.2^∘ N, just south of the northernmost site, which is not realistic. More reconstructions with good chronology are needed in the Central Arctic to determine which model behaviour is more realistic, and in particular if the summer ice-free Arctic simulated by the HadGEM3 model alone, among the 16 models, is possible. This would be key in assessing ESMs used for future projections with respect to climates with much warmer summers than today. This means that it is all the more crucial that there appear to be a nearly linear relationship between the ESM simulations of summer sea ice for the near future (years 50 to 70 of transient 1pctCO2 simulations) and that simulated for the LIG: the models which respond strongly to the LIG forcing also respond strongly for the 1pctCO2 forcing. This implies inter-comparisons between the LIG simulation and simulations with a moderate CO₂ increase (during the transition to high CO₂ levels) may yield insight into likely 21st century Arctic sea ice changes, especially if we achieve a more extensive characterisation of LIG Arctic sea ice from marine cores.