Glacial inception through rapid ice area increase driven by albedo and vegetation feedbacks

Abstract

We present transient simulations of the last glacial inception using the Earth system model CLIMBER-X with dynamic vegetation, interactive ice sheets, and visco-elastic solid Earth responses. The simulations are initialized at the middle of the Eemian interglacial (125 kiloyears before present, $\mathrm{ka}$) and run until 100 $\mathrm{ka}$, driven by prescribed changes in Earth's orbital parameters and greenhouse gas concentrations from ice core data.

CLIMBER-X simulates a rapid increase in Northern Hemisphere ice sheet area through MIS5d, with ice sheets expanding over northern North America and Scandinavia, in broad agreement with proxy reconstructions. While most of the increase in ice sheet area occurs over a relatively short period between 119 and 117 $\mathrm{ka}$, the larger part of the increase in ice volume occurs afterwards with an almost constant ice sheet extent.

We show that the vegetation feedback plays a fundamental role in controlling the ice sheet expansion during the last glacial inception. In particular, with prescribed present-day vegetation the model simulates a global sea level drop of only $\sim$ 20 $\mathrm{m}$, compared with the $\sim$ 35 $\mathrm{m}$ decrease in sea level with dynamic vegetation response. The ice sheet and carbon cycle feedbacks play only a minor role during the ice sheet expansion phase prior to $\sim$ 115 $\mathrm{ka}$ but are important in limiting the deglaciation during the following phase characterized by increasing summer insolation.

The model results are sensitive to climate model biases and to the parameterization of snow albedo, while they show only a weak dependence on changes in the ice sheet model resolution and the acceleration factor used to speed up the climate component.

Overall, our simulations confirm and refine previous results showing that climate–vegetation–cryosphere feedbacks play a fundamental role in the transition from interglacial to glacial states characterizing Quaternary glacial cycles.

1 Introduction

Modelling of glacial inceptions, i.e. the relatively fast (compared to the duration of glacial cycles) transitions of the Earth system from an interglacial (no significant ice sheets over the northern continents, except for Greenland) to a glacial state, caused by lowering of the boreal summer insolation, is a crucial test for the Milankovitch theory of glacial cycles . Milutin Milankovitch was the first who used mathematical modelling to test the hypothesis that the lowering of boreal summer insolation associated with precession and obliquity cycles is sufficient to cause growth of northern ice sheets. In his “Cannon of insolation”, using a simple energy balance model and taking into account the ice–albedo feedback, Milankovitch estimated that a typical decrease in summer insolation would cause a lowering of the snow line by about 1 $\mathrm{km}$ and would cause a southward expansion of ice cover up to 55° N.

A set of transient and equilibrium experiments with the Earth system models of intermediate complexity (EMIC) CLIMBER-2 demonstrated that glacial inception could be interpreted as a bifurcation transition in the Earth system in the phase space of orbital forcing, defined as the maximum summer insolation at 65° N. This idea was explicit already in classical work by but was not tested before with geographically explicit climate–ice sheet models. traced both equilibrium states in summer insolation-phase space and found a wide range of insolations for which multistability of the climate cryosphere system exists. They also show that there is a second important dimension of the phase space, the atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration. This finding was confirmed using a completely different model in . A detailed analysis of the critical relationship between summer insolation and ${\mathrm{CO}}_{\mathrm{2}}$ was performed with the CLIMBER-2 model in , where a logarithmic relationship between critical ${\mathrm{CO}}_{\mathrm{2}}$ concentration and insolation was found. This result has been recently corroborated by the new CLIMBER-X model .

Since the most recent glacial inception, which began between 120 and 115 $\mathrm{ka}$, is, compared to previous glacial inceptions, relatively well covered by paleoclimate data and occurred during one of the deepest (the second lowest during the late Quaternary) boreal summer insolation minima, it is often used for testing Milankovitch theory. Until recently, complex climate models based on general circulation models (GCMs) did not include ice sheet components and thus could not simulate the evolution of ice sheets, which can be directly compared with existing paleoclimate data. So-called time slice simulations with fixed orbital parameters and greenhouse gas (GHG) concentrations corresponding to the insolation minimum (115 or 116 $\mathrm{ka}$) have been traditionally performed instead. In such simulations, an establishment of perennial snow cover over at least several model grid cells in northern North America and/or Eurasia was considered to be the criterion for a successful simulation of glacial inception. A number of such simulations have been performed during the past 4 decades e.g.. Similar time slice experiments using climate model output to force a standalone ice sheet model have also been performed e.g.. Overall, the results of these simulations in terms of simulated (or not simulated) areas of perennial snow cover and their locations have been diverse. In many simulations, perennial snow cover did not appear even under such extreme orbital forcing or appeared in locations not supported by paleoclimate reconstructions. This is not surprising, since climate models used for paleoclimate studies usually have a coarse spatial resolution, which has been shown to have a significant impact on simulations of glacial inception e.g.. No less important is the fact that climate models tend to have climate biases often comparable in magnitude with the simulated climate response to orbital forcing. Some studies e.g. demonstrated a significant influence of climate biases on the simulation of glacial inception. It was also proposed that the failures to simulate glacial inception in climate models can be related to the omission of some positive feedbacks, such as the climate–vegetation feedback . Several studies demonstrated that a weakening of the Atlantic Meridional Overturning Circulation (AMOC) during glacial inception e.g. could potentially affect ice growth in North America or Scandinavia . The role of GHGs was also investigated e.g., but since the last glacial inception began when concentrations of ${\mathrm{CO}}_{\mathrm{2}}$ and other GHGs were close to pre-industrial values, it is clear that climate–carbon cycle feedbacks can contribute to the rate of ice sheet growth but not to the initial ice growth.

So far, only a few attempts have been made to simulate the last glacial inception with ice sheet models asynchronously coupled to atmospheric GCMs but have generally failed to grow enough ice to match sea level reconstructions, even when using permanent coldest orbit (116 $\mathrm{ka}$). Many more simulations of the last glacial inception and even full glacial cycles have been performed with different types of simplified, geographically explicit climate–ice coupled models. In particular, 2D energy balance climate models were coupled with 2D  and later 3D thermomechanical ice sheet models to simulate the last glacial cycle. More recently, different EMICs were applied in transient simulations of the last glacial inception and the entire last glacial cycle (or even several glacial cycles) . While EMICs are much more computationally efficient than GCMs, which allows for the performance of many long-term (orbital timescale) simulations, they usually have even larger climate biases than GCMs. This often causes significant misplacement of simulated ice sheets during glacial inception compared to paleoclimate reconstructions , and it was shown that the correction of temperature biases significantly improves the agreement between simulated and reconstructed ice sheets evolution during glacial inception .

EMICs are computationally fast enough to allow large ensemble of simulations to be performed, exploring the dependence of simulated glacial inception on model parameters e.g. and its sensitivity to different factors e.g. and found that the transition from interglacial to glacial state is characterized by two timescales associated with different mechanisms. They have shown that the initial growth of ice sheets begins nearly simultaneously at several nucleation centres in high latitudes and high elevations of the Northern Hemisphere, and the lateral spreading occurs only through the slow motion of ice sheets, resulting in rather low initial rate of ice volume increase. However, at a certain time corresponding to a sufficiently low summer insolation, an abrupt expansion of the ice sheet area occurs through the positive albedo feedback. This leads to a rapid increase in accumulation area and causes a significant acceleration of the ice volume growth rate. Over the past 20 years since the publication of , the phenomenon of rapid expansion of ice sheet accumulation area through albedo feedback mechanisms was either not found or not properly discussed by other authors, although it is likely that a similar mechanism operated in simulations described in . This is why it is essential to revisit this topic as new and improved modelling tools become available.

In this study we build on our experience in modelling glacial cycles with the CLIMBER-2 model and employ the newly developed Earth system model CLIMBER-X with interactive ice sheets, visco-elastic solid Earth response and dynamic vegetation to simulate the last glacial inception from 125 to 100 $\mathrm{ka}$. CLIMBER-X operates at a relatively high horizontal resolution of 5° for the climate model and 32 $\mathrm{km}$ for the ice sheet model and includes a rather detailed representation of the processes which have been found to be potentially important to simulate the interglacial–glacial transition, making it a powerful tool to explore mechanisms and feedbacks at play during the last glacial inception.

The paper first describes the model and in particular the ice sheet surface mass balance computation and the ice sheet coupling strategy (Sect. ), with more details provided in the Appendix. The experimental setup is then presented (Sect. ) followed by a description of the results of the model simulations (Sect. ) and finally discussions and conclusions.

2 Model description and setup

To simulate the last glacial inception we use the Earth system model CLIMBER-X , which is a major step forward in modelling capabilities compared to its predecessor CLIMBER-2, with a much higher spatial resolution both in climate and ice sheets and substantial improvements in all components: atmosphere, ocean, sea ice, vegetation, ice sheets, and surface mass balance scheme. The CLIMBER-X climate model component includes a semi-empirical, statistical dynamic atmosphere model, a 3D frictional geostrophic ocean model, a dynamic thermodynamic sea ice model, and a land surface and dynamic vegetation model, all with a horizontal resolution of 5° $\times$ 5°.

Additionally, the two ice sheet models SICOPOLIS and Yelmo are part of CLIMBER-X . In the present study we employ SICOPOLIS to model the Northern Hemisphere (NH) ice sheets since it has already been extensively applied to NH glaciation studies, particularly in the framework of the CLIMBER-2 model e.g., the predecessor of CLIMBER-X. In CLIMBER-X we make use of one of the latest versions of SICOPOLIS (v5.3; ), which is based on the shallow-ice approximation for grounded ice, the shallow-shelf approximation for floating ice and hybrid shallow-ice–shelfy-stream dynamics for ice streams . The enthalpy method of is used as thermodynamics solver. The ice sheet model is applied to the NH at a default horizontal resolution of 32 $\mathrm{km}$. The performance of the standalone SICOPOLIS model for present-day Greenland has been presented in . Since we concentrate on the NH, Antarctica is prescribed at its present-day state in this study based on the assumption that the Antarctic ice sheet contribution to global ice volume changes during the glacial inception is small. Different modelling studies confirmed that during this time period the Antarctic contribution to global sea level drop is only about 5 $\mathrm{m}$ e.g., about 10% of global sea level variations reconstructed for this period.

Figure 1

Schematic illustration of the coupling of the ice sheets to the other CLIMBER-X model components. The colours of the arrows indicate the coupling frequency, while the colours of the boxes represent the internal time step of the different model components. The abbreviations in the figure are as follows: SMB stands for surface mass balance, $T_{\text{ice}}$ stands for 10 $\mathrm{m}$ firn temperature, $Z_{\text{sur}}$ stands for surface elevation, $H_{\text{ice}}$ stands for ice thickness, $Z_{\text{bed}}$ stands for bedrock topography, $Z_{\text{base}}$ stands for ice sheet base elevation, BMB stands for basal mass balance, $T_{\text{ocn}}$ stands for seawater temperature, and $S_{\text{ocn}}$ stands for seawater salinity.

[Figure omitted. See PDF]

The coupling of the ice sheet model to the other CLIMBER-X components is schematically illustrated in Fig.  and further details of SICOPOLIS are presented in Appendix .

