The phase space of last glacial inception for the Northern Hemisphere from coupled ice and climate modelling

We present an ensemble of Last Glacial Inception (LGI) simulations for the Northern Hemisphere that largely captures inferred ice volume changes within proxy uncertainties. This ensemble was performed with LCice 1.0, a coupled ice sheet and climate model, varying parameters of both climate and ice sheet components, as well as the coupling between them. Certain characteristics of the spatio-temporal pattern of ice growth and subsequent retreat in both North America (NA) and Eurasia (EA) are sensitive to parameter changes, especially with respect to regional rates of ice growth and retreat. We find 5 that the initial inception of ice over NA and EA is best characterized by the nucleation of ice at high latitude and high elevation sites. Subsequent spreading and merger along with large-scale conversion of snow fields dominate in different sectors. The latter plays an important role in the merging of eastern and western ice regions in NA. The inception peak ice volume in the ensemble occurs approximately at 111 ka and therefore lags the summer 60°N insolation minimum by more than 3 kyr. Ice volumes consistently peak earlier over EA than NA. The inception peak in North 10 America is characterized by a merged Laurentide and Cordilleran ice sheet, with Davis Strait covered in ice in ∼80% of simulations. Ice also bridges Greenland and Iceland in all runs by 114 ka and therefore blocks Denmark Strait. This latter feature would thereby divert the East Greenland Current and Denmark Strait overflow and thereby potentially have a significant impact on ocean circulation. The Eurasian ice sheet at its inception peak varies across ensemble runs between a continuous ice sheet to multiple smaller ice caps. 15 In both continents, the colder high latitudes (Ellsmere and Svalbard) tend to grow ice through the entire simulation (to 102 ka), while lower latitudes lose ice after ∼110 ka. We find temperature decreases over the initial phases of the inception lead to the expansion of NA ice sheet area, and that subsequent precipitation increases contribute to its thickening. EA ice sheet area also expands with decreasing temperatures, but sea ice limits any increases in precipitation, leading to an earlier retreat away from the EA maximum ice sheet volume. 20 We also examine the extent to which the capture of both LGI ice growth and retreat constrains the coupled ice/climate model sensitivity to changing atmospheric pCO2. For a standard transient climate response experiment (1% increase in pCO2 until doubled), warming ranges between 0.6-2.0°C for our initial set of 500 simulations without LGI constraint. The warming is reduced to 0.7-1.4°C for the 55 member ensemble that captures both LGI ice growth and retreat. This therefore underlines the potential value of fully coupled ice/climate modelling of last glacial inception to constrain future climate change. 25 2 https://doi.org/10.5194/cp-2020-1 Preprint. Discussion started: 27 January 2020 c © Author(s) 2020. CC BY 4.0 License.

1 Introduction optimal balance in sensitivity tests between efficiency and proximity to shorter coupling time step solutions (Bahadory and Tarasov, 2018).
Fields passed from the ice sheet to the atmosphere include ice mask and surface elevation, the latter via one of the three included schemes (simple, envelope, and silhouette, the choice of which is under ensemble parameter control). The atmosphere to 155 ice coupling includes the monthly mean and standard deviation temperature and monthly mean precipitation, evaporation, wind direction and magnitude, and vertical temperature lapse-rate. LCice 1.0 uses an innovative scheme to downscale precipitation to the ice model grid that accounts for orographic forcing on the GSM grid resolution topography. Temperature downscaling uses the evolving vertical surface temperature gradient field of LOVECLIM. The coupler also includes a simple radiative cloud parameterization to compensate for the present-day prescribed radiative cloud cover of LOVECLIM. 160 In ice sheet-ocean interactions, the GSM determines the runoff routing, and passes freshwater fluxes to the ocean model, while the ocean model provides the GSM with vertical temperature profiles, required to calculate sub-shelf melt. Details of each component of the coupling and their influence are described in Bahadory and Tarasov (2018).
Given model limitations, there is no one best run in the ensemble. Instead, different runs have different features, each of which will likely have different patterns of misfits against inferred proxy records. In the following results, we crudely interpret 165 feature frequency in the ensemble to be a partial metric of feature likelihood, though this is far from a rigorous probabilistic analysis.

