Journal topic
Clim. Past, 14, 1441–1462, 2018
https://doi.org/10.5194/cp-14-1441-2018
Clim. Past, 14, 1441–1462, 2018
https://doi.org/10.5194/cp-14-1441-2018

Research article 18 Oct 2018

Research article | 18 Oct 2018

# The role of regional feedbacks in glacial inception on Baffin Island: the interaction of ice flow and meteorology

The role of regional feedbacks in glacial inception on Baffin Island: the interaction of ice flow and meteorology
Leah Birch1, Timothy Cronin2, and Eli Tziperman1 Leah Birch et al.
• 1School of Engineering and Applied Sciences and Department of Earth and Planetary Sciences, Harvard University, Cambridge, MA, USA
• 2Department of Earth, Atmosphere and Planetary Sciences, MIT, Cambridge, MA, USA

Correspondence: Leah Birch (lbirch@seas.harvard.edu)

Abstract

Over the past 0.8 million years, 100 kyr ice ages have dominated Earth's climate with geological evidence suggesting the last glacial inception began in the mountains of Baffin Island. Currently, state-of-the-art global climate models (GCMs) have difficulty simulating glacial inception, possibly due in part to their coarse horizontal resolution and the neglect of ice flow dynamics in some models. We attempt to address the role of regional feedbacks in the initial inception problem on Baffin Island by asynchronously coupling the Weather Research and Forecast (WRF) model, configured as a high-resolution inner domain over Baffin and an outer domain incorporating much of North America, to an ice flow model using the shallow ice approximation. The mass balance is calculated from WRF simulations and used to drive the ice model, which updates the ice extent and elevation, that then serve as inputs to the next WRF run. We drive the regional WRF configuration using atmospheric boundary conditions from 1986 that correspond to a relatively cold summer, and with 115 kya insolation. Initially, ice accumulates on mountain glaciers, driving downslope ice flow which expands the size of the ice caps. However, continued iterations of the atmosphere and ice models reveal a stagnation of the ice sheet on Baffin Island, driven by melting due to warmer temperatures at the margins of the ice caps. This warming is caused by changes in the regional circulation that are forced by elevation changes due to the ice growth. A stabilizing feedback between ice elevation and atmospheric circulation thus prevents full inception from occurring.

1 Introduction

For the past million years, Earth's climate has repeatedly fallen into ice ages with a periodicity of ∼100 000 years . How glacial inceptions occur is still a largely unresolved problem in climate science . The leading hypothesis for the cause of glacial inception is that orbital variations cause cooler summers leading to less melting, although previous studies found that additional amplifying feedbacks beyond the ice–albedo feedback are necessary to explain inception. Possible amplifiers of orbital forcing may be divided into large-scale non-local feedbacks, including ocean temperature and planetary-scale circulation changes, and local feedbacks associated with greater height and extent of mountain glaciers , regional circulation changes, clouds, and more. In this study, we seek to understand specifically the role of local ice–climate feedbacks in glacial inception using a regional climate model asynchronously coupled to a simple ice sheet model. Our regional model is centered over Baffin Island – the likely location of the last glacial inception .

Glacial inception of the Laurentide ice sheet is thought to have progressed fairly rapidly from Baffin Island in eastern Canada , with the associated drop in sea level estimated to be as large as 50 m in 10 kyrs . The Laurentide ice sheet has been identified as the primary initial cause of the drop in sea level , and some ice modeling studies have shown rapid inception within a few thousand years . However, this rapid growth required extreme precipitation increases or excessive climate cooling every thousand years, without any explicit explanation. A more modest sea level drop of 10–20 m in 10 kyrs has also been posited , which indicates a slower ice sheet growth rate and possibly ice growth relegated to Arctic Canadian islands .

Global climate models (GCMs) have been used to simulate glacial inception, focusing on large-scale feedbacks. They often find that snow does not accumulate in unglaciated areas or that glaciers do not spread into unglaciated areas. One reason GCMs have difficulty simulating the glacial inception is their coarse horizontal resolution, which lowers the topographic features dramatically . Snow then does not accumulate because the region is too warm due to artificially low elevations . Though topography is important , other aspects of the climate system may also be misrepresented in GCMs, like ocean feedbacks , dust , sea ice , vegetation , and meteorology . Increasing the model resolution has led to moderate improvements in the simulation of glacial inception , through the occurrence of areas of net snow accumulation, which can progressively expand as the climate cools. Given that GCMs had such difficulties in simulating glacial inception, our study focuses on the importance of regional feedbacks over Baffin Island.

While Milankovitch forcing explains some summer cooling, the large-scale circulation is an important driver of heat and moisture fluxes to inception areas. In the current climatology, large-scale cyclonic circulation over Baffin Bay is associated with cool summers , which may be all that is necessary for inception. However, and also noted increased storms and poleward moisture fluxes during inception. showed that, in addition to favorable circulation and temperature, snow could accumulate in new areas only after artificial large-scale cooling was prescribed. The following question remains: how does this cooling occur? One possible answer is the interaction of ice sheets and the atmosphere.

As the ice sheet grows, it influences the atmosphere due to albedo and surface elevation changes. The regional circulation can also be diverted by the changing topography. Anticyclones have been noted to occur over mountains (Smith1979), but ultimately the effects of mountains on the regional climate and large-scale circulation depend on the mean overlying flow and mountain size . Looking at large-scale ice growth in the context of glacial inception, argued using a simple model that slight changes in topography can influence the stationary waves, significantly changing the circulation and leading to colder temperatures. The topography and geographic location may be critical to the glacial inception process , but previous studies did not consider the effect of ice-growth-induced regional topography changes on Baffin Island.

To capture the interaction of ice and the atmosphere, climate models of intermediate complexity have been combined with ice sheet models. This type of coupled ice–climate model allows thousands of years of ice growth due to changes in insolation and CO2 to be simulated. Often, however, these models glaciate Alaska and western Canada very rapidly , which is not consistent with the geological record showing inception occurred first over Baffin Island in eastern Canada . Ice sheet models and GCM studies have a similar problem of not agreeing with the geological record. Combining GCMs and ice sheet modeling may allow for progress in understanding ice–atmosphere interactions , but the smoothed topography and other biases in these coarse-resolution climate models often need to be artificially corrected to obtain an accurate regional climate. Furthermore, mountain ice caps of the Canadian Arctic archipelago in the present climate are too small to be modeled by large-scale ice models – hence our use of a higher-resolution regional ice model.

The interaction of ice and atmosphere is a critical component of glacial inception, often left out of climate models. For instance, focused on a single ice cap on Baffin Island using a high-resolution cloud-permitting regional climate model without an ice flow model. That study found that ice would accumulate on the mountain glaciers of Baffin Island. The 4 km resolution of that study captured both the ablation and accumulation zones of mountain glaciers. By integrating the mass balance from the top of the ice cap downslope, implicitly represented the effects of ice flow. As a result, the Penny Ice Cap was deduced to grow outside of its present-day bounds due to decreased insolation and a sequence of cold summers and wet winters. Consistent with previous GCM studies, no snow accumulated in unglaciated areas.

Thus, the purpose of the current paper is to study the role of regional feedbacks in the inception over the major ice caps of Baffin Island, including the role of ice flow. By prescribing boundary conditions from the modern climate, we do not consider any role of non-local ice sheet growth in modifying the planetary wave pattern, which could aid or hinder inception on Baffin Island. This is consistent both with other studies and with our intention to study the effect of regional feedbacks. Similarly, we neglect the role of sea surface temperature (SST) changes, as there is debate over warm or cold oceans during inception . We choose to keep modern SSTs and identify whether inception may be described as a local process over Baffin Island, driven by Milankovitch forcing. Regional ice–albedo and height–mass balance feedbacks are expected to dominate at least some of the planetary-scale feedbacks at the very beginning of the inception process. We drive the model with boundary conditions corresponding to wet and cool summers from reanalysis, ideal for inception but not extremely uncommon during present-day climate. Such conditions may be the result of SST forcing that is not explicitly simulated, so in that sense we do include some SST effects. Our regional model configuration could also be driven by the boundary conditions from a GCM simulation of glacial inception, which would incorporate some planetary-scale feedbacks. However, given the inability of such models to simulate the inception and often the lack of substantial circulation changes over Baffin Island , we choose not to take that route and focus solely on regional feedbacks.

We investigate the feedbacks between a growing ice sheet and the atmosphere at higher spatial resolution than can be simulated with GCMs, using an shallow ice approximation ice flow model asynchronously coupled with the Weather Research and Forecasting model Skamarock et al. (WRF; 2008). We find that the change to ice sheet topography induces an anticyclonic atmospheric response over Baffin Island. This anticyclone causes a flow of warm air over the Baffin Island and Hudson Bay area, increasing ablation, stagnating further ice sheet growth, and thus indicating a negative feedback on inception at the regional scale.

The atmospheric and ice models used for the asynchronous coupling experiment are presented in Sect. 2. We present the results in Sect. 3, beginning with the mass balances and resulting ice thickness of our atmosphere–ice iterations. We then investigate the regional climate and local circulation impacts of expanding ice cover and increasing ice elevation. Sensitivity simulations exploring the importance of topography and insolation are also presented. Finally, the sensitivity of ice flow to ice viscosity and mass balance are examined. We conclude with a discussion in Sect. 4, including the limitations of this study and its implications.

2 Methods

Our goal in this study is to better understand the interaction between changing ice cover and the regional climate of Baffin Island and the feedback on the regional circulation over Canada. We therefore use a nested configuration of the WRF model with the outer domain covering all of North America with a 100 km resolution and the inner domain covering Baffin Island with a 20 km resolution, encompassing both the Penny and Barnes ice caps. To better understand the growth of these mountain glaciers during inception, we must simulate ice flow in the inner 20 km resolution domain, but WRF does not contain an ice flow model, which means the ice expansion must be calculated offline. Based off of , we developed a simple ice flow model that reduces the 3-D Stokes flow to two dimensions using the shallow ice approximation. The basal sliding velocity is assumed to be 0. This approximation results in a non-linear diffusion equation that represents ice flow by evolving ice thickness H in time based on horizontal gradients of surface elevation H*, along with surface mass balance G:

$\begin{array}{}\text{(1)}& \frac{\partial H}{\partial t}=\mathrm{\nabla }\cdot D\mathrm{\nabla }{H}^{*}+G.\end{array}$

In Eq. (1), the yearly surface mass balance G is calculated from WRF simulations and assumes that the density of ice is 900 kg m−3. The surface elevation H* and ice thickness H are linked by ${H}^{*}=B+H$, where B is the basal topography. The non-linear shear thinning of ice flow is represented by a diffusivity, D, which depends on both the ice thickness and surface elevation gradients:

$\begin{array}{}\text{(2)}& D=A{H}^{m+\mathrm{2}}{\left[{\left(\frac{\partial {H}^{*}}{\partial x}\right)}^{\mathrm{2}}+{\left(\frac{\partial {H}^{*}}{\partial y}\right)}^{\mathrm{2}}\right]}^{\left(m-\mathrm{1}\right)/\mathrm{2}}.\end{array}$

The parameter m for all model simulations is set to 3, which is a typical value for realistic ice simulations . The coefficient A in Eq. (2) is the creep parameter that influences the viscosity of the ice based on the temperature of the ice; i.e., the colder ice flows more slowly, while warmer ice flows more quickly. We use a creep parameter $A=\mathrm{3.5}×{\mathrm{10}}^{-\mathrm{25}}$ Pa−3 s−1 from . More details on this choice of creep parameter can be found in Appendix A, which also contains the details of the numerical scheme used to integrate this ice model in time.

