Simulating the mid-Holocene, last interglacial and mid-Pliocene climate with EC-Earth3-LR

As global warming is proceeding due to rising greenhouse gas concentrations, the Earth system moves towards climate states that challenge adaptation. Past Earth system states are offering possible modelling systems for the global warming of the coming decades. These include the climate of the mid-Pliocene (∼ 3 Ma), the last interglacial (∼ 129–116 ka) and the mid-Holocene (∼ 6 ka). The simulations for these past warm periods are the key experiments in the Paleoclimate Model Intercomparison Project (PMIP) phase 4, contributing to phase 6 of the Coupled Model Intercomparison Project (CMIP6). Paleoclimate modelling has long been regarded as a robust out-of-sample test bed of the climate models used to project future climate changes. Here, we document the model setup for PMIP4 experiments with EC-Earth3-LR and present the large-scale features from the simulations for the mid-Holocene, the last interglacial and the mid-Pliocene. Using the pre-industrial climate as a reference state, we show global temperature changes, large-scale Hadley circulation and Walker circulation, polar warming, global monsoons and the climate variability modes – El Niño–Southern Oscillation (ENSO), the Pacific Decadal Oscillation (PDO) and the Atlantic Multidecadal Oscillation (AMO). EC-Earth3-LR simulates reasonable climate responses during past warm periods, as shown in the other PMIP4-CMIP6 model ensemble. The systematic comparison of these climate changes in past three warm periods in an individual model demonstrates the model's ability to capture the climate response under different climate forcings, providing potential implications for confidence in future projections with the EC-Earth model.

1 Introduction

Looking back on Earth's history, the warm climate that we are now experiencing is no exception. Several warm periods offer possible geological analogues for the future: the mid-Pliocene warm period (3.264–3.025 million years ago), the last interglacial (129 000–116 000 years ago) and the mid-Holocene (6000 years ago). The mid-Pliocene is the most recent period with atmospheric CO₂ compared to the present (approximately 400 ppmv) (Pagani et al., 2010), with average annual surface temperatures of roughly 1.8 to 3.6 ^∘C warmer than pre-industrial temperatures, reduced ice sheet size and increased sea levels (Haywood et al., 2013). Global mean annual temperatures during the last interglacial were approximately 0.8 ^∘C (maximum 1.3 ^∘C) warmer than pre-industrial temperatures (Fischer et al., 2018) and the amplified seasonality characterized the northern latitudes (Marcott et al., 2013). During the mid-Holocene period, temperatures were 0.7 ^∘C warmer than pre-industrial temperatures (Otto-Bliesner et al., 2017), with increased seasonal temperatures and strengthened Northern Hemisphere (NH) monsoon (Marcott et al., 2013). The simulations for these warm periods are the key experiments in the Paleoclimate Model Intercomparison Project phase 4, contributing to phase 6 of the Coupled Model Intercomparison Project (PMIP4-CMIP6) (Kageyama et al., 2018). These past warm periods provide essential insights into the climate processes that act in the context of higher CO₂ or stronger insolation. The mid-Pliocene experiment is designed to elucidate the climate system's long-term response to a concentration of atmospheric CO₂, close to present at 400 ppm (Haywood et al., 2016). Simulations in the mid-Holocene and last interglacial periods provide an opportunity to examine the impact of two different radiative forcing changes on the climate when other forcings were relatively similar to those present (Otto-Bliesner et al., 2017). Although the simulations for these past warm periods are not the perfect analogue for future projections, each experiment has distinct differences in the climate system's forcing and initial condition. Available paleorecords provide a means of assessing whether we are correctly capturing the climate response to different forcings in the models being used for future climate change projections. These simulations can be used to evaluate the response of ocean circulation, Arctic sea ice, climate variability patterns (e.g., El Niño–Southern Oscillation), the global hydrological cycle and regional monsoon systems to elevated concentrations of atmospheric CO₂ or the redistribution of solar energy around the globe. The evaluation will be a direct out-of-sample test of state-of-the-art models' reliability to simulate climate change, particularly climate warming.

EC-Earth3-LR is one of the models that contribute to PMIP4-CMIP6. For an individual climate model, the PMIP4 experiments' contribution provides an opportunity to test the climate response to different forcing and an overview of the past climate change in one model system. All the participating modelling groups follow the PMIP4 protocols described in Kageyama et al. (2018). The detailed justifications of the experimental protocols for each period are described separately in specific PMIP4 experiment design papers: Otto-Bliesner et al. (2017) for the mid-Holocene (midHolocene) and last interglacial at 127 ka (lig127k) experiments and Haywood et al. (2016) for the mid-Pliocene experiment (midPliocene). These past periods are characterized by changes in boundary conditions such as greenhouse gas concentrations (e.g., midPliocene) and orbital parameters (e.g., midHolocene and lig127k). The large-scale features from the PMIP4 multi-model ensemble have been evaluated for respective PMIP4 experiments for midHolocene (Brierley et al., 2020), lig127k (Otto-Bliesner et al., 2021) and midPliocene (Haywood et al., 2020a).

In the present paper, we document the experimental setup and the model performance for the pre-industrial, the mid-Holocene, the last interglacial and the mid-Pliocene with EC-Earth3-LR. We expect that a systematic comparison across different past warm climates in one individual model can improve our understanding of global warming through different forcings and provide a more reliable scientific base to evaluate the future climate projections.

In the following presentation, we use piControl, midHolocene, lig127k and midPliocene to refer to the official PMIP-CMIP6 experiment names (except for midPliocene-eoi400 for which we use a short version: midPliocene). We use the abbreviations PI, MH, LIG and MPlio to refer to these four experiments in the illustrations.

2 EC-Earth model description

2.1 EC-Earth model components and the configuration

EC-Earth is a fully coupled Earth system model that integrates several state-of-the-art components in the climate system including atmosphere, ocean, sea ice, land and biosphere. It is developed by the European consortium of 27 research institutions to date and widely used in various studies on climate change (Hazeleger et al., 2010; Hazeleger et al., 2012). The atmospheric component of EC-Earth3 is the Integrated Forecasting System (IFS) model of the European Centre for Medium-Range Weather Forecasts (ECMWF), with cycle 36r4. The revised Tiled ECMWF Scheme for Surface Exchanges over Land incorporating land surface hydrology (HTESSEL) model is used for land surface use (Balsamo et al., 2009) as an integrated part of IFS. The ocean component is the Nucleus for the European Modelling of the Ocean (NEMO) (Madec, 2008). It uses a tripolar grid of poles over northern North America, Siberia and Antarctica. NEMO includes the Louvain-la-Neuve Sea-ice model version 3 (LIM3) (Vancoppenolle et al., 2012), a dynamic thermodynamic sea-ice model with five ice thickness categories.

