A Comparison of the CMIP6 midHolocene and lig127k Simulations in CESM2

Results are presented and compared for the Community Earth System Model version 2 (CESM2) simulations of the middle Holocene (MH, 6 ka) and Last Interglacial (LIG, 127 ka). These simulations are designated as Tier 1 experiments (midHolocene and lig127k) for the Coupled Model Intercomparison Project phase 6 (CMIP6) and the Paleoclimate Modeling Intercomparison Project phase 4 (PMIP4). They use the low‐top, standard 1° version of CESM2 contributing to CMIP6 DECK, historical, and future projection simulations, and to other modeling intercomparison projects. The midHolocene and lig127k provide the opportunity to examine the responses in CESM2 to the orbitally induced changes in the seasonal and latitudinal distribution of insolation. The insolation anomalies result in summer warming over the Northern Hemisphere continents, reduced Arctic summer minimum sea ice, and increased areal extent of the North African monsoon. The Arctic remains warm throughout the year. These changes are greater in the lig127k than midHolocene simulation. Other notable changes are reduction of the Niño3.4 variability and Drake Passage transport and a small increase in the Atlantic Meridional Overturning Circulation from the piControl to midHolocene to lig127k simulation. Comparisons to paleo‐data and to simulations from previous model versions are discussed. Possible reasons for mismatches with the paleo‐observations are proposed, including missing processes in CESM2, simplifications in the CMIP6 protocols for these experiments, and dating and calibration uncertainties in the data reconstructions.