The ice model requires the specification of initial surface elevation (H*), initial ice thickness (H), basal topography (B) that is time-independent, and mass balance (G) as an input for each time step. The ice thickness (H) is adapted from the present-day ice cover on Baffin Island from the NASA Operation IceBridge mission , while H* and G are taken from our WRF simulations. The grid size and time step are also user specified and chosen to be 20 km (same as WRF inner domain) and 1 year, respectively. With these atmospheric and ice models, we perform iterations, running WRF, using the mass balance calculated by WRF to force the ice model, updating the ice thickness and extent calculated by the ice model, and then running WRF again with the newly calculated ice cover as an input. This amounts to an asynchronous coupling of the two models, where a WRF run followed by the ice model is considered an iteration.

WRF , version 3.7, is a compressible, non-hydrostatic model, which has been shown to provide a reasonable simulation of mountainous and Arctic regions . Reanalysis data from the European Centre for Medium-Range Weather Forecasts Dee et al. (ECMWF; 2011) are used for setting the meteorological boundary and initial conditions, including sea surface temperatures and sea ice extent. also noted that the interim reanalysis (ERA-Interim) is a suitable choice of boundary conditions for WRF in the Arctic.

We use the WRF single-moment six-class microphysics scheme, which allows for ice, snow, and graupel processes . The community Noah land surface model with multiparameterization options (Noah-MP) is used for the land surface model. It includes four layers of snow and allows for refreezing, which is an important process on the Penny Ice Cap according to historical trends . Noah-MP caps the snow water content at 2000 kg m−2 by default, but we remove this limit. The radiation scheme is the Rapid Radiative Transfer Model for GCM applications Iacono et al. (RRTMG; 2008), which has been found to have minimal bias in WRF . CO2 is set to 290 ppm for a glacial inception scenario . The boundary layer is simulated with the Mellor–Yamada–Janjic scheme Janjic (MYJ; 1994). Lateral boundary conditions are prescribed at the edges of the outer domain, yet no nudging to observations is used in the domain interior, as we want to allow the atmosphere to respond to land surface changes. Such spectral nudging is often used in regional models to eliminate artificial circulation patterns that develop in spite of the prescribed boundary conditions . Comparing a present-day control run to ECMWF, we find (Fig. S1 in the Supplement) non-negligible deviations in geopotential height. Yet we show below that the anomalous circulation that develops in response to ice growth in our simulations is likely unrelated to these anomalies, as they develop gradually in response to the growing ice height. Note that these physics are the same as those used in our previous paper , except now we use the Grell convective scheme , which has been shown to give good results in the present day and the Last Glacial Maximum (LGM) Arctic climates.

We use realistic present-day land types from the WRF database as initial conditions for the ice extent on Baffin Island. Over these land grid points classified as ice, a snow water equivalent (SWE) of 5000 kg m−2 is prescribed – sufficient for maintaining snow cover for the duration of our runs and allowing us to capture a representative yearly mass balance for the chosen climatology. Greenland is initialized with present-day ice cover and elevation based on studies finding little influence of Greenland on the build up of the Laurentide ice sheet .

We first run a present-day meteorology simulation forced by ECMWF boundary conditions, representative of average and cold atmospheric conditions, also simulated by . With the above WRF configuration, we find a reasonable reproduction of the Arctic atmospheric state when we look at the full fields (Fig. S1). While there are deviations in geopotential and temperature between WRF and ECMWF, the differences in geopotential have been seen in other Arctic studies . The temperature bias occurs from differences in the snow cover between the models. For instance, ECMWF has been noted to have too large of a snow-covered area on coasts , promoting cooling on Baffin Island which is mostly coast at the reanalysis resolution. Also, June is the cause of the large temperature bias because ECMWF still has snow on the ground, associated with delayed snowmelt , while the snow in WRF melts earlier. July and August only differ in temperature by about 1 C, when the snow cover in ECMWF and WRF are more similar. With the warm temperature bias, WRF will likely not cause an artificial inception, which making it more suitable than models with a cold bias. Furthermore, ECMWF is known to have a cold bias and a significant low geopotential height in the central Arctic of the same magnitude as the deviation we see in Fig. S1. Therefore, with WRF having the opposite biases, the Arctic climate simulated is not abnormal. Our WRF simulation matches observations of precipitation on the coast of Baffin Island, important for the surface mass balance, of about 300 mm per year . In any case, the warm bias in our WRF control run does not prevent the circulation from bringing relatively cold and wet conditions to Baffin Island . Previous works also proceeded to look at climate sensitivity by comparing model experiments in spite of biases in the control run , as we do here, and we therefore do not believe the biases would have a significant impact on the sensitivity results here.

For our main set of glacial inception simulations, the incoming solar radiation is set to values consistent with 115 kya; the process of modifying the insolation in WRF based on the orbital parameters is detailed in the Appendix of . The conditions most conducive to snow accumulation are from October 1985 to October 1986. These cold meteorology boundary conditions caused the present-day ice cover to expand when combined with 115 kya insolation in the study of . Figure S3 in the Supplement shows cumulative mass balance for present-day average conditions, present-day cold meteorology conditions, and cold meteorology plus 115 kya orbital forcing. The combination of moderately cool temperatures and reduced summer insolation leads to net snow accumulation, hinting that inception might be viable given sufficiently positive regional feedbacks. A longer averaging of the forcing, as obtained from a climate model, would have some advantages. However, given the persistent difficulties in simulating glacial inception using such GCMs, our goal here is to minimize the introduction of biases from global climate simulations. We therefore prefer to choose the necessary cold and wet meteorology from the observed modern climate (reanalysis). We cannot explain how these ideal circulation patterns occur during inception, but a future study with a GCM might carefully explore this. We do note that forcing with a single-year meteorology may bias the ice growth patterns if that year is not representative of longer-term cold conditions.

## 2.1 Asynchronous coupling procedure

We begin our WRF–ice flow coupled model experiment with a WRF simulation that uses meteorology from October 1985 to October 1986 (Iteration 1). WRF is run for 5 years, but we only analyze the last 4 years for an average climatology. Accumulation occurs over areas initialized with ice cover at high elevations, and net mass balance is positive, allowing ice growth (Fig. 1a). The Iteration 1 simulation agrees with our previous 4 km simulation of the Penny Ice Cap : Milankovitch forcing and cold meteorology cause the positive mass balance needed for the ice cap to expand. This mass balance is then adapted for use in the ice model.

Figure 1Ice thickness and mass balance for ice sheet initialization: (a) ice thickness was adapted from data collected by the NASA Operation IceBridge mission with extent set by land use type in WRF. (b) Annual mass balance from the WRF Iteration 1 run. (c) Scatter plot of mass balances vs. elevation with a binned regression visualization

For the first iteration of the ice model, initial ice thickness is adapted from present-day ice cover from the NASA Operation IceBridge Mission (Fig. 1a). We use the surface elevation from the 20 km resolution Iteration 1 WRF domain. The mass balance needed to force the ice model is calculated from the WRF simulation. Figure 1b depicts the mass balance results from the Iteration 1 WRF simulation with cold meteorology, 115 kya insolation, and present-day ice cover. The mass balance over the ice caps is generally positive, except on the edges of the ice cap where there is net ablation. Over the rest of the domain, WRF cannot calculate a meaningful annual mass balance because all snowfall received melts off during the summer. Thus, to have a mass balance in ice-free grid cells, we characterize the mass balance as a function of elevation (G(H*)). This characterization involves sorting the domain by surface elevation and forming bins every 100 m, and allowing the mass balance during the ice model run to be updated at each grid point as the surface elevation and ice thickness change every time step. The average the mass balance in each bin is depicted by the purple crosses in Fig. 1c. Note the plateau in the mass balance at high elevations, which is expected from the elevation desert effect.

For our asynchronous coupling experiment, the ice model is run for 500 years, as any longer has been shown to inhibit ice growth artificially . We find 500 years sufficient to grow ice while maintaining computationally efficiency. After 500 years, the ice model has calculated a new ice extent and surface elevation (Fig. 2a) with which to force WRF. We import the surface elevation directly into the 20 km resolution WRF domain and use the ice extent to update the land type, i.e., the expansion of the ice cap. Much like our Iteration 1 simulation, over areas classified as ice, we prescribe a snow depth of 5000 kg m−2. The surface elevation and land type of the 100 km resolution parent domain are adjusted through interpolation of the 20 km resolution surface elevation and land type. Soil temperatures from the previous WRF run are used to initialize the next WRF iteration. In our simulations, soil temperatures take a year to equilibrate, leading us to agree with that present-day soil temperatures are too warm. Again using 115 kya insolation and cold meteorology, WRF is run for 5 years to characterize the mass balance on this new ice cap. This new run is the start of the next iteration. Thus, each iteration begins with a WRF run and ends after the 500-year ice model integration.

3 Results

In the context of glacial inception, expanding ice cover interacts with the atmosphere, causing changes to both the regional climate and atmospheric circulation near Baffin Island. We first examine the progression of the WRF–ice model iterations, looking at the expansion of ice cover and the resultant changes in mass balance. Then, in Sect. 3.1, the changes in mass balance are explained by exploring changes in temperature, ablation, precipitation, and cloud cover. In Sect. 3.2, we take a closer look at the regional circulation. Additional sensitivity simulations with the ice sheet model further clarify climate factors that could lead to Baffin Island glaciation in Sect. 3.5.

Figure 2Ice growth due to WRF–ice model iterations: (a) the ice thickness from the ice sheet model and (b) the mass balance from WRF. The average mass balance over ice-covered areas for each iteration is included in the right column.

The progression of WRF–ice model iterations (5 years of WRF and 500 years of the ice model) is laid out in Fig. 2. The ice thickness calculated by the ice model and then used in WRF is depicted in the left column, and the resulting annual mass balance from WRF in the right column. Also included in the right column is the spatial average of mass balance on the ice cap (defined as all grid cells with ice land cover type). Starting with the top left, we show the ice thickness resulting from running the ice model for 500 years using the Iteration 1 WRF mass balance (Fig. 1b, c). This new ice thickness, ice extent, and surface elevation computed by the ice model are used as inputs to the next WRF simulation. Then, WRF is run for 5 years resulting in the mass balance presented in the right column. After the second iteration of WRF, we bin the mass balance by elevation (as seen in Fig. 1c) and force the ice model with this new updated mass balance. A new ice cover is computed (Fig. 2a, second row) and used as an input in WRF. Subsequent iterations pictured follow the same procedure (going down Fig. 2, left to right).

Figure 3Characterizing annual mass balance vs. surface elevation. The domain is sorted into 100 m bins, and within each bin the mass balance is averaged. Elevations falling within the range of the bin have the associated mass balance.

The second iteration has the greatest ice extent, which is due to the large positive Iteration 1 mass balance. In the third iteration, the ice retreats from its margins, though there is still growth in the highest parts of the domain. Additional iterations show that the ice cap coalesces and increases in thickness to nearly 2000 m but expands slowly into new areas. After 5 kyrs, the ice accumulation on Baffin Island would only lower global sea level by 0.6 m. Although the average mass balance does generally increase with time, the rest of our analysis focuses on why the ice does not expand quickly.