A physically based surface energy and mass balance scheme (SEMIX, Surface Energy and Mass balance Interface for CLIMBER-X) is used to interface the ice sheet model with the atmosphere. In the ice sheet modelling community it is still common to use the simple positive-degree-day (PDD) model e.g. for the computation of surface melt. However, this simplified approach computes melt only from near-surface air temperature, with potentially important implications on the simulated response of the surface mass balance to orbital forcing and climate change e.g.. To address some of the limitations of the PDD scheme, so-called insolation–temperature methods (ITMs) have been developed that explicitly account for absorbed surface solar radiation . Recently, energy balance schemes or land surface models have also been applied to compute the surface mass balance of ice sheets (e.g. ). Similarly to these models, SEMIX also explicitly resolves the surface energy balance and is largely based on the ideas in , but with notable modifications and improvements. The climate model fields needed by SEMIX are first bilinearly interpolated from the coarse-resolution climate model grid ($\mathrm{5}\mathit{°}\times \mathrm{5}\mathit{°}$) onto the high-resolution ice sheet model grid, where SEMIX operates. Subsequently, temperature, humidity, and radiation fields are elevation corrected. We do not account for the additional orographic effect on precipitation. Such an effect is crudely accounted for by the atmospheric component of CLIMBER-X, which has a much higher resolution compared to CLIMBER-2, where a parameterization of the orographically forced precipitation was included . Snow albedo depends on snow grain size and the concentration of dust in snow, while bare ice sheet albedo transitions between firn-like values, clean-ice values, and dirty-ice values depending on the age of ice and the time it has been exposed to the atmosphere. CLIMBER-X includes an interactive dust cycle , and the dust concentration in snow is computed from the simulated dust deposition flux. The surface energy balance equation is then solved following a similar approach as over land and sea ice as described in and , and the surface mass balance is finally derived accounting for snowfall, rainfall, evaporation or sublimation, melt, refreezing, and runoff. The performance of SEMIX for Greenland in the fully coupled CLIMBER-X configuration has already been presented in . In order to compute the annual surface mass balance, SEMIX is called every 10 years over a full year with a time step of 1 d. As long as the forcing is slow enough, a higher frequency of SEMIX calls is not needed to accurately capture the climatological evolution of the surface mass balance because CLIMBER-X does not resolve weather and internal inter-annual variability. SEMIX is described in detail in Appendix .

Figure 2

(a) CLIMBER-X summer (JJA) near-surface air temperature bias relative to ERA5 reanalysis for the time period 1981–2010. (b) Difference in simulated summer (JJA) near-surface air temperature between 116 $\mathrm{ka}$ and pre-industrial values resulting only from the different orbital parameters.

[Figure omitted. See PDF]

The surface mass balance of ice sheets is largely controlled by air temperature and precipitation and relatively small climate model biases can result in large changes in simulated ice sheets e.g.. In accordance with Milankovitch theory, ice sheets are particularly sensitive to summer temperatures, which control snow and ice melt. Thus, even small summer temperature biases, which are common even in state-of-the-art climate models e.g., can therefore preclude a realistic simulation of ice sheet growth and decay. In CLIMBER-X, present-day summer temperature biases over NH land are generally within a few degrees (Fig. a). However, in the area around Hudson Bay and the Labrador Peninsula, simulated summer temperatures are up to $\sim$ 5 $\mathrm{°}\mathrm{C}$ too warm. Biases of this order of magnitude are comparable to the temperature changes induced by changing orbital configuration (Fig. b) and can consequently strongly affect the simulated ice sheet distribution. Therefore, similarly to , in SEMIX we implemented a temperature bias correction over northern North America as shown in Fig. . The bias correction is applied throughout the year as a constant offset in the computation of the surface energy balance. The reasoning behind the choice of applying summer temperature bias correction throughout the year is guided by the fact that, following Milankovitch theory, ice sheet surface mass balance is largely determined by ablation during summer, which is highly sensitive to temperature. We have actually also tested the use of bias correction applying monthly mean temperature differences, but the model results indicated absolutely negligible differences compared to using the mean summer bias. We also assume that the bias is a persistent feature of the model also under different boundary conditions and therefore apply the same bias correction at all times. Since the focus of our study is on the last glacial inception and the boundary conditions at $\sim$ 120 $\mathrm{ka}$ (just before the expected onset of ice sheet growth) are similar to pre-industrial values in terms of GHG concentrations and orbital parameters, we do not expect significant differences in temperature biases compared to the present day. The assumption on stationarity of the temperature biases is therefore, at least during the initial ice growth phase, well justified. Furthermore, since air temperature is closely related to the downwelling long-wave radiation at the surface, which affects the surface energy balance, we additionally correct the downwelling long-wave radiation using a simple quadratic relation derived from ERA5 reanalysis data. For more details, see Appendix .

As soon as an ice sheet gets in contact with the ocean its mass balance is also influenced by melting and freezing at the base of the floating ice shelves. At present, basal melt is important mostly for the Antarctic ice sheet, which is in extensive contact with the ocean, and for the interactions of the Greenland ice streams with the fjords. However, little is known about the role played by basal melt for ice sheets in the NH during the last glacial cycle. The basal melt is strongly influenced by the ocean circulation on the continental shelf and by the turbulent plume dynamics generated at the ice–water interface by the density difference between the fresh meltwater and the surrounding saline seawater e.g.. These processes can not be properly represented at the typically applied resolutions of global ocean circulation models, and several simple parameterizations have been put forward to represent the basal melt process in ice sheet models based on large-scale ambient properties of seawater and on the ice shelf base topography . Simple box models have also been developed to describe the ocean circulation in the cavities below the ice shelves , but they are specifically designed to be applied to Antarctica, and an extension to the NH is not straightforward. In CLIMBER-X we have therefore implemented the simple and general basal melt parameterizations of and , both of which rely on the difference between the ambient ocean temperature simulated by the ocean model and the freezing point temperature of seawater at the ice shelf base. The linear model of has been used in the simulations presented in this work. A more general 2D plume model also including the Coriolis effect that is general enough to be applied to any ice shelf geometry has been recently developed and will be considered for a future implementation in CLIMBER-X. See Appendix  for more details on the basal melt schemes in CLIMBER-X.

The presence of ice sheets also affects the simulated freshwater fluxes into the ocean. Freshwater produced by melt occurring at the surface of the ice sheets is routed to the ocean following the steepest surface gradient, and the same runoff directions are used also to route the water originating from basal melt of grounded ice to the ocean. The freshwater flux from basal melt of floating ice and ice shelf calving is applied locally to the ocean model. Both liquid water runoff and calving fluxes are distributed uniformly over several coastal grid cells. For calving fluxes the latent heat needed to melt the ice is also accounted for in the ocean model. CLIMBER-X so far does not include an iceberg model to explicitly simulate the fate of ice originating from ice shelf calving.

The glacial isostatic adjustment (GIA), which controls regional sea level change and surface displacements, is computed by the viscoelastic solid Earth model VILMA . In VILMA, mass conservation between ice and ocean water is considered, and the loading effect of the ice and gravity-consistent redistribution of water in the ocean is solved applying the sea level equation including the rotational feedback . VILMA takes the ice sheet loading as input and computes the relative sea level (RSL) changes. The sea-level equation at the surface is solved on a Gauss–Legendre grid of 512 $\times$ 256 grid points (n128), which is consistent with the spectral resolution in spherical harmonics up to a Legendre degree or order of 170, corresponding to a wavelength of $\sim$ 120 $\mathrm{km}$. In this study we use the 3D viscosity structure from , which is based on a 3D tomography model of the upper mantle seismic velocity structure that is converted into viscosity variations considering further constraints by geodynamics and lithospheric structure. In CLIMBER-X, the solid Earth dynamics determines the bedrock elevation through simulated changes in relative sea level. The relative sea level computed by VILMA is substracted as an anomaly from the reference high-resolution ($\approx$ 10 $\mathrm{arcmin}$) bedrock elevation from . This is done in order to preserve the small-scale structure, which is important both for runoff routing and for surface mass balance. Surface elevation, land–sea mask, and runoff routing directions are then updated accordingly. The surface elevation is then aggregated to the coarse climate model resolution, and the land and ocean grid cell fractions on the climate model grid are then derived. The climate model in CLIMBER-X can deal with land–sea mask changes, a feature that is rather exceptional for contemporary Earth system models e.g.. VILMA is called every 10 years, implying that the ice sheet model gets an update of the bedrock elevation field every 10 years and that the topography, the land–sea mask and the runoff directions are also updated every 10 years. Experiments using a higher coupling frequency of 1 year showed negligible differences. In a recent study coupling VILMA with the PISM ice sheet model showed a coupling time step of 100 years as being sufficient.

3 Experimental design

We first performed a transient Holocene simulation from 10 $\mathrm{ka}$ to present day (2000 $\mathrm{CE}$) to make sure that no glacial inception is simulated at present. Since for this purpose we are mainly interested in the growth of ice sheets outside of Greenland, the Holocene simulation was performed with a prescribed present-day Greenland ice sheet. In this simulation, the model is driven by prescribed changes in the orbital configuration and atmospheric concentrations of ${\mathrm{CO}}_{\mathrm{2}}$, ${\mathrm{CH}}_{\mathrm{4}}$ and ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ from ice core data and historical observations . Note that NH summer insolation is continuously decreasing over the Holocene, although not enough to cause a glacial inception at present . At the end of the transient Holocene simulation ice cover is consistently simulated only over parts of Baffin Island and Ellesmere Island, in the southeast of Iceland, on Svalbard, over part of Novaya Zemlya, and in a small area over the Rocky Mountains (Fig. ), in good agreement with glacier inventories .

Figure 3

Present-day (2000 CE) ice sheet extent and thickness from the reference transient Holocene CLIMBER-X simulation. The Greenland ice sheet is prescribed in this experiment.

[Figure omitted. See PDF]

Table 1

List of last glacial inception experiments. $T$ offset indicates the uniform and constant temperature offset added to the simulated surface air temperature in the SEMIX model when computing the surface energy and mass balance. Geo indicates whether the topography and land–ocean–ice sheet mask are updated in the simulations or not. Snow albedo offset indicates the uniform and constant offset in snow albedo applied in SEMIX.

Experiment	Ice	$T$ bias	$T$ offset	Vegetation	Geo	GHGs	Dust	Snow	Ice sheet	Acceleration
	model	corr	($\mathrm{°}\mathrm{C}$)				scaling	albedo	resolution
								offset	($\mathrm{km}$)
Ref	on	on	0	dyn	on	var	1	0	32	1
NoIce	off	–	0	dyn	on	var	1	0	–	1
NoIceFixVeg	off	–	0	fix PI	on	var	1	0	–	1
FixVeg	on	on	0	fix PI	on	var	1	0	32	1
FixGeo	on	on	0	dyn	off	var	1	0	32	1
FixGHG	on	on	0	dyn	on	fix 125 $\mathrm{ka}$	1	0	32	1
NoBiasCorr	on	off	0	dyn	on	var	1	0	32	1
Tm1	on	on	$-$1	dyn	on	var	1	0	32	1
Tp1	on	on	+1	dyn	on	var	1	0	32	1
Dust05	on	on	0	dyn	on	var	0.5	0	32	1
Dust2	on	on	0	dyn	on	var	2	0	32	1
$\mathit{\alpha }-$	on	on	0	dyn	on	var	1	$-$0.025	32	1
$\mathit{\alpha }+$	on	on	0	dyn	on	var	1	+0.025	32	1
Res16	on	on	0	dyn	on	var	1	0	16	1
Res64	on	on	0	dyn	on	var	1	0	64	1
Acc2	on	on	0	dyn	on	var	1	0	32	2
Acc5	on	on	0	dyn	on	var	1	0	32	5
Acc10	on	on	0	dyn	on	var	1	0	32	10
Acc20	on	on	0	dyn	on	var	1	0	32	20
Acc50	on	on	0	dyn	on	var	1	0	32	50
Acc100	on	on	0	dyn	on	var	1	0	32	100