Introduction
The modeling of paleoclimate in the Coupled Model Intercomparison Project (CMIP) and the Paleoclimate Modeling Intercomparison Project (PMIP) has long been used to understand and explain past environmental change (Braconnot et al., 2007(Braconnot et al., , 2012Joussaume et al., 1999) and as out-of-sample tests of models being used for projections of future climate (Folland et al., 1990;Jansen et al., 2007;Masson-Delmotte et al., 2013;McAvaney et al., 2001).
The middle Holocene (MH, 6 ka) has been a target of PMIP since its first phase using atmosphere only models because of evidence from data of warmer Northern Hemisphere (NH) summers and enhanced NH monsoons (COHMAP Members, 1988;Harrison et al., 1998) in response to enhanced summer insolation associated with the shift of perihelion to near the boreal autumn equinox as compared to near the boreal winter solstice at present (Berger & Loutre, 1991;Otto-Bliesner et al., 2017). A wide variety of archives have been analyzed since then to reconstruct additional climate variables, including lake levels (Coe & Harrison, 2002), growing degree days (Bartlein et al., 2011), biomass burning (Power et al., 2008), ENSO variability (Emile-Geay et al., 2016) and sea ice (Stein, Fahl, Schade, et al., 2017). The midHolocene is considered Tier 1 within CMIP6 and one of two paleoclimate experiments, with the Last Glacial Maximum (LGM), that serves as an entry card to PMIP4 along with a CMIP6 DECK preindustrial (PI) experiment as a paired control (Eyring et al., 2016;Kageyama et al., 2018).
For the first time, the Last Interglacial (LIG, 127 ka) has been included as a Tier 1 CMIP6-PMIP4 experiment to allow the assessment of model simulations for that warm period in the IPCC AR6 (Intergovernmental Panel on Climate Change Sixth Assessment Report; Otto-Bliesner et al., 2017). Previously, LIG simulations were made with older and/or lower-resolution versions of the models than were used for future projections (Lunt et al., 2013). Differences in these LIG simulations could also have arisen because of differences in the experimental designs: the orbital forcing and GHG concentrations varied among the simulations. Several reasons were given at the 2015 Cambridge PAGES QUIGS workshop for identifying a 127 ka time slice for the CMIP6-PMIP4 LIG experiment: Northern Hemisphere seasonal insolation anomalies are large with perihelion near boreal summer solstice, little or no remnants of the North American and Eurasian ice sheets, and sufficient time elapsed to minimize the imprint of the earlier deglaciation and Heinrich 11 (H11) meltwater event given dating uncertainties (Marino et al., 2015;Otto-Bliesner et al., 2017). For CMIP6-PMIP4, a common protocol was established for participating modeling groups to run with the same model and at the same resolution as their DECK simulations if participating in CMIP6 or their preindustrial control simulations if only participating in PMIP4 (Otto-Bliesner et al., 2017).
The family of CESM models (CSM1, CCSM3, CCSM4) have simulated the MH, starting with atmosphere-only simulations with fixed present-day sea surface temperatures (CSM1, Otto-Bliesner, 1999), low-resolution coupled atmosphere-ocean simulations (CCSM3, Otto-Bliesner et al., 2006), and coupled atmosphere-ocean simulations in CMIP5 with CCSM4 run at the same resolution (FV1x1) as the preindustrial, historical and future projections. CCSM4 experiments were also completed for the Last Glacial Maximum (LGM, Brady et al., 2013), Last Millennium (LM, Landrum et al., 2013), Last Interglacial, and mid-Pliocene Warm Period (mPWP, Rosenbloom et al., 2013) to understand the model response of the climate system to different climate forcings for documented states very different from the present and historical.
In this paper, we provide a short summary of the model components in Section 2 (for more details see Danabasoglu et al., 2020) and the experimental design in Section 3. In Section 4 we provide an analysis of the responses of CESM2 to the orbitally induced changes in insolation on the atmosphere, ocean, sea ice, and land. We compare and contrast the responses between the midHolocene and lig127k with their differences in the phasing and magnitudes of the seasonal and latitudinal distribution of insolation anomalies in terms of mean climate and climate variability. Comparison to previous modeling results from the family of CESM models and proxy reconstructions is included in Section 5.

Model Description
The CESM2 is a coupled model consisting of components for the atmosphere (CAM6, WACCM6), ocean (POP2), sea ice (CICE5), land (CLM5), and land ice (CISM2). We briefly summarize the component models relevant to our midHolocene and lig127k simulations. More details on each of these components, their coupling, and the science improvements can be found in Danabasoglu et al. (2020).
The low-top atmosphere model, the Community Atmosphere Model version 6 (CAM6) is used for the CESM2 midHolocene and lig127k simulations. CAM6 has a nominal 1°(1.25°in longitude and 0.9°in latitude) horizontal resolution. There are 32 vertical levels with the top at 2.26 hPa (about 40 km): 7 levels below 850 hPa, 10 levels between 850 and 200 hPa, and 15 levels above 200 hPa. It includes the four-mode Modal Aerosol Model (MAM4; Liu et al., 2016). A two-moment cloud microphysics scheme (MG2, Gettelman & Morrison, 2015) explicitly models the cloud-aerosols interactions, allowing for the inclusion of the indirect effects of aerosols. The Cloud Layers Unified by Binormals (CLUBB) is now used for representation of shallow convection, boundary layer, and grid-scale condensation (Bogenschutz et al., 2013). The orographic drag parameterizations have also been updated.
The ocean model is the Parallel Ocean Program version 2 (POP2), the same as in CCSM4 but with improvements in some physical parameterizations, including an estuary balance model for better handling of runoff from the land model and increased mesoscale diffusivities at depth to improve the representation of passive tracers. The ocean model also simulates ocean biogeochemistry (MARBL), including multiple nutrient co-limitation (N, P, Si, Fe) and three explicit phytoplankton (diatoms, diazotrophs, pico/nano phytoplankton), one implicit group (calcifiers), and one zooplankton group. The horizontal resolution is a nominal 1°w ith a uniform resolution of 1.125°in the zonal direction. The North Pole is displaced into Greenland. The meridional resolution is variable: 0.27°at the equator; increasing monotonically to 0.53°at 32°S, then constant farther south; 0.38°in the northwestern Atlantic Ocean and~0.64°in the northwestern Pacific Ocean. There are 60 vertical levels with a maximum depth of 5,500 m. The upper 16 vertical levels have a uniform thickness of 10 m. Below~160 m depth, the level thicknesses increase monotonically to~3,500 m depth. The deepest 2000 m have a nearly uniform thickness of~250 m extending to the bottom at~5,500 m. The sea ice model is CICE version 5 (Hunke et al., 2015), now incorporating mushy-layer thermodynamics and an updated melt pond scheme. It uses the same horizontal grid as the ocean model.
The land model is the Community Land Model version 5 (CLM5), with extensive improvements to the model's hydrological and ecological parameterizations. Its horizontal grid is the same as the atmospheric model. Atmospheric dust is mobilized from the land. Factors for determining soil erodibility and dust emission include the wind friction speed, vegetation cover, and soil moisture (Mahowald et al., 2006;Zender et al., 2003). All other aerosol emissions in our lig127k and midHolocene simulations are prescribed as in the piControl simulation as given by the CMIP6 datasets. A new river model (Model for Scale Adaptive River Transport, MOSART) is included, which simulates streamflow, time-varying channel velocities, channel water depth, and channel surface water variations.
The land-ice model is the Community Ice Sheet Model version 2.1 (CISM2.1), based on the Glimmer-CISM model, and includes a three-dimensional thermodynamical model and high-order velocity solver. It currently runs for Greenland on a 4-km grid. For the midHolocene and lig127k simulations, CISM2.1 is set to NO-EVOLVE, i.e. a fixed Greenland ice sheet. The Greenland ice sheet Surface Mass Balance (SMB) and surface temperature are downscaled from the 10-virtual elevation classes calculated in the land model to the 4-km CISM grid, allowing for future use in ice sheet-only simulations (Lipscomb et al., 2013).

Experimental Design and Methods
The CESM2 (CAM6) midHolocene and lig127k simulations use the CESM2.1.0 released version. We adopt the forcings and boundary conditions (Table 1) following the protocols of the Tier 1 experiments (midHolocene and lig127k) for CMIP6 and PMIP4, as defined in Otto-Bliesner et al. (2017), with two exceptions: vegetation and dust. Vegetation is prescribed as potential vegetation rather than preindustrial vegetation. The CESM2 interactive dust parameterization is modified to account for the influence on dust emissions exerted by the climate-induced changes in vegetation during the middle Holocene and Last Interglacial not included in our simulations with prescribed vegetation. These simulations are compared to the CESM2 CMIP6 piControl (see Eyring et al., 2016, andDanabasoglu et al., 2020, for details on the forcings, boundary conditions, and piControl simulation).
The Earth's orbital parameters (eccentricity, longitude of perihelion, and obliquity) are prescribed following Berger and Loutre (1991). These parameters affect the seasonal and latitudinal distribution and magnitude of solar energy received at the top of the atmosphere and, in the case of obliquity, the annual mean insolation at any given latitude. At 6 ka, the eccentricity is similar to PI but perihelion occurs near the boreal autumn equinox as compared to near the boreal winter solstice at PI. The orbit at 127 ka is characterized by larger eccentricity than at PI, and in addition perihelion occurs close to the boreal summer solstice. These orbital configurations result in large positive NH summer insolation anomalies at 6 ka and 127 ka as compared to PI (Figure 1). The anomalies are considerably greater at 127 ka than at 6 ka, associated with greater eccentricity and boreal summer solstice at perihelion at 127 ka.
Climate models generally adopt the modern calendar when outputting monthly averages. That is, the month lengths are fixed at their modern definitions. The altered orbit at 6 ka and 127 ka, as compared to PI, affects the Earth's transit speed along different parts of its orbit. So, for example, the number of days between the NH summer solstice and the NH fall equinox is about 4 days and 7 days shorter than at PI for 6 ka and 127 ka, respectively (see Figure 2 in Otto-Bliesner et al., 2017). Conversely, the number of days between the NH winter solstice and the NH spring equinox get progressively longer from PI to 6 ka to 127 ka. The net effect is that NH winters (December, January, February, DJF) were longer and summers (June, July, August, JJA) shorter from an insolation perspective at 127 ka and 6 ka than at PI. Bartlein and Shafer (2019) illustrate that the impacts on temperature and precipitation anomaly patterns related to changes in the length of months is greatest from September to December for 6 ka, and all but January to May for 127 ka. The PaleoCalAdjust v1.0 of Bartlein and Shafer (2019) accounts for the impact of these changes in the length of months at 6 ka and 127 ka. In the analyses of CESM2 results, we use PaleoCalAdjust v1.0 to convert the outputted monthly averages to 'celestial' calendar months.
Vegetation cover is prescribed in the midHolocene, lig127k, and piControl simulations, though the leaf area index and vegetation height are prognostic and can be affected by the changed climate. The land cover in the piControl simulation merges present-day satellite land cover with the LUH2 data for primary and secondary forest and nonforest land units, five crop groups, managed pasture, rangeland, and urban areas at 1850 CE (Lawrence et al., 2019). To remove anthropogenic influences on the land surface for the midHolocene and lig127k simulations, a new "potential vegetation" distribution was derived from natural vegetation land units, such that historical crop and urban areas of CESM2 are replaced with suitable natural vegetation types reflecting the distribution of vegetation in the absence of human activity under preindustrial environmental conditions. Potential vegetation distributions were generated using the LUMIP Land Use Harmonization version 2 (LUH2) data (Lawrence et al., 2016) starting with the land use distribution for the year of 850 CE. Any agriculture present in that year was replaced by the nearest neighboring natural land unit, producing the potential vegetation LUH2 land unit distribution. The largest changes in agricultural land use compared to PI are in India, Europe and China (Figure 2 and Figure S1). The urban land cover was set to zero for all land points in keeping with the lack of significant built infrastructure at 6 ka and 127 ka.
In CESM2, soil erodibility maps provide a scaling factor on dust emission fluxes. It is a spatially varying parameter summarizing the differences in susceptibility to erosion related to, e.g. soil textures and geomorphology (Mahowald et al., 2006;Zender et al., 2003). The CESM2 midHolocene and lig127k simulations do not include changes in vegetation associated with the simulated changes in precipitation and temperature during these time intervals (i.e., we use prescribed near-modern, pre-anthropogenic vegetation). To account for the influence on dust emissions exerted by changes in vegetation , we performed offline simulations with the coupled biogeography and biogeochemistry equilibrium model (BIOME4) (Kaplan et al., 2003) to determine midHolocene vegetation, and use these results to determine the location of potential dust source areas, as well as the area of a grid box available for dust generation, assuming that the nonvegetated fractions of arid areas can lead to dust emissions (Mahowald et al., 2006). We then incorporated this as an additional information layer of the soil erodibility map, at the same time removing the dependence of dust mobilization on TLAI (Total Leaf Area Index) from the main code, which would use the preindustrial vegetation instead, as described in Albani et al. (2015). A second step involves
further refining, at the continental level, the soil erodibility map. A set of scale factors for the soil erodibility map is defined for different macro areas (i.e., continents), based on achieving the best fit of dust deposition from preliminary midHolocene simulations with CESM2 to an observational benchmark, i.e., dust mass accumulation rates from paleodust archives (Albani et al., 2015). The midHolocene soil erodibility map was also used in the lig127k simulation. Note that a change in the coarse mode aerosol definition (from a geometric sigma of 1.6 to 1.2) between CAM5 (Albani et al., 2014) and the CAM6 used here, results in a more negative dust direct radiative effect, no longer matching available observations. Further analyses investigating this aspect in more detail would be useful.
Records from Antarctica and Greenland ice cores indicate 6 ka and 127 ka concentrations of the well-mixed GHGs: CO 2 , CH 4 , and N 2 O comparable to, albeit somewhat lower than PI (Table 1). All other forcings and boundary conditions-solar constant, ozone, volcanic aerosols, topography and bathymetry, and ice sheets -in the midHolocene and lig127k simulations remain the same as in the piControl simulation.

Results
4.1. Atmosphere 4.1.1. Surface Air Temperature Both 6 ka and 127 ka show large NH positive insolation anomalies during boreal summer as compared to PI ( Figure 1). The positive NH summer insolation anomalies at 6 ka force warming greater than 2°C over the NH continents in June-July-August (JJA) in the midHolocene simulation, with some regions at mid-latitudes over Eurasia experiencing surface air temperature increases of more than 4°C (Figure 3(f)). JJA warming over the NH continents in the lig127k simulation, associated with the much larger positive NH summer insolation anomalies, exhibit surface temperature increases approximately doubled as compared to 6 ka, to more than 4°C to over 8°C (Figure 3(e)). In the lig127k simulation, the response of surface air temperatures is more than just a NH response. The nature of the orbital parameters, with the positive insolation anomalies from April to September extending into the Southern Hemisphere (SH), results in considerable warming over South America, Australia, and southern Africa as well as in the Atlantic and Indian Ocean sectors of the Southern Ocean in the lig127k simulation. These responses are also reflected in JJA global and NH extratropical land averages (Table 2).
At 6 ka and 127 ka, negative insolation anomalies occur during December-January-February (DJF) as compared to PI (Figure 1). The continental and much of the oceanic regions in both hemispheres cool in response to the insolation forcing, with greater cooling in the lig127k simulation than the midHolocene simulation ( Figure 3(c), 3(d)). The exceptions are the Arctic, high-latitude North Atlantic, and Atlantic and Indian Ocean sectors of the Southern Ocean, which remain warmer than PI. Positive feedbacks with the cryosphere (sea ice and snow cover) and ocean provide the memory that allows simulated high-latitude warming, in DJF as well as annually (Figure 3(a), 3(b)).
While the largest anomalies in surface air temperature appear seasonally, there is a consistent pattern of weak annual mean warming at high latitudes and cooling across the tropics and subtropics notable in both the lig127k and midHolocene. This pattern is suggestive of the obliquity-driven changes to the annual mean insolation which brings about a weak increase of up to 4 W m −2 at high latitudes in both, and a decrease of about 1 W m −2 in midHolocene (about 0.5 W m −2 in lig127k) in the tropics. Interestingly, the increase in obliquity is larger in midHolocene than in lig127k, yielding slightly greater annual mean insolation anomalies, yet the annual mean surface temperature anomalies are weaker in the midHolocene in comparison, especially over the high northern latitudes. Also, the midHolocene has a greater decrease in greenhouse gas forcing than lig127k. This suggests other processes remain important in driving the resulting changes to the annual mean surface temperature such as ocean and cryosphere memory. Arctic surface ocean gains heat in response to sea ice loss in summer due to the greater open ocean, and due to higher ocean heat capacity (memory), the warmer ocean surface temperatures are retained into the fall and winter seasons, resulting in thinner winter sea ice. The thinner winter sea ice is consistent with a weaker insulating effect of sea ice, allowing increased air-sea heat exchange and hence increased winter air temperatures (Serreze & Barry, 2011).

Precipitation
The numerous improvements incorporated into the atmospheric component of CESM2 have had a positive impact on the simulation of precipitation, including reductions in the biases of the Intertropical Convergence Zone (ITCZ), over the tropical oceans, and the Amazon, Australia, and western U.S. regions (Danabasoglu et al., 2020). The seasonal and latitudinal structure of insolation anomalies at 6 ka and 127 ka promote enhanced NH monsoons and reduced SH monsoons in the midHolocene and lig127k simulations ( Figure 4), with little change in annual mean, global precipitation (Table 2). Over North Africa, JJA precipitation increases and extends northward into the Sahara as well as southward in central Africa as compared to the piControl, with these changes greater in the lig127k than the midHolocene. The Asian monsoon is enhanced in both interglacial simulations as compared to the piControl, with again greater increases in precipitation in the lig127k simulation than the midHolocene simulation. The lig127k simulation shows a band of enhanced JJA precipitation over the islands of Malaysia, Indonesia, and Papua New Guinea, also present in the midHolocene simulation but not statistically significant. Overall tropical (20°S-20°N) land, JJA precipitation increases by 8% in the midHolocene and 27% in the lig127k as compared to the piControl.
The SH monsoons and associated austral summer precipitation over southern Africa, Australia, and South America are reduced in both interglacial simulations as compared to the piControl, with greater decreases in precipitation in the lig127k than the midHolocene.
Over the oceans, increased precipitation is associated with the winter (DJF) mid-latitude storm tracks over the North Pacific and North Atlantic Oceans compared to the piControl, most prominently in the lig127k but also in the midHolocene, and decreased precipitation associated with the winter (JJA) mid-latitude storm tracks over the South Pacific, South Atlantic, and southern Indian Oceans compared to the piControl ( Figure 4). Annually, the ITCZ over the tropical oceans shifts northward with the interglacial forcings.

Dust Aerosol Optical Depths
As explained in Section 3, the midHolocene and lig127k simulations share the same soil erodibility map for dust (Otto-Bliesner et al., 2017). In this section we focus on the Holocene, for which observational constraints on the dust cycle are available (Albani et al., 2015). Simulated dust deposition (wet + dry) for the midHolocene captures fairly well the observational constraints represented by dust mass accumulation rates from paleodust archives integrated for the 6 ka time slice (5-7 ka) ( Figure 5(a); Albani et al., 2015).
Dust Aerosol Optical Depth (AOD) anomalies between midHolocene and piControl show a significant reduction in the North African dust plume over the continent that extends into the Atlantic Ocean ( Figure 5(b)).
In terms of global budgets, total dust deposition (or emissions) amounts to 2,200 Tg yr −1 for the piControl and 1,824 Tg yr −1 for the midHolocene cases, respectively. Compared with previous estimates using an older version of the CESM model (Albani et al., 2016), we have lower emissions in both climates and a more accentuated (17% vs 9%) decrease in dust emissions between midHolocene and piControl. Figure 5. Dust in the midHolocene and piControl simulations. a) Comparison of simulated total (wet + dry) dust deposition (contours) against observational constraints from paleodust archives (circles) for the midHolocene (Albani et al., 2015), using the same color scale. b) Difference in simulated dust aerosol optical depth between the midHolocene and piControl cases. c) Fraction of aerosol optical depth due to dust with respect to all-aerosols, in the piControl simulation.
At the same time, in our simulations with CESM2, dust AOD was reduced from 0.030 in piControl to 0.022 in the midHolocene case, accounting for 19% and 25% of the total AOD, respectively. In particular, we note how dust dominates AOD over the entire arid belt extending from North Africa to East Asia ( Figure 5(c)). In the absence of more specific diagnostic variables for dust, AOD gives a first-order idea of the potential importance of dust changes in mediating Direct Radiative Effects (DRE). We should note that dust interacts not only with solar radiation, but also with thermal radiation.