Atmosphere–land and ocean–sea-ice components are coupled through OASIS (Ocean, Atmosphere, Sea Ice, Soil) (Valcke and Morel, 2006; Craig et al., 2017). The remapping of runoff from the atmospheric grid points to the runoff areas on the ocean grid has been re-implemented to be independent of grid resolution. This was achieved by introducing an auxiliary model component and relying on the interpolation routines provided by the OASIS coupling. In the EC-Earth3-LR configuration, the time step of IFS at T159L62 resolution is 60 min except for radiation, which is solved only every 3 h (but updated with cloud cover information in between two full computations). The time step for NEMO at ORCA1L75 resolution is 45 min; we use the same time step for the ocean and for sea ice. The coupling between IFS and NEMO is 3 h.

The EC-Earth3 Earth system model version also includes the Tracer Model version 5 (TM5) atmospheric chemistry component, the Lund–Potsdam–Jena General Ecosystem Simulator (LPJ-GUESS) dynamic vegetation model and the Pelagic Interactions Scheme for Carbon and Ecosystem Studies (PISCES) ocean biogeochemistry component. The EC-Earth3 model contributes to CMIP6 in several configurations by using different resolutions or including different components. A detailed description of EC-Earth3 and its contribution to CMIP6 is documented in Döscher et al. (2021). For example, regarding the resolution, the standard version EC-Earth3 with atmospheric resolution T255L91 (horizontal ∼ 80 km) and ocean resolution ORCA1L75 (horizontal 1^∘) is used in most MIPs (e.g., CMIP, DCPP, LS3MIP, PAMIP, RFMIP, ScenarioMIP, VolMIP, CORDEX, DynVarMIP, SIMIP). The high-resolution model version with atmospheric resolution T511L91 (horizontal ∼ 25 km) and ocean resolution ORCA025L75 (horizontal 0.25^∘) is used for DCPP and HighResMIP decadal predictions.

The PMIP4 experiments require long spin-up and equilibrium simulations which are computationally demanding. We thus use model configurations without atmospheric chemistry and ocean biogeochemistry at low resolution, with atmospheric resolution T159L62 (horizontal ∼ 125 km) and ocean resolution ORCA1L75. For piControl, midHolocene, lig127k and midPliocene, we have used EC-Earth3-LR, a configuration of EC-Earth3 with only atmosphere, ocean and sea-ice components at low resolution (LR). The coupled vegetation in LPJ-GUESS is used for the past1000 simulation with the model configuration named EC-Earth3-veg-LR. Despite using the lowest resolution among the EC-Earth3 CMIP6 MIPs configuration, EC-Earth3-LR has a relatively higher resolution than the other PMIP4 models (Brierley et al., 2020).

2.2 Implemented physics that are needed for paleosimulations

Two critical physical processes, i.e., orbital forcing and physics over the ice sheet, crucial for paleoclimate simulations, are implemented in EC-Earth3. Below, we describe the implementation of these two physics in EC-Earth3.

2.2.1 Orbital forcing

Previous versions of EC-Earth treated the orbital forcing as constants for the present-day climate in IFS by following the International Astronomical Union's recommendations. These formulas are not valid for dates too far away from 1 January 2000 (if it is more than one century). For the paleoclimate simulations, the orbital forcing is calculated in the IFS atmosphere component (Berger, 1978). According to the model time, the new orbital parameterization uses the Julian calendar to calculate the cosine of the solar zenith angle. The Julian calendar has two types of years: a standard year of 365 d and a leap year of 366 d. To follow the Gregorian calendar used in IFS, which contains leap years, the orbital parameter calculation follows a simple cycle of three normal years and one leap year, giving an average year length of 365.25 d. The annual and diurnal cycle of solar insolation is thus represented with a repeatable solar year of precisely 365.25 d and with a mean solar day of precisely 24 h, respectively. The year must fall within the range of 1950 ± 10⁶. In the EC-Earth model run script, there are two ways to set the orbital forcing: either change the three orbital parameters to known values corresponding to the specific year (defaults are for 1950) or specify a number for the year. To determine the year, the exact year number applies for the AD year; e.g., for the year 1850 AD, we set “ifs_orb_iyear = 1850”. For the year before AD, the distance to the reference year (1950) is given. For example, the mid-Holocene, which is 6000 years before the present day, can be specified as −4950 (obtained from −6000 + 1950). For the equilibrium simulations, the orbital forcing is fixed for each year, and we set “ifs_orb_mode = fixed”.

2.2.2 Albedo parameterization of snow on ice sheets

The ice-sheet-albedo calculation is highly parameterized in the earlier version of EC-Earth by taking a constant value for areas with thick perennial snow cover. This is an important reason for why the Greenland ice sheet surface mass balance is poorly resolved in the model. The surface scheme assumes the ice sheet as 10 m of perennial snow, and it is in thermal contact with the underlying soil. The albedo and snow density are fixed at 0.8 and 300 kg m³, respectively. The fixed, high-snow-albedo value over the ice sheet makes it difficult for the snow melting due to lack of albedo feedback. This process becomes crucial, especially for past periods with large ice-sheet extents, such as Last Glacial Maximum. Therefore, we introduce a new albedo parameterization of snow on the ice sheet, where snow melting becomes possible and occurs often. This varying snow albedo scheme is also active for the configuration of EC-Earth3 coupled with an interactive ice sheet model for the Greenland ice sheet, i.e., EC-Earth3-GrIS (Döscher et al., 2021). A comparison among different albedo parameterizations on the Greenland ice sheet is summarized in Helsen et al. (2017) using EC-Earth3. The time-varying snow albedo better represents the ice and/or snow albedo feedback and is beneficial for both cold and warm climates. From EC-Earth3, the snow is not allowed to accumulate for more than 9 m of water equivalent. All excessive snow falls on the ice sheet is immediately redistributed into the ocean and will not accumulate on the ice sheet. This treatment has an advantage of energy and mass conservation, and the calving ice cools the ocean along the coast, which reduces the warm bias in the Southern Ocean.

The two abovementioned physics are included in the CMIP6 version of EC-Earth3, with on and off options. For the PMIP4 experiments, we have applied these new physical implementations. However, other CMIP6 experiments, e.g., future projections, have not used these new implementations.

3 The experiment design and setup with EC-Earth3-LR

3.1 PMIP4-CMIP6 protocols on climate forcings and boundary conditions

We follow the PMIP4 protocol and implement the climate forcings and boundary conditions for each PMIP4 experiment as documented in Kageyama et al. (2018). The climate forcings implemented in the four experiments are listed in Table 1.

Table 1 PMIP4-CMIP6 protocols on forcings and boundary conditions.