The transient last glacial inception simulations are started from the Eemian interglacial at 125 $\mathrm{ka}$ and run until 100 $\mathrm{ka}$. This time interval is characterized first by a decrease in NH summer insolation, with a minimum reached at $\sim$ 116 $\mathrm{ka}$ followed by an increase in insolation until $\sim$ 105 $\mathrm{ka}$ (Fig. a). The model is driven by prescribed changes in the orbital parameters and atmospheric concentrations of ${\mathrm{CO}}_{\mathrm{2}}$, ${\mathrm{CH}}_{\mathrm{4}}$ and ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ from ice core data (Fig. b). The initial conditions of the model runs correspond to the climate model in equilibrium with 125 $\mathrm{ka}$ boundary conditions. We therefore make the reasonable assumption that climate was in quasi-equilibrium at 125 $\mathrm{ka}$. Because the Greenland ice sheet was likely not in equilibrium with climate at this time, for practical reasons we choose to start from the Greenland ice sheet prescribed from present-day observations with a uniform ice temperature of $-$10 $\mathrm{°}\mathrm{C}$. This is justified because Greenland is not the focus of our study and plays a negligible role in the glacial inception over the NH continents. With this general setup we performed different simulations to test the importance of different processes, namely (i) vegetation feedback, (ii) ice sheet feedback, and (iii) role of ${\mathrm{CO}}_{\mathrm{2}}$ variations, and additional experiments to explore the sensitivity of the simulated last glacial inception to (i) the temperature bias correction in SEMIX, (ii) snow albedo parameterization, (iii) ice sheet model resolution, and (iv) climate acceleration factor. The full set of experiments is listed in Table .

Figure 4

Transient last glacial inception simulation with reference model parameters. (a) Maximum summer insolation at 65 $\mathrm{°}\mathrm{N}$ , (b) equivalent ${\mathrm{CO}}_{\mathrm{2}}$ concentration for radiation (including the radiative effect of ${\mathrm{CO}}_{\mathrm{2}}$, ${\mathrm{CH}}_{\mathrm{4}}$, and ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$; ). (c) Simulated ice sheet area. (d) Simulated global relative sea level change compared to reconstructions . The blue shading indicates the 1 standard deviation uncertainty range.

[Figure omitted. See PDF]

4 Modelling results

4.1 Reference transient simulation of the last glacial inception

In the transient simulation of the last glacial inception, until $\sim$ 120 $\mathrm{ka}$ only minor changes in ice sheet area (Fig. c) and sea level (Fig. d) are produced by the model. This is followed by a rapid increase in ice sheet area by $\sim$ 10 million square kilometres between $\sim$ 119 and $\sim$ 117 $\mathrm{ka}$ (Fig. c), driven by the further decrease in NH summer insolation (Fig. a). Afterwards, ice area stabilizes, while ice volume continues to grow and sea level drops to $\sim$ $-$35 $\mathrm{m}$ by $\sim$ 110 $\mathrm{ka}$. After $\sim$ 110 $\mathrm{ka}$ both ice area and volume decrease again following an increase in summer insolation and despite a global cooling induced by a pronounced decrease in the atmospheric concentration of greenhouse gases, mainly ${\mathrm{CO}}_{\mathrm{2}}$ (Fig. b). There is an overall good agreement between simulated and reconstructed sea level in terms of timing, while the amplitude of the changes is somewhat underestimated (Fig. d). Part of this discrepancy could be related to the missing contribution from Antarctica, which is prescribed at its present-day state in our simulations and could have contributed $\sim$ 5 $\mathrm{m}$ to the sea level drop e.g..

A two-stage relation between ice sheet area and ice sheet volume in the first part of the glacial inception simulation is clearly shown in Fig. . Initially, a substantial increase in ice area occurs at a nearly constant ice volume, followed by a large increase in ice volume while the ice area remains almost constant. During the deglaciation phase no two-stage behaviour in the area–volume relation is found.

Figure 5

Relation between NH ice sheet area changes and ice sheet volume changes for the glacial inception simulation with the reference model version. Symbols indicate the area and volume every 50 years. The colour scale shows the progression of time.

[Figure omitted. See PDF]

Figure 6

Ice sheet extent and thickness at different points in time for the reference last glacial inception simulation.

[Figure omitted. See PDF]

At 121 $\mathrm{ka}$ the simulated ice sheet extent is comparable to the present day but with some loss of ice in southern Greenland (Fig. a). The ice then starts to nucleate at high altitudes in the Arctic islands and over the mountain ranges of Scandinavia and the Ellesmere and Baffin islands (Fig. b) before rapidly expanding over the Canadian Arctic Archipelago and Scandinavia, mainly due to large-scale thickening of snowfields (Fig. c). This is in line with what has been found by . Starting from the Arctic islands, ice also spreads into the Barents and Kara seas and the Cordilleran ice sheet is established (Fig. c). The ice then generally grows thicker and slowly expands to reach its maximum extent at $\sim$ 115 $\mathrm{ka}$ (Fig. d). Between 115 and 110 $\mathrm{ka}$, ice sheets mainly grow thicker with little change in ice extent (Fig. e). Towards the end of the simulation at $\sim$ 100 $\mathrm{ka}$, about two-thirds of the maximum “glacial” ice is melted (Fig. f).

The simulated ice sheet cover at 110 $\mathrm{ka}$ compares reasonably well with the reconstructions of and for marine isotope stage (MIS) 5d (Fig. ), also considering the relatively large associated uncertainty range (shown by the best and maximum reconstructed ice extent in Fig. ). Models and reconstructions show a good agreement in ice sheet cover over Fennoscandia, Iceland, and Greenland. Simulated ice cover is possibly underestimated over eastern North America, particularly over the Labrador region, although reconstructions are highly uncertain over this region. A Cordilleran ice sheet is formed in the model, while Alaska remains largely ice free. Contrary to what is indicated by the reconstructions mentioned above, the Cordilleran and Laurentide ice sheets merge in the model simulation. Little ice cover is simulated by the model over eastern Siberia, where reconstructions suggest the presence of small ice caps covering the mountain ranges.

Figure 7

Simulated ice sheet extent at 110 $\mathrm{ka}$ (MIS5d) compared to the best and maximum extent reconstructions from and .

[Figure omitted. See PDF]

The relative sea level represents the geoid height change relative to the displaced Earth surface. Its signature at 110 $\mathrm{ka}$, shown in Fig. , mainly represents the sum of the viscoelastic subsidence due to the loading of the ice sheet as geoid changes due to gravitational changes and due to the mass conservation between ice and ocean. Accordingly, the RSL rise in the loaded regions is dominated by the subsidence of the Earth surface due to the viscoelastic response of the solid Earth, whereas the drop over the oceans reflects the mass conservation between ice and ocean, meaning the water equivalent of the ice sheet mass. Furthermore, the ice sheets are surrounded by regions of lower RSL due the forebulge of the flexed elastic lithosphere. The bedrock deformation and depression is largest at the centres of the ice sheets, corresponding to a RSL that is up to 500 $\mathrm{m}$ higher than pre-industrial values, while the far-field RSL is $\sim$ 35 $\mathrm{m}$ lower than pre-industrial values (Fig. ). Alaska shows the highest RSL values, which can be attributed to the low viscosity structure of the solid Earth in this region allowing a faster rebound in comparison to the cratonic areas of NE Canada and Scandinavia (see viscosity structure in Fig.  in Appendix ).

Figure 8

Simulated relative sea level changes relative to pre-industrial values at 110 $\mathrm{ka}$. The magenta lines indicate the ice sheet extent.

[Figure omitted. See PDF]

In terms of changes in climate, from 125 to 110 $\mathrm{ka}$ the model shows a pronounced summer cooling at middle to high northern latitudes (Fig. c) as a direct response to the decrease in summer insolation (Fig. a), amplified by a strong albedo increase (Fig. b), which is driven by a combination of expanding sea ice, a southward shift of the treeline (Fig. ), and the establishment of land ice. Going from 121 to 115 $\mathrm{ka}$, summer temperatures undergo a substantial cooling over most of the NH, particularly over land (Fig. ).

Figure 9

Northern Hemisphere zonal mean differences in (a) maximum summer insolation, (b) summer (JJA) surface albedo, and (c) summer (JJA) surface air temperature as a function of time for the reference glacial inception simulation relative to the pre-industrial values.

[Figure omitted. See PDF]

Figure 10

Forest fraction at different points in time simulated by the model in the reference experiment. The magenta lines indicate the ice sheet extent.

[Figure omitted. See PDF]

Figure 11

Summer (JJA) temperature difference relative to the pre-industrial values at different times in the reference glacial inception simulations. Note the non-linear colour scale.

[Figure omitted. See PDF]

4.2 Role of vegetation, ice sheet, and carbon cycle feedbacks

The effect of the vegetation feedback on the expansion of ice sheets during glacial inception can be quantified by running a glacial inception simulation with vegetation prescribed at its pre-industrial state. Practically, fixed vegetation in our simulation means that the plant functional type fractions are not allowed to change and that the maximum leaf area index is fixed. However, for deciduous plants, the seasonality of the phenology will still be affected by the changing climatic conditions. The vegetation feedback plays a crucial role for glacial inception in our simulations. This is clearly illustrated by the much smaller increase in ice sheet area and volume in simulations where vegetation is prescribed at its equilibrium pre-industrial state, compared to the reference glacial inception simulation with interactive vegetation (Fig. ).

Figure 12

(a) Ice sheet area and (b) relative sea level changes for model simulations with fixed vegetation (blue), fixed present-day topography and land–ocean–ice sheet mask (red), and fixed 125 $\mathrm{ka}$ greenhouse gas concentrations (green) compared to the reference simulation (black).

[Figure omitted. See PDF]

Figure 13

(a–c) Simulated forest fraction at different times during the glacial inception simulation with dynamic vegetation but prescribed present-day ice sheets (simulation NoIce). Differences in (d–f) surface albedo and (g–i) surface air temperature in late spring–early summer (May–June–July) between simulations with prescribed present-day ice sheets and dynamic (NoIce) and fixed (NoIceFixVeg) vegetation at different times.

[Figure omitted. See PDF]

Figure 14

Simulated ice sheet extent and thickness at 115 $\mathrm{ka}$ for different experiments: (a) reference, (b) fixed pre-industrial vegetation, (c) fixed 125 $\mathrm{ka}$ GHGs concentrations, (d) no temperature bias correction over North America, (e) reduced snow albedo, and (f) increased snow albedo.

[Figure omitted. See PDF]

The effect of dynamic vegetation on climate can be isolated from two simulations, one with prescribed vegetation and one with dynamic vegetation, both without interactive ice sheets. The strong vegetation feedback is explained by a pronounced southward shift of the northern boreal treeline as a response to the gradual decrease in summer insolation at high northern latitudes during glacial inception (Fig. a–c). The expansion of tundra at the expense of taiga results in a substantial increase in surface albedo during the snow-covered season (Fig. d–f) with a particularly strong associated cooling in spring, which also extends into the summer months (Fig. g–i). This explains why the ice extent in the simulation with prescribed vegetation is generally reduced compared to the interactive vegetation case (Fig. a and b).

Our results on the important role played by the vegetation feedback for the initiation of NH glaciation are consistent with a number of previous studies that have shown that the vegetation response during the last glacial inception amplifies the orbitally induced summer cooling in high northern latitudes , thus favouring the growth of ice sheets .

