Transient climate simulations of the Holocene (version 1) – experimental design and boundary conditions

The Holocene, which started approximately 11.5 ka, is the latest interglacial period with several rapid climate changes with timescales, from decades to centuries, superimposed on the millennium-scale mean climate trend. Climate models provide useful tools to investigate the underlying dynamic mechanisms for the climate change during this well-studied time period. Thanks to the improvements in the climate model and computational power, transient simulation of the Holocene offers an opportunity to investigate the climate evolution in response to time-varying external forcings and feedbacks. Here, we present the design of a new set of transient experiments for the whole Holocene from 11.5 ka to the preindustrial period (1850; HT-11.5 ka) to investigate both the combined and separated effects of the main external forcing of orbital insolation, atmospheric greenhouse gas (GHG) concentrations, and ice sheets on the climate evolution over the Holocene. The HT-11.5 ka simulations are performed with a relatively high-resolution version of the comprehensive Earth system model CESM1.2.1 without acceleration, both fully and singly forced by time-varying boundary conditions of orbital configurations, atmospheric GHGs, and ice sheets. Preliminary simulation results show a slight decrease in the global annual mean surface air temperature from 11.5 to 7.5 ka due to both changes in orbital insolation and GHG concentrations, with an abrupt cooling at approximately 7.5 ka, which is followed by a continuous warming until the preindustrial period, mainly due to increased GHG concentrations. Both at global and zonal mean scales, the simulated annual and seasonal temperature changes at 6 ka lie within the range of the 14 PMIP4 model results and are overall stronger than their arithmetic mean results for the Middle Holocene simulations. Further analyses on the HT-11.5 ka transient simulation results will be covered by follow-up studies.

1 Introduction

1.1 Climate evolution over the Holocene

Since the end of the Younger Dryas cooling event at ∼11.5 ka, the global climate has entered the Holocene period. The Holocene is the latest interglacial period in the Quaternary glacial–interglacial cycles, which has experienced strong global environmental changes mainly due to the retreat of continental ice sheets, changes in atmospheric greenhouse gas (GHG) concentrations, and variations in the seasonal insolation associated with changes in the Earth's orbital parameters and solar variability (Mayewski et al., 2004). Although previous attempts to subdivide the Holocene still remain inconclusive, the Holocene can roughly be divided into the following three periods (Nesje and Dahl, 1993): the Early Holocene, with a high boreal summer insolation and a large remnant ice sheet in North America, the Holocene Hypsithermal, with a continuing high boreal summer insolation but a small remnant North American ice sheet, and the Late Holocene Neoglacial, with a declining boreal summer insolation (Wanner and Ritz, 2011; Wanner et al., 2011; Walker et al., 2012).

The Early Holocene, approximately spanning from 11.5 to 9 ka, represents the last transition from full glacial to interglacial conditions, which is characterized by a warming trend, mainly in the northern mid- and high-latitudes, as revealed from numerous reconstructions (e.g., Koerner and Fisher, 1990; Yuan et al., 2004; Wang et al., 2005; Hald et al., 2007; Vinther et al., 2008; Brooks et al., 2012; Shakun et al., 2012; Birks, 2015). In tropical regions, the change in precipitation was more prominent during the Early Holocene, among which the African Humid Period is one typical example (Gasse, 2001; Hoelzmann et al., 2004), in addition to a humidity climate in the Arabian Sea, the Bay of Bengal (Gupta et al., 2003; Staubwasser, 2006), and northern South America (Haug et al., 2001), which was suggested to be fundamentally affected by precession (e.g., Renssen et al., 2004). However, a drier climate was inferred from multi-evidence in central and southern South America (Maslin and Burns, 2000; Seltzer et al., 2000; Núñez et al., 2002) and westerlies-dominated North Asia, whereas monsoon-influenced East Asia was substantially wetter but with an asynchronous peak of precipitation in northeastern, northwestern, and northern China (An et al., 2000; Chen et al., 2008, 2019; Zhang et al., 2011).

The Holocene Hypsithermal, covering from about 9 to 5–6 ka, generally refers to the interval of warmth associated with the peak Holocene temperature, which is also called the Holocene Thermal Maximum, the Holocene Climatic Optimum, Altithermal, or Megathermal (Deevey and Flint, 1957; Wanner et al., 2008; Ljungqvist, 2011). The timing and spatial pattern of this period in different regions have long attracted the attention of researchers. Relative to the preindustrial period, the reconstructed Holocene maximum temperature appeared during 10–8 ka in the North Atlantic and the adjacent polar regions, which occurred during 7–5 ka in northern Europe and northwestern North America, and registered between Early and Middle Holocene at the northern mid- and high-latitudes, e.g., during 8.5–3.0 ka in China (Shi et al., 1993), and appeared to be time-transgressive across the western hemisphere of the Arctic (Kaufman et al., 2004), suggesting that these local warm periods were very likely not globally synchronous (Jansen et al., 2007). Different from the Northern Hemisphere, the southern mid- and high-latitudes experienced a drop in temperature during the Middle and Late Holocene, following the Early Holocene warming period (Shevenell et al., 2011). Concerning the moisture conditions, paleoenvironmental data consistently showed that the mid-continental extra-tropics were wetter, with cooler summers than today, during the Middle Holocene (Harrison et al., 2015).