Sea Ice and Ice Sheets 4.2.1. Northern and Southern Hemisphere Sea Ice
The annual cycles of NH sea ice area show reduced areal coverage during boreal summer, but little change during boreal winter, in both the midHolocene and lig127k simulations as compared to the piControl ( Figure 6(a)). The September minimum is 33% less in the midHolocene and 75% less in the lig127k. NH sea ice area remains below that of the piControl from July to November in both the midHolocene and lig127k simulations. It should be noted that the piControl has a thin ice bias in CESM2(CAM6) from more melting during the boreal summer months resulting in a low summer ice area as compared to satellite observations (Danabasoglu et al., 2020).
The location of the climatological mean Arctic sea-ice edge in August-September (AS) denoted by the extent (fractional ice area >15%) retreats from the piControl to the midHolocene to the lig127k ( Figure S2). This corresponds to the markedly thinner sea ice across the NH polar region in both the midHolocene and lig127k in comparison to the piControl consistent with the enhanced seasonal cycle of insolation in the NH. Negative AS sea ice thickness anomalies in the central Arctic exceed~0.5 m in the midHolocene and~1 m in the lig127k with greater thinning found towards Greenland and the North American coast in both simulations.
The midHolocene annual cycle of SH sea ice area shows little change from the piControl (Figure 6(b)). The lig127k has a dampening of seasonal cycle of SH sea ice area with somewhat more area in austral summer and less area in austral winter than the piControl. The climatological mean February-March (FM) sea ice extent (fractional ice area >15%) is also roughly similar in the SH in all simulations, except for a somewhat greater equatorward extent in the Pacific sector of the Southern Ocean in the lig127k in FM ( Figure S2). The simulated responses of climatological mean FM sea ice thickness are rather weak in the SH in comparison to the response found in the NH however. The lig127k has a somewhat thicker sea ice distribution in the Pacific sector of the Southern Ocean, with thinning mostly elsewhere, except east of the Antarctic Peninsula. The midHolocene shows a greater thinning response in the Weddell Sea in comparison to the response found in the lig127k simulation.