The effect of ice sheets on glacial inception can be quantified with a simulation in which the land–ocean–ice sheet mask and the topography are prescribed at their pre-industrial state. This is equivalent to disabling the back-coupling of the ice sheets to the climate model and therefore suppressing the ice sheet feedback on climate. The model results indicate that this ice sheet feedback plays only a minor role during the ice growth phase until $\sim$ 115 $\mathrm{ka}$ (Fig. ). This is explained by the ice expansion being driven by a rapid increase of perennially snow-covered area rather than by a slow lateral expansion of ice sheets. However, the higher albedo of ice compared to ice-free land plays an important role in slowing down the ice sheet melt during the ice retreat phase following the rising summer insolation after 115 $\mathrm{ka}$ (Fig. ).

The role of the ${\mathrm{CO}}_{\mathrm{2}}$ variations during glacial inception can be estimated from a simulation where the atmospheric concentrations of the greenhouse gases are kept constant at their 125 $\mathrm{ka}$ values. Since the GHG concentrations show only small variations until $\sim$ 115 $\mathrm{ka}$ (Fig. b), it is not surprising that the GHG forcing plays only a minor role during the first phase of glacial inception (Fig. ). Hence, the simulated ice sheets in the experiment with prescribed constant GHGs are very similar to those in the reference simulation at 115 $\mathrm{ka}$ (Fig. a and c). However, the decrease in the equivalent ${\mathrm{CO}}_{\mathrm{2}}$ concentration after $\sim$ 115 $\mathrm{ka}$ is important for slowing down the ice sheet melt and limit deglaciation (Fig. ).

Figure 15

Simulated (a) ice sheet area and (b) relative sea level for different uniform temperature offsets in SEMIX (experiments Tp1 and Tm1 in Table ) and for the simulation without temperature bias correction (NoBiasCorr).

[Figure omitted. See PDF]

4.3 Sensitivity to climate model biases

It is known that ice sheets can be highly sensitive to relatively small temperature changes. For instance, it has been shown that the bifurcation point for the complete melting of the Greenland ice sheet could be at only a few degrees above pre-industrial values . We therefore decided to apply different uniform temperature offsets in the surface energy and mass balance model and use these for sensitivity tests. This is justified because the global mean surface air temperature is (i) quite uncertain and (ii) different state-of-the-art climate models produce very different global mean temperatures . Our simulations show that the last glacial inception is sensitive to relatively small temperature perturbations in the surface mass balance model (Fig. ). In particular, the difference in simulated sea level decrease between the experiment with a uniform cooling of 1 $\mathrm{°}\mathrm{C}$ and the experiment with a uniform warming of 1 $\mathrm{°}\mathrm{C}$ in SEMIX is $\sim$ 35 $\mathrm{m}$ (Fig. ).

Without the temperature correction over North America (Fig. ), the simulated ice sheet distribution over this continent is shifted from the east to the west, as expected (Fig. a,d). Little ice is simulated in the area around the Hudson Bay, while ice extends further in the northwest. A similar east–west displacement of the ice distribution has also been found in other models that share temperature biases similar to the CLIMBER-X ones over North America e.g..

4.4 Sensitivity to snow albedo parameterization

The albedo of snow is a function of snow grain size, with smaller grain sizes resulting in higher albedo e.g.. Fresh dry snow has a generally small grain size, but the grain size tends to increase as the snow undergoes metamorphism processes and in particular as melting occurs. Snow albedo is also affected by the presence of light-absorbing impurities, such as mineral dust or soot particles e.g.. Since the albedo of snow is generally high, small relative changes in snow albedo will result in large relative changes in co-albedo, which is the relevant quantity determining the amount of absorbed solar radiation at the surface. It is therefore expected that uncertainties in the parameterization of the albedo of snow will result in substantial differences in the surface energy and mass balance e.g.. To explore this we performed simulations in which we added a constant offset of $\pm$ 0.025 to the albedo of snow and experiments with half and double the dust deposition rate in SEMIX. These perturbations to the snow albedo result in a simulated sea level decrease during MIS5d that varies by more than 40 $\mathrm{m}$ (Fig. ). This is remarkable considering that the snow albedo changes introduced are of the order of only a few percent. In terms of spatial pattern, the simulated ice sheet is more extended and twice as large by volume if the snow albedo is higher and less extended if the snow albedo is lower than in the reference run (Fig. e and f).

Figure 16

(a) Ice sheet area and (b) relative sea level changes for model simulations with scaled dust deposition fields and offsets in snow albedo.

[Figure omitted. See PDF]

Figure 17

Simulated (a) ice sheet area and (b) relative sea level when different climate model acceleration factors are applied.

[Figure omitted. See PDF]

4.5 Sensitivity to climate model acceleration

Several efforts are ongoing to simulate the last glacial cycle with state-of-the-art Earth system models based on general circulation models e.g.. Computational time is a strong constraint for these models and acceleration techniques e.g. are one possible way to alleviate this problem. These are based on the fact that typical timescales of the atmosphere, ocean, and vegetation are much shorter than orbital timescales. This allows for the to artificial acceleration of the external forcings, in this case orbital parameters and GHGs concentration, which, considering that the atmosphere and ocean are also the most computationally expensive parts of an Earth system model, results in an effective speed-up of the model. The ice sheet model cannot be accelerated in the same way. However, this is not problematic because ice sheet models are comparatively less computationally demanding.

Figure 18

Simulated ice sheet extent and thickness at 115 $\mathrm{ka}$ for simulations with different climate acceleration factors: (a) reference run without climate acceleration, (b) acceleration factor of 10, and (c) acceleration factor of 50.

[Figure omitted. See PDF]

To assess the applicability of the acceleration technique, we perform additionally several experiments with different acceleration rates to explore how sensitive the ice sheet evolution during glacial inception is to the applied acceleration factor. Our results show only a relatively weak sensitivity of the ice sheet evolution to climate acceleration for acceleration factors up to $\sim$ 10 (Fig. ), confirming earlier results from the CLIMBER-2 model . Much larger accelerations factors do not allow for a proper representation of the positive climate feedbacks at work during glacial inception, resulting in reduced simulated ice extent and volume (Fig. ).

4.6 Sensitivity to ice sheet model resolution

We also tested the dependence of the simulated glacial inception on the resolution of the ice sheet model. A higher-resolution ice sheet model will result in a better preservation of the fine-scale topographic structure, and since the surface mass balance is strongly dependent on surface elevation, it is expected that better resolving mountain peaks would facilitate the formation of ice caps. It is unclear, however, whether this would facilitate the formation of large-scale ice sheets or not, also because better resolving mountains also implies that deep valleys are better resolved, which would inhibit the spreading of ice from isolated ice caps to larger-scale ice sheets.

Our simulations show only a weak dependence of the model results on the resolution of the ice sheet model in the tested range between 16 and 64 $\mathrm{km}$ (Fig. ). However, some local ice caps that are formed with a high-resolution ice sheet model are not resolved when the resolution is decreased (Fig. ).

Figure 19

Simulated (a) ice sheet area and (b) relative sea level for different resolutions of the ice sheet model.

[Figure omitted. See PDF]

Figure 20

Simulated ice sheet extent and thickness at 115 $\mathrm{ka}$ for simulations with different horizontal resolutions of the ice sheet model: (a) 16 $\mathrm{km}$, (b) 32 $\mathrm{km}$ (Ref), and (c) 64 $\mathrm{km}$.

[Figure omitted. See PDF]

5 Discussion and conclusions

We have presented the results of a set of transient simulations of the last glacial inception with the CLIMBER-X Earth system model, which includes an ice sheet model and a model for the solid Earth response to changes in ice sheet loading. This paper also describes the ice sheet coupling with the atmosphere and ocean, which will serve as a reference for future studies using the model with interactive ice sheets.

We have shown that, as a response to the decreasing summer insolation at high northern latitudes, the model simulates a rapid expansion of ice sheets over northern North America and Scandinavia between $\sim$ 120 and $\sim$ 116 $\mathrm{ka}$, which is driven mainly by large-scale snowfield thickening. This result is fully consistent with the concept of glacial inception as a bifurcation in the climate system caused by a strong albedo feedback . The rapid expansion of the ice sheet area was followed by ice volume growth that corresponds to a global sea level drop of nearly 40 $\mathrm{m}$ at 110 $\mathrm{ka}$, which is in reasonable agreement with paleoclimate reconstructions. Compared to the CLIMBER-2 results from , as expected the CLIMBER-X model results are in better agreement with available paleoclimate reconstructions in terms of simulated ice volume and spatial distribution of ice sheets. The abrupt increase of the North American ice sheet extent simulated in CLIMBER-X is nearly half that shown in , and this increase happens in two steps. Less ice is formed in eastern Siberia and Alaska in CLIMBER-X than in CLIMBER-2. At the same time, despite substantial differences in model formulations and spatial resolution, the results of the simulations with CLIMBER-X confirm the main findings presented in . The abrupt increase in the North American ice sheet area occurs in CLIMBER-X approximately at the same time and in the same area as in and through the same mechanisms. Our results confirm the critically important role of vegetation feedback demonstrated in . The albedo feedback associated with an increase in snow-covered area, sea ice extent, and the southward retreat of the boreal forest plays a crucial role in the rapid ice area expansion in our simulations. The vegetation feedback alone increases the maximum simulated ice sheet area by 50 $\mathrm{\%}$.

Our new modelling results show a strong sensitivity of simulated ice sheet evolution to the parameterization of clean snow albedo and to the impact of impurities on snow albedo. The ice sheet feedback and the variations in GHGs concentrations are of minor importance during the ice growth phase associated with decreasing summer insolation prior to $\sim$ 115 $\mathrm{ka}$ but are fundamental in maintaining the system in a glacial state during the subsequent period of increase in summer insolation, resulting in an only partial deglaciation.

The results of model simulations demonstrate that the reduction of climate biases (too high summer air temperatures over eastern North America) leads to significant improvements in the simulated spatial extent of the North American ice sheet, as also shown by .

The model results are not very sensitive to climate acceleration up to a factor $\sim$ 10, confirming earlier findings by and . Assuming that this finding also holds for more complex Earth system models, a climate acceleration factor of 10 would allow these models to run transient glacial inception simulations in a reasonable time using less computational resources. The resolution of the ice sheet model only marginally affects the model results, at least in the tested range between 16 and 64 $\mathrm{km}$. This is because a large-scale ice expansion over relatively flat terrain is the dominant mechanism leading to glacial inception in our model, while ice caps, which can be captured only if the topography is highly resolved, have only a very localized effect and are therefore not of fundamental importance.

The glacial inception simulations presented here are a first step towards simulating the full last glacial cycle with CLIMBER-X.

Appendix A SICOPOLIS ice sheet model

For the inclusion into CLIMBER-X, the original SICOPOLIS code has been restructured and organized into Fortran 90-derived types to allow running several instances of the model at the same time, one for each ice sheet domain. Although in the present study we employ a single ice sheet domain for the NH, CLIMBER-X can be set up to run with an arbitrary number of ice sheet model domains, with potentially different resolutions, simply by specifying a list of domain names in a namelist.

As already mentioned above, SICOPOLIS is based on the shallow ice approximation for grounded ice, the shallow shelf approximation for floating ice, and hybrid shallow-ice–shelfy-stream dynamics for ice streams , and the enthalpy method of is used as thermodynamics solver.

In this study we use a Weertman-type sliding law to relate basal drag (${\mathit{\tau }}_{\mathrm{b}}$) to basal velocity ($u_{\mathrm{b}}$):

