Glacial Inception on Baffin Island: The Interaction of Ice Flow andMeteorology

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 initial inception problem on Baffin Island by asynchronously coupling the Weather Research and Forecast model (WRF), configured as a high resolution inner domain over Baffin and an outer 5 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 10 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 larger-scale 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.

, which is not consistent with the geological record showing inception occurred first in eastern Canada (Clark et al., 1993). Ice sheet models and GCM studies (Dong and Valdes, 1995;Meissner et al., 2003) have a similar problem of not agreeing with the geological record. Combining GCMs and ice sheet modeling may allow for 5 progress in understanding ice-atmosphere interactions (Herrington and Poulsen, 2011;Gregory et al., 2012;Löfverström et al., 2014), 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 most GCMs.
The interaction of the ice and the atmosphere is critical to inception. Birch et al. (2017) used a high-resolution cloud-10 permitting regional climate model and found that ice would accumulate on the mountain glaciers of Baffin Island. The 4 km resolution of that study allowed us to capture both the ablation and accumulation zones of mountain glaciers. By integrating the mass balance from the top of the ice cap down slope, Birch et al. (2017) 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 15 purpose of the current paper is to explicitly include ice flow from the growing mountain glaciers, allowing us to investigate the feedbacks between a growing ice sheet and the atmosphere at higher spatial resolution than can be simulated with GCMs. For this purpose, we use a simple ice model based on the shallow ice approximation, asynchronously coupled with the Weather Research and Forecasting model (WRF, Skamarock et al., 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 Baffin Island and Our goal in this study is to better understand the interaction between changing ice cover and the regional climate of Baffin Island and the large-scale circulation over North America. We therefore use a nested configuration of the Weather Research and Forecasting Model (WRF Skamarock et al., 2008) 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 20km 5 resolution domain, but WRF does not contain an ice flow model, which means the ice expansion must be calculated offline.
Based off of Oerlemans (1981), we developed a simple ice flow model that reduces the 3D Stokes flow to two dimensions using the shallow ice approximation. The basal sliding velocity is assumed to be 0. This approximation results in a nonlinear 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: In Equation 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 nonlinear shear-thinning of ice flow is represented by a diffusivity, D, which depends on both the ice thickness and surface elevation gradients: 15 D = AH m+2 ∂H * ∂x 2 + ∂H * ∂y 2 (m−1)/2 bias in WRF (Cassano et al., 2011). CO 2 is set to 290 ppm for a glacial inception scenario (Petit et al., 1999;Vettoretti and Peltier, 2004). The boundary layer is simulated with the Mellor-Yamada-Janic scheme (MYJ, Janjic, 1994). We do not nudge to observations in any way other than prescribed lateral boundary conditions at the edges of our outer domain. Note that these physics are the same as those used in our previous paper (Birch et al., 2017), except now we use the Grell convective scheme (Grell and Dévényi, 2002), which has been shown to give good results in present day (Hines et al., 2011) and LGM (Bromwich 15 et al., 2005) Arctic climates.
With this WRF configuration we find a reasonable reproduction of the ECMWF atmospheric state when we simulate presentday climate. Furthermore, precipitation on the coast of Baffin Island matches observations at that point of about 300 mm per year (Zdanowicz et al., 2012). 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 (Kubatzki et al., 2006). 20 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 Appendix B. We use realistic present day land type from the WRF database as initial conditions for the ice extent. 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.

Asynchronous Coupling Procedure
We begin our WRF-ice flow coupled model experiment with the WRF control simulation, which uses meteorology from October 1985 to October 1986. The circulation during this year brought the coldest and wettest extremes to Baffin Island (Bromwich et al., 2002) and causes the present day ice cover to expand when combined with 115 kya insolation (Birch et al., 2017). WRF is run for 5 years, but we only analyze the last 4 years for an average climatology. Accumulation occurs over 30 areas initialized with ice cover at high elevations, and net mass balance is positive, allowing ice growth (Fig. 1a). The control simulation agrees with our previous 4 km simulation of the Penny Ice Cap (Birch et al., 2017): cold meteorology causes the positive mass balance needed for the ice cap to expand. This mass balance is then adapted for use in the ice model. 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 20km resolution control 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 control 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. We characterize 5 the mass balance as a function of elevation (G(H * )), which allows the mass balance at each grid point to be updated as the surface elevation and ice thickness change every time step. This characterization involves sorting the domain by surface elevation and forming bins every 100 m. 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, which was judged appropriate by Herrington 10 and Poulsen (2011). After 500 years, the ice model has calculated a new ice extent and surface elevation ( Fig. 2a)  In the context of glacial inception, expanding ice cover interacts with the atmosphere, causing changes to both the regional climate of Baffin Island and the large-scale circulation. 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 Section 3.1, the changes in mass balance are explained by exploring changes in temperature, ablation, precipitation, and cloud cover. In section 3.2, we take a closer look at the large-scale circulation. Additional sensitivity simulations with the ice sheet model further clarify climate 5 factors that could lead to Baffin Island glaciation in Section 3.5.
The progression of WRF-ice model iterations (5 years of WRF and 500 years of the ice model) is laid out in Figure 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 10 the ice model for 500 years using the control WRF mass balance. This new ice thickness, ice extent, and surface elevation computed by the ice model are used as inputs to the next WRF simulations. Then WRF is run for 5 years with the mass balance presented in the right column. After the second iteration of WRF, we bin the mass balance by elevation, and force the ice model with the new updated mass balance. A new ice cover is computed (left row, second column) and used as an input in WRF. Subsequent iterations pictured follow the same procedure.

15
The second iteration has the greatest ice extent, which is due to the large positive control 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 coalesceses 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. 20 The binned mass balance used in the ice model is plotted in Figure 3 with annual mass balance on the y-axis and surface elevation on the x-axis. The control WRF run is the "+" blue line with the highest 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 25 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 ten 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 illustrate desertification, or the mass balance becoming less positive over time as the ice cap increases in height and its top surface flattens.

30
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 equations 1 and 2 suggest that a 0.5 m/yr 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 Figure 2.

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 control run and subsequent iterations are pictured 5 in Figure 4. In the control 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 10 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 Figure 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.

15
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 variables 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 control simulation. In general, the high elevation parts of the region experience no melting during the summer, due to their sub-freezing     temperatures. As the ice cover expands, the highest values of melting are located near the edges of the ice cap. As expected, 20 melting decreases as the surface elevation increases.    Finally, we are interested in the role that clouds may play in glacial inception. Jochum et al. (2012) and Birch et al. (2017) 10 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, Jochum et al. (2012) noted that less low clouds form during glacial inception, while Birch et al. (2017) found that cloud properties change such that relative to present day, less shortwave is blocked from the surface. Our objective here is to examine this feedback during the integration of our asynchronously coupled model. In Figure 7, the net CRF in the control simulation is negative during the summer because the clouds block shortwave radiation 15 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 control run: the control 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.

20
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 larger scale geopotential height field, temperature, and winds in the outer domain. First, the geopotential height and winds at 500 mbar for the control simulation during the summer are depicted in Fig. 8. The wind field over Baffin Island generally comes 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 Figs. 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 5 Baffin Island alter the circulation, and we will examine this hypothesis below. The 2m temperature response in the outer domain ( Figure 9) is consistent with the changes in geopotential height. The control 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 5 warm southerly advection by the anomalous anticyclone.
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 topography in the coarser outer domain being generally lower than the 20km resolution 10 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.

The of Impact Baffin Island Ice Elevation on Large-Scale 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 15 elevation the same as the control 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.
In Figure 10, we examine the temperature of the outer domain for 4 experiments. Again, the control simulation summer 20 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", in c 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. 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 25 two simulations ("No H" and "H Only") is that if we sum the difference between their response and the control simulation, the result is very close to the difference between iteration two and the control run, as seen in Fig. 11. This indicates that the responses to the ice topography and ice land type (albedo and other surface properties) are linear.

Glacial Inception with Changing Orbital Parameters
The WRF-ice model iterations performed previously used the minimum insolation from 115kya, which is useful for diagnosing 30 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 Birch et al. (2017) 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 12 presents the mass balances calculated from WRF and used in each iteration of the ice model. Note the control simulation and 2nd iteration plotted here are the same as presented in Figure 3 because the change in insolation occurs at 114 kya, which is the third iteration in Fig. 12. The elevation desert effect is again observed and the mass balances become more 5 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. 12).

