Retreat and Regrowth of the Greenland Ice Sheet During the Last Interglacial as Simulated by the CESM2‐CISM2 Coupled Climate–Ice Sheet Model

During the Last Interglacial, approximately 129 to 116 ka (thousand years ago), the Arctic summer climate was warmer than the present, and the Greenland Ice Sheet retreated to a smaller extent than its current state. Previous model‐derived and geological reconstruction estimates of the sea‐level contribution of the Greenland Ice Sheet during the Last Interglacial vary widely. Here, we conduct a transient climate simulation from 127 to 119 ka using the Community Earth System Model (CESM2), which includes a dynamic ice sheet component (the Community Ice Sheet Model, CISM2) that is interactively coupled to the atmosphere, land, ocean, and sea ice components. Vegetation distribution is updated every 500 years based on biomes simulated using a monthly climatology to force the BIOME4 equilibrium vegetation model. Results show a substantial retreat of the Greenland Ice Sheet, reaching a minimum extent at 121.9 ka, equivalent to a 3.0 m rise in sea level relative to the present day, followed by gradual regrowth. In contrast, a companion simulation employing static vegetation based on pre‐industrial conditions shows a much smaller ice‐sheet retreat, highlighting the importance of the changes in high‐latitude vegetation distribution for amplifying the ice‐sheet response.