A1$${\mathit{\tau }}_{\mathrm{b}}={\displaystyle \frac{\mathrm{1}}{c_{\mathrm{s}}^{\mathrm{1}/p}}}{N}^{q/p}{u}_{\mathrm{b}}^{\mathrm{1}/p},$$ with $p=\mathrm{3}$ and $q=\mathrm{2}$. $N$ is the reduced basal pressure computed as the difference between the ice overburden pressure and the water pressure if the base of the ice sheet is below sea level: A2$$N=P_{\mathrm{i}}-P_{\mathrm{w}}={\mathit{\rho }}_{\mathrm{i}}g{h}_{\mathrm{i}}-{\mathit{\rho }}_{\mathrm{w}}gmax(\mathrm{0},z_{\text{bed}}),$$ where ${\mathit{\rho }}_{\mathrm{i}}$ and ${\mathit{\rho }}_{\mathrm{w}}$ are the ice and water densities, $g$ is the acceleration due to gravity, $h_{\mathrm{i}}$ is ice thickness, and $z_{\text{bed}}$ is the bedrock depth below sea level. The $c_{\mathrm{s}}$ sliding parameter depends on the assumed sediment fraction in each grid cell: A3$$c_{\mathrm{s}}=(\mathrm{1}-f_{\text{sed}})c_{\mathrm{s}}^{\text{rock}}+f_{\text{sed}}{c}_{\mathrm{s}}^{\text{sed}}.$$

The sediment fraction is taken to be a linear function of sediment thickness between two critical values $h_{\text{sed}}^{\text{min}}$ and $h_{\text{sed}}^{\text{max}}$: A4$$f_{\text{sed}}={\displaystyle \frac{h_{\text{sed}}-h_{\text{sed}}^{\text{min}}}{h_{\text{sed}}^{\text{max}}-h_{\text{sed}}^{\text{min}}}},$$ where $f_{\text{sed}}$ is computed based on sediment thickness data from and is shown in Fig. . The value for bedrock is set to $c_{\mathrm{s}}^{\text{rock}}$ $=$ 25 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{Pa}}^{-\mathrm{1}}$ following the tuning for the Greenland ice sheet in , while the value over thick sediments is simply taken as 10 times this value, i.e. $c_{\mathrm{s}}^{\text{sed}}$ $=$ 250 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{Pa}}^{-\mathrm{1}}$.

Figure A1

Sediment fraction used in the computation of the basal sliding coefficient $c_{\mathrm{s}}$.

[Figure omitted. See PDF]

A broader investigation of the sliding law will be performed in forthcoming papers discussing the simulation of glacial termination, when the sliding is expected to be more important.

Iceberg calving at the marine-terminating ice front is parameterized using a simple thickness-based threshold method in which the calving rate is computed as follows: A5$$C={\displaystyle \frac{h_{\text{ice}}^{\mathrm{c}}-h_{\text{ice}}}{{\mathit{\tau }}_{\mathrm{c}}}}.$$

Table A1

List of SICOPOLIS model parameters.

Parameter	Value	Description
$c_{\mathrm{s}}^{\text{rock}}$	25 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{Pa}}^{-\mathrm{1}}$	sliding coefficient for rocks
$c_{\mathrm{s}}^{\text{sed}}$	250 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{Pa}}^{-\mathrm{1}}$	sliding coefficient for sediments
$h_{\text{sed}}^{\text{min}}$	10 $\mathrm{m}$	sediment thickness below which bare rock is assumed
$h_{\text{sed}}^{\text{max}}$	500 $\mathrm{m}$	sediment thickness above which full sediment cover is assumed
$z_{\text{bed}}^{\text{shallow}}$	200 $\mathrm{m}$	shallow-water depth for calving
$z_{\text{bed}}^{\text{deep}}$	1000 $\mathrm{m}$	deep-water depth for calving
$h_{\text{ice}}^{\mathrm{c},\text{shallow}}$	50 $\mathrm{m}$	critical calving thickness in shallow waters
$h_{\text{ice}}^{\mathrm{c},\text{deep}}$	500 $\mathrm{m}$	critical calving thickness in deep waters
${\mathit{\tau }}_{\mathrm{c}}$	10 years	timescale for calving

It is applied only where the ice thickness, $h_{\text{ice}}$, is lower than the critical thickness for calving, $h_{\text{ice}}^{\mathrm{c}}$, also in all neighbouring ice points. The critical thickness for calving varies spatially and increases linearly with the depth of bedrock: A6$$h_{\text{ice}}^{\mathrm{c}}=h_{\text{ice}}^{\mathrm{c},\text{deep}}+\left(z_{\text{bed}}^{\text{fil}}-z_{\text{bed}}^{\text{deep}}\right)\cdot {\displaystyle \frac{h_{\text{ice}}^{\mathrm{c},\text{deep}}-h_{\text{ice}}^{\mathrm{c},\text{shallow}}}{z_{\text{bed}}^{\text{deep}}-z_{\text{bed}}^{\text{shallow}}}},$$ where $z_{\text{bed}}^{\text{fil}}$ is the bedrock elevation filtered with a Gaussian filter with a 100 $\mathrm{km}$ radius and $z_{\text{bed}}^{\text{shallow}}$, $z_{\text{bed}}^{\text{deep}}$, $h_{\text{ice}}^{\mathrm{c},\text{shallow}}$, and $h_{\text{ice}}^{\mathrm{c},\text{deep}}$ are model parameters (Table ). Although this is only a very crude representation, it is arguably more reasonable than assuming a uniform threshold that does not account for the different ocean dynamics in the narrow channels of the Canadian Arctic Archipelago and the deep open ocean. Lower values of $h_{\text{ice}}^{\mathrm{c}}$ allow for a more rapid expansion of floating ice shelves, thereby also affecting the rate of shelf ice spreading. However, since floating ice is usually thin, this has only a minor impact on the simulated ice volume during glacial inception.

The global geothermal heat flux product of , substituted with data from over Greenland, is applied as bottom boundary condition for the bedrock temperature equation at a depth of 2 $\mathrm{km}$ below the land surface in SICOPOLIS.

For all simulations presented in this study we used an annual time step for the thermodynamic part and a half-yearly time step for the dynamics.

Appendix B SEMIX surface energy and mass balance interface

SEMIX is an adaptation of the surface energy and mass balance interface (SEMI) to CLIMBER-X. Its purpose is to determine the surface boundary conditions, namely surface mass balance and surface ice temperature, for the ice sheet model. In order to do that SEMIX has to bridge the gap in resolution between the climate model and the ice sheet model, which is achieved through a downscaling of the climate variables (Appendix ). The surface energy and mass balance equations are then solved on the high-resolution ice sheet grid as outlined in Appendix . Since multiple ice sheet domains are allowed in CLIMBER-X, similar to the ice sheet model SEMIX can also run in several separate instances according to the defined ice sheet domains. SEMIX is called every 10 years over a full year with a time step of 1 $\mathrm{d}$.

B1 Downscaling of climate variables

SEMIX is driven by climate fields computed by the atmospheric model SESAM . Since SESAM and SEMIX operate on very different grids, the first step required in the SEMIX coupling is the mapping of the climate fields from the coarse resolution, regular lat–long SESAM grid onto the Cartesian coordinates on the stereographic plane where the ice sheet model, and consequently SEMIX, operates. The mapping is done by simple bilinear interpolation using the four closest neighbouring grid points. A list of the climate model variables needed by SEMIX is given in Table . In the following we will denote the climate variables mapped onto the ice sheet grid with a superscript $\mathrm{i}$.

Table B1

List climate model variables needed by SEMIX.

Variable	Description	Unit
$z_{\mathrm{s}}$	grid cell mean elevation	$\mathrm{m}$
${\mathit{\sigma }}_{\mathrm{z}}$	standard deviation of sub-grid surface elevation over land	$\mathrm{m}$
$T_{\text{atm}}$	atmospheric temperature	$\mathrm{K}$
$r_{\text{atm}}$	atmospheric relative humidity
$P$	total precipitation rate	$\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$
$U$	surface wind speed	${\mathrm{ms}}^{-\mathrm{1}}$
$D$	dust deposition rate	$\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$
$f_{\text{cld}}$	cloud cover fraction
${\mathrm{SW}}_{\mathrm{TOA}}^{\downarrow }$	downward short-wave radiation at TOA	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$
$\mathit{\mu }$	daily mean cosine of solar zenith angle	$\mathrm{rad}$
${\mathrm{SW}}_{\mathrm{VIS},\mathrm{dir}}^{\downarrow }$	downward short-wave visible radiation at surface, clear sky	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$
${\mathrm{SW}}_{\mathrm{NIR},\mathrm{dir}}^{\downarrow }$	downward short-wave near-infrared radiation at surface, clear sky	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$
${\mathrm{SW}}_{\mathrm{VIS},\mathrm{dif}}^{\downarrow }$	downward short-wave visible radiation at surface, cloudy sky	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$
${\mathrm{SW}}_{\mathrm{NIR},\mathrm{dif}}^{\downarrow }$	downward short-wave near-infrared radiation at surface, cloudy sky	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$
${\mathit{\alpha }}_{\text{VIS,dir}}$	surface albedo for visible radiation at surface, clear sky
${\mathit{\alpha }}_{\text{NIR,dir}}$	surface albedo for near-infrared radiation at surface, clear sky
${\mathit{\alpha }}_{\text{VIS,dif}}$	surface albedo for visible radiation at surface, cloudy sky
${\mathit{\alpha }}_{\text{NIR,dif}}$	surface albedo for near-infrared radiation at surface, cloudy sky
$\frac{\partial {\mathrm{SW}}^{\downarrow }}{\partial \mathit{\alpha }}{|}_{\text{VIS,dir}}$	partial derivative of ${\mathrm{SW}}_{\mathrm{VIS},\mathrm{dir}}^{\downarrow }$ wrt ${\mathit{\alpha }}_{\text{VIS,dir}}$	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$
$\frac{\partial {\mathrm{SW}}^{\downarrow }}{\partial \mathit{\alpha }}{|}_{\text{NIR,dir}}$	partial derivative of ${\mathrm{SW}}_{\mathrm{NIR},\mathrm{dir}}^{\downarrow }$ wrt ${\mathit{\alpha }}_{\text{NIR,dir}}$	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$
$\frac{\partial {\mathrm{SW}}^{\downarrow }}{\partial \mathit{\alpha }}{|}_{\text{VIS,dif}}$	partial derivative of ${\mathrm{SW}}_{\mathrm{VIS},\mathrm{dif}}^{\downarrow }$ wrt ${\mathit{\alpha }}_{\text{VIS,dif}}$	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$
$\frac{\partial {\mathrm{SW}}^{\downarrow }}{\partial \mathit{\alpha }}{|}_{\text{NIR,dif}}$	partial derivative of ${\mathrm{SW}}_{\mathrm{NIR},\mathrm{dif}}^{\downarrow }$ wrt ${\mathit{\alpha }}_{\text{NIR,dif}}$	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$
$\frac{\partial {\mathrm{SW}}_{\mathrm{NIR},\mathrm{dir}}^{\downarrow }}{\partial z}$	partial derivative of ${\mathrm{SW}}_{\mathrm{NIR},\mathrm{dir}}^{\downarrow }$ wrt height	$\mathrm{W}{\mathrm{m}}^{-\mathrm{3}}$
$\frac{\partial {\mathrm{SW}}_{\mathrm{NIR},\mathrm{dif}}^{\downarrow }}{\partial z}$	partial derivative of ${\mathrm{SW}}_{\mathrm{NIR},\mathrm{dif}}^{\downarrow }$ wrt height	$\mathrm{W}{\mathrm{m}}^{-\mathrm{3}}$
${\mathrm{LW}}^{\downarrow }$	downward long-wave radiation at the surface	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$
$\frac{\partial {\mathrm{LW}}^{\downarrow }}{\partial z}$	partial derivative of ${\mathrm{LW}}^{\downarrow }$ wrt height	$\mathrm{W}{\mathrm{m}}^{-\mathrm{3}}$
$T_{\text{2m,JJA}}^{\text{bias}}$	near-surface summer air temperature bias at present-day	$\mathrm{K}$