For all the four experiments, the solar constant is 1360.747 W m⁻², fixed for piControl at the mean value for the first two solar cycles of the historical period (1850–1871) (Eyring et al., 2016). The vernal equinox is fixed at noon on 21 March. Aerosols, including dust and volcanic aerosols, are set as the piControl value. Datasets for midPliocene from the PRISM4 reconstructions are downloaded from PlioMIP2 website: https://geology.er.usgs.gov/egpsc/prism/7.2_pliomip2_data.html (last access: 22 February 2021).

In agreement with the PMIP4 protocol, vegetation is prescribed as present-day climatology that is constant in time. The present-day climatological vegetation, based on ECMWF ERA-Interim (Dee et al., 2011) and specified as albedos and leaf area index (LAI) of the Moderate Resolution Imaging Spectroradiometer (MODIS), is used in the piControl, midHolocene and lig127k simulations.

3.2 Boundary condition implementation for midHolocene and lig127k

In Tier 1 of PMIP4 for the midHolocene and the lig127k experiments, the changes in boundary conditions are orbital parameters and greenhouse gas concentrations (Table 1). As described earlier, the setting for orbital forcing and greenhouse gas in EC-Earth3-LR is straightforward. In the model run script, we set “Orb_year = −4950” for midHolocene and “Orb_year = –125 050” for lig127k; the computed orbital parameters are the same as those listed in Table 1.

The changes in orbital parameters, particularly the higher obliquity in midHolocene and lig127k, lead to large changes in the insolation's seasonal and latitudinal distribution, with stronger changes in lig127k than in midHolocene (Fig. 1). The global change in annual mean insolation is negligible. The Northern Hemisphere (NH) receives 20 W m⁻² more insolation during boreal summer in June–August in midHolocene and 48 W m⁻² in May–July in lig127k, while in boreal wintertime it receives 14 W m⁻² less in February–March in midHolocene and 36 W m⁻² less during October–November in lig127k. The maximum changes occur in the high latitudes; i.e., the Arctic region (60–90^∘ N) receives 31 W m⁻² more in July in midHolocene and 70 W m⁻² more in June in lig127k. The Southern Hemisphere (SH) receives 23 and 56 W m⁻² more in September and October in midHolocene and lig127k. These large insolation changes in high latitudes cause polar amplification (Sect. 4.4). The positive insolation anomalies in summer months lead to strong land–sea contrast and enhance the monsoons (Sect. 4.5).

The greenhouse gas concentration is provided by CMIP6; for piControl, the values are taken from the values at the year 1850 in the historical forcing files. An easy way for other PMIP4 periods to do this is to create a new file by changing those values in the year 1850 to the values specified in Table 1.

Figure 1 Changes in the latitudinal and seasonal distribution of insolation (W m⁻²) as compared to piControl: (a) midHolocene and (b) lig127k. The modern calendar is used, with the vernal equinox on 21 March at noon.

3.3 Boundary condition implementation for midPliocene

In the PlioMIP2 protocol, the mid-Pliocene period is centred at an interglacial peak (Marine Isotope Stage KM5c) at around 3.025 Ma, which has similar orbital forcing as the present day. Since the differences between using the 3.025 Ma orbit and the present-day orbit are minimal (de Boer et al., 2017), the orbital parameters are kept the same as in piControl. For the PlioMIP2 Tier 1 experiment, the CO₂ concentration in the atmosphere is set at 400 ppm. CH₄ and N₂O are specified as identical to those in piControl experiment. We have used the enhanced boundary condition according to PRISM4 reconstruction (Dowsett et al., 2016), which considers the change in dynamic topography associated with mantle flow and glacial isostatic adjustment due to Piacenzian ice loading. The topography changes include the closure of the Bering Strait and the straits through the Canadian Arctic Archipelago. We implement the topography, land–sea mask, bathymetry and ice sheet from PRISM4 to the horizontal grid in the EC-Earth3-LR atmosphere and ocean (Fig. 2). The global distributions of soil and biome are modified to match the midPliocene land–sea mask and ice-sheet reconstruction. Due to the change of land–sea mask, there is an emergence of new land areas, and rivers are routed to the nearest ocean point.

Figure 2 Prescribed boundary conditions for mid-Pliocene adapted to EC-Earth resolution: (a) topography and bathymetry (in metres), (b) land–sea mask and ice sheet, (c) high vegetation and (d) low vegetation.

Following the PlioMIP2 protocol, we prescribe the vegetation from PRISM3 reconstruction data (Salzmann et al., 2008). The reconstructed vegetation is a nine-type mega-biome map classified after Harrison and Prentice (2003). In the EC-Earth3 HTESSEL land model, the vegetation type classification is based on the Biosphere-Atmosphere Transfer Scheme (BATS) model (Yang and Dickinson, 1996), which contains 20 biome classifications. Each grid point is characterized by a maximum of two dominant vegetation types, high vegetation and low vegetation, and their area fractions. The sum of the high and low vegetation coverage is always between 0 and 1. Therefore, a translation from mega-biome to HTESSEL classification is performed. We first regridded the mega-biome reconstruction data to the EC-Earth3 TL159 horizontal grid. The translation is then done grid by grid with the correspondence between the two biomes (Table 2). We show the final vegetation types as four low vegetation types and four high vegetation types in HTESSEL (Fig. 2). Due to the present translation scheme, each grid will have only one vegetation type, representing either high or low vegetation. See Table 2 for details.

Table 2 The translation/correspondence between reconstructed mega-biome to EC-Earth HTESSEL vegetation type. H and L refer to the distinction between high (H) and low (L) vegetation.

3.4 Initial conditions, spin-up and production

A piControl spin-up with EC-Earth3-LR configuration is carried out during the model tuning process and contains minor changes in the simulation process. We apply the initial conditions for piControl, midHolocene and lig127k from the end of the piControl spin-up phase (which ran for approximately 1000 years). Our piControl differs from the version of spin-up piControl by applying the new implementation of orbital forcing and land–ice physics. Therefore, it takes around another 200 years of spin-up phase for the piControl simulation to reach equilibrium. The criteria for an equilibrium state are measured by the global mean surface temperature trend ($<\pm \mathrm{0}$.05 ^∘C per century) and a stable Atlantic meridional overturning circulation (AMOC) (Kageyama et al., 2018). Due to the changed boundary conditions, the midHolocene and lig127k simulations reach the equilibrium state after about 200 years of spin-up. The midPliocene experiment starts from the initial condition of a Levitus ocean state, and we perform a spin-up for 1400 years until it reaches the equilibrium state. The change in the land–sea mask in midPliocene requires appropriate initial fields for the changed grid box. These fields are all modified by the model and are expected to reach equilibrium with the given boundary conditions.