Mass Balance 100m Bins
During the summer on Baffin Island, there is again cooling from the presence of ice, but the unglaciated parts of Baffin Island 10 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 control simulation.
The increases in insolation causes continental-scale warming, and when combined with topography changes, the warming on 15 Baffin Island can be over 5 K (Figure 13). Thus, the realistic insolation changes cause further warming during inception, which makes the occurrence of inception even more difficult to explain and simulate.  Figure 13. The surface temperature in the control and differences for successive iterations.  First, we identify if we have simulated a mass balance sufficient for a successful inception. From Figure 3, we note that the mass balance in the control simulation is the most positive at high elevations and least negative at low elevations. All   As noted above, we choose the creep parameter (A(T ) = A(−10 • C)) based on Cuffey and Paterson (2010). 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(-15 • C) respectively.
For these simulations, the mass balance used is from the control 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( -15 • ) 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 5 that a warmer temperature and thus larger A would allow the ice to spread more easily, leading to larger ice cover. Figure 14 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(−10 • C) causes the elevation to be 200 m higher than A(0 • C).
More calving occurs with A(0 • C) ice reaches the coast sooner. We also examine the impact of temperatures at about −15 • C ( Fig. 14c), which  determined to be the basal temperature of the initial ice sheet growth. As expected, 10 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 control mass balance to a larger area, outside of Baffin Island and 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 10m (Creveling et al., 2017). 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 115kya with a regional high resolution configuration of WRF, coupled asynchronously with a simple ice flow model. A net positive mass balance occurs in the 1st WRF iteration simulation driven by 115 kya insolation, cold present day (1986) boundary conditions, and present day ice cover, which could cover all of Baffin Island in an ice sheet after 10 kyr if that mass balance was maintained. However, 10 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 large-scale 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.
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 1st iterations of the coupled 15 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 (Chen and Lin, 2005;Cook and Held, 1988). Thus, a different meteorology choice might lead to sustained accumulation. Bromwich et al. (2002) identified exceptionally cold conditions on Baffin Island and their varying large scale circulation patterns. Generally, 20 the cold summer circulation flows from the northwest or north-north-west. 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 a elevations 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 25 meteorology at a 4 km resolution was in approximate equilibrium, leaning towards slight expansion (Birch et al., 2017). Net ablation occurs in the 20 km resolution because of the lower topography. It is possible that during inception that 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 we have a large sample of "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 30 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 (Roe and Lindzen, 2001b;Liakka et al., 2012). In particular, small ice sheets cause a wave response that enhances accumulation, while stationary waves induce less cooling as ice sheets increase in size (Liakka and Nilsson, 2010). 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, Herrington and Poulsen (2011) found the Hudson Bay warmed in their coupled atmosphere-ice simulations due to anticyclonic flow. Gregory et al. (2012) also found warming on the margins from anticylones 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 5 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, 10 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 would likely decrease and sea ice extent grow larger (Meissner et al., 2003). Allowing the SST to cool has been shown to induce cooling over land, even if the atmosphere is warm (Yoshimori et al., 2002), and Kageyama and Valdes (2000) noted that storm tracks change with colder SSTs. The amount of sea ice can impact accumulation over the ice sheet (Yin and Battisti, 2001;Kubatzki et al., 2006) due 15 to decreasing atmospheric moisture and evaporation from ocean, according to the temperature precipitation feedback (Källén et al., 1979;Ghil, 1994;Tziperman, 2000, 2001;Li et al., 2005). Some studies argued that the cooling during inception did not cause sea ice to extend outside of Baffin Bay until 80kya (Löfverström et al., 2014), meaning the Atlantic remained ice free, which was conducive to high precipitation. Goñi et al. (2005) also argued that SSTs remained warm until 105 kya, when substantial ice sheets had formed. Both Khodri et al. (2001) and Wang and Mysak (2002) found that including 20 ocean feedbacks caused increased moisture fluxes conducive to inception, while Jackson and Broccoli (2003) and Kaspar et al. (2007) 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. In either case, using present day ocean SSTs is likely inaccurate, but more work is needed to understand the impact of the ocean.
We also fix the CO 2 concentration to 290 ppm, although CO 2 decreased during inception. This decrease has been shown 25 to influence the accumulation of northern hemisphere ice sheets (Ganopolski and Calov, 2011;Vettoretti and Peltier, 2004), accounting for up to 50% more ice (Bonelli et al., 2009). We also do not allow for dynamic vegetation, but on Baffin Island this should not matter. Most of these studies dealing with dynamic vegetation refer to the treeline moving southward during glacial inception (Harvey, 1989;Goñi et al., 2005), enhancing ice growth, but on Baffin Island the vegetation is mostly tundra and shrubs, which are not an impediment to snow accumulation as trees would be (Otterman et al., 1984). 30 We choose to drive our model with 1986 boundary conditions based on the anomalous increase in precipitation and cool temperatures over Baffin Island. The large-scale flow induced by these boundary conditions in conjunction with ice growth causes the negative feedback discussed above, but it is possible that a complete reorganization of the pattern of Northern Hemisphere stationary waves could occur as part of the inception process, which we cannot capture with our regional model begun growing at the same time though the Laurentide increased in volume more rapidly (Kleman et al., 2013;Ganopolski et al., 2010). The European ice sheet may have influenced the formation of the Laurentide (Beghin et al., 2014), but such a feedback is missing in our regional model. Roe and Lindzen (2001b) noted that changes in the stationary wave pattern would affect the downstream climate, "the possiblity of a complicated set of interactions in which each of the two major Northern Hemisphere ice sheet influences evolution of the other". 5 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 (Andrews and Mahaffy, 1976;Marshall and Clarke, 1999). While many of these ice modes do not match geological records because they 10 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 control 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 is ∼ 400 m below the peaks found in realistic Baffin Island 15 topography. Our simulation with average meteorology demonstrates that resolution changes can 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 (Birch et al., 2017).
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 20 a year is recommended (Yang et al., 1995;Cosgrove et al., 2003). We tested running the model for longer than 10 years, but the simulations did not deviate much from years 2-5 of our simulations. Deser et al. (2014) showed that multiple ensemble members of an atmospheric model simulation can give different results on a time scale up to 50 years, but this may not be an issue here because we are driving our regional model with boundary conditions that constrain the synoptic and planetary-scale flows (Denis et al., 2002).