Results
The LCice 1.0 ensemble reproduces the reconstructed pattern of rapid ice sheet volume growth and retreat during the LGI in 55 of the 500 runs. The total Northern Hemisphere ice volume averaged over the ensemble of 55 runs is plotted in figure 1. No 170 single ensemble parameter determines which runs meet the filter condition (not shown).
The maximum ice volume achieved by the LCice 1.0 ensemble during inception is lower than that inferred by Lisiecki and Raymo (2005), but within the collective uncertainty of the three reconstructions presented here (Waelbroeck et al., 2002;Siddall et al., 2003;Lisiecki and Raymo, 2005). The ensemble mean maximum ice volume is about 5 m in sea level equivalent (SLE) short of the Red Sea record (Siddall et al., 2003, dark purple in figure 1). This under-estimation is likely due in part to the 175 absence of any contribution from the Antarctic ice sheet (and perhaps Patagonian and Tibetan ice caps). This under-estimation is also consistent with the fact that the simulated ice sheet volumes never reach the peak rate of ice growth indicated by any of the sea level reconstructions.
The timing of when the LCice 1.0 simulations achieve their maximum inception ice sheet volume is bounded by the three proxy-based reconstructions shown in figure 1. All but the Greenland ice sheets reach their maximum LGI ice volumes at 180 least 3 kyr after the 60°N summer insolation minimum (orange line in figure 1). The earliest retreat occurs in the Red Sea reconstruction. This reconstruction suggests a faster decrease in pre-stadial sea level compared to that of the other three records, and its timing of the sea level minimum and subsequent sea level rise is slightly advanced of the LCice ensemble mean. The LCice maximum ice sheet volume occurs approximately midway between the timing of minimum insolation at 60°N and A second test of the representativeness of these simulations for the LGI is made between temperature changes from a glaciological inversion of the GRIP ice core δ 18 O record (Dansgaard et al., 1993;Tarasov and Peltier, 2003) and annual-mean temperatures calculated from the model grid cell containing its location. The ensemble mean 2m temperature anomaly relative to 119 ka follows the general trend of GRIP reconstructed temperatures in figure 2 until ∼ 112 ka. Individual runs have higher 190 decadal to centennial scale variance than that of the GRIP record. However, the large millennial scale variability of the GRIP record inversion is not captured by the simulations. The ensemble-mean annual temperatures from the GRIP site subsequently diverge from reconstructed temperatures after approximately 111 ka. At this time, simulated temperatures increase at the 8 https://doi.org/10.5194/cp-2020-1 Preprint. Discussion started: 27 January 2020 c Author(s) 2020. CC BY 4.0 License. Figure 2. Annual mean 2 m temperature anomaly relative to 119 ka for the GRIP ice-core (green) (Dansgaard et al., 1993;Tarasov and Peltier, 2003), ensemble mean (thick black), and three individual runs (gray lines). The orange line depicts the timing of insolation changes at 60°N. GRIP site following insolation changes, whereas there is no evidence of a similar increase in the GRIP record temperature inversion. Instead, reconstructed GRIP temperatures exhibit multi-millennial timescale oscillations around stable, stadial (cold 195 state) temperatures. It is unclear what mechanism would sustain stadial temperatures over central Greenland under increasing insolation, especially since the simulations consistently predict that strong warming should result. It may be this discrepancy reflects in part a lack of accounting for at least two standard sources of uncertainty in water isotope to temperature inversions: changes in the moisture source region and changes in the seasonal distribution of precipitation.

Glacial inception phase-space 200
Having established that LCice 1.0 is able to capture both the ice sheet growth and retreat phases of the LGI, we explore the pattern(s) of the ice growth and retreat across ensemble members. We start by analyzing the spatial patterns of EA and NA ice sheets at two diagnostic time intervals: first, the early stage of ice build up, and second, during the peak of the inception around 9 https://doi.org/10.5194/cp-2020-1 Preprint. Discussion started: 27 January 2020 c Author(s) 2020. CC BY 4.0 License. 112 ka. Next, we explore the consistency of ice and climate evolution between these two intervals and during the subsequent retreat phase. Despite having different start times (due to different calendar start years between 122 ka and 119 ka and spinup lengths varying between 3 to 5 kyr), all simulations start growing ice in the first 100 years of simulation (figure 3.a). Therefore, we analyze the spatial patterns of the first appearance of ice in the first 1000 years of simulation, rather than aggregating simulations according to a common calendar year.

210
In NA, all runs have extensive glaciation over Ellesmere and eastern Devon Island after 100 years of transient simulation ( figure 3). Subsequently, ice starts to spread through the Arctic archipelago and Baffin Bay sector of Baffin Island. This is in agreement with past suggestions that the first ice nucleation in NA occurs over the Canadian Archipelago with further growth, merger, and then expansion to southern and western regions (Weertman, 1964). This result is also consistent with the ongoing presence of extensive glaciers and small ice caps in this region.