A second step involves downscaling of the climate fields to account for, e.g. the difference in surface elevation between the coarse climate model and the high-resolution ice sheet grid. A constant temperature lapse rate is used to adjust the atmospheric temperature to the actual surface elevation: B1$$T_{\text{atm}}=T_{\text{atm}}^{\mathrm{i}}+\mathrm{\Gamma }\left(z_{\mathrm{s}}-z_{\mathrm{s}}^{\mathrm{i}}\right).$$

Following and we use a value of $\mathrm{\Gamma }$ $=$ $-$5 $\mathrm{K}{\mathrm{km}}^{-\mathrm{1}}$ for the lapse rate. The near-surface air temperature is computed as the average of atmospheric temperature and skin temperature: B2$$T_{\text{2m}}=\mathrm{0.5}\cdot (T_{\text{atm}}-T_{\text{skin}}).$$

$T_{\text{2m}}$ is then used to separate total precipitation into rain and snow, with the fraction of precipitation falling as snow computed as follows: B3$$f_{\text{snow}}=\left\{\begin{array}{l}\mathrm{1},T_{\text{2m}}\le (T_{\mathrm{0}}-\mathrm{5})\\ \mathrm{0.1}\cdot (T_{\mathrm{0}}+\mathrm{5}-T_{\text{2m}}),(T_{\mathrm{0}}-\mathrm{5})<T_{\text{2m}}<(T_{\mathrm{0}}+\mathrm{5})\\ \mathrm{0},T_{\text{2m}}\ge (T_{\mathrm{0}}+\mathrm{5})\end{array}\right.,$$ with $T_{\mathrm{0}}$ $=$ 273.15 $\mathrm{K}$. Snowfall and rainfall rate are then derived from total precipitation as follows: $$\begin{array}{}\text{B4}& {\displaystyle }& {\displaystyle }{P}_{\text{snow}}=f_{\text{snow}}\cdot P^{\mathrm{i}},\text{B5}& {\displaystyle }& {\displaystyle }{P}_{\text{rain}}=(\mathrm{1}-f_{\text{snow}})\cdot P^{\mathrm{i}}.\end{array}$$

Near-surface air-specific humidity is computed from near-surface relative humidity and specific humidity at saturation as in the climate component of CLIMBER-X : B6$$q_{\text{2m}}=r_{\text{2m}}\cdot q_{\text{sat}}(T_{\text{2m}},p_{\mathrm{s}}),$$ with $r_{\text{2m}}=(r_{\text{atm}}+r_{\text{skin}})/\mathrm{2}$, $r_{\text{skin}}=q_{\text{atm}}/q_{\text{sat}}(T_{\text{skin}},p_{\mathrm{s}})$, and $q_{\text{atm}}=r_{\text{atm}}^{\mathrm{i}}\cdot q_{\text{sat}}(T_{\text{atm}},p_{\mathrm{s}})$.

The downward short-wave radiation fluxes at the surface on the ice sheet model grid are computed from the interpolated fluxes and adjusted using the partial derivatives of the radiation fields with respect to surface albedo and surface elevation. $$\begin{array}{}\text{B7}& {\displaystyle }\begin{array}{rl}{{\mathrm{SW}}^{\downarrow }}_{\text{VIS,dir}}& ={\mathrm{SW}}_{\mathrm{VIS},\mathrm{dir}}^{\downarrow ,\mathrm{i}}\\ & +{\left.{\displaystyle \frac{\partial {\mathrm{SW}}^{\downarrow }}{\partial \mathit{\alpha }}}\right|}_{\text{VIS,dir}}^{\mathrm{i}}\left({\mathit{\alpha }}_{\text{VIS,dir}}-{\mathit{\alpha }}_{\text{VIS,dir}}^{\mathrm{i}}\right)\end{array}\text{B8}& {\displaystyle }\begin{array}{rl}{\mathrm{SW}}_{\mathrm{NIR},\mathrm{dir}}^{\downarrow }& ={\mathrm{SW}}_{\mathrm{NIR},\mathrm{dir}}^{\downarrow ,\mathrm{i}}\\ & +{\left.{\displaystyle \frac{\partial {\mathrm{SW}}^{\downarrow }}{\partial \mathit{\alpha }}}\right|}_{\text{NIR,dir}}^{\mathrm{i}}\left({\mathit{\alpha }}_{\text{NIR,dir}}-{\mathit{\alpha }}_{\text{NIR,dir}}^{\mathrm{i}}\right)\\ & +{{\displaystyle \frac{\partial {\mathrm{SW}}_{\mathrm{NIR},\mathrm{dir}}^{\downarrow }}{\partial z}}}^{\mathrm{i}}\left(z_{\mathrm{s}}-z_{\mathrm{s}}^{\mathrm{i}}\right)\end{array}\text{B9}& {\displaystyle }\begin{array}{rl}{{\mathrm{SW}}^{\downarrow }}_{\text{VIS,dif}}& ={\mathrm{SW}}_{\mathrm{VIS},\mathrm{dif}}^{\downarrow ,\mathrm{i}}\\ & +{\left.{\displaystyle \frac{\partial {\mathrm{SW}}^{\downarrow }}{\partial \mathit{\alpha }}}\right|}_{\text{VIS,dif}}^{\mathrm{i}}\left({\mathit{\alpha }}_{\text{VIS,dif}}-{\mathit{\alpha }}_{\text{VIS,dif}}^{\mathrm{i}}\right)\end{array}\end{array}$$

B10$$\begin{array}{rl}{\mathrm{SW}}_{\mathrm{NIR},\mathrm{dif}}^{\downarrow }& ={\mathrm{SW}}_{\mathrm{NIR},\mathrm{dif}}^{\downarrow ,\mathrm{i}}\\ & +{\left.{\displaystyle \frac{\partial {\mathrm{SW}}^{\downarrow }}{\partial \mathit{\alpha }}}\right|}_{\text{NIR,dif}}^{\mathrm{i}}\left({\mathit{\alpha }}_{\text{NIR,dif}}-{\mathit{\alpha }}_{\text{NIR,dif}}^{\mathrm{i}}\right)\\ & +{{\displaystyle \frac{\partial {\mathrm{SW}}_{\mathrm{NIR},\mathrm{dif}}^{\downarrow }}{\partial z}}}^{\mathrm{i}}\left(z_{\mathrm{s}}-z_{\mathrm{s}}^{\mathrm{i}}\right)\end{array}$$

The first correction term is needed because the downward short-wave radiation flux at the surface depends itself on the surface albedo through its influence on multiple scattering between surface and atmosphere. A larger surface albedo will generally result in a larger downward short-wave radiation flux. The elevation correction term for the near-infrared component arises from the absorption of short-wave radiation by the atmosphere in the near-infrared band, which introduces the elevation dependence of the radiative flux. Finally, the net short-wave flux absorbed at the surface, which is the quantity entering the surface energy balance equation, is derived as a weighted average of clear-sky (direct radiation) and cloudy-sky (diffuse radiation) radiative fluxes using cloud cover fraction. B11$$\begin{array}{rl}{\mathrm{SW}}_{\mathrm{net}}& =\left(\mathrm{1}-f_{\text{cld}}^{\mathrm{i}}\right)\cdot ({{\mathrm{SW}}^{\downarrow }}_{\text{VIS,dir}}\cdot \left(\mathrm{1}-{\mathit{\alpha }}_{\text{VIS,dir}}\right)\\ & +{\mathrm{SW}}_{\mathrm{NIR},\mathrm{dir}}^{\downarrow }\cdot \left(\mathrm{1}-{\mathit{\alpha }}_{\text{NIR,dir}}\right))\\ & +f_{\text{cld}}^{\mathrm{i}}\cdot ({{\mathrm{SW}}^{\downarrow }}_{\text{VIS,dif}}\cdot \left(\mathrm{1}-{\mathit{\alpha }}_{\text{VIS,dif}}\right)\\ & +{\mathrm{SW}}_{\mathrm{NIR},\mathrm{dif}}^{\downarrow }\cdot \left(\mathrm{1}-{\mathit{\alpha }}_{\text{NIR,dif}}\right))\end{array}$$

A similar approach is also applied for the downscaling of the downward long-wave radiation at the surface using the partial derivative of the flux with respect to elevation. B12$${\mathrm{LW}}^{\downarrow }={\mathrm{LW}}^{\downarrow ,\mathrm{i}}+{{\displaystyle \frac{\partial {\mathrm{LW}}^{\downarrow }}{\partial z}}}^{\mathrm{i}}\left(z_{\mathrm{s}}-z_{\mathrm{s}}^{\mathrm{i}}\right).$$

B2 Temperature bias correction

As described in Sect. , a temperature bias correction is applied only over North America. The bias correction is applied to the atmospheric temperature and Eq. () therefore becomes

B13$$T_{\text{atm}}=T_{\text{atm}}^{\mathrm{i}}+\mathrm{\Gamma }\left(z_{\mathrm{s}}-z_{\mathrm{s}}^{\mathrm{i}}\right)-T_{\text{2m,JJA}}^{{\text{bias}}^{\mathrm{i}}},$$ where $T_{\text{2m,JJA}}^{{\text{bias}}^{\mathrm{i}}}$ is the mean summer (June–July–August) temperature bias in CLIMBER-X relative to the ERA5 reanalysis over the time period from 1981 to 2010 and is shown in Fig. . The same temperature bias field is applied at all time steps throughout the year.

Figure B1

Temperature correction applied in SEMIX.

[Figure omitted. See PDF]

Furthermore, since the downwelling long-wave radiation at the surface, which also affects the surface energy balance, is closely related to near-surface air temperature, we also correct the downwelling long-wave radiation for the temperature bias using a simple quadratic relation derived from ERA5 reanalysis data (Fig.). Equation () is then modified to B14$$\begin{array}{rl}{\mathrm{LW}}^{\downarrow }& ={\mathrm{LW}}^{\downarrow ,\mathrm{i}}+{{\displaystyle \frac{\partial {\mathrm{LW}}^{\downarrow }}{\partial z}}}^{\mathrm{i}}\left(z_{\mathrm{s}}-z_{\mathrm{s}}^{\mathrm{i}}\right)\\ & -(\mathrm{0.034}\cdot T_{\text{2m}}-\mathrm{5})\cdot T_{\text{2m,JJA}}^{{\text{bias}}^{\mathrm{i}}}.\end{array}$$

Figure B2

Scatterplot of downwelling long-wave radiation at the surface versus near-surface air temperature in the ERA5 reanalysis . The black line indicates the best quadratic fit to the data and the inset shows the corresponding quadratic equation, which is used for bias correction in the model.

[Figure omitted. See PDF]

Present-day annual precipitation simulated by CLIMBER-X over the region where ice sheets were growing during glacial inception is in reasonable agreement with observations (Fig.), and typical biases do not exceed 200 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. At the same time, the effect of 1 $\mathrm{°}\mathrm{C}$ summer temperature bias on annual snowmelt can be estimated using the classical positive degree day approach to be 200–500 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. The low bound corresponds to a melt season duration of 2 months and parameter $\mathit{\alpha }$ $=$ 3 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}\mathrm{°}{\mathrm{C}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$, while the upper bound corresponds to the melt season duration of 3 months and $\mathit{\alpha }$ $=$ 5. Since simulated summer temperature biases are about 3–5 $\mathrm{°}\mathrm{C}$, temperature biases are much more important than precipitation biases. This is why we only correct the temperature and do not apply any corrections to precipitation.

Figure B3

Simulated present-day annual precipitation in CLIMBER-X (a) compared to ERA-Interim reanalysis  (b). Panel (c) shows the precipitation bias in CLIMBER-X.

[Figure omitted. See PDF]

B3 Surface energy and mass balance of snow and ice

The surface energy balance in SEMIX is written as follows:

B15$${\mathrm{SW}}_{\mathrm{net}}+{\mathrm{LW}}^{\downarrow }-{\mathrm{LW}}^{\uparrow }-\mathrm{SH}-\mathrm{LE}-G=\mathrm{0},$$ where ${\mathrm{SW}}_{\mathrm{net}}$ is the net short-wave radiation absorbed at the surface, ${\mathrm{LW}}^{\downarrow }$ and ${\mathrm{LW}}^{\uparrow }$ are the downwelling and upwelling long-wave radiation at the surface, $\mathrm{SH}$ is the sensible heat flux, $\mathrm{LE}$ is the latent heat flux, and $G$ the heat flux into the snow and ice. Equation () is then solved for the skin temperature, $T_{\text{skin}}$, using the formulations for the energy fluxes described next.

The surface albedo values used to compute the net short-wave at the surface in Eq. are defined as follows: $$\begin{array}{}\text{B16}& {\displaystyle }{\mathit{\alpha }}_{\text{VIS,dir}}& {\displaystyle }=f_{\text{snow}}\cdot {\mathit{\alpha }}_{\text{VIS,dir}}^{\text{snow}}+(\mathrm{1}-f_{\text{snow}})\cdot {\mathit{\alpha }}^{\text{bg}},\text{B17}& {\displaystyle }{\mathit{\alpha }}_{\text{NIR,dir}}& {\displaystyle }=f_{\text{snow}}\cdot {\mathit{\alpha }}_{\text{NIR,dir}}^{\text{snow}}+(\mathrm{1}-f_{\text{snow}})\cdot {\mathit{\alpha }}^{\text{bg}},\text{B18}& {\displaystyle }{\mathit{\alpha }}_{\text{VIS,dif}}& {\displaystyle }=f_{\text{snow}}\cdot {\mathit{\alpha }}_{\text{VIS,dif}}^{\text{snow}}+(\mathrm{1}-f_{\text{snow}})\cdot {\mathit{\alpha }}^{\text{bg}},\text{B19}& {\displaystyle }{\mathit{\alpha }}_{\text{NIR,dif}}& {\displaystyle }=f_{\text{snow}}\cdot {\mathit{\alpha }}_{\text{NIR,dif}}^{\text{snow}}+(\mathrm{1}-f_{\text{snow}})\cdot {\mathit{\alpha }}^{\text{bg}},\end{array}$$ where the snow cover fraction $f_{\text{snow}}$ is determined also considering the effect of topographic roughness following and : B20$$f_{\text{snow}}=\tanh \left({\displaystyle \frac{h_{\text{snow}}}{\mathrm{10}\cdot z_{\mathrm{0}}}}\right){\displaystyle \frac{h_{\text{snow}}}{h_{\text{snow}}+\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}\cdot {\mathit{\sigma }}_{\mathrm{z}}^{\mathrm{i}}}},$$ where $h_{\text{snow}}$ is snow thickness and $z_{\mathrm{0}}$ is the surface roughness length. The snow albedo scheme is the same as used in the climate component of CLIMBER-X and includes a dependence on snow grain size and dust and soot concentration following . The background albedo is a weighted average between a constant bare-soil albedo and variable ice albedo: B21$${\mathit{\alpha }}_{\text{bg}}=f_{\text{ice}}\cdot {\mathit{\alpha }}^{\text{ice}}+(\mathrm{1}-f_{\text{ice}})\cdot {\mathit{\alpha }}^{\text{soil}}.$$