Greenland Ice Sheet
The Greenland ice sheet Surface Mass Balance (SMB) for the 1961-1990 mean of the CESM2(CAM6) historical simulation agrees well with the regional atmospheric climate model RACMO2 simulation, albeit with a positive SMB bias in southern Greenland associated with excessive snowfall (Danabasoglu et al., 2020;van Kampenhout et al., 2020). Associated with warmer climates of the LIG and MH, enhanced precipitation with storm tracks results in more positive SMB in southern and southeastern Greenland, particularly so in the lig127k, as compared to the piControl (Figure 7). Areas of negative SMB (more ablation than accumulation annually) expand along edges of the Greenland ice sheet, particularly in western Greenland and also northwest Greenland, in the midHolocene and lig127k simulations. At high elevations in central Greenland, the small, positive SMB in the piControl remains unchanged in the midHolocene and lig127k. Annual temperatures are warmer in the midHolocene and lig127k simulations but remain below freezing year-round. The extent of negative SMB is especially pronounced in the lig127k simulation, suggesting the Greenland ice sheet would lose mass and contribute to the average global mean sea level rise if CESM2 was coupled to CISM2 dynamically.

Ocean 4.3.1. Barotropic Streamfunction
The interglacial simulations show changes in the annual mean barotropic transport streamfunction. Globally the largest differences in the barotropic ocean transport compared to the piControl are found in the Southern Ocean and the North Atlantic (Figure 8). In the Southern Ocean, the Antarctic Circumpolar Current (ACC) weakens in both interglacials, but to a greater degree in lig127k as compared to midHolocene. The transport through the Drake Passage is reduced by~22.0 Sv in the lig127k and~13.6 Sv in the midHolocene (Table 2). In addition, the dipole pattern in the annual mean barotropic transport anomalies in the Pacific sector of the Southern Ocean suggests a modest strengthening of the barotropic transport within the South Pacific subtropical gyre and the ACC at the equatorward edge in comparison to the piControl. This dipole pattern is similar but weaker in midHolocene as compared to lig127k. The response of ACC is consistent with a small shift northward of the southern westerlies in the Pacific sector, resulting in a weakening of the zonal wind stress at the latitude of the Drake Passage (not shown).
In the North Atlantic (Figure 8), the annual mean barotropic streamfunction associated with the circulation within the subtropical gyre system strengthens in both orbitally driven cases relative to the piControl, with a greater strengthening found in the lig127k. The strengthening occurs north of the Florida Straits, and equatorward of~38°N, with a decrease in the barotropic streamfunction found at the northern edge of the subtropical gyre. In the subpolar North Atlantic, both orbitally driven cases show an increase in the cyclonic circulation of the subpolar gyre, with the lig127k exhibiting a greater anomaly. These anomalies in the barotropic streamfunction in the North Atlantic are consistent with a strengthening of the westerlies over the N. Atlantic, which increases the magnitude of the wind stress curl in the vicinity of the subpolar gyre and western boundary of the subtropical gyre.

Meridional Overturning Circulation and Deep Vertical Mixing
The Atlantic meridional overturning circulation (AMOC) appears to be weakly sensitive to changes in the combined forcing from the seasonal and latitudinal redistribution of insolation and the small decrease in atmospheric greenhouse gases ( Figure 9, Table 2). A weak increase in the strength of the meridional overturning circulation cell associated with the formation of the model North Atlantic Deep Water (NADW) is found north of~30°N in the upper thermocline and is reflected in the weak increase in the maximum AMOC (Table 2). In comparison to the piControl simulation, both paleoclimate simulations exhibit a weak increase in the northward circulation of southern-sourced bottom water in the deepest levels of the Atlantic basin (Table 2, Figure 9). The piControl shows a somewhat weaker circulation of Antarctic Bottom Water (AABW) than reported for the 1990-2014 period in the historical simulation in Danabasoglu et al. (2020).
The deepest late winter mixing takes place in the North Atlantic in association with the formation of model NADW. In association with the increase in AMOC is a change in the late winter deep vertical mixing in the subpolar gyre region. In the midHolocene and lig127k simulations, late winter deep vertical mixing decreases somewhat (by~150 m) in the mouth of the Labrador Sea, where deep later winter mixing reaches depths of 1,600 m in piControl, and intensifies greatly by~500 m deeper off the east coast of Greenland and directly south of Greenland within the subpolar gyre ( Figure S3). In the Nordic Seas south of Spitzbergen ( Figure S3) the late winter mixing decreases in lig127k and midHolocene. The pattern of mixed layer depth anomalies relative to the piControl is similar in both but more amplified in the lig127k compared to the midHolocene.
In the SH, the late winter vertical mixed layer depths associated with the formation of relatively fresh and cold Antarctic Intermediate Water, extending west and east of the southernmost tip of South America, increase in the midHolocene and lig127k, with larger increases found in the lig127k ( Figure S3). A modest reduction in the late winter mixed layer depth is noted however, over the band of relatively well mixed upper water associated with the lighter type of intermediate mode water, known as the Subantarctic Mode Water, which extends west from the western Pacific into the Indian Ocean. Improvements in ENSO characteristics (period, amplitude, seasonal cycle) in CCSM4 as compared to CCSM3 are retained in CESM2 (Danabasoglu et al., 2020). Comparable to observations, the Niño3.4 SST index in the piControl simulation has maximum power at periods of 3-7 years, though with the Niño3.4 SST standard deviation overestimated by~30% (Figure 10(f)). The annual cycle of the monthly standard deviation peaks in November-December-January with a minimum in June-July (Figure 10(e)). The Niño3.4 index has reduced power going from the piControl to midHolocene to lig127k, with a 20% reduction in the Niño3.4 standard deviation in the lig127k as compared to the piControl ( Table 2). The lig127k also exhibits a much flatter annual cycle of monthly standard deviations with the minimum shifted to August-September (Figure 10(a)).
The seasonal cycle of SST along the equator in the Pacific shows a strong sensitivity to the orbitally-driven changes in the seasonal cycle of insolation in the tropics ( Figure S4). The anomalies in the seasonal cycle are largest in the eastern to central Pacific where the equatorial thermocline is shallow and coupled ocean-atmosphere dynamics play an important role in the annual cycle of equatorial SST, and more muted west of the dateline where the thermocline is deep. The anomalies in the annual SST cycle with respect to the piControl are the same order of magnitude as the amplitude of the annual cycle. The anomalies in the seasonal cycle of SST along the equator in the central and eastern Pacific with respect to the piControl, are generally in phase with the tropical insolation anomalies as shown in Figure 1. Relative to the seasonal cycle in the piControl in both paleoclimate simulations, central to eastern Pacific seasonal warm anomalies are largest when the insolation anomalies are most positive (i.e., JJA for lig127k and JAS for midHolocene) and the greatest cooling relative to the piControl's seasonal cycle occurs when insolation anomalies are most negative (i.e., JFM for the midHolocene and DJF for lig127k).
The eastern Pacific seasonal warming in lig127k relative to the annual mean ( Figure S4) is shifted later by a few months in phasing and warms with greater amplitude than the piControl. The large negative insolation anomalies during the late boreal winter/early spring in lig127k tend to mute the start of the warming in the eastern equatorial Pacific, then the seasonal equatorial warming is enhanced by the large positive insolation anomalies starting in late spring and peaking in boreal summer. This tendency shifts the peak seasonal warming by a couple of months. The equatorial insolation anomalies are weaker and shifted later in midHolocene than in lig127k. The negative insolation anomaly in midHolocene peaks in early boreal spring, muting the spring warming, while the positive insolation anomaly starts later and extends into the boreal fall, muting the cooling tendency by strong equatorial upwelling. The lig127k shows greater cooling from November to March in comparison to the piControl and the midHolocene when the large negative lig127k insolation anomalies tend to enhance the tendency of seasonal cooling by equatorial upwelling.
In the central and eastern Pacific, the amplitude of the seasonal cycle of SST in the lig127k is larger than in the piControl, whereas the seasonal SST amplitude is noticeably weaker in the midHolocene, particularly along the eastern boundary of South America, and in the amplitude of seasonal cooling in boreal fall. By itself, the opposite responses of the magnitude of the seasonal cycle do not appear to explain the weaker ENSO amplitude in both paleoclimate simulations as compared to the piControl. Rather, the colder annual mean SST compared to the piControl, in addition to colder tropical SSTs in the boreal cold season, may act to moderate the amplitude of the ENSO events which tend to peak at the end of the year.