215
By 1000 years, more than 20% of runs have extensive ice over the Pacific Cordillera down to 48°N. Northwestern Alaska remains ice free for the first 1000 years in all runs as does the non-Cordilleran sector of NA below 61°N.
To get a more detailed sense of what glacial inception might look like, it is worth examining ice sheet evolution for one of the best fitting runs (to sea level proxies). By 119 ka, most of NA above 65°N has ice cover, though much of it with surface elevation less than 500 masl (figure 4). The Canadian Cordilleran at this time has near complete ice coverage with all surface 220 elevation above 1000 masl.
10 https://doi.org/10.5194/cp-2020-1 Preprint. Discussion started: 27 January 2020 c Author(s) 2020. CC BY 4.0 License. To capture the maximum in ice volume for EA and NA during the LGI, we consider time slices for 114 ka, 112 ka, 110 ka and 108 ka in figures 5 and 6. We aggregate our simulation results according to their boundary condition years rather than their simulation years.
At 114 ka, the Cordilleran is completely ice covered in all runs down to approximately 45°N. Central NA ice extends to approximately 55°N until a sharp northward turn of the southern ice margin over James Bay extending to the east (figure 5).

235
The Greenland and Iceland ice sheets are bridged by ice across Denmark Strait in all runs by 114 ka (with most runs having grounded ice right across the Strait). Also, Alaska is almost fully ice-covered in all of the simulations, while Labrador and eastern NA remains ice-free, likely due to warm model biases in this region.
The main differences in peak LGI NA ice extent between ensemble members occur: at the northwestern Alaskan ice margin (40% of ensemble runs cover Bering Strait at 114 ka), at the southern margin, and over Davis Strait. For the latter, approximately 240 80% of simulations create an ice bridge connecting the Laurentide and Greenland ice sheets across the Strait. This ice bridge generally starts out from a merger of opposing ice shelves. For some (but not all) ensemble runs, it can also ground right across the Strait and therefore isolate Baffin Bay from the Labrador Sea.
After the stadial peak in NA ice volume, the main variation between ensemble members appears in the rate of ice retreat. Initially, while the south-eastern ice margin rapidly retreats to higher latitudes in simulations with smaller ice sheets, simulations 245 with larger ice sheets show little change in ice extent. This difference in behaviour leads to the largest difference in ice extent over Hudson Bay at 110 ka, when the entire area is covered by approximately 20% of the simulations and 30% are ice-free in this region. By 108 ka, the Laurentide and Cordilleran ice sheets are separated in only 10% of the simulations, fewer than 20% of runs simulate a connected Greenland-Iceland ice sheet, and the ice bridge across Davis Strait remains in fewer than 10% of runs.

250
A key feature from the sample best run snapshots (figure 4) is the continuous slow thickening of Ellesmere Island ice right through to 105 ka. Thus, limited snow accumulation appears to be the major controlling climate factor for this region during LGI. The ice dome north of Hudson Bay also only attains it maximum elevation at 107 ka.
Similar to the early phases of the inception, ice extent over EA is more variable between ensemble members around the stadial peak time (116 ka to 112 ka) compared to NA (figure 6). The maximum area of 100% cross-ensemble continental 255 ice cover occurs at 116 ka, with a significant reduction by 114 ka. Fewer than 10% of runs increase their southern ice extent through to 112 ka. Scotland exhibits some ice cover in the majority of runs, but the North Sea remains ice-free.

North American ice sheet
In all NA regions in figure 8 except NA El , ice volume increases to a maximum sometime between 112 ka and 109 ka and then decreases. In NA El , the coldest region of NA, ice volume increases throughout the LGI in most simulations.
Generally, the ice sheet growth phase for each sector is more consistent between runs than its retreat phase. In sector NA Bf (figure 8b), ∼10% of simulations lose between 1 and 1.5 m SLE of ice between 112 and 107 ka and maintain a constant ice 270 volume afterwards. The rest of the runs show a range of behaviours, from almost no ice loss to 80% loss. In contrast, in NA Qb , the most southern and warmest sector, the maximum ice volume varies between almost zero to more than 1 m SLE, and no simulation sustains ice cover through to 102 ka. The NA Rc region spans the widest range of latitudes, but it also contains some of the highest-elevation sites of NA. It shows both strong ice growth and a wide range of ice loss scenarios over the LGI.
Notably, ice develops over western NA (NA Rc ) at the same time as it is growing in the east.