The Late Holocene Neoglacial or Neoglaciation, lasting from about 5–6 ka to the preindustrial period, features an obvious increase in the number of glacial advances worldwide, mainly due to decreasing solar forcing in the boreal summer (Solomina et al., 2015). The global mean surface temperature during this period displayed a consistent cooling trend, as indicated from multi-proxy-based reconstructions (Wanner et al., 2008; Ljungqvist, 2011; Marcott et al., 2013; Kaufman et al., 2020a, 2020b; Zhang et al., 2022), which was suggested to be mainly due to the influence of orbital forcing together with the melting ice sheet in North America. In contrast, the Late Holocene warming trend was also supported by simulations (e.g., Liu et al., 2014, 2018) and by several reconstructions from subfossil pollen (Marsicek et al., 2018) and stalagmite records during the winter season (Baker et al., 2017), which was mainly attributed to rising atmospheric greenhouse gases (GHGs) and winter insolation. Therefore, it still remains inconclusive as to whether the Late Holocene was cooler or warmer. To reconcile this so-called Holocene temperature conundrum, further efforts should be taken to improve, both in simulations, by including solar irradiance forcing, land use changes, dust effects, and stratospheric chemistry and dynamics, and in reconstructions with higher spatial and temporal resolutions, a balanced distribution between marine and terrestrial and between the Northern Hemisphere and Southern Hemisphere data, and the representing proxies from both summer and winter seasons (Bova et al., 2021; Wanner, 2021; Zhang et al., 2022), as well as in the model–data assimilation approaches (Osman et al., 2021).

Superimposed on the millennium-scale mean climate trend during the above Holocene periods are several rapid cooling events on the timescales from decades to centuries (Wanner et al., 2011; Walker et al., 2012). Most of those events are characterized by cooling at northern mid- and high-latitudes, aridity at low-latitude monsoon regions, and major atmospheric circulation changes (Mayewski et al., 2004), which were generally inferred to relate to the freshwater outbreak caused by ice melting in the North Atlantic before 8.2 ka, the solar variability superimposed on long-term changes in insolation during the Middle and Late Holocene, and explosive volcanic eruptions and fluctuations in the thermohaline circulation (Mayewski et al., 2004; Magny et al., 2007; Wanner et al., 2011).

1.2 Transient modeling of the Holocene

Transient modeling of the Holocene offers an opportunity to examine the climate evolution in response to time-varying external forcings and feedbacks. Combined with proxy data, it is the most useful tool to comprehensively understand the causes for past complex and interrelating climate processes. Since the 1990s, the Paleoclimate Modeling Intercomparison Project (PMIP) has served to coordinate paleoclimate simulations and model–data comparisons (Joussaume and Taylor, 1995; Braconnot et al., 2007; Taylor et al., 2012), and it now enters its fourth phase (PMIP4; Kageyama et al., 2018), contributing to the current sixth phase of the Coupled Model Intercomparison Project (CMIP6; Eyring et al., 2016) and the Intergovernmental Panel on Climate Change (IPCC) Sixth Assessment Report (IPCC, 2022). The major focuses in PMIP4 include not only equilibrium (or time slice) simulations to investigate the impact of changes mainly in orbital forcing, e.g., during the Middle Holocene (6 ka), but also transient simulations of the last deglaciation from 21 to 9 ka with time-varying orbital forcing, GHGs, ice sheets, and other geographical changes (Ivanovic et al., 2016) and of the Holocene from 6 to 0 ka with temporal changes in orbital forcing and GHGs (Otto-Bliesner et al., 2017).

Up to now, a dozen transient simulations covering the Holocene have been carried out using global climate models with different complexities. For example, FOAM, a fully coupled global atmosphere–ocean general circulation model with a low resolution of $\mathrm{4}{}^{\circ }\times \mathrm{7.5}{}^{\circ }$ (latitude × longitude), was used to run a transient simulation for the past 280 ka to examine the evolutionary response of global monsoons to orbital forcing, which is accelerated by a factor of 100 (Kutzbach et al., 2008). FAMOUS (Smith et al., 2008), a low-resolution ($\mathrm{5}{}^{\circ }\times \mathrm{7.5}{}^{\circ }$) version of the HadCM3 (Gordon et al., 2000) coupled model, was used to conduct a number of accelerated transient experiments (by a factor of 10) for the last 120 ka to investigate the climate of the last glacial cycle (Smith and Gregory, 2012). A set of transient simulations of the last deglaciation for the past 26 ka was performed with a low-resolution ($\sim \mathrm{3.75}{}^{\circ }\times \mathrm{3.75}{}^{\circ }$) Earth system model MPI-ESM1.2 to explore differences in climate response to different underlying ice sheet reconstructions and methods of meltwater distribution (Kapsch et al., 2020), as suggested within the PMIP4 deglaciation protocol (Ivanovic et al., 2016). The Transient Climate Evolution of the past 21 ka (TraCE-21 ka; He, 2011), run with the low-resolution ($\sim \mathrm{3.75}{}^{\circ }\times \mathrm{3.75}{}^{\circ }$) National Center for Atmospheric Research (NCAR) Community Climate System Model (CCSM) version 3 (CCSM3), is widely used for examining the climate evolution since the last glacial maximum at global and regional scales (e.g., Liu et al., 2014, 2021; Wu et al., 2021). A fully coupled atmosphere–ocean–sea ice model ECBilt-CLIO, with a horizontal resolution of $\sim \mathrm{5.6}{}^{\circ }\times \mathrm{5.5}{}^{\circ }$, was also used for transient simulations of the last deglaciation (last 21 ka) to investigate the effect of using accelerated (by a factor of 10) boundary conditions and of uncertainties in the initial state (Timm and Timmermann, 2007). Additionally coupled with a vegetation model, ECBilt-CLIO-VECODE3 was performed to investigate the global characterization of the Holocene thermal maximum for the last 9 ka (Renssen et al., 2012). A later-developed version, LOVECLIM, additionally coupled with carbon cycle and terrestrial ice sheets (Goosse et al., 2010), was employed to examine the Early Holocene (11.5–7 ka) warming (Zhang et al., 2016), the Holocene northern temperature and hydroclimate trends (Zhang et al., 2018, 2020), and the obliquity and CO₂ effects on southern climate during the past 408 ka with an acceleration factor of 5 (Timmermann et al., 2014).