Figure 3 The evolution during the 200-year spin-up period for the piControl simulation using annual mean data. (a) Global mean energy balance at the top of the atmosphere (W m⁻²), (b) global mean energy balance at the surface (W m⁻²), (c) NH and SH average surface air temperature (^∘C), (d) NH and SH average sea-ice extent (10⁶ km²), (e) global mean deep-ocean temperature averaged below 2500 m (^∘C) and (f) global mean deep-ocean salinity averaged below 2500 m (psu). The 200-year mean value is indicated with a dashed red line.

Figure 4 Evolution of annual mean globally averaged SST (^∘C, left panel) and annual mean maximum AMOC (Sv, right panel) in spin-up period for the (a, b) piControl, (c, d) midHolocene, (e, f) lig127k and (g, h) midPliocene simulations. The 200-year mean is indicated with a dashed red line.

The energy balance at the top of the atmosphere (TOA) and the surface are stable in all the simulations (Fig. 3a–b and Figs. S1–S3 in the Supplement). We expect a closure of energy in the atmosphere–ocean system for an equilibrium state, and all our simulations show approximately −0.5 W m⁻² imbalance (TOA energy minus surface). Reducing this energy imbalance is a primary target in the model tuning process under the pre-industrial climate condition, and the tuned parameters are not allowed to change in any other past or future scenarios according to the CMIP6 protocol. The energy imbalance differs for different resolution tunings. For example, in the EC-Earth3 standard resolution (T255L91), the energy imbalance is on the order of 0.25 W m⁻² (Döscher et al., 2021), and in the high-resolution T511L91 tuning version, the imbalance is 0.9 W m⁻² (Shiyu Wang, personal communication, 2020). The 0.5 W m⁻² energy imbalance in EC-Earth3-LR is in the pre-industrial control simulation range in 25 CMIP5 climate models (Hobbs et al., 2016). One reason may be the short spin-up, as the deep-ocean temperature and salinity have shown a larger trend (Figs. 3 and S1–S3) than that of the surface temperature (Fig. 3). The deep-ocean temperature has exhibited an increasing trend in all four simulations. The deep-ocean salinity shows a decreasing trend in piControl (Fig. 3) and midHolocene (Fig. S1) but a rising trend in lig127K (Fig. S2) and midPliocene (Fig. S3). Considering that the climates are warmer in the latter two conditions, both sea-ice melting in the polar regions and enhanced oceanic evaporation can influence the surface salinity in different ways. More diagnostics in ocean dynamics are needed to understand these in-depth ocean features in different climate scenarios. No apparent trend is seen in the sea-ice extent in the Arctic and Antarctic (Figs. 3 and S1–S3). The Arctic sea ice shows a similar low-frequency variability, as seen in AMOC (right panel in Fig. 3), indicating the association between the Arctic sea ice and AMOC. The dramatic overall sea-ice melting is seen in midPliocene (Fig. S3) in both the Arctic and Antarctic. In the other three experiments, the Antarctic sea-ice extent remains similar to that of piControl. These mean climate features in sea ice in the spin-up period remain in the production simulation, which is performed from the end of each spin-up. We present the large-scale features from 200 years of production simulations in Sect. 4.

We present the time evolution of the global mean sea surface temperature (SST) and the maximum AMOC in Fig. 4 for the four simulations. To be comparable, we show the 200-year spin-up period for all the four experiments. The trends in global mean sea surface temperature in the four spin-ups have met the criteria of equilibrium. In the EC-Earth3-LR simulations, the maximum AMOC mostly appears at 25^∘ N latitude and occasionally shifts northward to 45^∘ N. The AMOC is about 17.5 Sv in the piControl simulation and becomes stronger in the midHolocene, lig127k and midPliocene simulations. The multidecadal variability is apparent in both AMOC and SST in the piControl, midHolocene and lig127k simulations. The midPliocene simulation indicates the warmest period and shows less low-frequency variability.

3.5 Post-processing of the model output for PMIP4-CMIP6

Following the CMOR tables (Climate Model Output Rewriter, output format under all CMIP standards), we use the ece2cmor3 software package for post-processing model output. The ece2cmor3 package is a Python package developed by the EC-Earth community to convert the EC-Earth3 model output to a CMIP6-compliant format. It uses Climate Data Operators (CDO) (Schulzweida, 2019) and the CMOR library bindings to select variables and vertical levels, perform time average (or take the daily extreme), interpolate spectral and grid-point atmospheric fields to a regular Gaussian grid and calculate variables by an arithmetic combination of the original model fields. The ece2cmor3 package uses the Program for Climate Model Diagnosis and Intercomparison (PCMDI) CMOR library to produce NetCDF files with the appropriate format and metadata. All the data presented in this paper are publicly available on the CMIP6 database at https://esgf-node.llnl.gov/projects/cmip6/ (last access: 22 February 2021).

It is worth noting that particular caution should be paid when downloading the piControl data from ESGF data nodes. There are several piControl simulations by the EC-Earth3 model available in the CMIP6 database. It is appropriate to select the one with the model name “EC-Earth3-LR” as the other PMIP4 experiments' reference. The other piControl experiments by EC-Earth3 differ in the model configuration or model resolution.

4 Climate responses in the midHolocene, lig127k and midPliocene simulations

Here, we present the changes in the large-scale features of the past three warm climates simulated by EC-Earth3-LR, i.e., the two interglacial simulations midHolocene and lig127k, as well as the high-CO₂ featured midPliocene, with the comparison to the piControl reference simulation.

4.1 Large-scale climate response

We summarize the global annual mean values for climate anomalies from the three experiments to piControl in Table 3. The global mean temperature from the midHolocene simulation shows a slightly cooling, placing the simulation towards the upper end of the PMIP4 midHolocene ensemble (Brierley et al., 2020). The global mean precipitation anomaly is mild in the midHolocene simulation, with about 5 %–7 % more precipitation than in piControl over the land during the boreal summer and autumn, following the enhanced Northern Hemisphere monsoon (Sect. 4.5). The more considerable changes in orbital forcing in the lig127k result in 0.4 ^∘C warmings globally, and this warming is amplified over land and in the Arctic. In lig127k, the land precipitation increased 16 % in boreal summer, indicating a more vigorous monsoon than in the midHolocene simulation. The midPliocene simulation shows significant global warming reaching 4.8 ^∘C, which is among the warmest simulations with the highest Arctic amplification in the PlioMIP2 ensemble (Haywood et al., 2020b; de Nooijer et al., 2020). The precipitation increases by more than 10 % globally, and the land precipitation increased more than 20 % in the boreal summer and autumn seasons in the midPliocene simulation.