The binned mass balance used in the ice model is plotted in Fig. 3 with annual mass balance on the y axis and surface elevation on the x axis. The Iteration 1 WRF run is the “+” blue line with the most positive mass balance, particularly at the lowest and highest elevations. In subsequent WRF runs, the mass balance up to 1000 m becomes more negative at low elevations, which we find surprising. One might expect the equilibrium line to shift downslope as a result of regional cooling from the initial expansion of the ice cover, leading to a more positive mass balance. However, plotting the binned mass balance for all iterations shows a general trend of increasingly more negative mass balances at low elevations, with the equilibrium line altitude shifting upwards by several hundred meters over the course of 10 iterations. As a result, the ice caps coalesce and thicken at high elevations but stagnate and fail to expand horizontally. The mass balance at highest elevations from all iterations also illustrates desertification or the mass balance becoming less positive over time as the ice cap increases in height and its top surface flattens.

Lateral expansion is slow because the lower terrain on Baffin Island is not strongly sloped. Thus, ice must build up to a substantial height before the spread of ice outweighs ablation at the ice edge. For example, with flat basal topography, calculations based on Eqs. (1) and (2) suggest that a 0.5 m yr1 of ice thickness increase due to ice flow in a grid cell with no ice requires an adjacent grid cell ice thickness around 500 m. This approximate ice thickness on the margins occurs in our simulations as depicted in Fig. 2.

Figure 4The summer 2 m temperature in select asynchronous coupling iterations: (a) the temperature in the Iteration 1 WRF run. Differences in temperature between Iteration 1 and the (b) second and (c) last iterations. Note that “bin” refers to the mass balance binning procedure used in each iteration.

## 3.1 Mechanisms of changes in mass balance

A closer examination of temperature, melting, and precipitation can help us how the ice cap growth causes the changes to mass balance in the WRF simulations. We begin with the 2 m temperature on Baffin Island during the summer (June, July, August) when melting occurs. The differences in the temperature between the Iteration 1 run and subsequent iterations are pictured in Fig. 4. In the Iteration 1 simulation, we note that most of the ice-covered land remains below freezing during the summer (blue). At the edge of the ice cap, temperatures generally remain near zero (white), and the low-elevation tundra is above freezing (red). The largest expansion of ice cover occurs in the second iteration as noted above, which dramatically cools much of the region mostly due to the presence of snow and the albedo change. Also, note that the upper part of the ice cap also cools, likely due to increasing elevation, since 500 years of ice growth increased the surface elevation by 200 m. In the third iteration, the ice has retreated somewhat, which makes the low elevations warmer than in the second iteration. The integrated mass balance becomes more positive as seen in Fig. 2. As iterations of WRF and the ice model continue, we find that the area covered in ice cools, but the low-elevation tundra area surrounding the ice cap warms. The lapse rate effect due to increasing surface elevation progressively cools the ice cap, and below we explore why the low-elevation tundra warms in successive iterations.

Figure 5Melting in WRF in select iterations: (a) melting in the Iteration 1 WRF run. Differences in melting between Iteration 1 and the (b) second and (c) last iterations.

The mass balance is made up of ablation and precipitation. When precipitation is larger than ablation, the ice cap grows. We are interested in how both metrics change over the iterations of the coupled WRF–ice sheet models. We start by examining the melting over the summer, and how it relates to the temperature discussed above. Figure 5a depicts the yearly melt in the Iteration 1 simulation. In general, the high-elevation parts of the region experience no melting during the summer due to their subfreezing temperatures. As the ice cover expands, the highest values of melting are located near the edges of the ice cap. As expected, melting decreases as the surface elevation increases.

Figure 6Annual precipitation in WRF in select iterations: (a) the precipitation in the Iteration 1 WRF run. Precipitation differences between Iteration 1 and the (b) second and (c) last iterations.

Next, we look at the precipitation component of the mass balance. In this region, precipitation increases with elevation due to the orography (Fig. 6). For glacial inception, continued snowfall is critical in addition to low ablation. Figure 6b and c depict the differences in precipitation between the Iteration 1 simulation and ice model iterations. We see an elevation desert effect that reduces precipitation by 200 mm year−1 or ∼20 % of the annual precipitation. The highest elevations receive less precipitation as the iterations go on. However, at the edges of the ice cap, precipitation increases, because the change in elevation and surface slope triggers orographic precipitation there. It is interesting that precipitation increases with ice growth in the Barnes region of the ice cap, which may be a necessary condition for glacial inception to occur. The accumulation increases on the western side of Baffin Island, which is also the windward side, and as noted by , a characteristic of the inception process. This general increase in precipitation near the margins agrees with the maps of the mass balance (Fig. 2) indicating the greatest mass balances do not occur on the peaks of the island but at about 1200 m elevation.

Figure 7The summer cloud radiative forcing (CRF) in select iterations: (a) the CRF in the Iteration 1 WRF run. Differences in CRF between Iteration 1 and the (b) second and (c) last iterations.

Finally, we are interested in the role that clouds may play in glacial inception. and found clouds to be negative feedback during glacial inception. When insolation is changed from present day to that appropriate for the time of glacial inception at 115 kya, noted that less low clouds form during glacial inception, while found that cloud properties change such that relative to present day, less shortwave radiation is blocked from the surface. Our objective here is to examine this feedback during the integration of our asynchronously coupled model. In Fig. 7, the net cloud radiative forcing (CRF) in the Iteration 1 simulation is negative during the summer because the clouds block shortwave radiation from reaching the surface. However, as the ice cover expands, the clouds block less shortwave radiation, illustrated by the positive net shortwave CRF difference between later iterations and the Iteration 1 run: the Iteration 1 simulation net CRF is more negative in all parts of the domain. The albedo of increasing snow cover accounts for less negative CRF over newly glaciated areas, but the change to the shortwave CRF in the rest of the domain is due to changes in cloud cover and properties.

Figure 8Changes in 500 mbar geopotential height in the parent domain in select iterations with wind vectors. Wind speeds are 6 m s−1 over Baffin Island. (a) The 500 mbar geopotential height for the Iteration 1 simulation with average large-scale flow during the summer and (b–c) 500 mbar geopotential height differences during the summer between Iteration 1 and Iterations 2 and 10. Wind anomalies of 1 m s−1 are overlaid.

## 3.2 Interaction of Baffin Island inception and atmospheric circulation

We find the significant warming over Baffin Island to be surprising, as one might expect the expanded ice cover over the peaks of Baffin Island to lead to cooling (which would amplify the positive mass balance and lead to a larger ice sheet). We now analyze the warming by looking at the regional geopotential height field, temperature, and winds in the outer domain. First, the geopotential height and winds at 500 mbar for the Iteration 1 simulation during the summer are depicted in Fig. 8. The winds over Baffin Island generally come from the northwest with a low pressure trough covering the whole island. A weak cyclonic flow also exists over the island and Baffin Bay. Next, we examine the subsequent iterations and their anomalies in the winds and geopotential height at 500 mbar in Fig. 8b–c. The pressure increases south of Baffin Island, while it decreases over Baffin Bay to the north, and an anomalous anticyclone forms over the Hudson Bay region. It seems that the larger ice caps on Baffin Island alter the circulation, and we will examine this hypothesis below.

Figure 9Summer 2 m temperature in the parent domain in select iterations: (a) the 2 m temperature for Iteration 1 and (b–c) 2 m temperature differences during the summer between Iteration 1 and iterations 2 and 10.

The 2 m temperature response in the outer domain (Fig. 9) is consistent with the changes in geopotential height. The Iteration 1 simulation is depicted in the upper left with subsequent panels showing differences in the following iterations. On the ice caps of Baffin Island, the temperature generally grows colder because the elevation has increased, and newly glaciated areas cool due to both the presence of ice and the elevation increase. As iterations progress, the non-glaciated points of the island warm, which is also seen in the inner domain. Additionally, warming occurs outside of Baffin Island in the Hudson Bay region. These areas of warming, coinciding with geopotential height anomalies, appear to be caused by the circulation response, with warm southerly advection by the anomalous anticyclone.

Figure 10Scatter plots of summer (JJA) maximum (a) temperature anomaly and (b) geopotential height anomaly over Baffin Island vs. maximum ice elevation in each iteration.

This warming appears to be a negative feedback on ice sheet growth, and it explains the mass balances becoming more negative at low elevations on Baffin Island through subsequent iterations (Fig. 3). With a region so close to freezing during the summer, even a degree warming can prevent perennial snow from occurring during the critical inception period. As a caveat, we note that with the elevation in the coarser outer domain being generally lower than the 20 km resolution domain, the temperatures on Baffin Island are much warmer than is likely realistic. For instance, we do not see any summer temperatures significantly below freezing as found in the inner domain with more realistic topography.

It is known that regional models can develop artificial anomalous circulations , and it is important to verify that the negative feedback identified here is not strongly affected by such artifacts. Figure 10 shows two scatter plots of the temperature anomaly and the geopotential anomaly as function of ice height for Iterations 2–10. The anomalous circulation and heating develop progressively with the ice growth, suggesting that the anomalous circulation is driven by the changing ice topography and is likely not an artifact. We further investigate the effect of topography in the following section. The geopotential anomaly increases with ice height, with the exception of the last iteration when the ice sheet topography becomes more of a plateau. Given that the temperature anomaly is still strong in this last iteration, the feedback is still valid through this last iteration. The response we find is consistent with previous studies which noted that the size of the mountain is critical to the orographic response and slight changes to topography can cause significant changes . We speculate that the physical mechanism underlying the topography-induced anticyclone might be related to mixing of shallow-water potential vorticity (PV) by eddies near a localized high point in a weak background flow. Lateral eddy diffusion near a high point would stir high-PV air (due to lower thicknesses) away from the high terrain. The anticyclone is created by the influx of lower-PV air from the surrounding environment. The anticyclone would theoretically increase in strength as the surface height increased, as seen in Fig. 10. This assumes eddy mixing to be unchanged, which may explain the weakened anomaly in Iteration 10. The anomaly seen here differs from the more commonly considered atmospheric situation where advection of PV by the mean flow dominates horizontal mixing. A more thorough consideration of this simplified problem, including eddy mixing, mean flow, and planetary vorticity gradients, is left as a subject for future work.

## 3.3 The impact of Baffin Island ice elevation on regional circulation

Next, we diagnose the roles of ice and topography on circulation changes by setting up sensitivity simulations based on the second iteration, which has changes in ice extent and surface elevation (Fig. 2). For our “no H” scenario, we keep the surface elevation the same as the Iteration 1 simulation but use the same ice extent as in the second iteration. For an “H only” scenario, we remove all ice, changing the land type classification to tundra. We then set the surface elevation to the ice height in the second iteration. These two additional simulations were run for 5 years with cold meteorology and 115 kya insolation, and the analysis focuses on the last 4 years.

Figure 11Changes in 2 m summer temperature in the parent domain for topography and ice sensitivity simulations. (a) The 2 m temperature for the Iteration 1 simulation averaged during the summer (JJA). (b) Differences during the summer between Iterations 1 and 2. (c) Differences between the Iteration 1 and “no H”. (d) Differences between Iteration 1 and “H only”.