Concentrating more on the Holocene, transient simulations for the past 12 ka, run with an acceleration factor of 10 with the low-resolution ($\sim \mathrm{3.75}{}^{\circ }\times \mathrm{3.75}{}^{\circ }$) Community Earth System Model (CESM) version 1.0.3 (CESM1.0.3), showed a cooling trend during 5–0.15 ka by considering the volcanic forcing (Wan et al., 2020) and an increasing trend in the northern mid-latitude precipitation during 7–0 ka in response to the orbital forcing (Sun et al., 2020). A transient simulation of the last 8 ka, performed with the high-resolution ($\mathrm{1.875}{}^{\circ }\times \mathrm{1.875}{}^{\circ }$) MPI-ESM1.2, additionally forced by solar variability and volcanic forcing, further revealed that both global warming and cooling modes coexisted during the Holocene (Bader et al., 2020), which also captured the global vegetation transitions and the time-transgressive end of the African Humid Period, as revealed in proxy data (Dallmeyer et al., 2020, 2021). Using a low-resolution ($\sim \mathrm{3.75}{}^{\circ }\times \mathrm{3.75}{}^{\circ }$) atmosphere–ocean–sea ice coupled model, KCM, with an acceleration factor of 10, a transient climate simulation of the past 9.5 ka was conducted to examine the summer precipitation across the Asian monsoon region in response to the orbital forcing (Jin et al., 2014). Besides, two versions of the IPSL model were used to run the transient simulations of the last 6 ka, indicating larger impacts of different model versions and experimental setups on the mean state of climate than on the associated interannual-to-decadal variability (Braconnot et al., 2019a, b; Crétat et al., 2020). In addition, the four Holocene transient simulations run with MPI-ESM1.2, two versions of the IPSL model, and AWI-ESM2 (Sidorenko et al., 2019) were compared to support the orbital forcing of the El Niño–Southern Oscillation (ENSO) amplitude since the Middle Holocene (Carré et al., 2021) and indicate the importance of the application of calendar transformation in the analysis of climate simulations when we perform multimodel comparisons (Shi et al., 2022).

1.3 A new set of Holocene transient experiments

Taken together, among the abovementioned Holocene transient simulations, only seven of them covered the whole Holocene of the past 11.5 ka, i.e., simulations with LOVECLIM for the past 408 ka, FOAM for the past 284 ka, FAMOUS for the past 120 ka, MPI-ESM1.2 for the past 26 ka, CCSM3 and ECBilt-CLIO for the past 21 ka, and CESM1.0.3 for the past 12 ka. They were performed with climate system models or Earth system models with intermediate complexity, all of which have relatively low spatial resolutions varying from $\sim \mathrm{5.6}{}^{\circ }\times \mathrm{5.5}{}^{\circ }$ to $\sim \mathrm{3.75}{}^{\circ }\times \mathrm{3.75}{}^{\circ }$ and/or use acceleration techniques with a factor from 5 to 100 to reduce computational resources. Note that, on the one hand, the model resolution determines the level of detail in its representation of the physical processes related to the topography, heat transport, sea ice cover, etc., which may have large effects on the Holocene climate evolution (Zhang et al., 2018). On the other hand, acceleration in transient simulations can also have a significant impact on the local climate history, for example, which will lead to damped and delayed temperature response to the boundary conditions in the North Atlantic (Timm and Timmermann, 2007). In addition, most models used for the Holocene transient simulations are fully forced by time-varying boundary conditions of orbital configurations, GHGs, and ice sheets, some of which are additionally forced by meltwater flux, land use/land cover, solar variability, or volcanic forcing, making it difficult to isolate their relative contribution to the past climate evolution.

In this context, it is urgent to carry out a new set of unaccelerated transient experiments focusing on the whole Holocene, since 11.5 ka, using the comprehensive Earth system model with relatively high spatial resolutions, considering both the full-forcing and single-forcing effects of several of the most important boundary conditions, including the orbital parameters, GHG concentrations, and ice sheets, on the Holocene climate evolution. This set of Holocene transient simulations for the last 11.5 ka will also fill in the gap between the PMIP4 transient Core experiments for the last deglaciation from 21 to 9 ka (Ivanovic et al., 2016) and the PMIP4 Tier 3 transient simulations for the Holocene from 6 to 0 ka (Otto-Bliesner et al., 2017).

The aim of this paper is to outline the model description (in Sect. 2), experimental setup, and boundary conditions (in Sect. 3) for the new Holocene transient climate experiments of the last 11.5 ka (henceforth HT-11.5 ka for brevity), with some preliminary results being provided in Sect. 4. Further analyses on these simulation results will be covered by follow-up studies.

2 Model description