Model-Model and Model-Data Comparisons
To inform our understanding of the climate responses to the orbital forcing at 6 ka and 127 ka, comparing the CESM2 midHolocene and lig127k results to the latest proxy reconstructions is instructive. Also, model differences are considerable between the older versions of the models (CCSM3 and CCSM4) and CESM2, Figure 10. The monthly standard deviations (°C) and power spectra (°C 2 /cycles mo −1 ) of Niño3.4 SST in the lig127k, midHolocene, and piControl simulations. The mean annual cycle is removed prior to every calculation by subtracting the long-term monthly means. The best-fit first-order Markov red noise spectrum (red curve) and its 95% (blue curve) and 99% (green curve) confidence bounds are shown for the power spectra. The observed power spectra for 1920-2018 from the HadiSST dataset is shown in light gray in panel f).
which are associated with significant improvements made in the physics of the component models. Simulations with these three versions were conducted using the nominal 1°resolution displaced pole grid for sea ice and ocean components and either the FV0.9x1.25 horizontal grid in the atmospheric component (CESM2 and CCSM4) or the T42 (~2.8°resolution) grid in the CCSM3. CCSM4 was used in CMIP5 and PMIP3, whereas CCSM3 was used in PMIP2. Though there were small variations between the different PMIP and CMIP experimental protocols, the differences are not large enough to explain major differences between the three model versions for a given experiment.

Temperature
The JJA zonal-average surface temperature anomalies are compared for the three versions of the model (Figure 11). For the midHolocene, the JJA zonal-average surface temperature anomalies from~40°S to~20°N are negative and of similar magnitude across the three model versions. All have maximum warming in the NH mid-latitudes of approximately 1°C. Differences in the surface temperature anomalies at polar latitudes in both hemispheres are related to simulation of sea ice thickness and extent in the respective preindustrial simulations, as discussed in Section 5.4.
Notable for the lig127k simulations are that the JJA zonal-average surface temperature anomalies are positive at all latitudes except~25°N in CCSM4 and CESM2. These positive anomalies peak also similarly at 3-4°C in the NH mid-latitudes. As for the MH, differences in the surface temperature anomalies at polar latitudes in both hemispheres are related to simulation of sea ice thickness and extent in the respective preindustrial simulations.
The CESM2 midHolocene annual surface air temperature anomalies can be compared to the 6 ka paleo-temperature records from a new global compilation of quality-controlled, published, temperature-sensitive proxy records for the Holocene (Figure 12(a)) (Kaufman et al., 2020). Several regional features are in general agreement between the proxy records and the midHolocene simulation. These include east-west patterns of positive/negative 6 ka surface air temperature anomalies as compared to PI in eastern/western U.S., eastern China/Tibet, and western/eastern North Atlantic. CESM2 and the proxy reconstruction also indicate warmer surface temperatures over South Africa and Brazil at 6 ka as compared to PI. CESM2 simulates no change in annual temperatures at the locations of the East Antarctica ice cores in agreement with the 6 ka reconstruction, but also no change at the locations of the Greenland summit ice cores where the reconstruction suggests modest positive surface temperature anomalies at 6 ka. Over Europe and Canada, the temperature proxies indicate substantial warming nearby to substantial cooling making a credible model-data comparison difficult. Similarly, over North Africa, the only one record available does not allow an assessment on whether the simulated annual cooling is reasonable. The latter is consistent with the enhanced monsoon precipitation and cloud cover simulated in the midHolocene simulation.
Much less data are available to compare to the lig127k simulation of annual surface temperatures (see Otto-Bliesner et al., 2020 for discussion and tables of reconstructions). Similar east-west contrasts of simulated temperature anomalies over the U. S, China, and the North Atlantic seen in the midHolocene simulation are also present in the lig127k simulation (Figure 12(b)). As well, the lig127k has positive surface temperature anomalies over South Africa and Brazil that cannot yet be verified from proxy records. Over Europe, the temperature proxies show a generally consistent picture of annual warming in the lig127k simulation. The polar ice cores and the lig127k simulation suggest higher surface temperature at 127 ka than PI, though CESM2 underestimates this warming. Note that the protocol for the lig127k simulation adopts modern ice sheet extents and elevations for the Greenland and Antarctic ice sheets. The proxy-reconstructed positive surface temperature anomalies in the coastal upwelling regions off California, as expected are not seen in the lig127k simulation. Models such as CESM2 do not have high enough resolution to realistically simulate coastal upwelling.
The 6 ka proxy reconstruction indicates only modest changes in surface temperatures during the boreal summer over North America, Nordic Seas, Europe, western Eurasia and far northern Africa in regions where the midHolocene simulation predicts warmer surface temperatures (Figure 13(a), 13(c)). It should be noted that the lack of spatial coherence in the European temperature reconstructions complicates comparisons with the CESM2 outputs. Both time periods show patterns of greater summer warming over the NH continents than the North Atlantic, but the lig127k simulation has significantly greater positive surface temperature anomalies. Notable, is the good agreement of the simulated and proxy records over Baffin Island and the central North Atlantic for 127 ka (Figure 13(b)). On the other hand, proxy records south of Greenland and in the Nordic Seas indicate colder SSTs at 127 ka than PI (Figure 13(d)). This model-data mismatch could be associated with meltwater from potentially remnant ice sheets over Canada and Scandinavia (Barlow et al., 2018), which the lig127k simulation did not incorporate.

Precipitation
All three model versions (CCSM3, CCSM4, CESM2) successfully capture the large-scale observed modern features of precipitation over North Africa, with annual precipitation peaking at equatorial latitudes and