The ice cover fraction is computed from ice sheet thickness, $h_{\text{ice}}$, and topographic roughness as: B22$$f_{\text{ice}}=\tanh \left({\displaystyle \frac{h_{\text{ice}}}{{\mathit{\sigma }}_{\mathrm{z}}^{\mathrm{i}}}}\right).$$

For ice-free grid cells next to the ice sheet margin, where $h_{\text{ice}}=\mathrm{0}$, an ice thickness is computed instead as an average over the 3 $\times$ 3 grid cells neighbourhood and will therefore by definition be $>\mathrm{0}$. The ice albedo of newly formed ice is set to a constant value representative for firn, ${\mathit{\alpha }}^{\text{firn}}=\mathrm{0.7}$. As soon as ice starts to melt, the albedo gradually decreases towards a clean-ice albedo, assuming that the clean-ice albedo is reached when the firn layer is melted: B23$${\displaystyle \frac{\mathrm{\Delta }{\mathit{\alpha }}^{\text{ice}}}{\mathrm{\Delta }t}}=-{\displaystyle \frac{{\mathit{\alpha }}^{\text{ice}}-{\mathit{\alpha }}_{\text{clean}}^{\text{ice}}}{h_{\text{firn}}{\mathit{\rho }}_{\text{firn}}/M_{\text{ice}}}},$$ where $M_{\text{ice}}$ is the rate of ice melt, $h_{\text{firn}}$ is the constant thickness of the firn layer, and ${\mathit{\rho }}_{\text{firn}}$ is the density of firn, which is also considered to be constant (Table ). As the ice becomes snow free and therefore exposed to the deposition of dust and other wind-borne material that darken the ice, the albedo is assumed to decrease further towards a dirty-ice albedo, ${\mathit{\alpha }}_{\text{dirty}}^{\text{ice}}=\mathrm{0.4}$, with a timescale of 100 years. Whenever the annual surface mass balance is positive, indicating accumulation of snow and consequently formation of a firn layer, the ice albedo is reset to the albedo of firn (${\mathit{\alpha }}^{\text{ice}}={\mathit{\alpha }}^{\text{firn}}$). A few considerations are appropriate here on the representation of ice albedo in ice-free grid cells at the ice sheet margin. If the surface mass balance is positive in these points, the ice (or background) albedo is not very important because it implies that not all snow is melted during summer and therefore the background does not become exposed. On the other hand, if the surface mass balance is negative, the background albedo in the ice-free margin points should reflect the properties of the ice that could flow into these points from neighbouring grid cells through horizontal ice flow. In this case the ice albedo will determine how negative the mass balance is through its control on the absorbed radiation when all snow is melted. How negative the surface mass balance is will eventually determine whether the ice flowing into these grid cells will be melted completely or not and therefore plays a role in the velocity at which ice sheet cover can expand laterally. Based on these considerations the ice albedo in the grid cells at the ice sheet margin is computed as the average of the ice albedo of the neighbouring ice points.

Table B2

List of SEMIX model parameters.

Parameter	Value	Description
$\mathrm{\Gamma }$	5 $\mathrm{K}{\mathrm{km}}^{-\mathrm{1}}$	lapse rate of surfacetemperature
${\mathit{\alpha }}^{\text{soil}}$	0.2	bare-soil albedo
${\mathit{\alpha }}_{\text{firn}}^{\text{ice}}$	0.7	firn albedo
${\mathit{\alpha }}_{\text{clean}}^{\text{ice}}$	0.55	clean-ice albedo
${\mathit{\alpha }}_{\text{dirty}}^{\text{ice}}$	0.4	dirty-ice albedo
$h_{\text{firn}}$	100 $\mathrm{m}$	thickness of firn layer
${\mathit{\rho }}_{\text{firn}}$	500 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$	firn density
${\mathit{\rho }}_{\text{snow}}$	250 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$	snow density
${\mathit{\varphi }}_{\text{snow}}$	0.7	snow porosity
$C_{\mathrm{i}}$	2110 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$	specific heat capacity of ice
$L_{\mathrm{f}}$	3.34 $\times$ 10${}^{\mathrm{5}}$ $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}$	latent heat of fusion of ice
$f_{\text{rfz,snow}}^{\text{max}}$	0.5	maximum fraction ofrefreezing in snow layer
$h_{\text{snow}}^{\text{crit}}$	4 $\mathrm{m}$	critical snow thickness forrefreezing

The surface emitted long-wave radiation is given by the Stefan–Boltzmann law: B24$${\mathrm{LW}}^{\uparrow }=\mathit{ϵ}\mathit{\sigma }{T}_{\text{skin}}^{\mathrm{4}},$$ with $\mathit{ϵ}$ the surface emissivity and $\mathit{\sigma }$ the Stefan–Boltzmann constant.

The sensible heat flux is computed from the temperature gradient between the skin and near-surface air using the bulk aerodynamic formula: B25$$\mathrm{SH}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}{c}_{\mathrm{p}}}{r_{\text{aer}}}}(T_{\text{skin}}-T_{\text{2m}}),$$ where ${\mathit{\rho }}_{\mathrm{a}}$ is air density, $c_{\mathrm{p}}$ is the specific heat of air, and $r_{\text{aer}}$ is the aerodynamic resistance. Similarly, the latent heat flux over sea ice is expressed in terms of the specific humidity gradient between the surface and near-surface air: B26$$\mathrm{LE}=L{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}}{r_{\text{aer}}}}\left(q_{\text{sat}}\right(T_{\text{skin}})-q_{\text{2m}}).$$ $L$ is the latent heat of sublimation and $q_{\text{sat}}$ is the specific humidity at saturation. The aerodynamic resistance is computed from the wind speed, surface exchange coefficient, and bulk Richardson number following .

The conductive heat flux into the snow/ice ($G$) is computed as follows: B27$$G={\displaystyle \frac{\mathit{\lambda }}{\mathrm{\Delta }{z}_{\mathrm{1}}}}(T_{\text{skin}}-T_{\text{s,1}}),$$ where $\mathit{\lambda }$ is the heat conductivity of snow and ice and $\mathrm{\Delta }{z}_{\mathrm{1}}$ is the distance between the surface and the middle of the snow layer or uppermost ice layer with temperature $T_{\text{s,1}}$, depending whether a snow layer is present or not.

The prognostic terms in $T_{\text{skin}}$ in the formulation of the surface energy fluxes are then linearized using Taylor series expansion assuming that the temperature at the new time step, $T_{\text{skin},n+\mathrm{1}}=T_{\text{skin},n}+\mathrm{\Delta }{T}_{\text{skin}}$ with $\mathrm{\Delta }{T}_{\text{skin}}\ll T_{\text{skin}}$: $$\begin{array}{}\text{B28}& {\displaystyle }{T}_{\text{skin},n+\mathrm{1}}^{\mathrm{4}}=T_{\text{skin},n}^{\mathrm{4}}+\mathrm{4}{T}_{\text{skin},n}^{\mathrm{3}}(T_{\text{skin},n+\mathrm{1}}-T_{\text{skin},n}),\text{B29}& {\displaystyle }\begin{array}{rl}{q}_{\text{sat}}\left(T_{\text{skin},n+\mathrm{1}}\right)& =q_{\text{sat}}\left(T_{\text{skin}}\right)+{\displaystyle \frac{\mathrm{d}{q}_{\text{sat}}}{\mathrm{d}{T}_{\text{skin}}}}{|}_{T_{\text{skin}}=T_{\text{skin},n}}\\ & \times (T_{\text{skin},n+\mathrm{1}}-T_{\text{skin},n}).\end{array}\end{array}$$