The model used to run the HT-11.5 ka simulation in this study is the NCAR CESM version 1.2.1 (CESM1.2.1), which was released on December 2013 and was a relatively new version when we started this simulation in 2016. CESM is a fully coupled global Earth system model that provides state-of-the-art simulations of the Earth's past, present, and future climate states, with an earlier version 1 (CESM1) participating in the CMIP phase 5 (CMIP5; Taylor et al., 2012) and the later version 2 (CESM2) participating in the CMIP6 (Eyring et al., 2016). CESM1.2.1 consists of four components and can be coupled in different configurations. The atmospheric component we used here is the NCAR Community Atmosphere Model version 4 (CAM4; Neale et al., 2010), also used in the CCSM version 4 (CCSM4; Gent et al., 2011), which employs a finite volume dynamical core at a horizontal resolution of approximately $\mathrm{1.9}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$ (latitude 96× longitude 144 grid points) with 26 vertical layers. It contains notable improvements in the deep convection parameterization, with inclusions of more realistic dilution effects through an explicit representation of entrainment and subgrid-scale convective momentum transports, in the cloud fraction parameterization with a freeze-drying process contributing to a greater consistency between polar cloud fraction and water condensate properties, and in the radiation interface and computational scalability (Neale et al., 2010). The land component is the Community Land Model version 4.0 (CLM4; Lawrence et al., 2011) with the carbon nitrogen–dynamic global vegetation model turned off, which is also used in CCSM4 and run on the same horizontal grid as the atmospheric component. The oceanic component is the Parallel Ocean Program version 2 (POP2; Danabasoglu et al., 2012), based on the POP version 2.1 of the Los Alamos National Laboratory (Smith et al., 2010) as utilized in CCSM4, which uses an approximately 1^∘ horizontal grid with 60 vertical levels. Compared to previous model versions, the physical improvements in POP2 include a near-surface eddy flux parameterization, an abyssal tidally mixing parameterization, an overflow parameterization to represent the Denmark Strait, Faroe Bank Channel, Weddell Sea, and Ross Sea overflows, a submesoscale mixing scheme, vertically varying thickness and isopycnal diffusivity coefficients, and modified and generally reduced horizontal viscosity coefficients (Hurrell et al., 2013). The sea ice component is the Los Alamos Sea Ice Model (or referred to as the Community Ice Code) version 4.0 (CICE 4.0; Hunke and Lipscomb, 2008) with the same horizontal grid as the ocean component. These four components are coupled by the version 7 coupler (CPL7; Craig et al., 2012). More details on the model are available online at https://www.cesm.ucar.edu/models/cesm1.2/ (last access: 20 February 2022) and in Hurrell et al. (2013). Further information on the model performance, scalability, and throughput for different model components distinguishing between input/output, ocean and atmosphere costs, etc., are documented in the CCSM4/CESM1 special collection of the Journal of Climate (see http://journals.ametsoc.org/collection/CCSM4-CESM1, last access: 5 May 2022), and additional model simulations for the past and present climates can be found, for instance, in Goldner et al. (2014), Song and Zhang (2018, 2019), and Park et al. (2019).

3 Experimental designs for the HT-11.5 ka simulations

The HT-11.5 ka simulations of the Holocene focus on the period from 11.5 ka to the preindustrial period (1850), which will be the supplement for the PMIP4 Core simulations of the last deglaciation spanning from 21 to 9 ka (Ivanovic et al., 2016). Note that 0 ka is 1950, and 1850 is equivalent to 0.1 ka. All simulations spin up from the Early Holocene at 11.5 ka, and all boundary conditions are available from 11.5 ka to the preindustrial period. The initialization state at 11.5 ka and boundary conditions since 11.5 ka, including Earth's orbital parameters, atmospheric GHG concentrations, and ice sheets and orography, are summarized in Table 1 and described below. Boundary condition data for the HT-11.5 ka transient experiments are provided on the PMIP4 wiki for the transient deglaciation (PMIP Last Deglaciation Working Group, 2016).

Table 1 Summary of model boundary conditions for the Early Holocene (11.5 ka) spin-up and Holocene transient (11.5 ka to 1850) experiments; see the text for details. Data are available from PMIP4 last deglaciation wiki (PMIP Last Deglaciation Working Group, 2016). Boundary condition group headings are in bold. Note: ppm is parts per million and ppb is parts per billion.

3.1 Early Holocene spin-up

The Early Holocene spin-up experiment for the HT-11.5 ka simulations is an equilibrium or time slice simulation with the boundary conditions of orbital parameters, GHGs, and ice sheets prescribed at 11.5 ka (Table 1). Orbital parameters of eccentricity, obliquity, and longitude of perihelion are set to 0.01957212, 24.17958, and 270.2132^∘, respectively, with the date of vernal equinox fixed at the noon on 21 March. The solar constant is fixed to 1365 W m⁻², which is consistent with the preindustrial value for the PMIP phase 3 (PMIP3) and CMIP phase 5 (CMIP5) experiments but higher than the value recommended for the transient deglaciation simulation (e.g., 1361.0±0.5 W m⁻²; Mamajek et al., 2015) and the PMIP4/CMIP6 experiments (e.g., 1360.8±0.5 W m⁻²; Matthes et al., 2017; Otto-Bliesner et al., 2017). The atmospheric GHGs are prescribed to 267.4 ppm (parts per million) for carbon dioxide (CO₂), 652.4 ppb (parts per billion) for methane (CH₄), 251.1 ppb for nitrous oxide (N₂O), 0 for chlorofluorocarbons (CFCs), and the PMIP3/CMIP5 preindustrial (1850) value for ozone (O₃). These values are compatible with the time-evolving boundary conditions for the HT-11.5 ka simulations and consistent with the latest ice core age model (AICC2012; Veres et al., 2013) and records (Loulergue et al., 2008; Schilt et al., 2010; Bereiter et al., 2015). Prescribed ice sheets and the associated topography use the 11.5 ka data from the ICE-6G_C (VM5a) reconstruction (henceforth ICE-6G_C; Argus et al., 2014; Peltier et al., 2015). The initial ocean condition is taken from an archived 500-year spin up (CCSM4; 2^∘ resolution in the atmosphere and 1^∘ resolution in the ocean and preindustrial control; case b40.1850.track1.2deg.003). All other boundary conditions, such as the coastlines, bathymetry, vegetation and land cover, and aerosols, are prescribed at the preindustrial levels.

Forced by the above boundary conditions, the Early Holocene spin up is continuously integrated for 1500 years, with the absolute value of the trend in global mean surface air and sea surface temperatures being less than 0.05 ^∘C and the Atlantic Meridional Overturning Circulation being stable (approximately 25.4 Sv) for the last 100 years, which means the model has reached a quasi-equilibrium state, as suggested by Kageyama et al. (2018). All Holocene transient experiments start from this equilibrium state at 11.5 ka and run transiently until the preindustrial period at 1850.

3.2 Holocene transient simulations

3.2.1 Orbital parameters, solar constant, and insolation anomalies

