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.

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 1PMIP4-CMIP6 protocols on forcings and boundary conditions.

Download Print Version | Download XLSX

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

Download

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 2Prescribed 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 2The 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.

Download Print Version | Download XLSX

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

Download

Figure 4Evolution 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.

Download

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

Download Print Version | Download XLSX

Figure 5Mean 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 6Mean 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 7Mean 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 8The 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⁻¹.

Download

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

Download

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

Download

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 11Sea-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 12Global 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 13Global monsoon area (a) and monsoon intensity (b) for the NH, SH and globally in the four simulations.

Download

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 14Simulated 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 15Simulated 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 16Simulated 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.