Table 3 Global mean temperature and precipitation anomalies with respect to piControl for annual mean (ANN), December–January–February (DJF), March–April–May (MAM), June–July–August (JJA) and September–October–November (SON) means.

Figure 5 Mean annual temperature (^∘C) in (a) difference between piControl and ERA20C (1901–1930 mean), and the anomaly in (b) midHolocene, (c) lig127k and (d) midPliocene from piControl.

For the climate changes in a global pattern, we show a comparison between the annual mean temperature and precipitation in piControl with the early century (1901–1930) mean in ERA20C reanalysis data (Figs. 5a and 6a). EC-Earth3-LR tends to overestimate the temperature in high latitudes, with cold biases in the NH and warm biases in the SH. This is a reasonably robust feature of most global climate models, particularly prominent in the EC-Earth model, and the cause is attributed to the simulated sea ice (Koenigk et al., 2013; Chapman and Walsh, 2007). The dry biases in NH and wet biases in SH are noticed when comparing the simulated PI precipitation with ERA20C, primarily in the equatorial ocean area (Fig. 6a). In the following analysis, we compare the anomalies in the three paleoclimate experiments to piControl; any changes in the temperature reflect the responses to the external forcing.

Figure 6 Mean annual precipitation (mm d⁻¹) in (a) difference between piControl and ERA20C (1901–1930 mean), and the anomaly in (b) midHolocene, (c) lig127k and (d) midPliocene from piControl.

During the two interglacial periods, the orbital forcing is characterized by the redistribution of the insolation over the globe in different seasons (Fig. 1). This does not affect too much of the global energy balance; therefore, global mean climate changes are insignificant (Table 3). The spatial patterns of mean annual temperature show warming in the Arctic region in the two interglacial periods relative to the pre-industrial, with a larger magnitude in lig127k than in midHolocene (Fig. 5b and c). The warming in the southern high latitudes shown in the PMIP4 midHolocene ensemble (Brierley et al., 2020) is not apparent in the EC-Earth3-LR simulation. The simulated midPliocene warming is much stronger than those of the two interglacial periods, dominated by the robust Arctic amplification in EC-Earth3-LR (de Nooijer et al., 2020). In midPliocene, the strong warming amplification appears in both polar regions (Fig. 5d), which contrasts the weak warming in some parts in the southern high latitudes in midHolocene and lig127K (Fig. 5b and c). The seasonal features of the temperature response in the polar regions are discussed more in detail in Sect. 4.4. In all three warm periods, the response in annual mean temperature is stronger over land than over the ocean (Fig. 5 and Table 3), which is a common feature in most PMIP4-CMIP6 simulations (Brierley et al., 2020; Otto-Bliesner et al., 2021; Haywood et al., 2020b). These notable changes in polar amplification and land–sea contrast modulate the temperature gradients and lead to changes in the large-scale circulation (Sect. 4.4) and global monsoons (Sect. 4.5). The most pronounced and robust changes in the global land precipitation occur in the Sahel and the subcontinent of India, indicating an enhanced western African monsoon and Indian monsoon in the midHolocene, lig127k and midPliocene simulations (Fig. 6). The increased East Asian monsoon appears in midPliocene but is not evident in midHolocene and lig127k. More analyses on global monsoon are presented in Sect. 4.5. Over the Indian Ocean, the precipitation increases over the western Indian Ocean. It decreases in the eastern Indian Ocean during midHolocene (Fig. 6b) and lig127k (Fig. 6c), opposite to that in midPliocene (Fig. 6d). There is a northward shift of the Intertropical Convergence Zone (ITCZ) in the Atlantic and Pacific rain belts in the three simulations, with magnitude increases in midHolocene, lig127k and midPliocene. Increased precipitation is also observed in northern high latitudes in the three warm periods associated with the strong Arctic warming amplification.

4.2 Data–model comparison

The model simulations for lig127k and midPliocene are assessed using two proxy reconstructions of sea surface temperatures, following the methods presented in Williams et al. (2020) and Haywood et al. (2020a). Due to the low spatial coverage of available SST reconstructions for midHolocene, this period was excluded from the data–model comparison. The lig127k simulation is assessed against the global marine compilation by Hoffman et al. (2017), consisting of 86 annually reconstructed multi-proxy records from foraminifera, radiolaria, coccolithophores, alkenones, $\mathrm{Mg}/\mathrm{Ca}$ and diatoms, and the values are taken from the provided 127 ka time slice. The midPliocene simulation is assessed against the SST dataset provided by Foley and Dowsett (2019) covering a 30 000-year interval around the Marine Isotope Stage (MIS) KM5c peak. The proxy records are reconstructed using the alkenone paleotemperature (Uk37) method. The anomalies for both lig127k and midPliocene are calculated relative to pre-industrial (1870–1899) Hadley Centre Sea Ice and Sea Surface Temperature (HadISST) mean SSTs (Rayner et al., 2003) for the same data points.

Figure 7 Mean annual sea surface temperature (^∘C) anomalies between (a) lig127k and (b) midPliocene, and piControl. The shadings are overlaid by the proxy-data anomalies from (a) Hoffman et al. (2017) and (b) Foley and Dowsett (2019), calculated relative to HadISST (1870–1899 mean). Proxy-data anomalies for the data points are described in Sect. 4.2 for lig127k and midPliocene, respectively.

In Fig. 7, we show the mean annual sea surface temperatures anomalies from the two model simulations, calculated with respect to piControl. The lig127k simulation (Fig. 7a) shows a warming of the Northern Hemisphere with the highest values in the Arctic region. In contrast, the Southern Hemisphere shows a widespread cooling, overall consistent with the PMIP4 ensemble (Otto-Bliesner et al., 2021). There is generally good agreement globally between the lig127k simulation and the reconstructed SSTs. The proxies indicate high-latitude warming in both hemispheres and cooling along the Equator, which is consistent with the overall pattern of the model simulation. However, the magnitude of warming shown by the proxies in the Southern Hemisphere suggests larger warming than is captured in EC-Earth3-LR. Over the North Atlantic, proxy reconstructions indicate warming right next to cooling, making a comparison to the simulation difficult, but this suggests that EC-Earth3-LR overestimates the homogeneity of the North Atlantic warming in lig127k and fails to capture both the cooling and its magnitude in the Norwegian Sea. However, this feature is consistent with the results of the PMIP4 ensemble, which also fails to capture this cooling (Otto-Bliesner et al., 2021). The magnitude of warming of the coasts in the Southern Hemisphere (e.g., west of South Africa and east of New Zealand) shows large discrepancies between the model simulation and proxy reconstructions in signs of change in magnitude especially. Simulated lig127k SST suggests temperature anomalies of up to 0.5 ^∘C, while the proxy dataset indicates changes exceeding 5 ^∘C. This underestimation of the warming of coastal upwelling regions is seen in the PMIP4 ensemble. It is suggested to result from a warm bias in the piControl simulations, leaving little room for warming and the low resolution in the PMIP4 ensemble models failing to capture these narrow regions (Otto-Bliesner et al., 2021).