The Earth's orbital parameters (eccentricity, obliquity, and longitude of perihelion) are time-evolving through the Holocene from 11.5 ka to the preindustrial period, following Berger (1978), which affects the seasonal and latitudinal distribution and magnitude of solar radiation received at the top of the atmosphere. Here, all orbital parameters are interpolated from the originally temporal resolution of 1000 years to 1 year to produce continuous orbital forcings as model inputs. More specifically, the eccentricity declined through the Holocene, indicating a gradual decrease in the seasonal difference of incoming insolation since 11.5 ka. The obliquity was maximal at 9 ka, when the insolation difference was largest between low and high latitudes. The perihelion occurred close to the boreal summer solstice, autumnal equinox, and winter solstice at 11.5, 6 ka, and the preindustrial period, respectively. Due to the above changes, the summer (June–July–August) insolation average in the Northern Hemisphere slightly increased from 465 W m⁻² at 11.5 ka to a maximum of 468 W m⁻² at about 9.3 ka, and then continuously declined until the preindustrial period. June insolation at 60^∘ N gradually decreased through the Holocene, whereas December insolation at 60^∘ S increased since 11.5 ka and was maximal around 2 ka (Fig. 1a). Relative to the preindustrial period, June insolation anomalies generally increased with latitude and decreased with time, with positive maximums of ∼50 W m⁻² occurring at northern high latitudes during the Early Holocene (Fig. 2).

Figure 1 Temporal evolution of the external forcings in the HT-11.5 ka simulations since the Early Holocene (11.5 ka), for (a) June insolation at 60^∘ N (solid line) and December insolation at 60^∘ S (dashed line; Berger, 1978), (b) atmospheric CO₂ concentration (recent composite of EPICA Dome C, Law Dome, and West Antarctic Ice Sheet Divide records, Antarctica; Bereiter et al., 2015), (c) atmospheric CH₄ concentration (EPICA Dome C, Antarctica; Loulergue et al., 2008), (d) Atmospheric N₂O concentration (Talos Dome, Antarctica; Schilt et al., 2010), and (e) area of the ice sheets in the Northern Hemisphere, according to the ICE-6G_C reconstruction (Argus et al., 2014; Peltier et al., 2015).

Figure 2 Latitude–time anomalies in the June insolation relative to the average for the last 1000 years, according to Berger (1978).

The solar constant is fixed to 1365 W m⁻² throughout the run, which is in line with the PMIP3 experiment setup (PMIP PI Working Group, 2010). As mentioned in Otto-Bliesner et al. (2017), this value is higher than the value for the solar constant used by the models in PMIP4 (1360.7 W m⁻²), which leads to a global increase in incoming solar radiation compared to the PMIP4 experiments.

3.2.2 Atmospheric GHGs

During the Holocene, CFCs is fixed at 0 and O₃ is set to the PMIP3/CMIP5 preindustrial value as the Early Holocene spin up. Other atmospheric GHGs are time-evolving through the Holocene and adjusted to the AICC2012 age model (Veres et al., 2013), with CO₂ following Bereiter et al. (2015), CH₄ following Loulergue et al. (2008), and N₂O following Schilt et al. (2010). The high-resolution atmospheric CO₂ concentrations are derived from a composite of EPICA Dome C (EDC; Monnin et al., 2001, 2004), the Antarctic composite ice core data built from the Law Dome (MacFarling Meure et al., 2006; Rubino et al., 2013), and the West Antarctic Ice Sheet Divide (Marcott et al., 2014), which show a slightly decreasing trend from 11.5 to 7 ka and an increasing trend afterward (Fig. 1b). The record of atmospheric CH₄ concentrations is built from the Antarctic EDC (Loulergue et al., 2008), with a decreasing trend from the Early to Middle Holocene followed by a rising trend (Fig. 1c). The low-resolution atmospheric concentrations of N₂O derived from the Talos Dome ice core (Schilt et al., 2010) show large variabilities during the Holocene, with relatively high values around 11 ka and during 6–0 ka but low levels between 11 and 6 ka (Fig. 1d). All the above GHG records are interpolated to the temporal resolution of 1 year to produce continuous forcings as model inputs.

Note that the chronologically accurate values of CO₂, CH₄, and N₂O from the above reconstructions slightly differ from the defined 6 ka concentrations in PMIP4. The nearest measured CO₂ concentration to 6 ka is 266.7 ppm from Bereiter et al. (2015), which is higher than the 6 ka value of 264.4 ppm in PMIP4 derived from the original EDC ice core (Monnin et al., 2001, 2004) after a smoothing spline (Otto-Bliesner et al., 2017). In contrast, both the nearest measured CH₄ and N₂O concentrations of 586 ppb (Loulergue et al., 2008) and 252.4 ppb (Schilt et al., 2010) produced on the AICC2012 chronology (Veres et al., 2013) are lower than the 6 ka values of 597 and 262 ppb in PMIP4 (Otto-Bliesner et al., 2017), respectively, on the original EDC1 chronology (Spahni et al., 2005).

3.2.3 Ice sheets and topography

The ice sheet extent and topography are time-evolving from 11.5 ka to the preindustrial period using the ICE-6G_C reconstruction (Argus et al., 2014; Peltier et al., 2015), which is based on glacial isostatic adjustment modeling and constrained by global positioning system measurements of the vertical motion of the crust, exposure age dating of ice thicknesses, relative sea level histories, and space-based gravity data by the satellite system (Peltier et al., 2015). During the Holocene, the most prominent change in ice sheets occurs in the Northern Hemisphere, especially in North America. According to the ICE-6G_C reconstruction, the North American ice sheet retreated gradually from 11.5 to 9 ka, with the Laurentide ice sheet holding ∼12.5 m global mean sea level rise equivalent ice volume at 10 ka (Matero et al., 2020), retreated sharply from 9 to 7 ka, and almost disappeared after 6 ka (Fig. 3) following the data provided by Dyke (2004). The Eurasian ice sheet has retreated since 11.5 ka and deglaciated after ∼10 ka, following Gyllencreutz et al. (2007), whereas the Greenland ice sheet was stable overall from 11.5 to 9 ka and slightly retreated at its marginal regions during 8–7 ka (Fig. 3). Averaged over the Northern Hemisphere, the ice area continuously declined from $\sim \mathrm{1.2}\times {\mathrm{10}}^{\mathrm{7}}$ km² at 11.5 ka to $\sim \mathrm{0.2}\times {\mathrm{10}}^{\mathrm{7}}$ km² during the Middle Holocene (∼7–5 ka) and has remained little changed since then (Fig. 1e). In contrast, there is little variation in the southern hemispheric ice sheets during the Holocene, with only a slight decrease in the ice extent from 11.5 to 10 ka around the Antarctica (Fig. 3). The orography evolves to be consistent with the ice sheet extent, with the elevation of ice sheet below 2000, 2500, and 3500 m over Eurasia, North America, and Greenland, respectively, and reaching 4000 m over Antarctica through the Holocene (Fig. 3). Besides, land surface properties are also adjusted for changes in ice extent. Where ice retreats, land surface is modified to the prescribed vegetation at the preindustrial period.