10.1029/2020PA003957
Paleoceanography and Paleoclimatology becoming negligible north of 20°N ( Figure S5). Although the maximum precipitation in the equatorial regions is reduced in CESM2 from that simulated by CCSM3 and CCSM4, all three versions still overestimate annual precipitation south of 10°N. The summer monsoon over North Africa intensifies and shifts northward in the three model versions, though, not sufficient to explain the proxy reconstruction. This is also true for all comparable simulations in PMIP3 and PMIP4 for the midHolocene , highlighting a shortcoming of current models in simulating the African Humid climate. The CESM2 midHolocene simulation has less northward extent of the monsoon than the comparable simulations in earlier versions of the model. As also noted in the multi-model ensembles  and consistent with proxy reconstructions (Scussolini et al., 2019), the areal extents and total water precipitated for the North African monsoon are enhanced in the 127 ka simulations in CCSM4 and CESM2 as compared to their corresponding simulations at 6 ka.
It should be noted that none of the three model versions include interactive vegetation. The climatic changes over North Africa to the orbital forcing alone at 6 ka and 127 ka want to extend tropical xerophytic shrubland into the Sahara (see next paragraph). Incorporating shrub/savanna vegetation over the Sahara in Holocene simulations has been shown to impact the albedo and evapotranspiration of the surface and water budgets, reducing model and data mismatches (Pausata et al., 2016;Thompson et al., 2019). A vegetated Sahara would also reduce the dust loading over the Sahara. Previous modeling studies have come to different conclusions on the effects of this dust reduction on middle Holocene Saharan rainfall, depending on differences in the optical properties of dust, on whether changes in both dust and vegetation were considered, and whether indirect dust effects were also simulated. With only direct aerosol effects, Pausata et al. (2016) found an increase in Saharan rainfall with an assumed decreased dust loading at 6 ka, while Thompson et al. (2019) found that including the indirect effects of dust dampened this increase. With only direct aerosol effects (less absorbing dust in line with recent observational constraints) but no changes in vegetation, Albani and Mahowald (2019) found a positive feedback of dust on monsoon precipitation, fading as dust level are reduced during the midHolocene; on the other hand, using more absorbing dust would yield the opposite results. A previous modeling study suggests that soil feedback might also be important for driving the African monsoon northward during interglacials (Levis et al., 2004). As concluded by Hopcroft and Valdes (2019), there is no smoking gun in resolving this long-standing model-data mismatch.
temperatures for NH continents (Figure 3) result in expansion of temperate grassland and temperate xerophytic shrubland in North America and Eurasia, with the larger JJA temperature increases for lig127k (Figure 3) leading to expansion of grassland and shrubland over a larger area of North America and Eurasia than for midHolocene (Figure 14a and 14b). Warmer than present JJA temperatures also result in expansion of cold evergreen needleleaf forest northward in northern North America and northern Eurasia for midHolocene and lig127k. In eastern China, the northern limit of warm temperate evergreen broadleaf and mixed forest contracts south to~30°N for midHolocene and~28°N for lig127k. The northern limit of the CESM2 midHolocene warm temperate evergreen broadleaf and mixed forest is similar to that in CCSM4 6 ka biome simulations (Lin et al., 2019), but not as far north as observed for 6 ka (Ni et al., 2010). Although changes in atmospheric CO 2 concentrations may also affect biome distributions (e.g., , the CO 2 concentrations used in the biome simulations for the piControl, midHolocene, and lig127k are similar ( Table 1).
The CESM2-enhanced NH monsoons (Figure 4) produce biome shifts for midHolocene and lig127k compared to piControl (Figure 14). In North Africa, the northward shift of summer monsoon precipitation ( Figure S5) leads to an increase in the northward extent of simulated temperate and tropical xerophytic shrubland vegetation in North Africa for both the midHolocene and lig127k compared to piControl (Figure 14f). The midHolocene shrubland expansion reaches~20°N, a similar extent as in CCSM3 mid-Holocene vegetation simulations (Herold et al., 2012). The lig127k shrubland expansion reaches~22°N across North Africa, with a strip of shrubland extending along the northwest edge of Africa to the Atlas Mountains. This northward expansion of vegetation in the Sahara for midHolocene and lig127k is consistent with interglacial vegetation reconstructions for the region Larrasoaña et al., 2013), although the simulated biomes do not extend as far north across the Sahara as indicated by the reconstructions (Figures 14e and S5).
In the SH, CESM2 austral summer (DJF), temperatures are generally cooler than the piControl simulation for both the midHolocene and lig127k (Figure 3). In South America, tropical savanna and temperate and tropical xerophytic shrubland for midHolocene and lig127k expand compared to piControl (Figure 14). The midHolocene shift to more xeric vegetation is in general agreement with observations ( Figure 14e; Marchant et al., 2009). In southern Africa, temperate xerophytic shrubland expands northward and eastward for midHolocene with expansion of desert for lig127k. Warm temperate evergreen broadleaf and mixed forest decreases in southeast Africa for midHolocene and lig127k, with the midHolocene decreases in general agreement with 6 ka observations (Figure 14b and 14e; Jolly et al., 1998).

Sea Ice
The three model versions (CCSM3, CCSM4, CESM2) differ significantly in their simulations of piControl, midHolocene, and lig127k Arctic sea ice thicknesses and concentrations of Arctic summer (August-September) minima ( Figure 15). In comparison to the earlier versions, CESM2 shows a greater response to the orbital-driven forcing of the MH and LIG simulations in Arctic summer minimum (August-September, AS) sea ice extent and concentrations. The CCSM3 MH and PMIP3 CCSM4 MH simulations show a much smaller change in Arctic summer minimum extent and area as compared to their preindustrial controls, whereas there is a noticeable reduction in summer minimum extent and area in the CESM2 midHolocene. The CESM2 lig127k shows a greater reduction in summer minimum extent and area as compared to the CCSM4 127 ka simulation.
The Arctic summer sea ice thickness shows more sensitivity to MH forcing than sea ice extent in CCSM4, but the reduction in thickness over the central Arctic appears to be less in CCSM4 (−28%) than occurs in CESM2 (−36%). The central Arctic minimum sea ice thickness in CCSM3 shows a similar large sensitivity to the PMIP2 MH forcing compared to the preindustrial control (−34%) as shown in the CESM2 midHolocene. The CCSM3 preindustrial control shows the thickest Arctic minimum summer sea ice. This may be because in CCSM3, a 'present day' simulation with present day (1990 CE) greenhouse gases and orbital parameters was used to tune the top of atmosphere imbalance. The CCSM3 preindustrial control was set up with 1780 CE greenhouse gases and modern orbital parameters which yielded a much thicker central Arctic sea ice in comparison as a result.
Paleo-sea ice distributions from marine cores collected in the Arctic Ocean, Nordic Seas, and northern North Atlantic are available for the Holocene and Last Interglacial based on dinocysts, foraminifers, ostracods, and biomarkers (Stein, Fahl, Gierz, et al., 2017;Stein, Fahl, Schade, et al., 2017;Kageyama et al., 2020), though with some differences in the interpretations based on the different indicators. Following the phytoplankton-IP 25 (PIP 25 ) index of Müller et al. (2011), which combines the environmental information carried by IP 25 and other phytoplankton biomarkers, Stein, Fahl, Schade, et al. (2017) interpret the PIP 25 index as indicating Figure 15. Northern Hemisphere Sea ice thickness (m) during August-September through multiple versions of CESM. Shown are the PI controls and 6 ka simulations for CCSM3, CCSM4, and CESM2 and the 127 ka simulations for CCSM4 and CESM2. Monthly data have been converted to the celestial calendar. Sea ice thickness is shown where sea ice grid box fractional area is greater than 15%. reduced sea ice cover (0.3 to 0.5), seasonal sea ice cover including ice-edge situation (0.5 to 0.7) and extended to perennial sea ice cover (>0.7). The synthesis documented in Kageyama et al. (2020) classifies sites in terms of number of minimum and maximum months with sea ice cover greater than 15%. Stein, Fahl, Gierz, et al. (2017) conclude from their analysis of four ocean sediment cores in the Arctic (see Figure S6 and Supp. Table S2 for locations) dated to the LIG that for the three central Arctic sites perennial sea ice cover was probably predominant. Sea ice conditions were much more variable along the northern Barents Sea north of Svalbard and with the minimum sea ice concentrations towards almost ice-free summers. This agrees with the synthesis presented in Kageyama et al. (2020) that the three sites in the central Arctic were sea ice covered for at least 9 months and possibly perennial sea ice covered, and that the site north of Svalbard could have been ice free up to 6 months. The CESM2 lig127k simulation reproduces well the extent as interpreted in Kageyama et al. (2020), with one sea ice-free month (September) at the more southern sites (PS2138-2 and PS2757-8) but perennially sea ice covered at the two sites poleward of 85 N ( Figure S7). The CCSM4 simulation retains sea ice perennially at all four sites. Interestingly, the CESM2 midHolocene has the southern sites at the sea-ice edge.
For four sites around the periphery of the Arctic Ocean (see Figure S6 and Supp. Table S2 for locations), the biomarker proxy records analyzed by Stein, Fahl, Schade, et al. (2017) show minimum sea ice extent during the early Holocene and reduced, seasonally sea ice at 6 ka. CESM2 midHolocene simulates extended sea ice-free (<15% concentration) conditions during the summer at the locations of the Chukchi and Fram Strait cores and much reduced September sea ice at the locations of the Siberian and Laptev cores ( Figure S8). The CCSM4 and CCSM3 6 ka simulations retain sea ice perennially at all four sites.
Previous studies suggest that the mean-ice state in the control climate can influence the magnitude and spatial distribution of warming at high latitudes in future projections (Holland & Bitz, 2003). Thinner Arctic sea ice is more susceptible to summer melting than thicker Arctic sea ice. Across the family of CESM models, this relationship is clear with a correlation of 0.86 (Figure 16a). Arctic sea ice thickness, averaged for model grid cells with at least 15% concentration, varies substantially in the PI simulations, ranging from 3.2 m in CCSM3, 2.8 m in CCSM4, and 1.8 m in CESM2. CESM2 has much greater reductions in August-September minimum Arctic sea ice area at 6 ka and 127 ka than either CCSM4 or CCSM3.
There is a clear linear relationship between less Arctic sea ice coverage in August-September and warmer August-September average high northern latitude air temperature across the CESM family of models and for these three climate states (Figure 16b). This relationship, with a linear correlation coefficient of −0.91, is consistent with both the largest controlling processes, the albedo effect of sea ice on air temperatures and the insulating effect of sea ice on air-sea exchange. Reduced sea ice coverage lowers the average albedo of the Arctic surface, warming air temperatures and allowing increased ocean heat uptake, while also reducing the insulating effect of sea ice which allows increased air-sea heat exchange.