The midPliocene simulation (Fig. 7b) shows general warming globally. A strong polar amplification appears in both hemispheres, consistent with the PlioMIP2 ensemble (Haywood et al., 2020a; de Nooijer et al., 2020) and the proxy dataset (Foley and Dowsett, 2019). The strong North Atlantic warming seen in proxies is captured well by the model simulation, both spatially and in its magnitude. However, EC-Earth3-LR overestimates the warming around the equatorial belt and to some degree also the sign of change, especially in the Central American region where proxy records suggest some degree of cooling during midPliocene. Most considerable discrepancies are found off the coast of California, with up to 8 ^∘C difference between model and proxy data, and in the Sea of Japan, where simulated SST shows a cooling, while the proxy data suggest warming.

4.3 Response in Hadley circulation and Walker circulation

The changes in global temperature pattern, shown above, could induce large atmospheric circulation anomalies, i.e., in the Hadley circulation and Walker circulation. Previous studies have demonstrated that the meridional surface air temperature gradient between low latitudes and mid-to-high latitudes determines the changes in Hadley circulation (Corvec and Fletcher, 2017; Burls and Fedorov, 2017). Indeed, the weakening of the Hadley circulation seen in the NH in the three experiments (Fig. 8a, c and e) corresponds to the decreased temperature gradient between tropics and mid-to-high latitudes (Fig. 5b–d). However, the thermal contrast is less pronounced in the Southern Hemisphere in the midHolocene and lig127k simulations. The two simulations even show slight cooling in the SH compared to the tropics (Fig. 5b and c), resulting in a slightly stronger meridional temperature gradient, and induce a stronger Hadley circulation (Fig. 8c). The similar results can be seen in Fig. 9a; all three simulations show a weaker Hadley circulation in both hemispheres, except lig127k and midHolocene in the Southern Hemisphere. All three warm period simulations show the Hadley circulation is shifted roughly polewards around 30^∘ N–20^∘ S, with respect to the PI simulation. In the lig127k and midPliocene simulations, we observe a poleward widening of the Hadley circulation in both hemispheres. This widening is less evident in the midHolocene simulation.

Figure 8 The changes in annual mean Hadley circulation (represented by zonal-mean meridional mass stream function in the left panels) and Walker circulation (represented by zonal mass stream function averaged between 10^∘ S–10^∘ N over Indo-Pacific region 60^∘ E–150^∘ W in the right panels). The shadings indicate the changes in the (a, b) midHolocene, (c, d) lig127k and (e, f) midPliocene simulations with respect to piControl. The contours represent the climatology from the piControl simulation (solid contours indicate positive values and dashed contours indicate negative values). Units: 10¹⁰ kg s⁻¹.

Figure 9 The meridional (a) and zonal (b) mean of vertical integrated meridional and zonal mass stream function in Fig. 8. Units: 10¹⁰ kg s⁻¹.

The Walker circulation in piControl is characterized as ascending in the Maritime Continent and western Pacific and descending in the eastern Pacific (contours in Fig. 8b, d and f), which is consistent with previous studies (Peixoto and Oort, 1992; Kamae et al., 2011; Bayr et al., 2014). Compared to piControl, all three warm period simulations show strong ascension over the tropical western Pacific, indicating a strengthening and westward shift of the Pacific Walker circulation cell (Fig. 8b, d and f). These features are further shown in Fig. 9b, and it is clear that all three simulations indicate a stronger and westward shift of Pacific Walker circulation. There is a more westward shift in the lig127k and midPliocene simulations than the midHolocene simulation, corresponding to less ascension over the Maritime Continent. The increasing magnitude of shift from piControl to midHolocene, lig127k and midPliocene can be associated with reduced El Niño–Southern Oscillation (ENSO) variability (Sect. 4.6).

4.4 Climate responses in the polar region

The poles' response to the different boundary conditions is of particular interest as the largest temperature anomalies are observed here (Fig. 5), which are known to associate with strong feedback at the poles such as ice-albedo and carbon cycle feedbacks (Masson-Delmotte et al., 2013). The amplification of global temperature anomalies in the Arctic is a near-universal feature of historical simulations in response to the increased greenhouse gas. It is also observed in the instrumental record (Serreze and Barry, 2011). Similar amplification of global temperature anomalies is not observed in the Antarctic, which may be attributed to the rapid removal of surface heat in the Southern Ocean, which, in turn, limits the ability of climate feedbacks to amplify the warming (Stroeve et al., 2007). Arctic amplification is prominent in the EC-Earth3-LR midPliocene simulation, with Arctic (60–90^∘ N) temperature anomalies exceeding global temperature anomalies by a factor of 2.4 (Table 3). The Arctic warming in the midHolocene and lig127k simulations is not as striking as in the midPliocene simulation. Given that the global mean temperature anomalies in the midHolocene and lig127k simulations are small (Table 3), it is inappropriate to apply an amplification ratio as discussed in Hind et al. (2016).

The minimum Arctic surface air temperature (SAT) warming is expected in the summer because of increased ocean heat uptake following reductions in sea ice, while maximum SAT warming is expected in the (boreal) autumn and winter following the release of this heat (Pithan and Mauritsen, 2014; Serreze et al., 2009; Yoshimori and Suzuki, 2019; Zheng et al., 2019). This expected seasonality of SAT is not present in the Arctic in the midHolocene and lig127k simulations (Fig. 10). There are peaks in Arctic SAT warming, amounting to +2.3 and +5.8 ^∘C in the midHolocene and lig127k simulations, respectively, in the boreal summer (Fig. 10a) due to positive insolation changes (Fig. 1). Negative insolation changes during the autumn may offset some of the expected autumn excess heat release. In the midPliocene simulation, Arctic SAT warming peaks in the autumn at +15.2 ^∘C, but there is no exact summer minimum (Fig. 10a). In the Antarctic, the SAT seasonality in the two interglacial simulations does not differ significantly from the piControl simulation. In the midPliocene simulation, Antarctic warming peaks in late austral winter and early spring (July–September), with maximum warming reaching +10.4 ^∘C in August.