Figure 3 Northern Hemisphere (top) and Southern Hemisphere (bottom) ice sheet elevation at 11.5, 10, 8, 6, and 0 ka (1950) for the ICE-6G_C reconstruction at 10 arcmin horizontal resolution; elevation is plotted where the fractional ice mask is more than 0.5 (Peltier et al., 2015).

The original ICE-6G_C data set is provided at 10 arcmin ($\sim \mathrm{0.17}{}^{\circ }$) horizontal resolution and at intervals of 500 years through the Holocene. Both included ice extent and topography data are interpolated to the time interval of 250 years and updated every 250 years in the transient run. Considering that coastlines, bathymetry, and land–sea distributions are little varied during the Holocene, they are fixed as the preindustrial states in all runs. In addition, the meltwater flux change is not considered here, which is one of the options provided for the PMIP4 transient deglaciation experiments. All other forcings, such as the vegetation, land surface, and dust parameters, are prescribed at the preindustrial levels.

3.2.4 Experimental setup

In total, four transient experiments are run in the HT-11.5 ka simulations, namely one all-forcing experiment and three single-forcing experiments. The orbital parameters, atmospheric GHGs, ice sheets and topography, as detailed above, are synchronously evolving with time in the all-forcing experiment (hereinafter HT-ALL), while the three single-forcing experiments are forced by only one of the above forcings, with the other two forcings fixed at 11.5 ka (hereinafter HT-ORB, HT-GHG, and HT-ICE, respectively; Table 2). All four experiments are transiently run from 11.5 ka to the preindustrial period at 1850 for 11.4 kyr. Monthly output data are provided for each simulation.

Table 2 Boundary conditions for the HT-11.5 ka experiments.

4 Preliminary simulation results

On the millennium timescale, the annual global mean surface (2 m) air temperature (GMST) in the HT-ALL simulation slightly decreases from 11.5 ka (∼13.20 ^∘C) to 7.5 ka (∼12.93 ^∘C), with an abrupt cooling at approximately 7.5 ka, and then continuously increases until the preindustrial period (13.57 ^∘C; Fig. 4). The warming trend since the Middle Holocene simulated by CESM1.2.1 in our HT-11.5 ka transient experiment is similar to that in the CCSM3, LOVECLIM, and FAMOUS transient simulations (Liu et al., 2014) but opposite to the cooling trend as derived from the multi-proxy reconstructions (Marcott et al., 2013; Kaufman et al., 2020a), which is the so-called Holocene temperature conundrum. The abrupt cooling at approximately 7.5 ka is also seen in the LOVECLIM transient experiment (Liu et al., 2014). In addition, the GMST evolution from Early to Middle Holocene in our CESM1.2.1 transient simulation is somewhat different from the warming trend both in reconstructions and other Holocene transient simulations (Marcott et al., 2013; Liu et al., 2014; Kaufman et al., 2020a), and the reasons for the discrepancies need for further investigation.

Figure 4 Global annual mean surface air temperature since 11.5 ka under the all-forcing simulation. The time series is smoothed by a 101-year moving average. The gray circles at 6 and 0.1 ka stand for the Middle Holocene and preindustrial results simulated by the arithmetic mean of the 14 PMIP4 models, with gray vertical error bars representing plus and minus 1 standard deviation for 14 models.