• A simulation from 127 to 119 ka with CESM2-CISM2 reveals retreat and incipient regrowth of the Greenland Ice Sheet for the Last Interglacial • The Greenland Ice Sheet reaches a minimum extent at 121.9 ka with a contribution of 3.0 m sea-level equivalent relative to present day • Simulated forest and tundra expansion in Canada and Greenland leads to greater warming and ice loss than with a pre-industrial land cover (Berger & Loutre, 1991) provided the primary forcing for higher summer temperatures over most of the northern hemisphere (NH). Feedbacks among the cryosphere, ocean, and atmosphere (Masson-Delmotte et al., 2013;Otto-Bliesner et al., 2013) allowed these higher surface temperatures to be maintained year-round in the Arctic and over Greenland (Otto-Bliesner et al., 2020). Although the circumstances are different, understanding the behavior, processes, and feedbacks in the Earth system during the LIG provides insights relevant to what we might expect during future global warming.
Global mean sea level during the LIG was likely 6-9 m higher than present (Dutton et al., 2015;Kopp et al., 2009Kopp et al., , 2013. The magnitude of peak sea level is still an active research area, with new insights being gained from glacial isostatic adjustment and dynamic topography modeling (Austermann et al., 2017;Whitehouse, 2018). By ∼127 ka, evidence from proximal marine records suggests little or no remnants of the North American and Eurasian ice sheets (Otto-Bliesner et al., 2017).
It is an open question as to what portion of sea-level rise was due to the retreat of the GrIS, as opposed to the collapse of parts of the Antarctic Ice Sheet (Clark et al., 2020) and thermal expansion of the oceans (McKay et al., 2011;Shackleton et al., 2020;Turney et al., 2020). In a synthesis of modeling and data studies, Dutton et al. (2015) estimated 0.6-3.5 m in sea-level rise due to contributions from Greenland.
The first ice-sheet modeling estimates of the GrIS contribution to peak LIG sea level used the shallow ice approximation (SIA; Hutter, 1983) and "index methods" to estimate the climate forcing based on reconstructed ice-core temperatures (Cuffey & Marshall, 2000;Huybrechts, 2002;Lhomme et al., 2005;Tarasov & Peltier, 2003). Snapshot simulations with general circulation models (GCMs) have provided a more direct attribution/connection of the influence of orbital forcing on the evolution of the LIG GrIS, though only using one-directional climate-ice sheet coupling (Born & Nisancioglu, 2012;Otto-Bliesner et al., 2006;Robinson et al., 2011;. These earlier ice-sheet model estimates used relatively coarse horizontal resolution (∼10-20 km) and computed the SMB using either a positive degree-day scheme (most) or an empirical scheme accounting for solar insolation (Robinson & Goelzer, 2014). Each of these studies developed their best estimates of GrIS melting by choosing simulations that were most consistent with independent evidence from ice cores. Coupled climate-ice sheet simulations have also been attempted, with limited success. Helsen et al. (2013) asynchronously coupled a global climate model to a regional climate model and an ice sheet model. Goelzer et al. (2016) used a coarse resolution (∼5° atmosphere, 3° ocean), coupled GCM-ice sheet model to simulate the Eemian GrIS evolution but required scaling of the temperature-anomaly forcing to prevent the complete collapse of the GrIS.
In the past several years, significant progress has been made in coupling ice sheet models within complex Earth system models. In particular, the current versions of the Community Earth System Model and the Community Ice Sheet Model, CESM2 and CISM2, have been coupled interactively (Muntjewerf et al., 2021). The aim of this study is to use CESM2 and CISM2 to conduct a transient, coupled global simulation of the LIG climate with a dynamic, evolving GrIS. This paper presents the results of this work as an improved modeling estimate of how the GrIS responded during this period. The sections below describe the modeling methods, results, ice-sheet behavior, and implications.

Methods
We conduct an 8000-year transient simulation over the interval from 127 to 119 ka starting from the Paleoclimate Modeling Intercomparison Project Phase 4 (PMIP4) and Coupled Model Intercomparison Project Phase 6 (CMIP6; Eyring et al., 2016) lig127k simulation with CESM2 (hereafter CMIP6 lig127k; Otto-Bliesner et al., 2020). The 127 ka time slice was chosen for CMIP6 to capture the large, positive NH seasonal insolation anomalies, with few or no remnants of the North American and Eurasian ice sheets, the GrIS similar to its current extent, and sufficient time to minimize the imprint of the previous deglaciation and the Heinrich-11 meltwater event (Marino et al., 2015). This time span to 119 ka captures the declining NH summer insolation anomalies and effects of the cooling climate on the GrIS.

Model Description
We use the Community Earth System Model version 2.1 (CESM2; http://www.cesm.ucar.edu; Danabasoglu et al., 2020), a global climate model, which consists of interactively coupled, prognostic component models of atmosphere, ocean, land, sea ice, and (optionally) land ice. The atmosphere model is the Community Atmosphere Model (CAM6); the ocean model is the Parallel Ocean Program (POP2); the land model is the Community Land Model (CLM5); the sea-ice model is CICE5, and the land-ice model is the Community Ice Sheet Model version 2.1 (CISM2). CISM2 is a higher-order thermo-mechanical ice-sheet model that solves equations for mass, momentum, and energy of ice sheets (Lipscomb et al., 2019). The various prognostic component models are connected through the Common Infrastructure for Modeling the Earth (CIME), which contains the coupling infrastructure, support scripts, data models, and utility libraries for running CESM2 as a single executable. The simulations are run on the Cheyenne supercomputer at the NCAR-Wyoming Supercomputing Center, with a typical cost of 3500 CPU-hr per simulated year, and a throughput of ∼27 CESM2 years per wall-clock day. For 8000 ice-sheet years accelerated by a factor of 5 (described below in Section 2.4), the simulation of 1600 CESM years costs ∼5.6 million CPU-hr.
For this study, CESM2 is run at the same nominal 1-degree resolution as the CMIP6 DECK (Diagnostic, Evaluation, and Characterization of Klima) and future scenario experiments (Eyring et al., 2016 Nowicki et al., 2016).  and Muntjewerf, Sellewold, et al. (2020) presented results of coupled, transient CESM2-CISM2 simulations for future climate scenarios; this study presents the first such simulation for a paleoclimate scenario. CISM2 is run on a fixed rectangular mesh at a resolution of 4 km, with a depth-integrated higher-order velocity solver (Goldberg, 2011), a pseudo-plastic basal sliding law (Aschwanden et al., 2016), and ice assumed to calve as soon as it becomes afloat. Isostatic rebound is computed using an elastic lithosphere-relaxing asthenosphere model with a relaxation time scale of 3000 years (Rutt et al., 2009). See Lipscomb et al. (2019) for more details on CISM2 dynamics. Muntjewerf et al. (2021) describe the CESM2-CISM2 coupling in detail. The ice-sheet surface mass balance (SMB) is computed in the land model, CLM5, with multiple elevation classes in each glaciated grid cell to better resolve sub-grid-scale topography. The SMB is accumulated and averaged on annual time scales and downscaled to the CISM2 grid, with linear interpolation between elevation classes in the vertical direction and bilinear interpolation in the horizontal. Normalization factors are applied to the accumulation and ablation to conserve the total water mass. CISM2 then evolves dynamically, and the new ice-sheet extent and thickness are returned to CLM5, where the land surface types (i.e., glaciated vs. vegetated land units) and surface elevation are updated annually to be consistent with CISM2. In CAM6, the surface elevation over Greenland is updated every 10-atmosphere years, using an offline tool (described in Lofverstrom et al., 2020) in a version of CIME with the capability for case-specific workflows. The ice sheet and ocean models are partly coupled in the form of solid (calving) and liquid (basal melting) freshwater fluxes sent from CISM2 to POP2. No floating ice shelves are in this simulation (hence no sub-ice-shelf melting), and the ocean geometry is fixed.

Initialization
Initializing ice sheet simulations with CMIP6-class models is critically important but challenging. For the present day, various strategies have been implemented, including spin-up from an earlier stage of ice-sheet history and optimization using satellite observations of velocity and ice thickness (Goelzer et al., 2018). Differences in initial states translate into large differences between models in their projections of the future sea-level contribution of the GrIS (Goelzer et al., 2020). Previous simulations of the response of the LIG GrIS have usually adopted either the "index method" using Antarctic and Greenland ice cores, or asynchronous coupling to climate model simulations, to approximate a temperature and precipitation history for Greenland over past glacial-interglacial cycles (Helsen et al., 2013;Quiquet et al., 2013;. In both cases, the ice model parameter settings were chosen to be most consistent with independent evidence from ice cores. Here, we adopt the initial state of the GrIS from the spun-up CESM2-CISM2 under pre-industrial climate forcing using the procedure described by Lofverstrom et al. (2020). This initial state is chosen for several reasons. First, the spin-up procedure accounts for the memory of the past climate during the last deglaciation and is geometrically and thermodynamically consistent with the simulated pre-industrial climate state. Second, the initial state is the same as used in future projections with CESM2 and CISM2   (Menviel et al., 2019) and the non-local bedrock signal from surrounding ice sheets (Bradley et al., 2018). Thus, we do not have validation benchmarks for a spin-up under the climate conditions leading up to the LIG.
The spun-up GrIS contains 8.2 m sea-level equivalent (SLE), about 12% more than the present-day ice sheet which contains 7.4 m SLE (Morlighem et al., 2017). Here and elsewhere, ice mass and volume are stated in units of SLE, assuming 1 mm SLE = 361.8 Gt ice, with a uniform ice density of 917 kg m −3 . Because of the greater mass, the bedrock beneath the GrIS is slightly depressed compared to the present day. At the end of the spinup, the ice sheet is close to equilibrium with the pre-industrial climate. At the start of our transient simulation, an initial rapid decline in ice-sheet area and volume is attributed to differences between pre-industrial and LIG forcing. This initial "shock" to the ice sheet disappears in the first 500 years during which low-lying marginal ice is melted.
The atmosphere, ocean, land, and sea ice components are initialized from the CMIP6 lig127k simulation with CESM2 (Otto-Bliesner et al., 2020). The potential vegetation surface boundary condition of the CESM2 lig127k simulation is replaced by vegetation biomes as simulated by BIOME4 (Kaplan et al., 2003) using monthly mean climate data from the CMIP6 lig127k simulation (see Section 2.5 for more details).

Last Interglacial Orbital Configuration and Greenhouse Gases
Changes in Earth's orbital elements over the interval from 127 to 119 ka were large. Perihelion advanced from early July to November, and obliquity decreased from 24.04° to 22.86° over the same period (Figures 1a-1c). At a time of relatively high eccentricity (∼0.04), these orbital changes resulted in a decrease in boreal summer (e.g., July) insolation at 65°N, for example, from 483.0 W m −2 at 127 ka (55.6 W m −2 greater than present), to 426.1 W m −2 at 119 ka (1.2 W m −2 lower than present), a decrease of about 13% ( Figure 1d). As a result of changes in the timing of perihelion and eccentricity, the length of individual months defined using an angular (or celestial) calendar also changed (Joussaume & Braconnot, 1997;Kutzbach & Gallimore, 1988), with boreal summer months increasing in length over the interval, and boreal autumn and early winter months decreasing in length (e.g., July defined using the angular calendar lengthened from 27.81 to 30.89 days, and November shortened from 32.47 to 28.31 days). The CESM2 output used as input to the vegetation model described below (Section 2.5) was summarized using a present-day "no leaps" calendar, so we apply month-length corrections to the input data for BIOME4 following the approach of Bartlein and Shafer (2019).
Atmospheric CO 2 concentration is held constant at 275 ppm from 127 to 119 ka, the value adopted by Otto-Bliesner et al. (2017, see their Table 1) for the CMIP6 lig127k experiment. In contrast, Köhler et al. (2017) suggested that the CO 2 concentration ranged between 268.01 and 279.04 ppm over the interval, with an average of 275.46 ppm, which would produce small variations in radiative forcing relative to the constant 275 ppm values.
Other boundary condition values are set using those of the lig127k experiment, except for vegetation cover (see Section 2.5 below).

Acceleration of Orbital and Ice-Sheet Model Parameters
To facilitate an 8000-year transient simulation, we accelerate the ice-sheet model by a factor of 5 relative to the broader Earth-system model. This acceleration is chosen based on computational costs and response times of the different model components. The atmosphere model is much more computationally expensive than the ice-sheet model, so it is advantageous to run multiple ice-sheet years for each atmosphere year. An acceleration factor of 5 for the ice-sheet model was selected after preliminary design simulations with factors ranging from 2 to 20, as well as tests to examine the response of the ocean to freshwater fluxes from the GrIS (ocean temperature, heat uptake, salinity, circulation patterns, and other diagnostic indicators). For every CESM2 year, the resulting climate is used to force CISM2 for 5 years, with fluxes from the final CISM2 year sent back to the coupler. This results in a violation of freshwater mass conservation, but the freshwater flux is relatively low by oceanographic standards: the maximum freshwater flux of 0.0184 Sverdrup (Sv; 1 Sv = 10 6 m 3 s −1 ) occurs near the beginning of the simulation when the climate is warmest, with a mean freshwater flux for the entire simulation of 0.0069 Sv. To put this magnitude into context, the total mean freshwater river runoff in the Arctic from 2000-2010 was ∼0.13 Sv, an order of magnitude higher, not including input from Greenland (Haine et al., 2015).
Changes in the orbital parameters over time (eccentricity, obliquity, and precession) are accelerated with the ice-sheet year, advancing 5 years for each CESM2 year. The atmosphere adjusts quickly to changes in the overall model state (typically within a few months), while the deep ocean responds on a time scale of hundreds or thousands of years. Therefore, if the acceleration is too great, the deep ocean state will not keep up with the changing climate (Varma et al., 2016). To complete the 8000-year transient ice-sheet and orbital simulation, this acceleration scheme requires a total of 1600 coupled CESM2 years.

Representation of Vegetation Changes
For other CESM2-only Quaternary applications (e.g., Otto-Bliesner et al., 2020), vegetation distribution in the land model is prescribed based on potential vegetation in the absence of anthropogenic influence (i.e., by removing the agricultural and urban land cover), under a pre-industrial climate. Observed data (e.g., pollen records) suggest that, during the LIG warm period, the distribution of global vegetation was different than present (Hoogakker et al., 2016;Tarasov et al., 2013). At NH polar latitudes, pollen and plant macrofossil evidence indicate that, except in Alaska and central Canada, boreal forest extended to the Arctic coast in some areas (Edwards et al., 2003;LIGA, 1991;Lozhkin & Anderson, 1995), with corresponding lower albedo and increased transpiration (Swann et al., 2010), resulting in enhanced Arctic warming. These vegetation changes and their feedbacks on the climate system are important to include when simulating the GrIS evolution during the LIG.
The land component of CESM2 (CLM5) does not include functional dynamic vegetation. Here we use BIOME4 version 4.2b2, an equilibrium vegetation model (Kaplan et al., 2003), to simulate global biomes at 500-year intervals for 127-119 ka during the LIG, and for the pre-industrial time period (approximately 1850 CE). Input climate data for BIOME4 consist of monthly climatologies of near-surface air temperature (°C), percent possible sunshine (%), total precipitation (mm), and absolute minimum temperature (°C), which were processed following the protocol of Harrison et al. (2014). CESM2 total cloud percent is converted to percent possible sunshine using geographically weighted regressions between CRU CL 1.0 (an observational dataset; New et al., 1999) percent clouds and CRU CL 2.0 (New et al., 2002) percent sunshine. Absolute minimum temperatures are calculated by applying the experiment minus control long-term mean differences in the mean temperature of the coldest month to present-day absolute minimum temperatures. For the LIG biome simulations, we use 20-year monthly climatologies from the end of each 100-year CESM2 simulation segment. With orbital and ice-sheet acceleration, these 20-year climatologies represent the last 100 orbital years of each 500-year orbital and ice-sheet interval ( Figure  S1 in Supporting Information S1). The CESM2 climate data are regridded to a 0.5-degree global grid resolution using the Climate Data Operators version 1.9.3 (http://mpimet.mpg.de/cdo) 'remapbil' bilinear interpolation function. The 0.5-degree data are calendar adjusted using PaleoCalAdjust v1.0 (Bartlein & Shafer, 2019).
To correct for biases in the 20-year monthly climatologies, long-term mean differences are calculated for each climate variable with respect to the CESM2 pre-industrial control (piControl) simulation. The long-term mean differences are then applied to CRU CL 2.0 data (30-year mean data for 1961-1990 regridded to the 0.5-degree grid; New et al., 2002). This approach captures the spatial variability of climate that strongly influences vegetation in topographically complex regions.
The CRU CL 2.0 air temperature data over the GrIS and northeast Canada ice caps are lapse-rate-adjusted to the bedrock elevations under the ice using the elevation differences between the ETOPO1 bed and ice elevation data (Amante & Eakins, 2009). In addition to the climate data, the BIOME4 input data include orbital eccentricity and obliquity for each 500-year interval (Bartlein & Shafer, 2019), and soil water-holding capacity and percolation rate (Harrison et al., 2014). We use pre-industrial soil properties for both the pre-industrial and LIG biome simulations. For grid cells deglaciated by the simulated GrIS retreat during the LIG, we extrapolate soil water-holding capacity and percolation rate to the deglaciated grid cells from nearby grid cells with pre-industrial soil data.
For the pre-industrial biome simulation, input climate data for BIOME4 are from the CESM2 CMIP6 DECK piControl simulation (Otto-Bliesner et al., 2020). The piControl temperature, precipitation, and cloud data are re-gridded using the same procedure as described above for the LIG climate data but are not calendar-adjusted.
The piControl climate data are bias-corrected using the CESM2 historical simulation as the base period for calculating long-term mean differences. Atmospheric CO 2 concentration for the pre-industrial biome simulation is 284.7 ppm.
The simulated biomes for each 500-year time step from 127-119 ka are regridded and mapped to corresponding percentages of plant functional types (PFTs) in CLM5 using a biome-to-PFT table developed for DeepMIP . For areas of Greenland above 60°N latitude and between 70°W and 10°W longitude, BI-OME4-simulated cold evergreen needleleaf forest (corresponding to 80% boreal needleleaf evergreen tree, 10% boreal broadleaf deciduous tree, and 10% boreal broadleaf deciduous shrub PFTs) is replaced with the low and high shrub tundra biome (100% boreal broadleaf deciduous shrub PFT) for deglaciated land grid cells that previously had 100% ice cover at 127 ka. This adjustment was made because deglaciated regions more realistically support tundra rather than forest immediately after deglaciation. The updated PFT land surface data are used in the simulation of the next 500 CISM2 years (100 CESM2 years). This offline coupling of biome distributions every 500 ice-sheet years allows for the incorporation of quasi-dynamic vegetation changes in the simulation. A schematic of coupling between the models is shown in Figure S1 in Supporting Information S1.

Results
We present results from our 127-119 ka transient simulation, focusing on the evolution and response of the GrIS in this interglacial period. The full model output contains additional global climate data documenting the transient climate of the LIG and the influence of a diminishing GrIS.

Temperature and Vegetation
The simulation period 127-119 ka was selected to begin with the highest Arctic summer temperature anomalies and to capture the influence of this initially warm climate and subsequent cooling while beginning after the Heinrich-11 event and at a time when the GrIS may have been similar to, or slightly larger than, its present-day size (Otto-Bliesner et al., 2017). By 121 ka, the summer temperature anomalies relative to the pre-industrial piControl climate decrease except over the deglaciated parts of Greenland, where these anomalies reflect the temperature differences between a vegetated land surface during the LIG and an ice-covered pre-industrial surface (Figures 2  and S4 in Supporting Information S1).
The simulated biomes change substantially over the course of the LIG (Figures 2 and S6 in Supporting Information S1), evolving in response to the changing climate. From 127-123 ka, the northern border of boreal forest (e.g., cold evergreen needleleaf forest and cold deciduous forest biomes) is simulated to extend farther north than its pre-industrial extent, before starting to retreat southward (Figures 2 and S7 in Supporting Information S1). Similarly, the southern border of boreal forest migrates poleward from 127-121 ka compared to pre-industrial but moves equatorward after 121 ka as the climate cools (Figures 2 and S7 in Supporting Information S1). The biome simulations indicate expansion and subsequent retreat of tree taxa from the Arctic coast in northern Canada and Russia, in general agreement with inferred LIG vegetation (Edwards et al., 2003;Hoogakker et al., 2016;Kienast et al., 2011;LIGA, 1991;Lozhkin & Anderson, 1995;Tarasov et al., 2013). In Greenland and northeast Canada, warmer LIG temperatures also simulate increases in forest and shrub biomes, implying the presence of woody taxa, which is generally supported by observations, although data are limited (CAPE-Last Interglacial Project Members, 2006;Crump et al., 2021;Fréchette et al., 2006).

Response of the Greenland Ice Sheet
During the LIG, the GrIS recedes significantly from its initial extent. Figure 3 shows ice thickness, surface mass balance, and surface velocity every 2000 model years (see Figures S12-S14 in Supporting Information S1 for 500-year snapshots). The ice sheet retreats under the warm Arctic climate during the first half of the LIG, but this retreat slows and reverses as the NH summer temperature anomalies decrease later in the simulation (Figure 2).
The GrIS exhibits a distinctive pattern of retreat. The most pronounced ablation and retreat occur along the western edge of the ice sheet, especially in the southwest, and along the northern margin. Meanwhile, the eastern extent remains relatively stable throughout the simulation, with ablation generally limited to the ice-sheet margin and a net accumulation in the southeast. This pattern of accumulation and ablation persists through the simulation, maintaining ice thickness in the east as the western ice margin retreats inland. The regional difference in ice accumulation and ablation leads to an eastward migration of the ice divide, as shown by the surface velocity (Figures 3 and S14 in Supporting Information S1). As the western portion of Greenland becomes increasingly deglaciated and exposed, a saddle forms between the larger northern part of the ice sheet and the southern dome. The saddle point becomes progressively thinner until the two masses separate around 122 ka. After the minimum ice volume and area are reached at 121.9 ka, the two domes reconnect by ∼120.5 ka ( Figure S14 in Supporting Information S1).
To further examine how the ice sheet evolves, Figure 4 shows timeseries of mass-balance fluxes over the simulation period. The total mass balance is the sum of three component fluxes: calving flux (iceberg discharge), basal mass balance (BMB) from ice melting at the bed, and surface mass balance (SMB) from accumulation and ablation. The negative total mass balance is initially dominated by the calving flux, but this decreases in magnitude as the ice-sheet margin retreats from the coast and low-lying marginal areas. The BMB component is small compared to the calving flux; while it does vary, its magnitude remains fairly constant relative to the other fluxes. The SMB is generally negative for the first ∼4,000 years of the simulation (127-123 ka), implying more ablation than accumulation, and thus a net mass loss. The interannual variability is large, however, and some years have a net positive SMB even during the warmest part of the simulation. The total mass balance is negative for the first ∼5,000 years (127-122 ka) of the simulation. Together, the calving flux and BMB compose the freshwater flux from CISM2 to the ocean. Liquid meltwater runoff is sent to the ocean from the land model. The maximum freshwater flux (0.0184 Sv) occurs early in the simulation when the ice sheet is most extensive and the Arctic summer insolation anomalies are highest, while the minimum freshwater flux (0.0039 Sv) occurs at 122.9 ka. The influence of increased freshwater from Greenland as the ice sheet retreats has a negligible effect on the Atlantic meridional overturning circulation.
The GrIS volume decreases almost linearly for the first 3000 years of the simulation (until 124 ka), at which point its collapse slows ( Figure 5). Ice-sheet area decreases rapidly in the first few hundred years of the simulation as thin areas of the ice sheet near the coast are lost, and more slowly thereafter. Figure 5 shows a time series of GrIS  Figure S1 in Supporting Information S1). (f-j) BIOME4 simulated biomes are shown on the 0.5-degree grid used for biome simulations. Greenland cold evergreen needleleaf forest is replaced with low and high shrub tundra as described in Section 2.5. (k-o) Simulated change in forest extent compared to simulated pre-industrial forest extent ( Figures S6 and S7 in Supporting Information S1). For each land grid cell, gray indicates a forest biome is simulated by BIOME4 for both the LIG and pre-industrial time periods, tan indicates a forest biome is simulated for the pre-industrial but not for the LIG time period, and dark green indicates a forest biome is simulated for the LIG but not for the pre-industrial time period. White land grid cells indicate that a non-forest biome is simulated for both the LIG and pre-industrial time periods. Greenland Ice Sheet extent is displayed for grid cells with ≥50% ice cover in the BIOME4 veg simulation. Figures S2-S7 in Supporting Information S1 show data for other time periods. We ran a companion simulation for 6000 years, from 127-121 ka, with prescribed pre-industrial vegetation cover assuming no anthropogenic influence (no-Anthro veg, the same vegetation used in the CMIP6 lig127k simulation; Otto-Bliesner et al., 2020). The evolving GrIS volume and area are shown for comparison in Figure 5. Apart from the vegetation distribution remaining constant, this simulation follows the same methods as described in Section 2 for the main simulation. In this simulation, the GrIS deglaciation is muted, and the ice sheet does not lose nearly as much volume and area as with time-evolving vegetation. Instead, the ice sheet reaches a minimum area at 124.3 ka and a minimum volume of 6.8 m SLE at 124.1 ka, after which it starts to regrow. The total ice loss at 124.1 ka is 1.4 m SLE, or 0.6 m SLE relative to the present day. The simulation was stopped at 121 ka after regrowth was well underway. The influence of vegetation on ice-sheet response is discussed further in Section 4.2.

Comparison to Geologic Data
Ice core and marine records allow an evaluation of our simulation of GrIS extent during the Last Interglacial. Six deep ice-coring projects have reached ice layers back to the LIG. Enriched-δ 18 O ice is preserved at the base of all deep boreholes (reaching bedrock), though records are either not complete or not stratigraphically intact. Coeval measurements of CH 4 concentrations in bubbles trapped in the ice allow dating when correlated to the longer Antarctic ice cores. Also, marine cores offshore contain LIG sediment sourced from Greenland (Colville et al., 2011;Hatfield et al., 2016).
Consistent with the GRIP and GISP2 ice-core reconstructions (Johnsen et al., 2001;Raynaud et al., 1997;Yau et al., 2016), the summit location in central Greenland remains glaciated through the LIG in our simulation (Figure 6). Ice also remains at the NGRIP ice core site located northwest of GISP and GRIP and at the deposition site of the NEEM ice core. Drilling of the NGRIP ice core recovered an undisturbed climate record from the LIG, but only to 123 ka at bedrock (NGRIP, 2004). The NEEM ice core (NEEM, 2013), although folded with disturbances of the stratigraphic layering, has been interpreted as ice remaining during the LIG. Modeling using a 3D ice-sheet flow model provides an estimate of the deposition site of NEEM ice, located about 205 km upstream of the coring site at 128 ka (NEEM, 2013;Yau et al., 2016). Data from bedrock under the GISP2 ice core and offshore sediment records indicate a stable ice sheet in East Greenland for most of the past 1 million years (Bierman et al., 2016;Schaefer et al., 2016).
Our modeled GrIS retreats from the site of the Camp Century ice core by ∼125 ka. Silty layers sitting just above bedrock were retrieved by the drilling team (Hansen & Langway, 1966). Johnsen et al. (2001) speculated that the ice retrieved at the bottom of this core may have been deposited during the regrowth of the GrIS at the Camp Century site, having melted away during the early LIG. This suggestion has been reinforced by thermoluminescence (TL) dating of sediments at the bottom of the core (Christ et al., 2021) which indicates ice-free conditions sometime prior to 45 ka. The southern dome of the GrIS in our simulation reaches its minimum extent at 121.9 ka, before regrowing and reattaching to the northern dome at 120.5 ka. Ice remains along the high elevations in southeast Greenland. This behavior agrees with interpretations of sediment cores drilled at ODP 646 and MD99-2227 from the Eirik Drift, off southern Greenland, which suggests that the southern GrIS was present during the LIG, though smaller than in the Holocene. The isotopic composition of retrieved sediments confirms that the ice sheet eroded all three Proterozoic terranes of southern Greenland throughout the LIG (Colville et al., 2011;Hatfield et al., 2016). Ice-rafted debris (IRD) indicates significant melting and calving in the early LIG, with much less IRD after ∼125 ka, though with melting throughout the LIG (Irvalı et al., 2020;Rohling et al., 2019). Pollen concentrations preserved in these sediment cores are much higher than during the Holocene, with abundant alder and boreal ferns, indicative of proximal source and shrub tundra and boreal forests over southern Greenland during the LIG (de Vernal & Hillaire-Marcel, 2008). Ice is simulated to remain (barely) at the Dye-3 ice core location at 122 ka, with much-reduced thickness. Similar to Camp Century, Johnsen et al. (2001) proposed that the silty ice just above bedrock in Dye-3 could be from the regrowth of the ice sheet in the latest stages of the LIG. This superimposed ice is typical of the early growth stages of an ice sheet (Koerner, 1989). Alternatively, ice may have flowed into the region from central Greenland or from a surviving but isolated southern dome (Lhomme et al., 2005). DNA analysis of the Dye-3 silty ice concluded that the basal ice at this location predates the LIG, though assumptions and uncertainties of the age estimates cannot rule out the possibility of a LIG age for the Dye-3 basal ice (Willerslev et al., 2007).

Sea-Level Rise Contribution From Greenland
In the transient simulation with changing vegetation distribution (BIOME4 veg), the peak sea-level contribution from the GrIS is 3.8 m at 121.9 ka. Since our initial ice sheet at 127 ka (about 8.2 m SLE) is larger than the present-day ice sheet (about 7.4 m. SLE; Morlighem et al., 2017), this contribution can also be interpreted as 3.0 m SLE relative to present day. The latter value falls within the range of earlier estimates of GrIS contribution during the LIG. Dutton et al. (2015) summarized several previous studies, with a general consensus of 0.6-3.5 m contribution from the GrIS. However, some subsequent work, based on ice-sheet modeling and climate reconstruction from Greenland Summit ice cores, suggests a higher contribution of 5.1 m (±1 m) (Yau et al., 2016). Earlier modeling work also suggests a plausible contribution of 4-5.5 m from the GrIS during the LIG (Cuffey & Marshall, 2000).
At its minimum extent at 121.9 ka, the northern and southwestern parts of the GrIS are deglaciated. Past modeling studies show significant differences in the extent of the northern retreat (Haywood et al., 2019). However, unlike many of these previous studies, GrIS evolution in CESM2-CISM2 is not constrained by Greenland ice core data (see Section 3.3). In particular, an assumption that ice remained at the Camp Century location in northwest Greenland (shown in Figure 6) throughout the LIG is not made (or required) here. Sedimentological evidence from the Eirik Drift off southwest Greenland indicates that the southern GrIS was present during the LIG, though smaller than in the Holocene, with a sea-level contribution of 1.6-2.2 m SLE (Colville et al., 2011).
The timing of the peak sea-level contribution remains debated in the literature. For instance, the transient modeling simulations summarized in the Intergovernmental Panel on Climate Change Fifth Assessment Report (IPCC AR5) suggest a peak contribution (minimum ice volume) between 124.0 and 121.0 ka (Masson-Delmotte et al., 2013), consistent with the more recent work of Goelzer et al. (2016) and our simulation. Some reconstructions suggest a double-peak sea-level high stand during the LIG (Rohling et al., 2019), while others do not (Barlow et al., 2018). Another area of debate is the relative timing of the contributions from Greenland and Antarctic Ice Sheets (Rohling et al., 2019).
The results presented here constitute a new plausible estimate of sea-level rise during the LIG from the GrIS. While this is a significant advance in coupled climate-ice-sheet modeling, a few limitations should be noted. First, CESM2 does not dynamically evolve Antarctica. Including both dynamic ice sheets coupled with the global climate may reveal additional important feedbacks in the Earth system (Fyke et al., 2018). Ocean-driven forcing of the GrIS is also not included (such as sub-shelf melting of fjord-terminating outlet glaciers). However, ocean salinity changes transiently to reflect net freshwater fluxes between different Earth-system reservoirs as the ice sheet retreats and advances. Furthermore, our simulation begins from an ice sheet spun up under pre-industrial (c. 1850) climate to a configuration slightly larger (8.2 m SLE) than the present day. The GrIS undergoes an initial shock for ∼500 years as it adjusts to the 127 ka climate. A sensitivity simulation, initialized from a smaller ice sheet (7.5 m SLE, close to the modern value of 7.4 m) after adjusting to the shock, shows a similar rate of ice loss (∼1 mm yr −1 SLE) beginning from 127 ka. This simulation, if continued to 122 ka or beyond, would likely yield a smaller minimum ice sheet than the main simulation. Accounting for local and non-local isostatic adjustments prior to 127 ka and longer transient simulations beginning from the penultimate glaciation, could further illuminate ice-sheet conditions at 127 ka and throughout the LIG.

Sensitivity to Vegetation Changes
The transient vegetation changes (Figures 2 and S6 and S7 in Supporting Information S1) upstream of the GrIS have a significant influence on its subsequent evolution, as shown by the contrasting ice-sheet behavior in simulations with (BIOME4 veg) and without (no-Anthro veg) changing vegetation ( Figure 5). Figures 7 and S10 in Supporting Information S1 show boreal summer (June-August) mean surface temperature differences between the no-Anthro veg and BIOME4 veg simulations. Over high-latitude NH land areas, temperature differences demonstrate the effect of different PFTs in the two simulations. Areas of large negative anomalies from 125 to 121 ka, such as in northern Asia and North America, indicate where the no-Anthro veg vegetation consists of boreal broadleaf deciduous shrub and C3 Arctic grass PFTs, while the BIOME4 veg vegetation is composed of boreal forest tree PFTs. Vegetation types with lower albedo (e.g., forest) and higher evapotranspiration rates can enhance Arctic warming (Bartlein et al., 2015;Schurgers et al., 2007;Swann et al., 2010). Similarly, large regions of positive anomalies over North America and Asia, particularly from 127 to 123 ka, show areas where the BIOME4 vegetation is composed of temperate grassland PFTs, producing cooler surface temperatures (Gallimore et al., 2005) than in the no-Anthro veg simulation, where the vegetation in these regions is composed of tree PFTs (Figures 7, S8-S11, and S15 in Supporting Information S1). These contrasting impacts on near-surface air temperature are particularly evident in spring and early summer (April through June). Although our approach of changing the vegetation distribution every 500 years is not equivalent to a fully coupled dynamic vegetation model, it is a step forward compared to paleoclimate simulations that hold vegetation distribution constant. Our  results show substantially different ice-sheet behavior with and without evolving vegetation, a point previously made by , and emphasize the need for improved representation of vegetation in land models (Lawrence et al., 2019).

Reversibility of Ice-Sheet Collapse
Our results suggest that even with a significant retreat of the GrIS, regrowth is possible when the climate changes to cooler summer conditions. For example, the simulation using constant pre-industrial vegetation without anthropogenic influence (no-Anthro veg) loses 16% of the initial ice-sheet area and 17% of the initial ice-sheet volume at its minimum at 124.1 ka, then regrows to 97% of its initial area and 91% initial volume by 121 ka ( Figure 5). In the simulation with time-evolving BIOME4 vegetation (BIOME4 veg), the ice sheet reaches a smaller minimum extent at 121.9 ka, having lost 45% of the initial ice area and 46% of the initial volume. This retreat is followed by readvance, with the ice sheet recovering to 64% of the initial area and 59% of the initial volume by 119 ka (Figure 5). Ridley et al. (2010) found thresholds of irreversible ice-sheet collapse after a loss of 80%-90% of the ice-sheet volume. Gregory et al. (2020) simulated a threshold of ∼4 m SLE for regrowth of the ice sheet under 20th-century conditions. Our LIG simulations do not reach these extreme thresholds of loss before regrowth begins.

Conclusions
In this study, we have presented results from a transient, coupled global climate and ice-sheet model simulation using the CMIP6 models CESM2 and CISM2, representing the evolution and response of the climate and GrIS during the Last Interglacial period from 127 to 119 ka. We accelerate the ice-sheet model and orbital forcing by 5 times relative to the atmosphere, land, sea ice, and ocean model components. Vegetation distribution is updated every 500 years using a novel approach of simulating biome distributions based on the evolving climate.
Our simulation suggests a significant retreat of the ice sheet, with a peak sea-level contribution at 121.9 ka of 3.8 m from the initial model state at 127 ka (or 3.0 m sea-level rise relative to the present-day GrIS). In contrast, a companion simulation with a prescribed pre-industrial vegetation distribution shows much less retreat of the ice sheet and reaches a minimum ice volume at 124.1 ka with a peak contribution of 1.4 m SLE from the initial state (or 0.6 m sea-level rise relative to the present-day GrIS). This difference in behavior shows the importance of vegetation-climate-ice sheet interactions and feedbacks for polar ice-sheet retreat.
We present this modeling effort as an advance using a coupled, high-complexity global Earth system model (CESM2) and higher-order ice-sheet model (CISM2), simulating the LIG GrIS with the same atmospheric and ocean model resolutions as used for CMIP6 future projections. The model output contains a wealth of data that can be used to explore and analyze different aspects of the simulated LIG climate. Our simulations assume an initial state for the GrIS at 127 ka, and our sensitivity simulation indicates how using a different initial state may affect ice sheet size and ice loss. Coupled simulations beginning earlier during the penultimate glaciation (i.e., before 130 ka) could help to better constrain the state of the ice sheet at 127 ka and its subsequent influence on sea level during the LIG. Follow-up simulations should strive to incorporate truly dynamic vegetation and a coupled, evolving Antarctic Ice Sheet.
As the Earth warms under anthropogenic climate change, the GrIS behavior presented here for the Last Interglacial can be used to inform projections of future ice-sheet retreat, illustrating a possible retreat pattern in a warmer climate. Our simulation also demonstrates how this ice sheet is able to reverse a retreat trajectory by regrowing over long time scales with a cooling climate.