In Fig. 11, we examine the temperature of the outer domain for four experiments. Again, the Iteration 1 simulation summer temperature and difference from the second iteration (having both ice and elevation changes) are shown, but we also include the “no H” and “H only” scenarios, in Fig. 11c and d, respectively. With “no H”, a significant regional cooling occurs, encompassing both Baffin Island and the Hudson Bay region, which is the pattern we had anticipated to cause rapid ice growth. However, with “H only”, warming on both Baffin Island and around Hudson Bay occurs. The geopotential height anomalies presented in Fig. S2 in the Supplement agree with the temperature patterns seen here. Thus, the cause of the negative feedback due to large-scale anticyclonic circulation and warming in our simulations is the change in topography. An interesting aspect of these last two simulations (“no H” and “H only”) is that if we sum the difference between their response and the Iteration 1 simulation, the result is very close to the difference between the Iteration 2 and the Iteration 1 runs, as seen in Fig. 12. This indicates that the responses to the ice topography and ice land type (albedo and other surface properties) are linear.

Figure 12(a) Differences during the summer 2 m temperature between the first and second iterations. (b) Sum of temperature differences from “no H” and “H only” scenarios, implying a linear response due to changes in ice cover and surface elevation

## 3.4 Glacial inception with changing orbital parameters

The WRF–ice model iterations performed previously used the minimum insolation from 115 kya, which is useful for diagnosing the relationship between the ice cover and atmosphere. However, that is obviously not realistic. Integrated insolation increased after its minimum at 115 kya, impacting temperature and glacial mass balance. We use the procedure detailed in the Appendix of to set the insolation in WRF for 115–111 kya. After each WRF run, the ice model is integrated for 500 years using the same mass balance procedure as detailed above. This set of simulations explores the mass balance that results from a combination of ice flow and insolation changes.

Figure 13 presents the mass balances calculated from WRF and used in each iteration of the ice model. Note the Iteration 1 simulation and second iteration plotted here are the same as presented in Fig. 3 because the change in insolation occurs at 114 kya, which is the third iteration in Fig. 13. The elevation desert effect is again observed, and the mass balances become more negative at low elevations, hindering lateral growth. Thus, these iterations progress in a similar manner as the case with constant insolation in that the lateral growth is slow though the height of the ice sheets does grow, but the high-elevation plateaus can be 200 m lower than the previous set of simulations with constant 115 kya. In this set of simulations with realistic insolation, the mass balances at low elevations are more negative, explaining the slower growth (Fig. 13).

During the summer on Baffin Island, there is again cooling from the presence of ice, but the unglaciated parts of Baffin Island warm considerably. As the summer insolation increases, the summer temperature of all of North America, not just the Hudson Bay region, increases as well. When we look at the precipitation in the domain, we observe a similar pattern to the previous set of simulations. Accumulation increases on the edges of the ice cap due to orographic precipitation. Then, at the top of the ice cap, the elevation desert effect kicks in. The ice cap continually receives less precipitation than in the Iteration 1 simulation. The increases in insolation cause continental-scale warming, and when combined with topography changes, the warming on Baffin Island can be over 5 K (Fig. 14). Thus, the realistic insolation changes cause further warming during inception, which makes the occurrence of inception even more difficult to explain and simulate.

Figure 13The mass balance binned by elevation used in simulations with realistic insolation from 115 to 111 kya, where the first and second iterations use 115 kya insolation, third and fourth use 114 kya insolation, and so forth to Iteration 10 with 111 kya.

## 3.5 Sensitivity analysis of the ice flow model

The combination of WRF and the ice flow model does not indicate ice sheets would grow to cover all of Baffin Island, which would likely be necessary to explain the sea-level drop found in the paleo-record. The change in ice surface elevation actually causes a negative feedback on inception by affecting the large-scale circulation: temperatures around the ice increase inhibiting ice expansion. This result suggests a closer examination of the ice model and its sensitivities. In particular, we vary the creep parameter A to see if rapid ice expansion is possible with different constant choices. We also explore how a successful inception – matching paleo-record estimates – could occur from our simulated mass balance function G(H*) and limiting the calving loss of ice in Hudson Bay. The ice model in these sensitivity experiments is now run for 10 kyrs, unlike our previous asynchronous simulations that ran for 5 kyrs total.

Figure 14The surface temperature in the Iteration 1 run, and differences for successive iterations with insolation changes: (a) 115kya, (b) 114 kya, and (c) 111 kya.

First, we identify if we have simulated a mass balance sufficient for a successful inception. From Fig. 3, we note that the mass balance in the Iteration 1 simulation is the most positive at high elevations and least negative at low elevations. All subsequent WRF runs simulate a negative mass balance at low surface elevations. Using the Iteration 1 mass balance, the ice model is integrated over 10 kyrs. Of course, this situation is not realistic as our set of simulations shows changes in mass balance due to circulation changes as the ice sheet grows. We prescribe the mass balance only on Baffin Island, which was the location of our inner domain, and this causes Baffin Island to be mostly covered by ice, corresponding to 1.75 m sea level drop over 10 kyrs. Ice model simulations forced with the mass balances from Iterations 2–10 never reach the western coast of Baffin Island in 10 kyrs.

Figure 15Sensitivity of the ice model to choice of constant creep parameter A(T), which is related to the temperature profile within the ice as detailed by . (a) A evaluated at 10 C or A(10 C); (b) A(0 C); (c) A(15 C).

As noted above, we choose the creep parameter (A(T)=A(−10C)) based on . We now test the sensitivity of the ice flow when A(T) is evaluated at warmer and colder temperatures, or A(0 C) and A(−15C), respectively. For these simulations, the mass balance used is from the Iteration 1 WRF simulations, and when run for 500 years, the ice model changes are minor. We confirm with WRF simulations that these slight changes in ice elevation from A(0 C) and A(−15C) do not significantly affect the mass balance in the second iteration; the ice still retreats. Thus, we want to examine how the creep parameter affects ice flow on longer timescales, which is why we integrate the ice model in time for 10 kyrs. We expect that a warmer temperature and thus larger A would allow the ice to spread more easily, leading to larger ice cover. Figure 15 shows that changing the creep parameter can increase or decrease the thickness of a Baffin Island ice cap by about 20 % but does not alter the ice extent much after 10 kyr. On average, A(−10C) causes the elevation to be 200 m higher than A(0 C). More calving occurs with A(0 C) ice reaching the coast sooner. We also examine the impact of temperatures at about −15C (Fig. 15c), which determined to be the basal temperature of the initial ice sheet growth. As expected, the ice flows more slowly and accumulates at high elevations. These sensitivities to the creep parameter, however, are small compared to the sensitivity of ice extent and thickness to the mass balance G(H*) as inferred from successive iterations of WRF.

We also question what would happen if we apply the Iteration 1 mass balance to a larger area, outside of Baffin Island and our chosen domain. Based on the temperature and precipitation patterns of the outer domain, the other side of Hudson Bay could have a similar mass balance, especially if the domain were resolved at a higher resolution. We do not include the Torngat Mountains on mainland Canada just to the south of Baffin Island, as this area appears to receive less precipitation and experience higher temperatures. This agrees with observations on the lack of ice today, indicating the possibility of later glacierization . The ice accumulation with the positive mass balance from our Iteration 1 WRF simulation on mainland Canada causes an equivalent sea level drop of 2.8 m in 10 kyrs.

Figure 16The result of running the ice model for 10 kyrs with the maximum mass balance, found in the first WRF iteration. (a) Allowing a positive mass balance on mainland Canada, surrounding all of Hudson Bay. (b) G(H*) extrapolated to ${H}^{*}<\mathrm{0}$ in Hudson Bay, which allows for ice flow instead of immediate calving of ice.

The previous simulations assume immediate calving over any ocean points, which makes sense for the kilometer-deep Baffin Bay. The topography on the eastern side of Baffin Island is steep, possibly preventing ice shelves from forming. However, Hudson Bay is only about 200 m deep. We turn off immediate calving for this shallow water body and instead prescribe a negative mass balance. We extend the regression G(H) linearly to the bathymetry of Hudson Bay. This allows for ice to flow slowly into a larger area surrounding Baffin Island, creating a fairly large ice sheet and corresponding to a sea level drop of 6.6 m. This last scenario gets close to more recent constraints of sea level drop from modeling simulations at around 100 kya of 10 m . This growth of the ice sheet corresponds a successful inception, though still on the slower side of observational constraints, suggesting that the negative feedback due to the response of the large-scale circulation to expanding ice cover is what prevents our asynchronously coupled model from achieving a successful inception.

4 Discussion

In this paper, we have attempted to simulate the last glacial inception beginning at 115 kya with a regional high-resolution configuration of WRF, coupled asynchronously with a simple ice flow model. A net positive mass balance occurs in the first WRF iteration simulation driven by 115 kya insolation, cold present-day (1986) boundary conditions, and present-day ice cover. If the mass balance from this first iteration was maintained, an ice sheet could cover all of Baffin Island after 10 kyr. However, updating the ice cover and surface elevation in WRF creates a negative feedback from changes in the large-scale circulation: summer temperatures over Baffin Island increase due to a regional anticyclonic circulation. This causes the mass balance over Baffin Island to become more negative at low elevations, stagnating further ice expansion in the ice flow model. This is reminiscent of the findings of and , who noted the continued growth of the ice sheet is slowed by changes in large-scale circulation that arise from atmosphere–topography interactions. While the studies of and found the negative feedback slowed ice growth mostly when the ice sheet was fairly developed, of a continental scale, the negative feedback analyzed here occurs upon inception, as soon as the ice caps over Baffin Island begin to grow. Previous studies could not resolve this feedback with their coarse-resolution GCMs, while our coupled ice sheet – regional atmospheric model enables us to examine it.

The flow over Baffin Island driven by the 1986 boundary conditions is consistently from the northwest during the summer, with weak speeds of about 6 m s−1 at 500 mbar. The increase in the ice surface elevation after the first iterations of the coupled model causes an anticyclonic response of about 1 m s−1 to develop to the south of Baffin Island, while a cyclone develops in the north, diverting the winds and bringing in warmer air from the south. This leads to net ablation and arrests the glaciation process. This circulation response is likely dependent on the initial flow (direction and speed) and topographic height . Thus, a different meteorology choice might lead to sustained accumulation. identified exceptionally cold conditions on Baffin Island and their varying large-scale circulation patterns. Generally, the cold summer circulation flows from the northwest or north–northwest. We examine the role of meteorology directly using boundary conditions that differ from those used in our cold meteorology simulations. For this purpose, WRF is driven using 1979 boundary conditions, which is an average year in regards to temperature and precipitation, and the flow over Baffin Island generally comes from the west. WRF simulates a negative mass balance for this year, and the ice sheet retreats to an elevation above 1400 m when the ice model is forced with this mass balance for 500 years. In contrast, the ice cap from our average 1979 meteorology at a 4 km resolution was in approximate equilibrium, leaning towards slight expansion . Net ablation occurs in the 20 km resolution because of the lower topography. It is possible that during inception the circulation pattern was different and allowed for both the initial ice expansion and continued cold conditions. WRF did not simulate such a pattern in this study, nor did the few “cold” meteorological conditions in the ERA-Interim reanalysis.