Equation () can then be solved explicitly for the skin temperature at the new time step, $T_{\text{skin},n+\mathrm{1}}$. If the skin temperature is above freezing the surface energy fluxes are diagnosed first with the skin temperature greater then 0 $\mathrm{°}\mathrm{C}$ and then with skin temperature set to 0 $\mathrm{°}\mathrm{C}$. The difference between the sum of the energy fluxes is then used to directly melt snow and/or ice, without the need for the snow or ice temperature to be at melting point.

The heat transfer in the snow–ice column is represented by a one-dimensional heat diffusion equation. A single snow layer is represented on top of five unevenly spaced vertical layers in the ice reaching a depth of 15 $\mathrm{m}$. The heat flux $G$ is applied as top boundary condition, while a no-flux condition is applied at the bottom of the ice column. If the temperature of the snow or ice layers is greater than 0 $\mathrm{°}\mathrm{C}$, the excess energy is used to melt snow or ice. Liquid water produced by snowmelt or added to the snow layer from rainfall can be refrozen in the snow layer or refreeze to form superimposed ice. The liquid water available for refreezing is B30$$F_{\text{avail}}=(M_{\text{snow}}+P_{\text{rain}})\mathrm{\Delta }t,$$ where $M_{\text{snow}}$ is snowmelt and $\mathrm{\Delta }t$ the time step in SEMIX (1 $\mathrm{d}$). The maximum amount of water that is allowed to refreeze is taken to be a linear function of the thickness of the snow layer and accounts for the amount of water that has already been refrozen since the start of the melt season: B31$$F_{\text{max}}={\mathit{\varphi }}_{\text{snow}}{h}_{\text{snow}}{\mathit{\rho }}_{\text{snow}}-\sum F,$$ where ${\mathit{\varphi }}_{\text{snow}}$ is snow porosity, ${\mathit{\rho }}_{\text{snow}}$ is snow density, and $F$ is refreezing. The maximum amount of water available for refreezing then becomes B32$$F_{\text{avail}}^{\text{max}}=min(F_{\text{avail}},F_{\text{max}}).$$

The amount of water that can potentially be refrozen depends on the “cold content” of the snow layer: B33$$F_{\text{pot}}={\displaystyle \frac{C_{\mathrm{i}}}{L_{\mathrm{f}}}}{h}_{\text{snow}}{\mathit{\rho }}_{\text{snow}}\cdot (T_{\mathrm{0}}-T_{\text{snow}}),$$ with $C_{\mathrm{i}}$ the specific heat capacity of ice, $L_{\mathrm{f}}$ the latent heat of fusion, and $T_{\text{snow}}$ the temperature of the snow layer. The actual refreezing rate is then simply the minimum between available and potential: B34$$F={\displaystyle \frac{min\left(F_{\text{avail}}^{\text{max}},F_{\text{pot}}\right)}{\mathrm{\Delta }t}}.$$

A fraction $f_{\text{rfz,snow}}$ is then assumed to refreeze within the snow layer, while the rest forms superimposed ice: B35$$f_{\text{rfz,snow}}=f_{\text{rfz,snow}}^{\text{max}}\cdot min\left(\mathrm{1},{\displaystyle \frac{h_{\text{snow}}}{h_{\text{snow}}^{\text{crit}}}}\right),$$ where $f_{\text{rfz,snow}}^{\text{max}}$ and $h_{\text{snow}}^{\text{crit}}$ are model parameters (Table ).

Finally, the surface mass balance and runoff are diagnosed as follows: $$\begin{array}{}\text{B36}& {\displaystyle }& {\displaystyle }\mathrm{SMB}=P_{\text{snow}}-E-M_{\text{snow}}-M_{\text{ice}}+F\text{B37}& {\displaystyle }& {\displaystyle }R=M_{\text{snow}}+M_{\text{ice}}+P_{\text{rain}}-F.\end{array}$$

These are then integrated over the whole year and passed to the ice sheet model. The surface ice temperature, which is needed by the ice sheet model as top boundary condition in the temperature equation, is computed as the annual average temperature of the deepest ice layer ($\sim$ 10 $\mathrm{m}$).

Appendix C IMO ice shelf basal melt model

In CLIMBER-X we have implemented the simple and general ice shelf basal melt parameterizations of and , both of which rely on the difference between the ambient water temperature derived from ocean or lake model and the freezing point temperature at the ice shelf base.

The basal mass balance in the linear model of is

C1$$\mathrm{BMB}=k_{\mathrm{1}}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}{C}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{i}}{L}_{\mathrm{f}}}}(T_{\mathrm{w}}-T_{\mathrm{f}}),$$ where $k_{\mathrm{1}}$ is a model parameter (Table ), ${\mathit{\rho }}_{\mathrm{w}}$ is a typical seawater density, $C_{\mathrm{w}}$ is the specific heat capacity of water, ${\mathit{\rho }}_{\mathrm{i}}$ is the density of ice, and $L_{\mathrm{f}}$ is the latent heat of fusion of ice. $T_{\mathrm{w}}$ is either the seawater temperature from the ocean model for marine-terminating margins or the lake water temperature for lake-terminating margins, horizontally extrapolated to the ice sheet model grid and vertically interpolated to the depth of the ice shelf base. $T_{\mathrm{f}}$ is the freezing point temperature at the base of the ice shelf : C2$$T_{\mathrm{f}}=\mathrm{0.0939}-\mathrm{0.057}\cdot S_{\mathrm{w}}-\mathrm{7.64}\times {\mathrm{10}}^{-\mathrm{4}}\cdot z_{\mathrm{b}},$$ where $S_{\mathrm{w}}$ is the water salinity, derived similarly to $T_{\mathrm{w}}$, assuming that the salinity of lakes is zero, and $z_{\mathrm{b}}$ is the depth of the ice shelf base below sea level.

The basal mass balance in the model of is similar but relies on a quadratic dependence on the temperature gradient: C3$$\mathrm{BMB}=k_{\mathrm{2}}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}{C}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{i}}{L}_{\mathrm{f}}}}\left|T_{\mathrm{w}}-T_{\mathrm{f}}\right|\cdot (T_{\mathrm{w}}-T_{\mathrm{f}}),$$ where $k_{\mathrm{2}}$ is a model parameter; all other terms in this equation have already been defined above.

Similarly to , in order to avoid unrealistic ice shelf expansion over the deep ocean, we apply an additional scaling of basal melt with the depth of the bedrock elevation: C4$$f_{\mathrm{BMB}}=\mathrm{1}+{\displaystyle \frac{max\left(\mathrm{0},z_{\text{bed}}-z_{\text{bed}}^{\text{ref}}\right)}{z_{\text{bed}}^{\text{ref}}}},$$ with $z_{\text{bed}}^{\text{ref}}$ being a model parameter (Table ).

IMO is called every month to resolve the seasonal cycle in ocean and lake temperature, but the coupling with the ice sheet model is done yearly.

Table C1

List of IMO model parameters.

Parameter	Value	Description
$k_{\mathrm{1}}$	1 $\times$ 10${}^{-\mathrm{4}}$	parameter for the model
$k_{\mathrm{2}}$	5 $\times$ 10${}^{-\mathrm{4}}$	parameter for the model
$z_{\text{bed}}^{\text{ref}}$	200 $\mathrm{m}$	critical bedrock depth for basal melt scaling

Appendix D VILMA solid Earth model

The solid Earth model VILMA solves the field equations of a viscoelastic incompressible and self-gravitating continuum in the spherical domain for the mantle and lithosphere of the Earth. The fluid core and loading processes at the surface are considered boundary conditions at the respective boundaries. Lateral changes in viscosity are considered and are following the formulation as an initial value problem discussed in : the field equations are solved in the spherical harmonic domain but integrated over time with an explicit time differencing scheme. In view of solving the momentum equation in the spectral solution, the code is rather efficient and in view of the time-differencing scheme the coupling with the Earth system model is straightforward . The loading is considered as a pressure field boundary condition applied at the surface, where mass conservation is considered solving the sea level equation at each time step. The time step interval is constrained by the minimum ratio of viscosity versus shear modulus (the Maxwell time). It is usually about 20 years for a standard Earth structure but has to be reduced to about 1 year for a 3D structure containing structural features like low viscous regions in tectonically active regions. For the 3D viscosity field in this study we chose the 3D model $v_{\mathrm{1}}{.0}_s\mathrm{16}$ from (Fig. ) as it shows the smallest misfit with observational data in independent GIA models. In the present work we set a minimum viscosity of 10${}^{\mathrm{19.5}}$, which allows a time step of 10 years to be used. Due to the fact that the time step is rather small, an iteration at each time step as discussed in can be neglected.

Figure D1

Viscosity at 200 $\mathrm{km}$ depth from as used in the model simulations.

[Figure omitted. See PDF]

The relative sea level determined by VILMA is used as the spatially variable correction of the considered reference topography, $h_{\text{topo}}\left(t_{\mathrm{0}}\right)$. In this way, changes in the topography are considered with respect to the changing geoid defined here as the mean sea level at the respective time step: D1$$h_{\text{topo}}\left(t\right)=h_{\text{topo}}\left(t_{\mathrm{0}}\right)-h_{\text{rsl}}(t-t_{\mathrm{0}}),h_{\text{rsl}}\left(t_{\mathrm{0}}\right)=\mathrm{0}.$$

This view is consistent with the natural definition of topography in the understanding of Carl Friedrich Gauss. It also means that all changes in elevation are expressed as measured relative to the mean sea level at that time.

Changes in topography due to variations in sea level are considered furthermore by updating the land–sea mask at each time step. This holds also for the conditions of floating versus grounded ice. In view of completing formulations e.g., the effect of rotational deformations was implemented recently in VILMA following the formulation of . Rotational deformations have to be considered, especially when discussing the effect of GIA on geodetic observables like GNSS, EOP, and gravity or when discussing future sea level changes e.g..

Code and data availability

The source code of the climate component of CLIMBER-X v1.0 as used in the simulations of this paper is archived on Zenodo (10.5281/zenodo.7898797, ). The output of the reference model simulation is available at 10.5281/zenodo.10643295 .

Author contributions

MW, AG, RC, and ST designed the study. RG provided the SICOPOLIS model code and RC, RG, and JB assisted in the implementation of SICOPOLIS into CLIMBER-X. RC and MW adapted and tested SICOPOLIS for the integration into CLIMBER-X. MW, RC, and AG developed SEMIX. VK and MB provided the VILMA code and contributed to its coupling. MW coupled, tested, and tuned the different model components. MW ran the simulations and prepared the figures. MW wrote the paper, with contributions from all co-authors.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Special issue statement

This article is part of the special issue “Icy landscapes of the past”. It is not associated with a conference.

Acknowledgements

The authors gratefully acknowledge the European Regional Development Fund (ERDF), the German Federal Ministry of Education and Research, and the Land Brandenburg for supporting this project by providing resources on the high-performance computer system at the Potsdam Institute for Climate Impact Research.

Financial support

Matteo Willeit and Meike Bagge are funded by the German climate modelling project PalMod supported by the German Federal Ministry of Education and Research (BMBF) as a Research for Sustainability initiative (FONA) (grant nos. 01LP1920B, 01LP1917D, 01LP1918A). Ralf Greve was supported by Japan Society for the Promotion of Science (JSPS) KAKENHI (grant nos. JP17H06104 and JP17H06323) and by the Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT) through the Arctic Challenge for Sustainability project ArCS II (grant no. JPMXD1420318865). The publication of this article was funded by the Open Access Fund of the Leibniz Association.

Review statement

This paper was edited by Pepijn Bakker and reviewed by two anonymous referees.