25
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 (Van den Berg et al., 2006). As noted above, lateral expansion only occurs when substantial ice thickness exists in an adjacent grid cell. However, we find that simulations with a 4km resolution yield similar ice flow patterns to the 20 km results we present above. We are 30 neglecting basal sliding that may develop as the ice base melts due to geothermal heat flux (Marshall and Clark, 2002). The initiating ice sheet was likely frozen at the bed (Marshall et al., 2000), which makes our choice appropriate. Finally, we are not accounting for glacial isostatic adjustment. Herrington and Poulsen (2011) found that the isostatic adjustment was not critical for glacial inception, though including the adjustment revealed a small positive feedback. Therefore, contrary to conventional 20 Clim. Past Discuss., https://doi.org/10.5194/cp-2018-5 Manuscript under review for journal Clim. Past Discussion started: 7 February 2018 c Author(s) 2018. CC BY 4.0 License. wisdom, these studies suggest that isostatic adjustment might actually favor glacial inception over Baffin Island, and may even 35 be critical to the process by reducing the strength of negative elevation-circulation feedbacks on ice growth (as in the "No H" sensitivity test).
Though we find negative feedbacks on the growth of ice sheets due to clouds and large-scale circulation preventing inception in our model, high resolution regional modeling still has potential to make progress on the inception problem. We have captured the mass balance on critical mountainous regions in our control run, which GCMs cannot simulate accurately. When combined 5 with ice flow, the control 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. We do find this equivalent sea level drop is closing in on more recent sea level estimates (Creveling et al., 2017). The temperature and precipitation in this scenario is not very different from present day, meaning maintaining this meteorology could be all that is necessary for inception. The warming feedback of the large-scale circulation could be 10 mitigated by changes in the pattern of atmospheric stationary waves. Using a circumpolar regional model is the logical next step because it can capture the mass balance of the nucleation points of both the Laurentide and Fennoscandian ice sheets, as well as exploring the possible impact of a smaller Greenland ice sheet upon exiting the last interglacial.
The surface of the ice is colder than the base, meaning that throughout a real ice sheet temperature, and thus A(T ), vary. 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 ) (Weertman, 1957;Goodman and Pierrehumbert, 2003) in 10 Equation A2: Since A(T (z)) varies by height, we need a vertical temperature profile of our ice sheet, which WRF does not provide. We can calculate one using the average temperature of the ice surface from WRF and the geothermal heat flux of 55 W m −1 over 15 Baffin Island, assuming a linear vertical temperature profile. Then we evaluate A(T ) for this temperature profile as found in Equation A and take the vertical average. Typically, this vertical average captures A(T ) at warmer temperatures at the base of ice sheet, since A(T ) does not vary linearly but exponentially. Thus, it can be an order of magnitude larger at the base of the temperature compared to the surface of ice sheet. For our simulations, we calculate A(T ) = 3.5 × 10 −25 Pa −3 s −1 , which is the same as A(−10 • C) from Cuffey and Paterson (2010). Incidentally, −10 • C is basal temperature of ice sheets at the beginning of 20 inception , meaning our procedure is capturing the approximate temperature at the base of the ice sheet, which is where flow is more likely with warmer temperatures. Thus, we feel confident in our choice of A for this simulation.

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 Figure A3, where n is the time step, i is the node in the x direction, and j is the node in the y direction. 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 Oerlemans (1981) 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