We find significant warming on Baffin Island and the surrounding Hudson Bay region in response to the initial ice growth due to anticyclonic flow. The formation of anticyclones over expanding ice sheets has been noted to allow cold northerly winds to flow over the ice sheet, enhancing inception . In particular, small ice sheets cause a wave response that enhances accumulation, while stationary waves induce less cooling as ice sheets increase in size . In contrast, the anticyclone in our simulations forms south of the ice growing on Baffin Island, bringing significant warming even over these relatively small mountain glaciers. Similar to our study, found the Hudson Bay warmed in their coupled atmosphere–ice simulations due to anticyclonic flow. also found warming on the margins from anticyclones and reduced cloudiness. In both of these studies, though, the warming only slowed growth; the regional cooling due to the growing ice sheet kept the ice advancing. We find the heat fluxes into the region by the developing anticyclone circulation to completely overpower any regional cooling from the snow–albedo feedback, stagnating the ice sheet growth. We do observe the development of katabatic winds in response to ice sheet growth in our model, possibly contributing to the drying out of the region and decreasing cloud cover, but these winds do not cool the area outside of the glaciers sufficiently to allow for inception.

We do not find a positive feedback that would lead to a more positive mass balance due to the initial ice growth. However, some of our assumptions and modeling choices could be preventing such a positive feedback. For instance, our WRF simulations prescribe the SSTs and sea ice to those of 1986, but during glacial inception the SST could decrease and sea ice extent grow larger . Allowing the SST to cool has been shown to induce cooling over land, even if the atmosphere is relatively warm , while noted that storm tracks change with colder SSTs in a manner conducive for inception. In contrast, and argued that the oceans change slower than the atmosphere, so the oceans may remain warmer and ice free at the beginning of inception. argued that SSTs remained warm until 105 kya, when substantial ice sheets had formed. noted that the northern Atlantic remained close to modern-day temperatures through 115 kya until glacial stage 5d around 110 kya. argued that observations support a warm ocean, similar to an interglacial, during inception, though their focus is on the Southern Hemisphere. Some past GCM studies yielding some degree of glacial inception have also found Atlantic Ocean temperatures to remain similar to the present day .

Oceans may also affect moisture availability, not just heat fluxes. The amount of sea ice can impact accumulation over the ice sheet due to decreasing atmospheric moisture and evaporation from ocean, according to the temperature–precipitation feedback . Some studies argued that the cooling during inception did not cause sea ice to extend outside of Baffin Bay until 80 kya , meaning the Atlantic remained ice free, which was conducive to high precipitation. Both and found that including ocean feedbacks caused increased moisture fluxes conducive to inception, while and noted that storm frequency increased. Including ocean feedbacks could lead to increased moisture fluxes and thus precipitation, or it could lead to cooler temperatures decreasing melting, which could be more important for inception. With present-day fluxes prescribed at the boundaries, and the focus of this study being regional feedbacks, we keep SSTs at modern values for consistency, and effectively investigate the “warm” interglacial temperature ocean scenario, which does have support in observations .

We fix the CO2 concentration to 290 ppm, although CO2 decreased during inception. This decrease has been shown to influence the accumulation of Northern Hemisphere ice sheets , accounting for up to 50 % more ice . We also do not allow for dynamic vegetation, but on Baffin Island the tree canopy does not exist. Most of these studies dealing with dynamic vegetation refer to the treeline moving southward during glacial inception , enhancing ice growth, but on Baffin Island the vegetation is mostly tundra and small shrubs, which are not an impediment to snow accumulation as trees would be . However, the absence of dynamic vegetation means we cannot explore the cooling effects of albedo changes outside of Baffin Island due to vegetation changes. The southward treeline trend and the subsequent albedo impact, though, have been shown to occur later in the inception process , justifying our neglect of this feedback.

We choose to drive our model with meteorology conducive to a glaciation inception, and force the model with anomalously high precipitation and cool temperatures over Baffin Island in 1986. We are assuming that the planetary waves' response to the growing Laurentide ice sheet (Wells1983) was not significant until the size of the Laurentide was substantial, as has been found in other modeling studies of glacial inception . The regional flow induced by these boundary conditions in conjunction with ice growth causes the negative feedback discussed above. The strength of this negative regional feedback that prevents inception from proceeding suggests that a reorganization of the planetary-scale Northern Hemisphere circulation could occur as part of the inception process and counter the regional negative feedback. We cannot capture such a change in the general circulation with our regional model and prescribed boundary conditions, indicating a study of non-local processes could complement this study of local Baffin Island feedbacks. For instance, the Laurentide ice sheet and Fennoscandian ice sheet in Europe are thought to have begun growing at the same time, though the Laurentide increased in volume more rapidly . The European ice sheet may have influenced the formation of the Laurentide , but such a feedback is missing in our regional model. noted that changes in the stationary wave pattern would affect the downstream climate.

We limit our study to Baffin Island based on geological records of it being the nucleation of the Laurentide, but regional modeling of other Canadian mountain ranges could help identify how easily other areas glaciated, contributing to sea level drop. Growth of ice at high elevations simultaneously throughout eastern Canada is not impossible, which could all coalesce into the Laurentide. For instance, various ice sheet models have shown ice beginning to accumulate all over Canada . While many of these ice modes do not match geological records because they simulate the growth of ice in the west, ice growth on many of the more eastern mountain ranges could perhaps impact the large-scale circulation in a manner conducive to inception.

Though we do find net accumulation in our Iteration 1 simulation using a 20 km resolution, there are inherent issues that come from this resolution. The topography of Baffin Island is very rough but is greatly smoothed by the 20 km resolution. More importantly, the highest elevation peaks with 20 km resolution are ∼400 m below the peaks found in realistic Baffin Island topography. Our simulation with average meteorology demonstrates that resolution increases can easily switch the mass balance from net accumulation to net ablation. For cold meteorology, the resulting mass balance is similar to what we found with a 4 km resolution .

To get an average mass balance to drive the ice sheet model, WRF is run for 5 years; we use the average mass balance from the last 4 years. We find that the land model required a 1-year spin-up time, although often a spin-up time of more than a year is recommended . We tested running the model for longer than 10 years, but the simulations did not deviate much from years 2–5 of our simulations. showed that multiple ensemble members of an atmospheric model simulation can give different results on a timescale up to 50 years, but this may not be an issue here because we are driving our regional model with ideal inception boundary conditions that constrain the synoptic and planetary-scale flows .

Using the shallow ice approximation neglects important ice dynamics; in particular, it neglects lateral stresses and does not allow for the development of ice streams, which cause rapid ice flow. It has also been suggested that a finer grid size in the ice model could make for more accurate flow, as coarser resolutions have been shown to halt ice flow . As noted above, lateral expansion only occurs when substantial ice thickness exists in an adjacent grid cell. However, we find that simulations with a 4 km resolution (Fig. S6 in the Supplement) yield similar ice flow patterns to the 20 km results we present here. We are neglecting basal sliding that could occur as the ice base melts due to geothermal heat flux , but the initiating ice sheet was likely frozen at the bed , which makes our choice appropriate. Finally, we are not accounting for glacial isostatic adjustment. found that the isostatic adjustment was not critical for glacial inception, though including the adjustment revealed a small positive feedback. Therefore, contrary to conventional wisdom, these studies suggest that isostatic adjustment might actually favor glacial inception over Baffin Island, and may even be critical to the process by reducing the strength of negative elevation–circulation feedbacks on ice growth (as in the “no H” sensitivity test). Our binning procedure and the resulting surface mass balance (SMB) forcing routine, while simple, are consistent with the way many simple ice sheet models are forced. More importantly, the results of this procedure, showing that ice elevation causes a negative warming feedback, should be robust regardless of how the ice elevation was calculated, as it is a result of the atmospheric model itself. We found that finer bin sizes did not significantly affect the expansion of ice. However, it would be good to test alternative SMB routines, such as using mean monthly air temperatures and precipitation from the previous WRF integrations, and downscale them to the surface elevation of ice grid points.

Though we find negative feedbacks on the growth of ice sheets due to clouds and circulation preventing inception in our model, high-resolution regional modeling still has the potential to make progress on the inception problem. We have captured the mass balance on critical mountainous regions in our Iteration 1 run, which GCMs cannot simulate accurately. When combined with ice flow, the Iteration 1 mass balance causes a significant portion of Baffin Island and the surrounding area to become ice covered. Though the equivalent sea level drop is only 6 m in 10 kyrs, the geological constraints on sea level drop from the Laurentide are not tight, and the calculated sea level drop may be consistent with more recent sea level estimates . However, this scenario neglects the ice–elevation feedback on the atmosphere that we found significant. Given the geological importance of Baffin Island on glacial inception, additional work is needed to understand how Baffin Island cooled against two warming influences: the regional anticyclone circulation and increasing insolation from 115 kya onward. The warming feedback due to local circulation changes that we find could be mitigated by changes in the pattern of large-scale atmospheric circulation, including stationary and planetary waves. Using a circumpolar regional model could potentially capture the mass balance of the nucleation points of both the Laurentide and Fennoscandian ice sheets at a sufficiently high resolution, allowing an investigation in the downstream interplay of these two ice sheets. Another interesting study could involve exploring the possible impact of a smaller Greenland ice sheet upon exiting the last interglacial. Going a step further to more non-local feedbacks (such as SST changes) influenced by changing orbital parameters, a next step would be simulating the last glacial inception with a GCM and using the results as boundary conditions in a regional model.

Data availability
Data availability.

The scripts used and the model output are available upon request to the corresponding author.

Appendix A: Shallow ice model set-up

The simple ice sheet model is based off of , where ice flow is represented by a non-linear diffusion equation and the diffusivity contains a constant term called the creep factor. This parameter determines the ice viscosity and depends on the temperature of the ice , meaning it can be denoted A(T) as illustrated in Eq. (A1):

$\begin{array}{}\text{(A1)}& A\left(T\right)={A}_{\mathrm{0}}{e}^{-Q/RT},\end{array}$

if $T<-\mathrm{10}$C, then ${A}_{\mathrm{0}}=\mathrm{3.61}×{\mathrm{10}}^{-\mathrm{13}}$ Pa−3 s−1 and $Q=\mathrm{60}×{\mathrm{10}}^{\mathrm{3}}$ J mol−1, but if $T\ge -\mathrm{10}$C, then ${A}_{\mathrm{0}}=\mathrm{1.734}×{\mathrm{10}}^{\mathrm{3}}$ Pa−3 s−1 and $Q=\mathrm{139}×{\mathrm{10}}^{\mathrm{3}}$ J mol−1.

The surface of the ice is colder than the base, meaning that throughout a real ice sheet, temperature, and thus A(T), varies. In fact, the basal temperature often reaches 0 C for thick enough ice sheets. For an ice model using the shallow ice approximation, the standard procedure is to take the vertical average of A(T) in Eq. (A2):

$\begin{array}{}\text{(A2)}& A\left(T\right)=\frac{\mathrm{1}}{h}A\left(T\left(z\right)\right)\mathrm{d}z.\end{array}$

Since A(T(z)) varies by height, we need a vertical temperature profile within the ice sheet. We calculate it using the average temperature of the ice surface from WRF and a geothermal heat flux of 55 W m−1, assuming a linear vertical temperature profile. We then evaluate A(T) for this temperature profile as found in Eq. (A1) and take the vertical average. Typically, this vertical average is influenced by A(T) at the warmer temperatures at the base of the ice sheet, since A(T) varies exponentially with temperature. We calculate $A\left(T\right)=\mathrm{3.5}×{\mathrm{10}}^{-\mathrm{25}}$ Pa−3 s−1, which is the same as A(−10C) from . We find that using Robin's solution instead of a linear profile confirms that the basal temperature is below freezing, consistent with our shallow ice model formulation. The basal temperature we use is not far from about −10C, also consistent with the basal ice temperature in other studies of glacial inception .