Dust Aerosols
As noted in Section 4.1.3, the midHolocene captures fairly well the observational constraints, represented by dust mass accumulation rates from paleodust archives integrated for the 6 ka time slice (5-7 ka) (Albani et al., 2015). The evaluation of dust deposition for the piControl case is more complicated, because of the scarcity and uncertainties of records encompassing this time period, as e.g. marine core tops are lost during the drilling process, and there may not be the temporal resolution required (Albani et al., 2016;Carslaw et al., 2017;Hooper & Marx, 2018;Mahowald et al., 2010). However, we can use dust mass accumulation rate estimates for the late Holocene as a loose constraint (Figure 17). In this case too, the overall fit of model results with data is reasonable for the most part, but when we compare the midHolocene/piControl ratio we see some significant deviations. Some of the discrepancies are probably due to the poorer fit of the model with data in the piControl case compared to the midHolocene case; in this respect we stress how the piControl soil erodibility map was not updated consistently with the midHolocene map, therefore without the constraint to match this middle to late Holocene ratio from the observations. In addition, we note how observations from the same region sometimes show opposing values of the midHolocene/piControl ratio (Figure 17e), which could be interpreted as either due to finer temporal and spatial scale features that cannot be captured by the model or by uncertainties in the observational estimates. In any case, we stress how our simulations capture the most relevant changes in the dust cycle, represented by the reduction of dust emissions from the Sahara (Figures 17e, 17f). Figure 17. Comparison of simulated total (wet + dry) dust deposition against observational constraints from paleodust archives. (a,b) midHolocene, (c,d) piControl, and (e,f) the ratio midHolocene/piControl. On the left panels, contour maps represent model output, and filled circles represent observations, plotted using the same color scale. Panels on right depict scatterplots for the corresponding variables and time period, highlighting with colors the geographical locations of the paleodust archives used for the comparison. Dust mass accumulation rates used as observational constraints were integrated from the original time series between 5 and 7 ka for the midHolocene, and between 1 and 3 ka for the piControl case (Albani et al., 2015). The same analysis repeated using a subset of those observational records, with samples reaching into the last millennium, depict a similar picture (not shown).
Previous work with CESM1 indicates a net perturbation of the global atmospheric radiative budget due to dust DRE (effective Direct Radiative Perturbation: eDRP) at the top of the atmosphere (TOA) of −0.06 W m −2 and -0.03 W m −2 for the piControl and midHolocene, respectively (Albani & Mahowald, 2019). Considering that these global budgets result from partial compensation of both positive and negative net perturbations at the TOA, the overall DRE perturbation to the atmospheric radiative budget was better expressed by the absolute values of|0.27|W m −2 (piControl) and|0.23|W m −2 (midHolocene), i.e. the global average of the absolute values of eDRP (Albani & Mahowald, 2019). The same study found that dust tends to enhance the West African monsoon, potentially acting as a negative feedback mechanism during the middle to late Holocene transition, i.e. as more humid conditions tend to suppress dust emissionstherefore reducing dust amplification of the monsoon.
In CESM2 a different aerosol scheme is used, i.e. a modal scheme with dust distributed in two modes and accounting for aerosol indirect effects, whereas the work previously done with CESM1 relied on a bulk aerosol scheme with no inclusion of indirect effects. Since diagnosing the impact of aerosol-cloud indirect effect requires additional simulations (Wang et al., 2011), which were not done here, we do not know the impact of the changes in dust on the liquid clouds and the resulting impact on radiative forcing, but these are likely to be small due to the small change in cloud condensation nuclei numbers from coarse mode dust particles (Mahowald et al., 2011).