During the preindustrial period, the annual GMST of 13.57 ^∘C in our HT-ALL simulation is very close to that of 13.62 ^∘C from the arithmetic mean of the 14 PMIP4 models, although with a spread among the individual models ranging from 12.59 ^∘C in IPSL-CM6A-LR to 14.98 ^∘C in MIROC-ES2L (Table 3). At 6 ka, the annual GMST is approximately 13.16 ^∘C, which resembles the value of 13.30 ^∘C from the Middle Holocene simulation in the 14 PMIP4 model arithmetic mean and is within the range from 12.20 ^∘C (IPSL-CM6A-LR) to 14.53 ^∘C (MIROC-ES2L) in the individual model simulations (Table 3). Relative to the average for 1750–1850, there is a global annual mean cooling of about 0.41 ^∘C at 6 ka in the HT-ALL simulation, which lies within the range of the results from the PMIP4 models. In the 14 PMIP4 models, the Middle Holocene cooling varied from 0.18 ^∘C in FGOALS-g3 to 0.49 ^∘C in FGOALS-f3-L with respect to the preindustrial period, with a standard deviation of 0.11 ^∘C among the models and an average of 0.32 ^∘C for the arithmetic mean of all models (Fig. 5; Table 3). The geographical distribution of the annual mean surface air temperature change at 6 ka in our HT-ALL simulation is further compared with that from the arithmetic mean of the 14 PMIP4 models (Fig. 6). Compared to the preindustrial period, there is an annual cooling worldwide, except for a slight warming in the high latitudes and over the western Europe, North Pacific, southern North America, and parts of tropical land at 6 ka in the HT-ALL simulation (Fig. 6a). The annual cooling is particularly strong (>1 ^∘C) over northern Africa and North America, southern Asia, and part of Southern Ocean. This large-scale distribution is similar to the spatial pattern of the change between the Middle Holocene and preindustrial simulations in the 14 PMIP4 model mean, although some differences exist in magnitude (Fig. 6a, b; Brierley et al., 2020). Quantitatively, the HT-ALL simulation shows an annual cooling of 0.18, 0.40, 0.61, 0.40, 0.31, and 0.20 ^∘C at 6 ka averaged for the zonal band of 60–90^∘ N, 30–60^∘ N, 0–30^∘ N, 30^∘S–0^∘, 60–30^∘ S, and 90–60^∘ S, respectively. All those values are within the spread of the 14 PMIP4 models, with the largest spread over high latitudes, but are slightly stronger than the cooling in their arithmetic mean by 0.05–0.14 ^∘C, especially in the Northern Hemisphere (Fig. 7a). Due to the large seasonal changes in insolation, there are notable seasonal differences in the surface air temperature change at 6 ka, with the December–January–February (DJF) mean change dominating the annual pattern. In response to the decreased insolation in DJF, a strong cooling occurred at 6 ka over the global land regions, especially over the northern continents excluding Greenland, and a weak cooling appeared over the global oceans except for a warming in the Arctic and Southern Ocean (Fig. 6c). In contrast, the increased June–July–August (JJA) insolation at 6 ka led to a notable warming over the northern mid- and high-latitude continents, a weak warming over the northern mid- and high-latitude oceans and southern continents, and a cooling in the other regions (Fig. 6e). The seasonal difference in the spatial distribution of the surface air temperature changes at 6 ka in the HT-ALL simulation resembles that for the Middle Holocene simulations in most PMIP4 models and their arithmetic mean (Fig. 6d, f; Brierley et al., 2020). In a quantitative manner, the HT-ALL simulation shows a consistent DJF cooling of 0.68 (0.31), 0.95 (0.22), and 0.83 (0.43) ^∘C averaged over the 30^∘ zonal band of northern (southern) high-, mid-, and low-latitudes, respectively (Fig. 7b). For JJA, it displays an average warming of 0.99 and 0.50 ^∘C over the northern high- and mid-latitudes but an average cooling varying from 0.41 to 0.66 ^∘C over the other four zonal bands (Fig. 7c). On average, there is a global mean cooling of 0.60 for DJF and 0.17 ^∘C for JJA (Table 3). All those zonal and global mean changes in seasonal temperatures are within the range of the Middle Holocene changes in the 14 PMIP4 models, which are consistent with the changes from the 14 PMIP4 model arithmetic mean in signs but overall larger than the latter in magnitudes (Fig. 7b and c and Table 3). More specifically, both PMIP4 model spreads and differences of zonal mean changes between the HT-ALL and 14 PMIP4 model mean simulations are largest at northern high latitudes in DJF and at high latitudes in both hemispheres in JJA. As a whole, both at global and zonal mean scales, the annual and seasonal mean changes in surface air temperature at 6 ka in our HT-ALL simulation lie within the range of the 14 PMIP4 model results and are overall stronger than their arithmetic means. In particular, the annual, DJF, and JJA mean GMSTs simulated by CESM1.2.1 in the HT-ALL simulation are lower than those by CESM2 in the PMIP4 both for 6 ka and preindustrial simulations (Table 3). However, the Middle Holocene annual and DJF cooling magnitudes simulated by CESM1.2.1 (0.41 and 0.60 ^∘C) are nearly double those in CESM2 (0.20 and 0.38 ^∘C), and the JJA cooling in the former (−0.17 ^∘C) is opposite to and 4 times the JJA warming in the latter (0.04 ^∘C; Table 3). There are also some differences, mainly in magnitudes, for the zonal mean changes in annual and seasonal temperatures at 6 ka between the two versions of CESM (Fig. 7).

Table 3 Global annual, DJF, and JJA mean surface air temperatures (units in ^∘C) for the preindustrial period, Middle Holocene, and the difference between the Middle Holocene and preindustrial periods in simulations from the 14 PMIP4 models, their arithmetic mean (MMM), and HT-ALL. In the HT-ALL simulation, the preindustrial and Middle Holocene refers to the average for 1750–1850 and 6 ka ± 0.05 ka, respectively.

Figure 5 Global annual mean surface air temperature anomaly since 10 ka under the all-forcing (HT-ALL – black) and the single-forcing (HT-ORB – blue; HT-ICE – green; HT-GHG – orange) simulations relative to the last 101 years of the all-forcing simulation. All of the time series are smoothed by a 101-year moving average. The triangle, square, cross, diamond, and circle symbols at 6 ka stand for the Middle Holocene minus preindustrial results simulated by the 14 PMIP4 models and their multimodel arithmetic mean (MMM).

Figure 6 Changes in (a, b) annual, (c, d) DJF, and (e, f) JJA mean surface air temperatures (units in ^∘C) (a, c, e) between the averages for 6 ka ± 0.05 ka and the last 101 years (i.e., 1750–1850) in the HT-ALL simulation and (b, d, f) between the Middle Holocene and preindustrial simulations in the arithmetic mean of the 14 PMIP4 models, as listed in Fig. 5. The zonal mean changes are shown on the right of each panel. The dotted areas in (b), (d), and (f) represent less than 70 % of the PMIP4 models agree with the sign of multimodel mean change.

Figure 7 Zonal averaged changes in (a) annual, (b) DJF, and (c) JJA mean surface air temperatures (units in ^∘C) between the Middle Holocene and preindustrial simulations in the 14 PMIP4 models and their arithmetic mean (MMM), as well as between the averages for 6 ka ± 0.05 ka and the last 101 years (i.e., 1750–1850) in the HT-ALL simulation. Red, blue, gray, black, green, and orange bars represent the changes averaged over the zonal bands of 60–90^∘ N, 30–60^∘ N, 0^∘–30^∘ N, 30^∘ S–0^∘, 60–30^∘ S, and 90–60^∘ S, respectively. The error bar for the multimodel mean change represents plus and minus 1 standard deviation for the 14 PMIP4 models.