## A1 Finite difference approximation

We can use the flux form for the finite difference approximation. We use forward Euler and second-order centered difference as seen in Fig. A3, where n is the time step, i is the node in the x direction, and j is the node in the y direction.

$\begin{array}{ll}\text{(A3)}& & \frac{{H}_{i,j}^{n+\mathrm{1}}-{H}_{i,j}^{n}}{\mathrm{\Delta }t}& \phantom{\rule{1em}{0ex}}=\frac{{D}_{i+\frac{\mathrm{1}}{\mathrm{2}},j}^{n}\left(\frac{{H}_{i+\mathrm{1},j}^{*n}-{H}_{i,j}^{*n}}{\mathrm{\Delta }{x}_{i}}\right)+{D}_{i-\frac{\mathrm{1}}{\mathrm{2}},j}^{n}\left(\frac{{H}_{i,j}^{*n}-{H}_{i-\mathrm{1},j}^{*n}}{\mathrm{\Delta }{x}_{i-\mathrm{1}}}\right)}{\frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{\Delta }{x}_{i}+\mathrm{\Delta }{x}_{i-\mathrm{1}}\right)}\\ & \phantom{\rule{1em}{0ex}}+\frac{{D}_{i,j+\frac{\mathrm{1}}{\mathrm{2}}}^{n}\left(\frac{{H}_{i,j+\mathrm{1}}^{*n}-{H}_{i,j}^{*n}}{\mathrm{\Delta }{y}_{j}}\right)+{D}_{i,j-\frac{\mathrm{1}}{\mathrm{2}}}^{n}\left(\frac{{H}_{i,j}^{*n}-{H}_{i,j-\mathrm{1}}^{*n}}{\mathrm{\Delta }{y}_{j-\mathrm{1}}}\right)}{\frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{\Delta }{y}_{j}+\mathrm{\Delta }{y}_{j-\mathrm{1}}\right)}+G\end{array}$

The numerical expression for the diffusivity is presented in Eq. (A7):

$\begin{array}{ll}\text{(A4)}& & {D}_{i+\frac{\mathrm{1}}{\mathrm{2}},j}^{n}=A{H}_{i+\frac{\mathrm{1}}{\mathrm{2}},j}^{m+\mathrm{2}}\left[{\left(\frac{{H}_{i+\mathrm{1},j}^{*}-{H}_{i,j}^{*}}{\mathrm{\Delta }{x}_{i}}\right)}^{\mathrm{2}}\right& \phantom{\rule{1em}{0ex}}{+{\left(\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{{H}_{i+\mathrm{1},j+\mathrm{1}}^{*}-{H}_{i+\mathrm{1},j-\mathrm{1}}^{*}}{\mathrm{\Delta }{y}_{j-\mathrm{1}}+\mathrm{\Delta }{y}_{j}}+\frac{{H}_{i,j+\mathrm{1}}^{*}-{H}_{i,j-\mathrm{1}}^{*}}{\mathrm{\Delta }{y}_{j-\mathrm{1}}+\mathrm{\Delta }{y}_{j}}\right)\right)}^{\mathrm{2}}]}^{\left(m-\mathrm{1}\right)/\mathrm{2}}\\ \text{(A5)}& & {D}_{i-\frac{\mathrm{1}}{\mathrm{2}},j}^{n}=A{H}_{i-\frac{\mathrm{1}}{\mathrm{2}},j}^{m+\mathrm{2}}\left[{\left(\frac{{H}_{i,j}^{*}-{H}_{i-\mathrm{1},j}^{*}}{\mathrm{\Delta }{x}_{i-\mathrm{1}}}\right)}^{\mathrm{2}}\right& \phantom{\rule{1em}{0ex}}{+{\left(\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{{H}_{i,j+\mathrm{1}}^{*}-{H}_{i,j-\mathrm{1}}^{*}}{\mathrm{\Delta }{y}_{j-\mathrm{1}}+\mathrm{\Delta }{y}_{j}}+\frac{{H}_{i-\mathrm{1},j+\mathrm{1}}^{*}-{H}_{i-\mathrm{1},j-\mathrm{1}}^{*}}{\mathrm{\Delta }{y}_{j-\mathrm{1}}+\mathrm{\Delta }{y}_{j}}\right)\right)}^{\mathrm{2}}]}^{\left(m-\mathrm{1}\right)/\mathrm{2}}\\ \text{(A6)}& & {D}_{i,j+\frac{\mathrm{1}}{\mathrm{2}}}^{n}=A{H}_{i,j+\frac{\mathrm{1}}{\mathrm{2}}}^{m+\mathrm{2}}\left[\left(\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{{H}_{i+\mathrm{1},j+\mathrm{1}}^{*}-{H}_{i-\mathrm{1},j+\mathrm{1}}^{*}}{\mathrm{\Delta }{x}_{i-\mathrm{1}}+\mathrm{\Delta }{x}_{i}}\right\right\right& {{\phantom{\rule{1em}{0ex}}+\frac{{H}_{i+\mathrm{1},j}^{*}-{H}_{i-\mathrm{1},j}^{*}}{\mathrm{\Delta }{x}_{i-\mathrm{1}}+\mathrm{\Delta }{x}_{i}}))}^{\mathrm{2}}+{\left(\frac{{H}_{i,j+\mathrm{1}}^{*}-{H}_{i,j}^{*}}{\mathrm{\Delta }{y}_{j}}\right)}^{\mathrm{2}}]}^{\left(m-\mathrm{1}\right)/\mathrm{2}}\\ \text{(A7)}& & {D}_{i,j-\frac{\mathrm{1}}{\mathrm{2}}}^{n}=A{H}_{i,j-\frac{\mathrm{1}}{\mathrm{2}}}^{m+\mathrm{2}}\left[\left(\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{{H}_{i+\mathrm{1},j+\mathrm{1}}^{*}-{H}_{i-\mathrm{1},j+\mathrm{1}}^{*}}{\mathrm{\Delta }{x}_{i-\mathrm{1}}+\mathrm{\Delta }{x}_{i}}\right\right\right& {{\phantom{\rule{1em}{0ex}}+\frac{{H}_{i+\mathrm{1},j}^{*}-{H}_{i-\mathrm{1},j}^{*}}{\mathrm{\Delta }{x}_{i-\mathrm{1}}+\mathrm{\Delta }{x}_{i}}))}^{\mathrm{2}}+{\left(\frac{{H}_{i,j}^{*}-{H}_{i,j-\mathrm{1}}^{*}}{\mathrm{\Delta }{y}_{j-\mathrm{1}}}\right)}^{\mathrm{2}}]}^{\left(m-\mathrm{1}\right)/\mathrm{2}}.\end{array}$

Finally, our use of flux form requires Eq. (A11).

$\begin{array}{}\text{(A8)}& {H}_{i+\frac{\mathrm{1}}{\mathrm{2}},j}^{m+\mathrm{2}}={\left(\frac{{H}_{i+\mathrm{1},j}+{H}_{i,j}}{\mathrm{2}}\right)}^{m+\mathrm{2}}\text{(A9)}& {H}_{i-\frac{\mathrm{1}}{\mathrm{2}},j}^{m+\mathrm{2}}={\left(\frac{{H}_{i,j}+{H}_{i-\mathrm{1},j}}{\mathrm{2}}\right)}^{m+\mathrm{2}}\text{(A10)}& {H}_{i,j+\frac{\mathrm{1}}{\mathrm{2}}}^{m+\mathrm{2}}={\left(\frac{{H}_{i,j+\mathrm{1}}+{H}_{i,j}}{\mathrm{2}}\right)}^{m+\mathrm{2}}\text{(A11)}& {H}_{i,j-\frac{\mathrm{1}}{\mathrm{2}}}^{m+\mathrm{2}}={\left(\frac{{H}_{i,j}+{H}_{i,j-\mathrm{1}}}{\mathrm{2}}\right)}^{m+\mathrm{2}}\end{array}$

This model uses a time step of 1 year. At each time step, a new ice thickness is calculated based the equations above. Then, with the new ice thickness and basal topography known, a new surface elevation is calculated (${H}^{*}=B+H$).

The edge and center of the ice sheet can cause D=0, which is unrealistic behavior. Thus, we set a minimum D, but as noted by and confirmed in our simulations, the ice model does not depend strongly on this choice.

We also characterize the mass balance according to elevation, where G is actually G(H*). Thus, the mass balance is also updated after each time step, since it depends on H*. For our main set of simulations, we characterized the mass balance using elevation bins.

Supplement
Supplement.

Author contributions
Author contributions.

LB performed the research work; all authors contributed to the planning, interpretation, and writing of manuscript.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

This work was funded by the Harvard University Climate Change Solutions Fund and by the NSF P2C2 program, grant OCE-1602864. Timothy Cronin was supported by NSF grant AGS-1740533. Eli Tziperman thanks the Weizmann Institute for its hospitality during parts of this work. We would like to acknowledge high-performance computing support from Yellowstone (ark:/85065/d7wd3xhc) provided by NCAR's Computational and Information Systems Laboratory, sponsored by the National Science Foundation. We would like to thank two anonymous reviewers and the handling editor, Prof. Marit-Solveig Seidenkrantz, for their guidance in refining this paper. We also thank Alexander Robel for helpful conversations on the ice model.

Edited by: Marit-Solveig Seidenkrantz
Reviewed by: two anonymous referees

References

Abe-Ouchi, A., Segawa, T., and Saito, F.: Climatic Conditions for modelling the Northern Hemisphere ice sheets throughout the ice age cycle, Clim. Past, 3, 423–438, https://doi.org/10.5194/cp-3-423-2007, 2007. a

Abe-Ouchi, A., Saito, F., Kawamura, K., Raymo, M. E., Okuno, J., Takahashi, K., and Blatter, H.: Insolation-driven 100,000-year glacial cycles and hysteresis of ice-sheet volume, Nature, 500, 190–193, 2013. a

Andrews, J. and Barry, R.: Glacial inception and disintegration during the last glaciation, Ann. Rev. Earth Pl. Sc., 6, 205–228, 1978. a

Andrews, J. and Mahaffy, M.: Growth rate of the Laurentide Ice Sheet and sea level lowering (with emphasis on the 115,000 BP sea level low), Quaternary Res., 6, 167–183, 1976. a, b

Andrews, J. and Miller, G.: Quarternary history of northern Cumberland peninsula, Baffin Island, NWT, Canada: Part IV: Maps of the present glaciation limits and lowest equilibrium line altitude for north and south Baffin Island, Arctic and Alpine Research, 45–59, 1972. a

Beghin, P., Charbit, S., Dumas, C., Kageyama, M., Roche, D. M., and Ritz, C.: Interdependence of the growth of the Northern Hemisphere ice sheets during the last glaciation: the role of atmospheric circulation, Clim. Past, 10, 345–358, https://doi.org/10.5194/cp-10-345-2014, 2014. a

Birch, L., Cronin, T., and Tziperman, E.: Glacial Inception on Baffin Island: The Role of Insolation, Meteorology, and Topography, J. Climate, 30, 4047–4064, 2017. a, b, c, d, e, f, g, h, i, j, k, l, m, n