Arctic SAT warming will lead to a reduction in the sea-ice extent. The sea-ice extent defines a region as either “ice-covered” or “not ice-covered”. In the model grid, for each grid cell, either the cell has ice (usually a value of “1”) or the cell has no ice (usually a value of “0”). Here, we apply a commonly used threshold of 15 % (such as used by National Snow and Ice Data Center; NSIDC) to determine the ice labelling, meaning that if the model grid cell has greater than 15 % ice concentration, the cell is labelled as “ice-covered”. Compared with piControl, the simulations show significant decreases in sea-ice extent (SIE), with anomalies in the Arctic peaking at −1.8 × 10⁶ km² in midHolocene, −5.2 × 10⁶ km² in lig127k and −11.8 × 10⁶ km² in midPliocene. SIE anomalies peak earliest in the year in the lig127k simulation, followed by midHolocene (Fig. 10c). In the midPliocene simulation, the Arctic experiences sea-ice-free conditions from August to October, and anomalies peak both in the (boreal) autumn and the spring. The different timing of the SIE anomaly peaks may be explained by a runaway sea-ice-albedo feedback (Feng et al., 2019) in which early summer reductions in SIE lead to more warming through decreased albedo and, in turn, even more significant reductions in SIE. Since the positive insolation changes in the Arctic occur earlier in the year in lig127k than in midHolocene, this feedback would be triggered earlier in the year in lig127k. The importance of this feedback depends on the annual amplitude of SIE. Since this is relatively small in midPliocene, it results in the weakest link between the seasonal SIE and SAT cycles. This may explain the absence of a distinct minimum in Arctic SAT anomalies in the boreal summer in the midPliocene simulation. The seasonal cycle of Antarctic SIE anomalies closely resembles the SAT cycle, and anomalies peak in September at −9.2 × 10⁶ km² in the midPliocene simulation. Anomalies are around zero throughout the year in the midHolocene and lig127k simulations, consistent with the limited changes in SAT observed at these latitudes throughout the year.

Figure 10 Seasonal changes of Arctic mean surface air temperature (a, b) and sea-ice extent (c, d) in the Arctic (a, c, 60–90^∘ N average) and Antarctic (b, d, 60–90^∘ S average) for each simulation.

The sea-ice edge in the months of maximum and minimum SIE in each simulation is depicted in Fig. 11. The maximum SIE in the Arctic is observed in March for each simulation (Fig. 11a). For the Arctic minimum SIE, observed in August, the sea-ice edge is limited to the Greenland Sea for the piControl and midHolocene simulations in the Atlantic domain, while it does not reach further than the Fram Strait in the lig127k simulation (Fig. 11c). In the Antarctic, the maximum SIE is in September (Fig. 11b). Sea ice is present along most of Antarctica's coast during September in the midPliocene simulation (Fig. 11b). The minimum SIE is observed in February. The sea-ice edge is mainly limited to the Weddell Sea for the piControl, midHolocene and lig127k simulations, while no sea ice is present in the midPliocene simulation (Fig. 11d). The midPliocene SIE anomalies are the lowest in the summer because low sea-ice concentrations during this season limit SIE decrease in the mid-Pliocene, while SIE can still decrease more in the PI simulation.

Figure 11 Sea-ice edge in the Arctic (a, c) and Antarctic (b, d) in the maximum (a, b) and minimum (c, d) SIE months in each simulation. Maximum sea-ice extents occur in March in the Arctic and in September in the Antarctic, while minimum sea-ice extents are seen in August in the Arctic and in February in the Antarctic.

4.5 Responses in global monsoon systems

We compare the global monsoon area in different simulations in Fig. 12. Here, the global monsoon area is defined as the regions where the annual range precipitation exceeds 2 mm d⁻¹ and the local summer precipitation exceeds 55 % of the annual precipitation (Liu et al., 2009). The annual range refers to the precipitation difference between MJJAS (May to September) and NDJFM (November to March) in the NH and difference between NDJFM and MJJAS in the SH. The precipitation averaged for the 5 months over the monsoon area represents the global monsoon intensity (Zhou et al., 2008).

There is little difference in land monsoon regions between piControl and midHolocene (Fig. 12b), but there is still an increased area globally by 9 % in midHolocene with the majority of this increase (85 %) located in the SH (Figs. 12b and 13a). The most significant difference is found in the southern Indian Ocean, where South African and Australian monsoon regions meet (Fig. S4). However, the SH monsoon intensity decreases in the same period (Fig. 13b). The expansion in this monsoon region is not linked to any visible change in summer temperature or sea level pressure (Fig. S5a and b). It is likely caused by the precipitation decrease over the eastern Indian Ocean during the SH winter being more extensive than the reduction during the SH summer, leading to stronger seasonal differences (Fig. S4a and b).

Figure 12 Global monsoon area in the (a) piControl, (b) midHolocene, (c) lig127k, (d) the midPliocene simulations, with the piControl global monsoon boundary as red lines in panels (b)–(d).

The last interglacial experiences northward expansions of NH monsoon regions (northern Africa, Asia and North America) compared to pre-industrial conditions. In contrast, the expansions in the SH monsoon regions are mainly located over the ocean (Fig. 12c). Although the monsoon area increases in both hemispheres, the intensity has only increased in the NH and decreased in the SH (Fig. 13b). MJJAS exhibits positive temperature anomalies and negative pressure anomalies over land in lig127k, favouring a strengthening of monsoons (Webster et al., 1998). NDJFM experiences a global decrease in temperature except for the Arctic region, positive pressure anomalies over land and negative pressure anomalies over the ocean, leading to a weakening of the negative land–ocean pressure gradients often linked to monsoonal flows. This result is consistent with the PMIP4-CMIP6 lig127k ensemble, showing that the insolation changes during lig127k lead to an intensification of boreal summer monsoons in the NH and a weakening of the austral summer monsoons in the SH (Otto-Bliesner et al., 2021).

Figure 13 Global monsoon area (a) and monsoon intensity (b) for the NH, SH and globally in the four simulations.

The global monsoon area increases by 8 % during the midPliocene simulation, which comes mainly from the NH (75 %), while the increase in monsoon intensity comes mostly from the SH (9 % in the SH compared to 3 % in the NH) (Fig. 13). Compared to the other two interglacial periods, midPliocene is the only period with a higher monsoon intensity in the SH than in the NH. All NH monsoon regions except North America experience a northward expansion in MJJAS (Fig. 12d), consistent with strengthened latitudinal and land–sea temperature contrast and deepening of the low-pressure areas over Eurasia, northern Africa and the North Atlantic (Fig. S4e). It is also consistent with the northward expansion of the Hadley circulation (Fig. 8e). The precipitation intensity increase found in SH monsoon regions (NDJFM) is not linked to any strengthened land–sea pressure gradients (Fig. S4f). It might result from both a slightly strengthened land–sea thermal gradient for South America and southern Africa and a general temperature increase leading to higher specific humidity and increased moisture transport to the monsoon regions.

4.6 Responses in the modes of climate variability