Figure 5 further displays the annual GMST evolution since 10 ka in the all-forcing and single-forcing simulations. All time series are the anomalies relative to the last 101 years (1750–1850) of the HT-ALL simulation. In response to the orbital forcing alone, the GMST slightly decreases during the past 10 ka, with an averaged cooling of approximately 0.2 ^∘C in the HT-ORB simulation. The GMST first decreases from 10 to 6.5 ka and then steadily increases until now in the HT-GHG simulation, which closely follows the evolution of the atmospheric concentration of CO₂ and is also modulated by that of CH₄ and N₂O. In the HT-ICE simulation, the GMST increases from 10 to 6 ka in response to the retreat of ice sheets in the Northern Hemisphere and then slightly decreases possibly due to the lag effect of ice sheet melting. Taken together, the weak cooling trend during the Early Holocene in the HT-ALL simulation is mainly driven by both changes in orbital insolation and GHG concentrations, while the retreat of ice sheets plays an opposite role. On the other hand, the strong warming trend since 7.5 ka is mainly induced by increased concentrations of GHGs, which largely results from that of CO₂, with a simultaneously positive contribution from the retreat of ice sheets during 7.5–6 ka and a negative contribution from the orbital insolation throughout the past 7.5 ka. It is worth mentioning that, relative to the average for 1750–1850 in the HT-ALL experiment, the annual GMST anomalies since 10 ka in the three single-forcing experiments do not sum to the anomaly in the all-forcing experiment in our HT-11.5 ka simulations (Fig. 5), which holds true for the single- and all-forcing experiments in the CCSM3, LOVECLIM, and FAMOUS transient simulations (Smith and Gregory, 2012; Liu et al., 2014) and indicates a level of nonlinearity in the climate system.

5 Summary

The Holocene presents a host of opportunities to study the climate evolution over the latest interglacial period, including the mean climate trend on the millennium scale and abrupt changes on the timescales from decades to centuries. The climate model provides a useful tool to investigate the underlying dynamic mechanism for the climate change during this well-reconstructed time period; in particular, the improvements of the model and computational power make it possible to run transient simulations to comprehensively understand the climate evolution in response to time-varying external forcings and feedbacks. Several modeling studies have begun this work, but most of them are performed with climate system models or Earth system models with low spatial resolutions and/or acceleration techniques to reduce computational resources. Therefore, we have carried out a new set of unaccelerated transient experiments covering the whole Holocene since 11.5 ka (HT-11.5 ka), using a comprehensive Earth system model with relatively high spatial resolutions, considering both the full-forcing and single-forcing effects of the main boundary conditions.

The HT-11.5 ka simulations are run with the NCAR CESM1.2.1 with an atmospheric resolution of approximately $\mathrm{1.9}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$ in the horizontal and 26 layers in the vertical and an oceanic resolution of approximately 1^∘ and 60 levels in the horizontal and vertical, respectively. All simulations are spin up from the Early Holocene at 11.5 ka and transiently run from 11.5 ka to the preindustrial period at 1850, with the boundary conditions of Earth's orbital parameters, atmospheric GHGs, and ice sheets and orography that have been provided on the PMIP4 wiki for the transient deglaciation (PMIP Last Deglaciation Working Group, 2016). In total, four transient experiments are run, namely one all-forcing experiment forced by synchronously time-evolving orbital parameters, atmospheric GHGs, and ice sheets and three single-forcing experiments forced by only one of the above forcings with the other two forcings fixed at 11.5 ka. The all-forcing simulation result shows a slight decrease in the annual GMST from 11.5 to 7.5 ka, due to both changes in orbital insolation and GHG concentrations, with an abrupt cooling at approximately 7.5 ka, which is followed by a continuous warming until the preindustrial period, mainly resulting from increased GHG concentrations. Relative to the preindustrial period, both at global and zonal mean scales, the simulated annual and seasonal mean changes in surface air temperature at 6 ka are within the range of values in the 14 PMIP4 model simulations and overall stronger than those from their arithmetic mean. The Late Holocene warming trend simulated by CESM1.2.1 in our all-forcing simulation is consistent with that in the CCSM3, LOVECLIM, and FAMOUS transient simulations, as shown in Liu et al. (2014), but inconsistent with the reconstructed cooling trend in Marcott et al. (2013) and Kaufman et al. (2020a).

The present HT-11.5 ka transient simulations for the whole Holocene, spanning from 11.5 ka to the preindustrial period, are designed to investigate both combined and separated effects of the main external forcing of orbital insolation, atmospheric GHGs, and ice sheets on the climate evolution over the Holocene, which will be the supplement for the PMIP4 transient Core experiments for the last deglaciation from 21 to 9 ka (Ivanovic et al., 2016) and the PMIP4 Tier 3 transient simulations for the Holocene from 6 to 0 ka (Otto-Bliesner et al., 2017). Although this new set of transient Holocene simulations are run with a high-resolution version of the comprehensive Earth system model, the meltwater flux change, solar irradiance and volcanic forcing, land use and vegetation changes, dust and aerosol effects, and stratospheric chemistry and dynamics are not fully considered here, which will have some effects on the results. As a first step, this study provides the detailed model description, experimental design, boundary conditions, and some preliminary results for the new set of HT-11.5 ka transient simulations. Further analyses on the simulation results, including the mean climate evolution and abrupt climate changes over the Holocene both at global and regional scales, the underlying dynamic mechanisms, and comparisons with other transient simulations covering the Holocene and multi-proxy-based reconstructions, will be carried out by a series of follow-up studies.