275
One pattern that emerges most strongly in NA Qb is that the runs with larger ice sheets tend to have delayed peak times. This is consistent with the observation in the previous section that runs with the largest NA ice sheet extent retreat more slowly than those with smaller ice sheets. In assessing the contributions of sea ice, the AMOC and the jet stream, summer sea ice has the strongest correlations with temperature and precipitation changes in EA. Late winter sea ice area shows no consistent pattern of change over this time period and is not related to ice sheet volumes in either NA or EA (see supplemental figure 1). However, its summer extent varies in correspondence with Northern Hemisphere temperatures: it peaks prior to the minimum in insolation at 60°N, remains 345 extended, and then decreases. The onset of major sea ice retreat at approximately 113 ka is in phase with a rapid acceleration of both NA and EA summer warming and annual-mean precipitation. Deciphering the causal relationships of this phasing requires future sensitivity experiments. However, one can infer that sea ice likely has a positive feedback role for both precipitation and temperature at this time.
Neither the AMOC (supplemental figure 3) nor the wintertime jet stream exhibit any clear consistent changes that coincide 350 with temperature and precipitation or ice sheet changes. In 80% of the runs, the AMOC gradually increases during the glacial inception to a maximum of 22 Sv around 108 ka. After this, it decreases once more to its initial values of 16 to 18Sv. In the remaining 20% of runs, the AMOC oscillates between two values.

Caveat about marine sectors
Ice sheet growth in marine sectors is found to be highly sensitive to the treatment of sub-shelf melt, even at high latitudes. This 370 is particularly evident in figure 6, where marine ice sheet margins are at times extended straight lines. These lines match the boundaries for different ocean temperature sectors in LCice, which propagates the vertical temperature profile from assigned upstream diagnostic sites to whole downstream ocean sectors for computing submarine ice melt (Bahadory and Tarasov, 2018).

Brief comparison to past geological inferences
To our knowledge, there is no community-based geologically-inferred MIS 5 ice margin reconstruction for NA. Aside from the issue of Alaska (and certain adjacent parts of the Yukon), our results are, within (large) age uncertainties, consistent with the till stratigraphy presented in Clark et al. (1993).
Eurasia also lacks a clear geologically-inferred LGI stadial extent. However, the geologically-inferred Early Weichselian 410 (MIS 5)ice extent maximum of Svendsen et al. (2004, nominal 90 ka in) generally encloses (and for much of the southern margin largely tracks) the 50% ensemble distribution (figure 6). The main regional exceptions are more extensive ice on the western coast of Svalbard and extensive marine ice on the western Norwegian coast. We leave it to members of the geological community to execute more detailed and up-to-date comparisons with our ensemble chronologies.

4.5
Is there a single very likely spatio-temporal pattern of LGI ice sheet evolution?

415
To partially characterize the range of the spatio-temporal patterns of ice sheet evolution in our ensemble, we consider the intersectorial relationships of maximum ice volume for each ensemble run (figure 13). The absence of correlation in maximum ice volumes for different sectors will indicate that that there are multiple temporal patterns of ice development in these regions.
For NA, the northern Arctic (NA El ) sector maximum ice volume has no obvious correlation with that of other sectors. This is consistent with the continual growth of ice throughout the simulations in this region. All other NA sectors display relatively 420 strong correlations aside from a threshold response for NA Qb relative to the Pacific Cordillera (NA Rc ).
For EA, again the most northern and continuously growing sector (EA Sv ) has relatively no correlation in maximum ice volume with other sectors (figure 13). The relatively northern and largely marine eastern sector (EA Kr ) has a strong correlation with the two Fennoscandian sectors for the 5 runs with maximum (EA Kr ) greater than 2.5 mSLE. For the other runs, the correlation is much weaker and with a much lower mean slope, perhaps indicative of a threshold in ocean temperatures 425 controlling subshelf melt and enabling ice calving.
There are no strong correlations between NA and EA regions (c.f. supplemental figure 4). There is moderate correlation between the Baffin Island sector NA Bf and western Fennoscandia EA W F , perhaps reflecting ocean circulation connections between Baffin Bay and the GIN Seas. More limited correlation exist between NA Bf and eastern Fennoscandia (EA EF and and between the western Cordillera (EA Rc ) and western Fennoscandia EA W F ). The only other possible relation of note is the 430 absence of large maximum ice volumes for the eastern Kara Sea region (EA Kr ) when ice volumes are near maximal for all NA sectors south of Ellesmere (with only 5 runs for this case, the relationship is tentative).
The only clear indication of a bifurcation in regional temporal evolution is the presence of both early and late timing of maximum regional ice volume for Fennoscandia (EA EF and EA W F ) for a range of regional maximum ice volumes ( figure   8). The extent to which possible associated bifurcations in sea ice extent and stationary atmospheric waves (described in the 435 results section) may play a role in this must await future analysis.