Bonelli, S., Charbit, S., Kageyama, M., Woillez, M.-N., Ramstein, G., Dumas, C., and Quiquet, A.: Investigating the evolution of major Northern Hemisphere ice sheets during the last glacial-interglacial cycle, Clim. Past, 5, 329–345, https://doi.org/10.5194/cp-5-329-2009, 2009. a, b, c

Bromwich, D. H., Toracinta, E. R., and Wang, S.-H.: Meteorological perspective on the initiation of the Laurentide Ice Sheet, Quatern. Int., 95, 113–124, 2002a. a, b, c, d

Bromwich, D. H., Wang, S.-H., and Monaghan, A. J.: ERA-40 representation of the Arctic atmospheric moisture budget, in: ECMWF Workshop on Reanalysis, ECMWF Re-Analysis Proj. Rep. Ser, Citeseer, 3, 287–298, 2002b. a

Bromwich, D. H., Toracinta, E. R., Oglesby, R. J., Fastook, J. L., and Hughes, T. J.: LGM summer climate on the southern margin of the Laurentide Ice Sheet: wet or dry?, J. Climate, 18, 3317–3338, 2005. a

Calov, R., Ganopolski, A., Petoukhov, V., Claussen, M., Brovkin, V., and Kubatzki, C.: Transient simulation of the last glacial inception. Part II: sensitivity and feedback analysis, Clim. Dynam., 24, 563–576, 2005. a

Cassano, J. J., Higgins, M. E., and Seefeldt, M. W.: Performance of the Weather Research and Forecasting model for month-long pan-Arctic simulations, Mon. Weather Rev., 139, 3469–3488, 2011. a, b, c

Chen, S.-H. and Lin, Y.-L.: Orographic effects on a conditionally unstable flow over an idealized three-dimensional mesoscale mountain, Meteorol. Atmos. Phys., 88, 1–21, 2005. a, b, c

Clark, P., Clague, J., Curry, B. B., Dreimanis, A., Hicock, S., Miller, G., Berger, G., Eyles, N., Lamothe, M., Miller, B., Mott, R. J., Oldale R. N., Stea, R. R., Szabo, J. P., Thorleifson, L. H., and Vincent, J. S.: Initiation and development of the Laurentide and Cordilleran ice sheets following the last interglaciation, Quaternary Sci. Rev., 12, 79–114, 1993. a, b, c

Cook, K. H. and Held, I. M.: Stationary waves of the ice age climate, J. Climate, 1, 807–819, 1988. a, b, c

Cortijo, E., Duplessy, J., Labeyrie, L., Leclaire, H., Duprat, J., and Van Wearing, T.: Eemian cooling in the Norwegian Sea and North Atlantic ocean preceding continental ice-sheet growth, Nature, 372, 446–449, 1994. a, b

Cosgrove, B. A., Lohmann, D., Mitchell, K. E., Houser, P. R., Wood, E. F., Schaake, J. C., Robock, A., Sheffield, J., Duan, Q., Luo, L., Wayne Higgins, R., Pinker, R. T., and Tarpley, J. D.: Land surface model spin-up behavior in the North American Land Data Assimilation System (NLDAS), J. Geophys. Res.-Atmos., 108, https://doi.org/10.1029/2002JD003316, 2003. a

Creveling, J. R., Mitrovica, J. X., Clark, P. U., Waelbroeck, C., and Pico, T.: Predicted bounds on peak global mean sea level during marine isotope stages 5a and 5c, Quaternary Sci. Rev., 163, 193–208, 2017. a, b, c

Cuffey, K. M. and Paterson, W. S. B.: The physics of glaciers, Academic Press, 2010. a, b, c, d, e, f, g

Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M., McNally, A. P., Monge-Sanz, B. M., Morcrette, J.-J., Park, B.-K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J.-N., and Vitart, F.: The ERA-Interim reanalysis: configuration and performance of the data assimilation system, Q. J. Roy. Meteor. Soc., 137, 553–597, https://doi.org/10.1002/qj.828, 2011. a

Denis, B., Laprise, R., Caya, D., and Côté, J.: Downscaling ability of one-way nested regional climate models: the Big-Brother Experiment, Clim. Dynam., 18, 627–646, 2002. a

de Noblet, N. I., Prentice, I. C., Joussaume, S., Texier, D., Botta, A., and Haxeltine, A.: Possible role of atmosphere-biosphere interactions in triggering the Last Glaciation, Geophys. Res. Lett., 23, 3191–3194, 1996. a

Deser, C., Phillips, A. S., Alexander, M. A., and Smoliak, B. V.: Projecting North American climate over the next 50 years: Uncertainty due to internal variability, J. Climate, 27, 2271–2296, 2014. a

Dong, B. and Valdes, P. J.: Sensitivity studies of Northern Hemisphere glaciation using an atmospheric general circulation model, J. Climate, 8, 2471–2496, 1995. a

Drusch, M., Vasiljevic, D., and Viterbo, P.: ECMWF's global snow analysis: Assessment and revision based on satellite observations, J. Appl. Meteorol., 43, 1282–1294, 2004. a

Dutra, E., Balsamo, G., Viterbo, P., Miranda, P. M., Beljaars, A., Schär, C., and Elder, K.: An improved snow scheme for the ECMWF land surface model: Description and offline validation, J. Hydrometeorol., 11, 899–916, 2010. a

Farrell, B. F. and Abbot, D. S.: A mechanism for dust-induced destabilization of glacial climates, Clim. Past, 8, 2061–2067, https://doi.org/10.5194/cp-8-2061-2012, 2012. a

Fréchette, B., Wolfe, A. P., Miller, G. H., Richard, P. J., and de Vernal, A.: Vegetation and climate of the last interglacial on Baffin Island, Arctic Canada, Palaeogeogr. Palaeocl., 236, 91–106, 2006. a

Ganopolski, A. and Calov, R.: The role of orbital forcing, carbon dioxide and regolith in 100 kyr glacial cycles, Clim. Past, 7, 1415–1425, https://doi.org/10.5194/cp-7-1415-2011, 2011. a

Ganopolski, A., Calov, R., and Claussen, M.: Simulation of the last glacial cycle with a coupled climate ice-sheet model of intermediate complexity, Clim. Past, 6, 229–244, https://doi.org/10.5194/cp-6-229-2010, 2010. a, b, c, d

Gardner, A. S. and Sharp, M.: Influence of the Arctic circumpolar vortex on the mass balance of Canadian high Arctic glaciers, J. Climate, 20, 4586–4598, 2007. a

Ghil, M.: Cryothermodynamics: The chaotic dynamics of paleoclimate, Physica D, 77, 130–159, 1994. a

Gildor, H. and Tziperman, E.: Sea ice as the glacial cycles' climate switch: Role of seasonal and orbital forcing, Paleoceanography, 15, 605–615, 2000. a, b, c

Gildor, H. and Tziperman, E.: A sea ice climate switch mechanism for the 100-kyr glacial cycles, J. Geophys. Res.-Ocean., 106, 9117–9133, 2001. a, b, c, d

Glisan, J. M. and Gutowski, W. J.: WRF summer extreme daily precipitation over the CORDEX Arctic, J. Geophys. Res.-Atmos., 119, 1720–1732, 2014. a

Glisan, J. M., Gutowski Jr, W. J., Cassano, J. J., and Higgins, M. E.: Effects of spectral nudging in WRF on Arctic temperature and precipitation simulations, J. Climate, 26, 3985–3999, 2013. a, b

Goñi, M. S., Loutre, M., Crucifix, M., Peyron, O., Santos, L., Duprat, J., Malaizé, B., Turon, J.-L., and Peypouquet, J.-P.: Increasing vegetation and climate gradient in Western Europe over the Last Glacial Inception (122–110 ka): data-model comparison, Earth Planet. Sc. Lett., 231, 111–130, 2005. a, b, c, d

Goodman, J. C. and Pierrehumbert, R. T.: Glacial flow of floating marine ice in “Snowball Earth”, J. Geophys. Res.-Ocean., 108, https://doi.org/10.1029/2002JC001471, 2003. a

Gregory, J. M., Browne, O. J. H., Payne, A. J., Ridley, J. K., and Rutt, I. C.: Modelling large-scale ice-sheet-climate interactions following glacial inception, Clim. Past, 8, 1565–1580, https://doi.org/10.5194/cp-8-1565-2012, 2012. a, b, c, d, e

Grell, G. A. and Dévényi, D.: A generalized approach to parameterizing convection combining ensemble and data assimilation techniques, Geophys. Res. Lett., 29, 1–38, 2002. a

Harvey, L. D.: Milankovitch forcing, vegetation feedback, and North Atlantic deep-water formation, J. Climate, 2, 800–815, 1989. a

Herrington, A. R. and Poulsen, C. J.: Terminating the Last Interglacial: The role of ice sheet–climate feedbacks in a GCM asynchronously coupled to an ice sheet model, J. Climate, 25, 1871–1882, 2011. a, b, c, d, e, f, g

Hines, K. M., Bromwich, D. H., Bai, L.-S., Barlage, M., and Slater, A. G.: Development and testing of Polar WRF. Part III: Arctic land, J. Climate, 24, 26–48, 2011. a

Hong, S.-Y. and Lim, J.-O. J.: The WRF single-moment 6-class microphysics scheme (WSM6), Asia-Pac. J. Atmos. Sci., 42, 129–151, 2006. a

Iacono, M. J., Delamere, J. S., Mlawer, E. J., Shephard, M. W., Clough, S. A., and Collins, W. D.: Radiative forcing by long-lived greenhouse gases: Calculations with the AER radiative transfer models, J. Geophys. Res.-Atmos., 113, https://doi.org/10.1029/2008JD009944, 2008. a

Ives, J.: Indications of Recent Extensive Glacierization in North-Central Baffin Island, NWT, J. Glaciol., 4, 197–205, 1962. a

Ives, J. D., Andrews, J. T., and Barry, R. G.: Growth and decay of the Laurentide Ice Sheet and comparisons with Fenno-Scandinavia, Naturwissenschaften, 62, 118–125, 1975. a

Jackson, C.: Sensitivity of stationary wave amplitude to regional changes in Laurentide ice sheet topography in single-layer models of the atmosphere, J. Geophys. Res.-Atmos., 105, 24443–24454, 2000. a

Jackson, C. and Broccoli, A.: Orbital forcing of Arctic climate: mechanisms of climate response and implications for continental glaciation, Clim. Dynam., 21, 539–557, 2003. a, b, c

Janjic, Z. I.: The step-mountain eta coordinate model: Further developments of the convection, viscous sublayer, and turbulence closure schemes, Mon. Weather Rev., 122, 927–945, 1994. a

Jochum, M., Jahn, A., Peacock, S., Bailey, D. A., Fasullo, J. T., Kay, J., Levis, S., and Otto-Bliesner, B.: True to Milankovitch: glacial inception in the new Community Climate System Model, J. Climate, 25, 2226–2239, 2012. a, b, c, d

Kageyama, M. and Valdes, P. J.: Impact of the North American ice-sheet orography on the Last Glacial Maximum eddies and snowfall, Geophys. Res. Lett., 27, 1515–1518, 2000. a

Kageyama, M., Charbit, S., Ritz, C., Khodri, M., and Ramstein, G.: Quantifying ice-sheet feedbacks during the last glacial inception, Geophys. Res. Lett., 31, https://doi.org/10.1029/2004GL021339, 2004. a