Summary and Discussion
The midHolocene and lig127k experiments with CESM2 in CMIP6 and PMIP4 provide an out-of-sample assessment of the future projection experiments with CESM2. In concert with paleo-data for these time periods, they allow evaluations of the responses and connections within CESM2 to different radiative forcings. Particular foci of the midHolocene and lig127k experiments are high-latitude amplification of temperatures and low-latitude monsoon precipitation. They are also important for understanding and explaining past environments recorded in the geologic record.
The primary forcing for the midHolocene and lig127k are the orbitally induced changes in the seasonal and latitudinal distribution of insolation. In particular, the orbital configurations resulted in positive NH summer insolation anomalies at 6 ka and 127 ka as compared to PI. The anomalies are considerably greater at 127 ka than at 6 ka, associated with greater eccentricity and boreal summer solstice at perihelion at 127 ka. Atmospheric greenhouse gas levels were similar to those of the preindustrial period, land ice likely only remained over Greenland and Antarctica, and the continental configurations were almost identical to modern.
In response to the positive NH summer insolation anomalies, both experiments simulate JJA warming at NH middle and high latitudes (particularly over North America, Eurasia, and Greenland), reduced Arctic summer sea ice, and enhanced NH summer monsoon precipitation. These responses are greater in the lig127k than midHolocene simulation. Notably, the Arctic goes almost ice free in September in the lig127k simulation with thin sea ice in the central Arctic, and Greenland has much more extensive ablation zones in its western and northern coastal zones suggesting potential retreat if not fixed in the simulation. The lig127k simulation simulates the strong summer warming reconstructed for the Baffin region. Both the areal extent and cumulative rainfall is greater for the North African monsoon in the lig127k than midHolocene simulation Otto-Bliesner et al., 2020).
The most consistent regional pictures in the annual temperature reconstructions for 6 ka and 127 ka are the positive temperature anomalies over Greenland for both time periods and at 127 ka for Antarctica. CESM2 simulates the patterns but not the magnitudes of the reconstructed anomalies. Several regional features are in general agreement between the proxy records and the midHolocene simulation. These include east-west patterns of positive/negative 6 ka surface air temperature anomalies as compared to PI in eastern/western U.S., eastern China/Tibet, and western/eastern North Atlantic. The lig127k simulation also displays similar patterns but little data are available to confirm.
Vegetation integrates the effects of temperature and precipitation. Although the CESM2 midHolocene and lig127k simulations do not include predictive vegetation, the global terrestrial biomes simulated by BIOME4 using CESM2 midHolocene and lig127k climate data can be used for comparison to observed global vegetation patterns. Warmer than present JJA temperatures result in expansion of cold evergreen needleleaf forest northward in northern North America and northern Eurasia for midHolocene and lig127k, in generally good agreement with pollen and macrofossil evidence. Xerophytic shrubland expands in northern Africa in the midHolocene, and even more so in the lig127k in response to the simulated greater area extent of cumulative rainfall of the North African monsoon. Comparison of the midHolocene simulated to the 6 ka observed biomes indicates that CESM2 underestimates this expansion.
Globally the largest differences in the barotropic ocean transport compared to the preindustrial control are found in the Southern Ocean and the North Atlantic, with a weaker Antarctic Circumpolar Current and strengthened North Atlantic subtropical and subpolar gyres in both paleoclimate states, but to a greater degree in the lig127k as compared to the midHolocene. The AMOC and AABW are only weakly sensitive to changes in the combined forcing from the seasonal and latitudinal redistribution of insolation. Improvements in ENSO characteristics (period, amplitude, seasonal cycle) in CCSM4 as compared to CCSM3 are retained in CESM2, though the Niño3.4 SST index in the piControl simulation is overestimated by~30%. The Niño3.4 index has reduced power going from the piControl to midHolocene to lig127k, with a 20% reduction in the Niño3.4 standard deviation in the lig127k as compared to the piControl and a much flatter annual cycle in the lig127k simulation. The seasonal cycle of Pacific equatorial SST strengthens in lig127k and weakens in midHolocene compared to the piControl. The seasonal anomalies in the eastern Pacific relative to the seasonal cycle of the piControl appear phased with tropical insolation anomalies. Colder annual mean SST across the tropics plus enhanced seasonal cooling in DJF with respect to the piControl, especially notable in lig127k, may be responsible for muting ENSO events which tend to peak at the end of the year.
Model differences are considerable between the older versions of the models (CCSM3 and CCSM4) and CESM2, which are associated with significant improvements made in the physics of the component models and resolutions of the simulations. The summer monsoon over North Africa intensifies and shifts northward in the three model versions, though, less so than indicated by observed biomes. The family of models agree on warmer summers at middle and high NH latitudes with this warming greater at 127 ka than 6 ka. Differences in the polar surface temperature anomalies are related to simulation of sea ice thickness and extent in the respective preindustrial simulations with CESM2 having an improved reproduction of these features as compared to observationally-based estimates (Danabasoglu et al., 2020).
The selection of only two intervals, midHolocene and lig127k, for CMIP6-PMIP4 Tier1 simulations was designed to maximize both the multi-model ensemble size and opportunities for model evaluation with data syntheses. Indeed, 13 midHolocene and 17 lig127k simulations were contributed by modeling centers to CMIP6. Uncertainties in the boundary and initial conditions for the MH and LIG as well as missing processes in CESM2, and in simulations by many other modeling groups, may be important for explaining mismatches with the data. These will be explored in future experiments with CESM2.
Recognizing this, PMIP proposed CMIP6 Tier 2 sensitivity experiments for the midHolocene and lig127k to explore the feedbacks between vegetation and climate: midHolocene-veg and lig127k-veg (Otto-Bliesner et al., 2017). Reconstructions indicate that xerophytic shrubland (Figures 14 and S5) extended into the Sahara. Previous modeling studies suggest that incorporating shrub/savanna vegetation over the Sahara as well as soil feedbacks in CCSM/CESM are important for the northward extension of the African monsoon during the Holocene (Levis et al., 2004;Thompson et al., 2019). Replacement of tundra with boreal forests over the NH high-latitude continents, as indicated by pollen and macrofossil data (Figure 14), amplifies high-latitude warming in CCSM/CESM simulations (Feng et al., 2020;Swann et al., 2010). In addition, a new optional feature in the CESM land model is the demographically structured dynamic vegetation model (Functionally Assembled Terrestrial Ecosystem Simulator), which will allow prediction of biome boundaries directly from plant physiological traits Lawrence et al., 2019).
Matches of the CESM2 lig127k simulation and data reconstructions could be affected by both the design of the CMIP6 experiment and uncertainties in the dating of the LIG records. The 127 ka time slice for the CMIP6-PMIP4 LIG experiment was chosen to account for dating uncertainties and the imprint of the previous deglaciation and Heinrich 11 (H11) meltwater event on temperature. Data reconstructions for 127 ka (Figures 12 and 13) suggest that the bipolar response of the climate system to melting of the NH ice sheets, with cold surface temperature anomalies in the North Atlantic and warm surface temperature anomalies over the Southern Ocean and Antarctica, may still be present at 127 ka. In addition, the lig127k simulation protocols included the Greenland and Antarctic ice sheets prescribed as modern. However, it is possible that the West Antarctic Ice Sheet retreated early in the LIG and that the Greenland ice sheet was reduced in extent as compared to present. The CMIP6 Tier 2 sensitivity experiments: lig127k-H11, lig127-ais, and lig127k-gris, will allow the effects of these uncertain protocols to be explored. In addition, a transient LIG simulation is providing an opportunity to examine the time-dependent evolution of the climate and Greenland ice sheet in CESM and CISM to the LIG orbital forcing.
The CESM community has been participating in about 20 CMIP6-endorsed MIPs, with results documented in the AGU CESM2 Virtual Special Issue (https://agupubs.onlinelibrary.wiley.com/doi/toc/10.1002/ (ISSN)1942-2466.CESM2). The CESM2 lig127k and midHolocene experiments have potential implications for confidence in future projection simulations with CESM2 (DeRepentigny et al., 2020;Meehl et al., 2020). More than half of the 17 CMIP6 lig127k simulations have a retreat of the Arctic minimum (August-September) ice edge similar to the average of the last two decades . Across these CMIP6 models, Kageyama et al. (2020) noted a nearly linear relationship between the simulations of Arctic summer sea ice in their 1pctCO2 simulation at time of doubling and their lig127k simulations. That is, the models which respond strongly to the LIG forcing also respond strongly for the 1pctCO2 forcing. CESM2 has a high equilibrium climate sensitivity (ECS) of 5.3°C (Gettelman et al., 2019); HadGEM3 similarly has a high ECS of 5.5°C. Both predict an almost ice-free or ice-free Arctic in their lig127k experiments and the predicted year of disappearance of September sea ice in the SSP8-8.5 scenario of 2038 and 2035, respectively . The CMIP6 and PMIP4 paleoclimate simulations with CESM2 combined with proxy reconstructions allow an assessment of whether this high ECS is plausible (Feng et al., 2020;Zhu et al., 2020).
With the availability of LIG ice and marine core records, the transient CESM2-CISM2 LIG simulation with an evolving Greenland ice sheet can inform the ISMIP6 future projection simulations with CESM2 and their predictions of sea level rise. The CMIP6 (and CESM2) paleoclimate simulations thus provide an opportunity to better constrain projections of future climate change.

Data Availability Statement
CESM2 is freely available at http://www.cesm.ucar.edu/models/cesm2/. Simulation data for CESM2 used in this study are freely available from the Earth System Grid Federation (ESGF) at esgf-node.llnl.gov/search/ cmip6 and CCSM4 from the Climate Data Gateway at NCAR at https://www.earthsystemgrid.org/.