Lastly, we examine the simulated climate variability modes on interannual and longer timescales associated with changes in sea surface temperature patterns, namely, ENSO, the Pacific Decadal Oscillation (PDO) and the Atlantic Multidecadal Oscillation (AMO). We define the interannual variability of monthly SST anomalies in the tropical Pacific basin as ENSO spatial pattern and the leading empirical orthogonal function (EOF1) time series as the ENSO index. The PDO is calculated by the leading EOF of monthly SST anomalies in the North Pacific and PC1 time series following (Mantua and Hare, 2002). Similarly, the AMO is estimated by the leading EOF of monthly SST anomalies in the North Atlantic (Guan and Nigam, 2009).

Figure 14 Simulated SST interannual variability (standard deviation of SST anomalies after removing seasonal cycle, unit: ^∘C) for all simulations. Black boxes are the Niño3.4 and Atlantic Niño regions (for the El Niño–Southern Oscillation and Atlantic Niño, respectively).

Figure 15 Simulated PDO patterns. The spatial patterns show the leading empirical orthogonal function (EOF) of SST anomalies over the North Pacific (after removing the global mean SST anomaly). Note that the EOF calculation is restricted to the North Pacific (outlined by the black box). The global pattern is the regression coefficient of the monthly SST anomalies at each location onto the normalized principal component (PC) time series (unit: ^∘C/^∘C). The time series shows the associated PC time series (blue line) and the 10-year running mean (red line).

Figure 16 Simulated AMO patterns. The spatial patterns and time series are derived using the same EOF and regression method as PDO but for the North Atlantic region (outlined by the black box).

The simulations show a decline of ENSO intensity from piControl to midHolocene, lig127k and midPliocene, which is universal across the equatorial central/eastern Pacific (Fig. 14). The Atlantic Niño variability is consistently weakened. The reduced ENSO variability during the mid-Holocene is consistent with previous studies, either from proxy reconstructions (see a synthesis in Emile-Geay et al., 2016) or model simulations (see a review by Lu et al., 2018). The larger insolation anomalies in lig127k compared to midHolocene (Fig. 1) contribute to the more robust weakening of ENSO variability through reduced Bjerknes feedbacks (Liu et al., 2014; Lu et al., 2019).

The simulated PDO patterns generally reproduce the positive SST anomalies comparably in the equatorial eastern Pacific and the northeast Pacific and negative anomalies in the central and western North Pacific (Fig. 15). This internal PDO variability in paleoclimate conditions is consistent with a previous mid-Holocene study using earlier-generation climate models (An and Park, 2013). However, different simulations show remarkable changes in the relative amplitude: in midHolocene and lig127k, the anomalies in the tropics are much larger and are almost equal to those in the northeast Pacific. In midPliocene, the overall anomalies are weaker and stretch westward to the west coast of the tropical Pacific, contributing to the westward shift of the Walker circulation (Fig. 8f).

The AMO patterns in the EC-Earth3-LR simulations are characterized by a basin-wide mode in the North Atlantic (Fig. 16). The low-frequency mode has also been simulated in an earlier study during the Holocene (Wei and Lohmann, 2012). These anomalies also have a horseshoe shape (Gastineau and Frankignoul, 2015) with larger subpolar and subtropical anomalies. The spatial patterns of the AMO are mostly unchanged in the piControl, lig127k and midPliocene simulations but seem more predominant in the midHolocene simulation.

The spectral features of the simulated ENSO, PDO and AMO are detected using the power spectrum method. The power spectra show that the dominant periodicities of ENSO, PDO and AMO lie in the 2- to 7-, 8- to 20- and 20- to 40-year timescales (Fig. S6); all are statistically significant, although the spectral peak differs slightly. These periodicities are consistent with present-day observations (see a review by Deser et al., 2010) and highlight the model's capability to reproduce the representative natural internal variability.

5 Conclusions

This paper documents the four PMIP4-CMIP6 experiments performed by EC-Earth3-LR by providing the information on the model, experimental setups and results on large-scale climate responses. The model outputs from the four simulations are published on ESGF for more applications by the research community. We expect that these data will be further used in studies through model–data comparison to understand the past climate change and climate variability, model–model and model–data comparisons to improve the climate model performance, or multi-period comparison to obtain an overview of the climate history.

We have followed PMIP4 protocol and adapted the climate forcing and boundary conditions in EC-Earth3-LR. For the two interglacial simulations, most of the boundary conditions are prescribed as the pre-industrial conditions. The boundary conditions such as vegetation, aerosols and dust forcings, are essential to sensitivity studies and should be appropriately prescribed during these climate conditions in the future simulations. We try to perform the spin-up runs as long as possible, but due to limited computation resources and time constraints, the deep ocean may not be fully at equilibrium in some of the simulations. The surface parameters have reached an acceptable equilibrium state as required by PMIP4.

EC-Earth3-LR simulates reasonable climate response during past warm periods, as shown in paleo-reconstructions and the other PMIP4-CMIP6 model ensemble. The different forcings cause these warm periods: the insolation change forces the two interglacial simulations, and the mid-Pliocene simulation is forced by the high CO₂ concentration comparable to today. The results show that the climate response in the latter forcing is dramatically more significant than the former ones, implying the ongoing global warming is critical in the context of long climate history.

Climate changes are amplified over the land in the higher latitudes, especially in the Arctic region. Arctic warming is apparent in all simulations with a strong link between the seasonal cycles of surface air temperature and sea ice and a strong link between the surface air temperature and insolation anomalies in the midHolocene and lig127k simulations. The climatic response in the Antarctic region is less prominent, and only small changes are observed in the two interglacial simulations.

The temperature gradient changes led to a poleward expansion of Hadley circulation in all three warm period simulations, with reduced intensity in the Northern Hemisphere. In the Southern Hemisphere, the Hadley circulation strengthened in midHolocene and lig127k but weakened in midPliocene. We also found a strengthening and westward shift of the Pacific Walker circulation in all three periods related to the ENSO changes.

The global monsoon area increased in midHolocene compared to piControl, with most of the increase located in the SH. However, the SH monsoon intensity decreased due to stronger seasonal precipitation differences. In lig127k, the global monsoon area expanded northward, and the monsoon intensity increased in the NH while it decreased in SH due to changes in the land–ocean thermal contrast in both hemispheres. For midPliocene, both the global monsoon area and intensity increased, where the NH mainly contributed to the area increase and the SH to the intensity increase.

EC-Earth3-LR reproduces reasonable spatial patterns and frequency features of several intrinsic climate variability modes such as ENSO, PDO and AMO. ENSO variability is weakened from the pre-industrial era towards midHolocene, lig127k and midPliocene, in response to intensified insolation changes in seasonality. The changes in PDO and AMO are more subtle and complex and will need further investigation.