Källén, E., Crafoord, C., and Ghil, M.: Free oscillations in a climate model with ice-sheet dynamics, J. Atmos. Sci., 36, 2292–2303, 1979. a

Kaspar, F., Spangehl, T., and Cubasch, U.: Northern hemisphere winter storm tracks of the Eemian interglacial and the last glacial inception, Clim. Past, 3, 181–192, https://doi.org/10.5194/cp-3-181-2007, 2007. a

Khodri, M., Leclainche, Y., Ramstein, G., Braconnot, P., Marti, O., and Cortijo, E.: Simulating the amplification of orbital forcing by ocean feedbacks in the last glaciation, Nature, 410, 570–574, 2001. a, b, c, d

Kilpeläinen, T., Vihma, T., and Ólafsson, H.: Modelling of spatial variability and topographic effects over Arctic fjords in Svalbard, Tellus A, 63, 223–237, 2011. a

Kleman, J., Fastook, J., and Stroeven, A. P.: Geologically and geomorphologically constrained numerical model of Laurentide Ice Sheet inception and build-up, Quatern. Int., 95, 87–98, 2002. a

Kleman, J., Fastook, J., Ebert, K., Nilsson, J., and Caballero, R.: Pre-LGM Northern Hemisphere ice sheet topography, Clim. Past, 9, 2365–2378, https://doi.org/10.5194/cp-9-2365-2013, 2013. a, b

Kubatzki, C., Claussen, M., Calov, R., and Ganopolski, A.: Sensitivity of the last glacial inception to initial and surface conditions, Clim. Dynam., 27, 333–344, 2006. a, b

Kurtz, N. T. and Farrell, S. L.: Large-scale surveys of snow depth on Arctic sea ice from Operation IceBridge, Geophys. Res. Lett., 38, https://doi.org/10.1029/2011GL049216, 2011. a, b

Lambeck, K. and Chappell, J.: Sea level change through the last glacial cycle, Science, 292, 679–686, 2001. a

Lee, W.-H. and North, G.: Small ice cap instability in the presence of fluctuations, Clim. Dynam., 11, 242–246, 1995. a

Lehner, F., Born, A., Raible, C. C., and Stocker, T. F.: Amplified inception of European Little Ice Age by sea ice–ocean–atmosphere feedbacks, J. Climate, 26, 7586–7602, 2013. a

Li, C., Battisti, D. S., Schrag, D. P., and Tziperman, E.: Abrupt climate shifts in Greenland due to displacements of the sea ice edge, Geophys. Res. Lett., 32, https://doi.org/10.1029/2005GL023492, 2005. a, b

Liakka, J. and Nilsson, J.: The impact of topographically forced stationary waves on local ice-sheet climate, J. Glaciol., 56, 534–544, 2010. a

Liakka, J., Nilsson, J., and Löfverström, M.: Interactions between stationary waves and ice sheets: linear versus nonlinear atmospheric response, Clim. Dynam., 38, 1249–1262, 2012. a

Löfverström, M., Caballero, R., Nilsson, J., and Kleman, J.: Evolution of the large-scale atmospheric circulation in response to changing ice sheets over the last glacial cycle, Clim. Past, 10, 1453–1471, https://doi.org/10.5194/cp-10-1453-2014, 2014. a, b

Marshall, S. and Clarke, G.: Ice sheet inception: subgrid hypsometric parameterization of mass balance in an ice sheet model, Clim. Dynam., 15, 533–550, 1999. a, b, c

Marshall, S. J. and Clark, P. U.: Basal temperature evolution of North American ice sheets and implications for the 100-kyr cycle, Geophys. Res. Lett., 29, 1–67, 2002. a

Marshall, S. J., Tarasov, L., Clarke, G. K., and Peltier, W. R.: Glaciological reconstruction of the Laurentide Ice Sheet: physical processes and modelling challenges, Can. J. Earth Sci., 37, 769–793, 2000. a

Marshall, S. J., James, T. S., and Clarke, G. K.: North American ice sheet reconstructions at the Last Glacial Maximum, Quaternary Sci. Rev., 21, 175–192, 2002. a, b

Meissner, K., Weaver, A., Matthews, H., and Cox, P.: The role of land surface dynamics in glacial inception: a study with the UVic Earth System Model, Clim. Dynam., 21, 515–537, 2003. a, b

Milankovitch, M.: Kanon der Erdebestrahlung und seine Anwendung auf das Eiszeitenproblem, Königlich Serbische Akademie, 1941. a

Niu, G.-Y., Yang, Z.-L., Mitchell, K. E., Chen, F., Ek, M. B., Barlage, M., Kumar, A., Manning, K., Niyogi, D., Rosero, E., Tewari, M., and Youlong, X.: The community Noah land surface model with multiparameterization options (Noah-MP): 1. Model description and evaluation with local-scale measurements, J. Geophys. Res.-Atmos., 116, https://doi.org/10.1029/2010JD015139, 2011. a

Oerlemans, J.: Some basic experiments with a vertically-integrated ice sheet model, Tellus, 33, 1–11, 1981. a, b, c

Oerlemans, J.: On glacial inception and orography, Quatern. Int., 95, 5–10, 2002. a, b

Otieno, F. O. and Bromwich, D. H.: Contribution of Atmospheric Circulation to Inception of the Laurentide Ice Sheet at 116 kyr BP*, J. Climate, 22, 39–57, 2009. a, b

Otieno, F. O., Bromwich, D. H., and Oglesby, R.: Atmospheric circulation anomalies due to 115 kyr BP climate forcing dominated by changes in the North Pacific Ocean, Clim. Dynam., 38, 815–835, 2012. a, b, c, d

Otterman, J., Chou, M., and Arking, A.: Effects of nontropical forest cover on climate, J. Clim. Appl. Meteorol., 23, 762–767, 1984. a

Petit, J.-R., Jouzel, J., Raynaud, D., Barkov, N. I., Barnola, J.-M., Basile, I., Bender, M., Chappellaz, J., Davis, M., Delaygue, G., Kotlyakov, V. M., Legrand, M., Lipenkov, V. Y., Lorius, C., Pepin, L., Ritz, C., Saltzman, E., and Stievenard, M.: Climate and atmospheric history of the past 420,000 years from the Vostok ice core, Antarctica, Nature, 399, 429–436, 1999. a, b

Porter, D. F., Cassano, J. J., and Serreze, M. C.: Local and large-scale atmospheric responses to reduced Arctic sea ice and ocean warming in the WRF model, J. Geophys. Res.-Atmos., 117, D11115, https://doi.org/10.1029/2011JD016969, 2012. a

Rind, D., Peteet, D., and Kukla, G.: Can Milankovitch orbital variations initiate the growth of ice sheets in a general circulation model?, J. Geophys. Res.-Atmos., 94, 12851–12871, 1989. a, b

Roe, G. and Lindzen, R.: A one-dimensional model for the interaction between continental-scale ice sheets and atmospheric stationary waves, Clim. Dynam., 17, 479–487, 2001a. a, b

Roe, G. H. and Lindzen, R. S.: The mutual interaction between continental-scale ice sheets and atmospheric stationary waves, J. Climate, 14, 1450–1465, 2001b. a, b, c

Ruddiman, W., McIntyre, A., Niebler-Hunt, V., and Durazzi, J.: Oceanic evidence for the mechanism of rapid northern hemisphere glaciation, Quaternary Res., 13, 33–64, 1980. a, b

Screen, J. A. and Simmonds, I.: Erroneous Arctic temperature trends in the ERA-40 reanalysis: A closer look, J. Climate, 24, 2620–2627, 2011. a

Skamarock, W., Klemp, J., Dudhia, J., Gill, D., Barker, D., Duda, M., Huang, X., Wang, W., and Powers, J.: A description of the Advanced Research 488 WRF Version 3, Tech. rep., NCAR Technical Note NCAR/TN-475+STR, 489, 113 pp., https://doi.org/10.5065/D68S4MVH, 2008. a, b, c

Smith, R.: The influence of the earth's rotation on mountain wave drag, J. Atmos. Sci., 36, 177–180, 1979. a

Stokes, W. L.: Another look at the ice age, Science, 122, 815–821, 1955. a, b

Van den Berg, J., Van de Wal, R., and Oerlemans, J.: Effects of spatial discretization in ice-sheet modelling using the shallow-ice approximation, J. Glaciol., 52, 89–98, 2006. a

Vavrus, S., Philippon-Berthier, G., Kutzbach, J. E., and Ruddiman, W. F.: The role of GCM resolution in simulating glacial inception, Holocene, 21, 819–830, 2011. a

Vettoretti, G. and Peltier, W.: Post-Eemian glacial inception. Part II: Elements of a cryospheric moisture pump, J. Climate, 16, 912–927, 2003. a

Vettoretti, G. and Peltier, W.: Sensitivity of glacial inception to orbital and greenhouse gas climate forcing, Quaternary Sci. Rev., 23, 499–519, 2004. a, b

Vimeux, F., Masson, V., Jouzel, J., Stievenard, M., and Petit, J.: Glacial–interglacial changes in ocean surface conditions in the Southern Hemisphere, Nature, 398, 410–413, 1999. a

Waelbroeck, C., Labeyrie, L., Michel, E., Duplessy, J. C., McManus, J., Lambeck, K., Balbon, E., and Labracherie, M.: Sea-level and deep water temperature changes derived from benthic foraminifera isotopic records, Quaternary Sci. Rev., 21, 295–305, 2002. a

Wang, Z. and Mysak, L. A.: Simulation of the last glacial inception and rapid ice sheet growth in the McGill Paleoclimate Model, Geophys. Res. Lett., 29, https://doi.org/10.1029/2002GL015120, 2002. a, b

Wang, Z., Cochelin, A.-S. B., Mysak, L. A., and Wang, Y.: Simulation of the last glacial inception with the green McGill Paleoclimate Model, Geophys. Res. Lett., 32, https://doi.org/10.1029/2005GL023047, 2005.  a

Weertman, J.: Deformation of floating ice shelves, J. Glaciol., 3, 38–42, 1957. a

Wells, G. L.: Late-glacial circulation over central North America revealed by aeolian features, in: Variations in the global water budget, Springer, 317–330, 1983. a

Williams, L. D.: An energy balance model of potential glacierization of northern Canada, Arctic Alpin. Res., 11, 443–456, 1979. a, b

Yang, Z.-L., Dickinson, R., Henderson-Sellers, A., and Pitman, A.: Preliminary study of spin-up processes in land surface models with the first stage data of Project for Intercomparison of Land Surface Parameterization Schemes Phase 1 (a), J. Geophys. Res.-Atmos., 100, 16553–16578, 1995. a

Yin, J. H. and Battisti, D. S.: The importance of tropical sea surface temperature patterns in simulations of Last Glacial Maximum climate, J. Climate, 14, 565–581, 2001. a

Yoshimori, M., Reader, M., Weaver, A., and McFarlane, N.: On the causes of glacial inception at 116 kaBP, Clim. Dynam., 18, 383–402, 2002. a

Zdanowicz, C., Smetny-Sowa, A., Fisher, D., Schaffer, N., Copland, L., Eley, J., and Dupont, F.: Summer melt rates on Penny Ice Cap, Baffin Island: Past and recent trends and implications for regional climate, J. Geophys. Res.-Earth, 117, F02006, https://doi.org/10.1029/2011JF002248, 2012. a, b