Geologic evidence suggests that the last glacial inception (115 kya) occurred within the mountains of Baffin Island. Global climate models (GCMs) have difficulty simulating this climate transition, likely because of their coarse horizontal resolution that smooths topography and necessitates the use of cumulus parameterizations. A regional configuration of the Weather Research and Forecasting (WRF) Model is used to simulate the small-scale topographic and cloud processes neglected by GCMs, and the sensitivity of the region to Milankovitch forcing, topography, and meteorology is tested. It is found that ice growth is possible with 115-kya insolation, realistic topography, and slightly colder-than-average meteorology, represented by specific years within the past three decades. The simulation with low GCM-like topography shows a negative surface mass balance, even with the relevant orbital parameter configuration, demonstrating the criticality of realistic topography. The downslope growth of the ice sheets is studied by looking at the sensitivity of the mass balance to initial snow cover prescribed beyond that of the present day. It is found that the snow-albedo feedback, via its effects on the mass balance, allows such larger snow cover to persist. Implications for GCM studies of glacial inception are discussed.
Since the mid-Pleistocene transition ~800 thousand years ago (kya), 100 000-yr glacial cycles have dominated Earth’s climate variability (Petit et al. 1999). The transition into an ice age from an interglacial period (similar to today’s climate) is known as glacial inception. The area surrounding Hudson Bay, including Baffin Island, Labrador–Ungava, and Keewatin, has been identified as a likely inception site based on sediment cores (Clark et al. 1993), arguments pertaining to the orography (Oerlemans 2002; Ives et al. 1975), and previous studies (Williams 1979; Kleman et al. 2002; Otieno and Bromwich 2009). According to Clark et al. (1993), the growth of an ice sheet in North America from Baffin Island westward occurred over 20 000 years. Ice-volume reconstructions (Petit et al. 1999) confirm the time scale of this glacial inception process. We are interested in exploring the conditions necessary for the initial high-elevation mountain glacier expansion into the low-lying surrounding area, which occurs on a much shorter time scale. Using a high-resolution regional atmospheric model, we focus on the roles of complex topography and clouds on glacial inception.
According to Milankovitch theory, the driver or pacemaker of glacial cycles is changes in Earth’s orbital parameters (Milankovitch 1941; Hays et al. 1976), perhaps by inducing wet winters and cool summers (Bromwich et al. 2002). A more updated view is that “pacemaking” is due to nonlinear phase locking of a self-sustained glacial cycle to the insolation forcing, which therefore sets the timing of inceptions and terminations (Hyde and Peltier 1987; Gildor and Tziperman 2000; Tziperman et al. 2006; Crucifix 2013). It has been hypothesized that weaker summer insolation (in particular weaker integrated insolation; Huybers 2006) leads to decreased melting and thus glacial inception. Historically, state-of-the-art global climate models (GCMs) have difficulty simulating the transition from interglacial to glaciated due to Milankovitch forcing alone (Rind et al. 1989).
A successful glacial inception simulation in a model without ice flow is thought to include the expansion of perennial snow into initially nonglaciated regions. Early GCM studies either failed to simulate perennial snow cover or the pattern of snow expansion did not match geological records (Dong and Valdes 1995; Meissner et al. 2003), accumulating snow outside of the Hudson Bay region. The introduction of SSTs from the time of the last glacial inception and ocean feedbacks to GCMs has been shown to lead to increased poleward moisture transport and snow accumulation over the Laurentide region (Khodri et al. 2001; Yoshimori et al. 2002). Vettoretti and Peltier (2004) considered CO2 and orbital parameters, obtaining snow growth northward of inception sites inferred from geologic records. Simpler models have also been used to understand the relative importance of other climate factors that could influence glacial inception. These factors include the effects of nonlinear response to Milankovitch forcing (Le Treut and Ghil 1983), the impact of sea ice on snow accumulation (Gildor and Tziperman 2000), vegetation feedbacks (Crucifix and Loutre 2002; Wang et al. 2005), dust (Lambert et al. 2008; Calov et al. 2005)—particularly its ability to switch off the hydrological cycle (Farrell and Abbot 2012)—and the role of coupled Milankovitch–meridional ocean circulation feedbacks (Timmermann et al. 2010). Otieno and Bromwich (2009) thoroughly examined and nicely summarized the inception problem using a land surface model, noting that insolation, wet winters, and cold summers were all important but insufficient for inception, which necessitated an artificial 4°C of cooling. Jackson and Broccoli (2003) also illustrated that increased precipitation in northeast Canada at about 115 kya is a distinct possibility owing to increased storm activity. Most recently Jochum et al. (2012) was able to simulate inception over the Baffin Island region by using a higher GCM resolution that allowed for a more accurate representation of the topography. The importance of using realistic topography was also noted by Marshall and Clarke (1999), Wang and Mysak (2002), and Vettoretti and Peltier (2003). Although the work of Jochum et al. (2012) was a success in many ways, the topography over Baffin Island was still at least 600 m below the observed height of the mountainous terrain in the region, and there was snow accumulation in areas that may not have been glaciated. Furthermore, the model they used has a substantial cold bias in its present climate state, of up to 5°C over Baffin Island in the annual mean—raising questions about whether they simulated inception for the right reasons.
Given the demonstrated critical role of accurate topography over the inception sites, our objective is to study the inception problem using a model that represents this topography as accurately as possible. We, therefore, use the Weather Research and Forecasting (WRF) Model in a high-resolution (4 km) regional configuration over Baffin Island, focusing specifically on the Penny Ice Cap, which is one of the two large ice caps on the island today remaining from the Laurentide Ice Sheet. Ice cores from the Penny Ice Cap (Fisher et al. 1998) illustrate that accumulation rates on the ice cap can range from 0.18 to 0.36 m yr−1, which is something our high-resolution study will allow us to examine. The use of a high-resolution model also allows us to address the known important role of cloud radiative forcing (CRF) in the Arctic and subarctic mountain glaciers. From measurements and modeling (Weller 1972; Shupe and Intrieri 2004; Zhang et al. 1996), we know that CRF has a dominant role in the mass balance in the Arctic. The CRF can range from ~40 to −75 W m−2 (Zhao and Garrett 2015), which translates to a positive CRF (warming) in the winter months and cooling during the summer of about −50 W m−2. Low cloud feedbacks have been suggested to partially offset the insolation decreases (Jochum et al. 2012), meaning they cause less cooling in a glacial inception scenario. These issues are not well addressed by GCMs because their coarse resolution necessitates the use of inaccurate and uncertain convection and cloud parameterizations. Our model can be run in a cloud-resolving configuration without convective parameterizations, allowing us to study the important interactions between topography, orographic precipitation, and cloud formation. Using this model we explore and analyze the impact of topography, clouds, orbital parameters, and meteorology on the melting and accumulation of mountain glaciers.
We view the process of glacial inception as consisting of at least three mechanisms, two of which we explore in this paper. First, ice caps may expand from their current locations in a way that may be very sensitive to climate but is fundamentally reversible and linear for small-amplitude changes. Second, snow cover may alter regional climate and surface mass balance through the snow-albedo feedback. This second mechanism may simply act as an amplifier but can also lead to instability and irreversible transitions in idealized models (Lee and North 1995), although such instability is poorly understood in the context of glacial inception. Our main focus is the first mechanism (reversible ice cap growth), but we also address the second (regional snow-albedo feedback). A third mechanism underlying glacial inception, not dealt with here, is the development of a large-scale ice sheet whose growing elevation leads to further climate feedbacks including both cooling, which reduces melting, as well as the elevation-desert effect, which reduces accumulation at large heights (Weertman 1976; Oerlemans 1989). This third mechanism, known as the height–mass balance feedback, can also lead to instability and irreversible transitions, but we cannot directly address it because our model simulations are run for only a few years without an ice sheet model for glacial flow.
Our main results show that an accurate representation of topography beyond what can currently be done in state-of-the-art GCMs is critical for correctly simulating the mass balance of mountain glaciers. We demonstrate and analyze the importance of orographic precipitation. The results indicate that clouds, in the presence of artificially lowered topography, respond by inducing cooling and reducing melting over the low topography. We, therefore, suggest that the results of GCMs—even if they find inception under the appropriate orbital forcing—need to be critically examined in view of this compensation between clouds and smoothed topography. Overall, with the correct topography, we find that a series of cold and wet years can lead to glacier expansion on Baffin Island, roughly doubling the areal extent of the ice cap. We find little evidence for irreversible snow-albedo feedback at the spatial scales of our simulation, but we do show that the snow-albedo feedback can amplify orbitally forced changes in the surface mass balance and lead to further growth of perennial snow.
a. WRF Model
We use the WRF Model (Skamarock et al. 2008), version 3.7, a fully compressible, nonhydrostatic model, which has been shown to reproduce reasonable values in the Arctic (Cassano et al. 2011) and over complex terrain (Kilpeläinen et al. 2011). Reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF; Dee et al. 2011) are used for setting the meteorological boundary and initial conditions, including sea surface temperatures and sea ice extent. Cassano et al. (2011) also noted that the ECMWF interim reanalysis (ERA-Interim) is a suitable choice of boundary conditions for WRF, leading to minimal bias in the Arctic, specifically in temperature. We set up a nested experiment with a 12-km-resolution outer domain (Fig. 1a). Technically, there is some concern that this region may be sensitive to our domain size choice, but a test simulation with 20-km resolution produces similar results, giving us confidence in our 12-km parent domain. The inner domain has a 4-km resolution (Fig. 1b), and the actual elevation of the Penny Ice Cap, exceeding 2000 m, is well resolved at a 4-km grid spacing. We use realistic land cover from the WRF database (note the extent of Baffin Island covered by ice land-use type is outlined in white in Fig. 1). Over these land grid points classified as ice, we prescribe a snow water equivalent (SWE) of 5000 kg m−2—sufficient for maintaining snow cover for the duration of our runs. We focus our analysis on changes over this initial ice cap.
We use the WRF single-moment 6-class microphysics scheme, which allows for ice, snow, and graupel processes (Hong and Lim 2006). Noah-MP (Niu et al. 2011) 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 (Zdanowicz et al. 2012). Noah-MP caps the snow water content at 2000 kg m−2 by default, and we removed this limit. The radiation scheme is the Rapid Radiative Transfer Model for GCM applications (RRTMG; Iacono et al. 2008), which has been found to have minimal bias in WRF (Cassano et al. 2011). In this scheme, we set the CO2 to 290 ppm for our glacial inception scenarios (Petit et al. 1999; Vettoretti and Peltier 2004). Since this is a high-resolution study at 4 km, we chose to allow WRF to be cloud resolving, turning off cumulus parameterization. The boundary layer is simulated with the Mellor–Yamada–Janić scheme (MYJ; Janjić 1994). We do not nudge to observations in any way other than prescribed lateral boundary conditions, as our goal is to simulate a possible atmosphere 115 kya. With this WRF configuration we find a reasonable reproduction of the ECMWF boundary conditions when we simulate present-day climate. Furthermore, precipitation on the coast matches observations at values of about 300 mm yr−1, and the equilibrium line altitude is about 1400 m in our model, also agreeing with observations (Zdanowicz et al. 2012).
To better understand the importance of insolation J, meteorology (playing a role via both temperature and precipitation rate and referred to together as T), and topography Z on glacial inception, we conduct six experiments (Table 1). Plotting the integrated insolation (Huybers 2006) over the past 140 000 years, which includes the last ice age and previous interglacial, indicates an insolation minimum at 115 kya (Fig. 2a). Ice and sediment core records (Petit et al. 1999; Clark et al. 1993) confirm that the last glacial inception began at this time. Therefore, we choose this orbital configuration as our base state. We also run the model for the insolation of 128 kya (time of maximum integrated insolation over the past 140 ka) and present day (which coincidentally is a mean value for the past 140 ka) insolation, denoted and , respectively (Fig. 2a). Our method for setting the chosen orbital configuration in WRF’s radiation module is detailed in the appendix.
We examine the sensitivity of the Penny Ice Cap region to heat and moisture fluxes by using present-day interannual variability to supply boundary conditions for the regional model. We do not comment on how such altered boundary conditions would arise in a past climate simulation but instead suggest that the present-day variability can be used as general proxy for regional warming or cooling—which could occur naturally or as a result of orbital or CO2 forcing. Therefore, using the ERA-Interim dataset spanning the past 30 years, we determine average, cold, and warm boundary conditions. We also confirm that these anomalies are not unique to ECMWF but are also consistent with NCEP reanalysis data (Kalnay et al. 1996). For our control run, we identify 1980–82 as the three consecutive years that are the least anomalous (within one standard deviation and close to 0 anomaly) regarding summer temperature and yearly snowfall over the Baffin Island region (Figs. 2b,c). We run this simulation from October 1979 to October 1982 with 115-kya insolation and realistic topography. Model runs representing cold and warm perturbation simulations use boundary conditions based on 3-yr periods of the most anomalous summer temperatures and yearly snowfall (with at least one year outside of one standard deviation). We identify 1986–88 as anomalously wet with cold summers and 2008–10 as anomalously dry with warm summers (Figs. 2b,c). These simulations with cold and warm boundary conditions use insolation from 115 kya and realistic topography. We confirm that the atmosphere (heat and moisture fluxes) is strongly constrained by the boundary conditions, and for a regional study, these fluxes are mostly independent of the outer domain size.
To test the role of topography in inception, we run a case with modified topography (run ) using the topography from the 1°-resolution Community Earth System Model (Hurrell et al. 2013) from CMIP5. Over the Penny Ice Cap, this configuration lowers the peak elevation to ~500 m (Fig. 1c), whereas the 4-km resolution is over 1900 m, just shy of the realistic 2000-m peak elevation (Fig. 1b).
We also test the sensitivity of the region to initial snow cover extent by imposing an initial ice cap above a threshold elevation, which is varied across experiments (Table 2). We impose a glacial inception scenario by using 115-kya insolation and run this set of experiments for one year with both average (October 1979–October 1980) and cold (October 1985–October 1986) boundary conditions. We vary the initial ice cap extent by changing the land type to snow/ice and setting the SWE to be 5000 kg m−2 for all grid points above the chosen threshold elevations of 800 (approximate Penny Ice Cap extent), 400, and 0 m (all land covered in snow). Additionally, we perform one simulation with average meteorology and no initialized ice cap, meaning tundra land type is prescribed everywhere.
In section 3a, we study the first step of glacial inception—reversible mountain glacier growth over Baffin Island—by considering the sensitivity of the mass balance to insolation, meteorology, and topography. In section 3b, we examine the second step of the inception process: the potential for regional snow-albedo feedback to amplify initial changes in glacier extent.
a. Sensitivity to insolation, temperature, and topography
Milankovitch theory explains glacial inception mostly through decreased summer insolation and therefore ablation A, yet total precipitation P matters too, of course. In our simulations, there is very little liquid precipitation, and when it does occur, most freezes and contributes to the snowpack (not shown). Thus, we find that year-round all types of precipitation over Baffin Island are important to the mass balance, and we define the net snow accumulation for each year as . The difference in net accumulation between our control simulation and all other experiments is denoted , which varies based on our chosen perturbations of insolation , temperature , and topography (see Table 1).
We begin our analysis of glacial inception by looking at average snow depth on the Penny Ice Cap (approximately 800 m and above) over model runs that are three years long (Fig. 3a). There is net melting for all simulations, except for the cold simulation (Table 1), which is 0.5°C colder and has about 7% more precipitation than the control simulation with average meteorology. In our control case with 115-kya insolation, the snow depth only decreases by 100 kg m−2 and appears to be approximately in equilibrium according to the area-integrated mass balance (see Fig. 5b). Figures 3b,c confirm that the simulation with cold meteorology has the least melting during the summer and the most precipitation over three years, leading to the only experiment with a positive net accumulation over the initialized ice cap. Our simulations with 115-kya insolation have shorter melting seasons than those using insolation from present day and 128 kya , as expected. The only exception is when the experiment is forced with warm boundary conditions and 115-kya insolation, which leads to the longest and strongest melting seasons causing a loss of 500 kg m−2 yr−1 on the ice cap. Thus, we conclude that meteorology, via its effects on temperature and accumulation, has the largest impact on mass balance in these experiments. The difference in summer temperature between the colder and warmer boundary conditions is about 1.5°C and leads to a melting difference of 500 kg m−2 yr−1 (Fig. 3a). The yearly response of the melting to the imposed warm and cold boundary condition temperature perturbations is ~300 kg m−2 °C−1.
In agreement with Milankovitch theory, we find that stronger insolation ( and ) causes ablation to increase, by 100 and 200 kg m−2 yr−1, respectively, compared to our control simulation with minimum insolation from 115 kya. However, we can also see the effect of local insolation changes on precipitation by looking at Fig. 3c. Our simulations with 115-kya (purple), present-day (green), and 128-kya insolation (orange) have slightly different precipitation rates. Thus, Milankovitch forcing causes local changes in precipitation of ~20 kg m−2 yr−1 with the 115-kya simulation having the least amount of precipitation and 128-kya having the most. For our glacial inception scenario, we find that 115-kya insolation causes cooler temperatures, less moisture in the atmosphere, and decreased precipitation. In particular, we find that the rainfall rates are reduced when insolation decreases, but snowfall is fairly consistent across the runs (not shown). This result contrasts with studies with coupled atmosphere–ocean models that find increased moisture flux and precipitation in the Arctic (Khodri et al. 2001; Jackson and Broccoli 2003) during glacial inception and disagrees with Yoshimori et al. (2002), who found increased total precipitation but decreased snowfall. Thus, we are interested in examining how large-scale circulation is impacted by insolation changes, what is necessary for our ideal (cold and wet) conditions, and if the increased moisture flux found in earlier studies can be counteracted by regional temperature decreases.
The importance of realistic topography Z is demonstrated in Fig. 3a, which shows that using a low, GCM-like topography causes snow loss equivalent to when the model is forced with anomalously warm weather (warm) or increased integrated insolation (128 kya). In only three years, lowered topography alone causes a decrease in snow depth of about 1000 kg m−2, but this decrease in snow depth is not from melting alone (Figs. 3b,c). The low-topography simulation experiences ~500 kg m−2 more melting and ~500 kg m−2 less precipitation than the control experiment, indicating that precipitation and melting play an equal role in the mass balance. Since the control experiment and the low topography experiment are forced with the same atmospheric boundary conditions, their differences in accumulation indicate that orographic precipitation may be important in this region. We confirm that convective precipitation is mostly inactive in this region, and we diagnose possible orographic enhancement of precipitation from the spatial pattern of accumulation in Fig. 4a. The highest rates (1000 kg m−2 yr−1) occur on the windward (southwestern) side of the Penny ice and lower (200 kg m−2 yr−1) on the lee (northeastern) side. In contrast, when the simulation uses smoothed GCM-like topography, the snowfall is mostly uniform over the domain owing to the absence of significant topography (Fig. 4d). The low-topography simulation also has no snowfall during the summer and reduced snowfall during most of the rest of the year (not shown). Thus, realistic orography enhances the amount and alters the type of precipitation, which is not captured by GCMs, whose topography is substantially lower owing to their coarse resolution.
The effect of topography on the mass balance is further explored using spatial averages (Fig. 4). Precipitation is enhanced by orography, but the topography also greatly affects ablation. In the lower-topography case (Fig. 4e), ablation is fairly uniform, while in the realistic topography case (Fig. 4b) the highest amount of melting is localized on the edges of the initialized snow cover at about 800 m, and there is no melting at the top of the Penny Ice Cap. The pattern of mass balance for the realistic case (Fig. 4c) suggests that the upper ice cap would grow each year with realistic topography, while mass balance for the low-topography simulation indicates that the ice cap would be completely eliminated (Fig. 4f).
We investigate ice cap growth further by examining elevation-sorted mass balance in Fig. 5. We rank each land grid cell by the surface elevation and aggregate the points into bins representing 5% of the area each (Fig. 5a). Note that the ice cap lies above 800 m, which is in the 65th height percentile. In Fig. 5b, plotting the annual change in snow depth identifies the accumulation zone (values above zero), the ablation zone (values below zero), and the equilibrium line (the transition from net accumulation to net melting). In Fig. 5b, only the low-topography simulation (black dotted line) shows net melting over the entire domain. Increasing the insolation causes the equilibrium line to shift to higher altitudes, as seen by comparing the control 115-kya case (purple) with the increased-insolation 128-kya case (orange). Similarly, the combination of cold meteorology and 115-kya insolation results in the lowest equilibrium line.
In an actual ice cap, accumulation at the higher elevations leads to ice flow to lower elevations where the ice is ablated, such that over the extent of the ice cap, total accumulation equals total ablation. Our model does not include ice flow explicitly, and we, therefore, represent its effects indirectly. For that purpose, we integrate the yearly mass balance on the ice cap (65th percentile of height and above) beginning at the upper part of the domain where and integrating downslope until the integrated mass balance equals zero, indicated by large colored solid dots in Fig. 5b. This location, where the integrated mass balance vanishes, indicates the anticipated location of the equilibrium ice cap edge. We perform this and future integrations only over ice-covered area, where the mass balance can actually be evaluated. This, of course, should be viewed as a crude representation of ice flow and the feedbacks that it would induce. The present-day experiment shows that the ice cap extent is at the 70th percentile, corresponding to a height of ~900 m. This suggests that if run with an ice flow model, the ice cap in this run would equilibrate just above the current height of the Penny Ice Cap at 800 m (65th percentile). The 128-kya and warm meteorology experiments would retreat further upslope to heights of 1100 (85th percentile) and 1300 m (95th percentile), respectively. The low-topography experiment has no accumulation zone; the entire ice cap is therefore expected to disappear. The integrated mass balance of the control (solid purple line in Fig. 5b) goes to zero at the edge of the ice cap, signifying that the ice cap edge would be maintained at 800 m with 115-kya insolation, realistic topography, and average meteorology. Finally, for the cold simulation (dashed blue line in Fig. 5b), the ice cap maintains a positive integrated mass balance as indicated by the blue arrow and is expected to grow—and perhaps develop into an ice sheet. The average grid cell of the initial ice cap in this simulation accumulates 0.2 m yr−1; over time, ice flow would allow the ice cap to expand downslope past its initial bounds. In the next subsection, we analyze the ice cap growth by considering simulations with different initial ice cap extents.
Our cloud-resolving and high-resolution model allows us to examine the role of clouds in glacial inception and provide an independent test of GCM studies based on parameterized clouds and convection. Jochum et al. (2012) found low clouds to be a negative feedback when it comes to Milankovitch forcing: weaker insolation leads to less low clouds and therefore to a weaker shortwave CRF. Clouds also lead to warming that may counteract the Milankovitch-induced cooling. Figure 6 shows the shortwave (SW), longwave (LW), and net (SW + LW) CRF sorted by elevation for the summer (June–July–August) during the first year of the model run. For present-day insolation (green line in Fig. 6c), the model CRF is −40 W m−2 at low elevations, comparing well with other studies on CRF in the Arctic (Zhao and Garrett 2015). The figure shows only four of the runs because interannual variation in CRF for the warm and cold simulations is too large to allow us to deduce any clear temperature dependence. Shortwave cloud forcing becomes less negative as insolation decreases (cf. orange and purple curves in Fig. 6a), and this effect is consistent over the three years we ran the model, which is why we only look at the first year. Next, we consider the interplay of clouds and topography. The low-topography simulation forms more low clouds over the entire domain, and these clouds have a higher liquid water path and more negative CRF. The net CRF for low topography only becomes similar to that of realistic topography in the lowest part of the domain, where the two simulations have similar cloud properties. We conclude that clouds act as a negative feedback not only in response to Milankovitch changes but also in response to topography variations. This suggests that GCM studies of the inception problem that do not represent the topography accurately may suffer CRF biases, which may allow the development of ice caps due to cloud response to the modified topography.
b. Sensitivity of mass balance to initial snow extent
The previous set of experiments involves an essentially linear and reversible sensitivity to various inception forcing factors, and from our above analysis, we infer the location of the equilibrium ice cap edge by integrating the mass balance to the point where its net vanishes. As discussed in the introduction, the second stage of the inception process involves the (potentially nonlinear and irreversible) snow-albedo feedback and the effect of larger areas of snow cover on the local climate, particularly the survival of perennial snow. This feedback is our focus here.
To examine this second inception mechanism, our next set of experiments explores the sensitivity of the mass balance to a specified snow cover that is beyond the present-day ice cap extent. We perform these runs for both average (1980) and cold (1986) meteorology boundary conditions—please note there are differences in mass balance between Fig. 5, which is the average of three years, and Fig. 7, which is only a yearlong simulation. As explained in section 3b, we run such experiments for initial snow cover that is specified above threshold elevations of 800, 400, and 0 m.
For each of the two boundary conditions, we find that the equilibrium line (transition from accumulation to ablation) occurs at nearly the same elevation, although increased snow cover does cause a slight downslope shift (cf. same-color curves in Fig. 7a). Note that the mass balance may be estimated only over height ranges that are initially snow covered, and, therefore, we plot it only over these elevations. The average boundary conditions have an equilibrium line at about 1000 m (80th-percentile elevation bin), and the cold boundary conditions have an equilibrium line at about 800 m (70th percentile). We again calculate the integrated mass balance over the ice cap in each experiment. We denote the location of the zero-integrated mass balance with a solid circle, and if the integrated mass balance remains positive down to the edge of initial ice extent, indicating growth, we use an arrow. For the 800-m-threshold experiment denoted by the dotted blue curve in Fig. 7a, the integrated mass balance does not vanish over the specified initial snow extent, indicating that the 800-m ice cap would grow. This is different from the 3-yr simulations in Fig. 5b, because we ran the model for only 1 yr, which has the least amount of melting.
With average meteorology, the ice cap is expected to expand to 700 m (60th percentile; purple circle in Fig. 7a) as identified by the 400-m and all-snow (0 m) simulations. This conclusion can only be reached by specifying larger snow cover than the present-day ice cap. Furthermore, the cold boundary conditions result in a doubling of ice cap areal extent (relative to average meteorology), with the ice cap edge at about 300 m (15th percentile; blue circle in Fig. 7a). Such meteorology (cold summers, larger precipitation) persisting over several years can, therefore, lead to ice growth over much of the area surrounding the Penny Ice Cap.
We also examine ablation A and precipitation P separately in Figs. 7b,c. As expected, melting decreases at higher elevations where temperature is colder, but we see a slight fanning out of the yearly melt rates downslope for different ice cap extents. For instance, when we lower the threshold elevation from 400 to 0 m (all snow), less melting occurs (cf. solid and dashed curves in Fig. 7b), and we will show next that this is because the temperature is lower in the all-snow experiment because of the snow-albedo feedback. The yearly precipitation does not differ greatly between the various ice cap extents, except above 1200 m (85th percentile), where the all-snow simulations have 50 kg m−2 yr−1 less snowfall, colder temperatures, and less moisture in the atmosphere. Figure 7c indicates that over the ice cap the cold boundary conditions cause about 300 kg m−2 yr−1 of precipitation more than the average boundary conditions. Note that this precipitation rate is different from the simulations run for three years (1986–88) because 1986 is the coldest and wettest of these three years. Thus, we find that the downward expansion of the ice cap edge from 800 to 300 m is largely dependent on precipitation resulting from the large-scale circulation prescribed at the boundaries.
To confirm that the differences in melting between the different initial snow-cover experiments are due to the snow-albedo feedback, we examine the summer temperature over the domain for average boundary conditions (Fig. 8). Temperatures are below freezing where ice is present above 800 m in the control case (Fig. 8a). When we change the initial snow cover threshold elevation to 400 m, the 2-m temperature below 800 m decreases since snow is now present, and below 400 m the temperature is only ~1°C colder. When snow covers the entire domain, the snow-albedo feedback causes colder temperatures than the 400-m simulation over the same elevations (note the darker blue pattern). The temperature above 800 m does not change as we lower the ice cap threshold elevation from 800 to 0 m, but when the simulation is initiated without snow cover, temperatures are up to 3°C warmer during the summer even though an ice cap is developing.
Thus, we identify the snow-albedo feedback to be significant in this region from decreased summer temperature and lessened melting with increased snow cover. For another perspective on the effect of the snow-albedo feedback, we consider a scatterplot of positive degree-days (PDD) against melting for the simulations with average meteorology (Fig. 9). We use these plots to determine how changes in net accumulation relate to changes in temperature and how topography impacts melt rates. The data points are color coded by gridcell elevation, where blue is the highest (and coldest) part of the domain and red is the lowest. Melting increases where there are more positive degree-days, and there are generally more positive degree-days at lower elevations. For an initial snow cover threshold of 800 m, we observe a rather linear relationship between melting and PDD (Fig. 9a). As the threshold elevation is decreased, the snow-albedo feedback decreases melting rates at a given elevation. For instance, the edge of the 800-m ice cap experiences ~1000 kg m−2 yr−1 of melt, but when the entire domain is covered with snow, only elevations below 300 m melt at that rate. We can observe this melt-rate decrease in the downward shift of blue and purple points as we move from Fig. 9a to Fig. 9c (800-m elevation threshold to 0 m).
These same scatterplots show two additional interesting points. First, in all simulations, the highest elevations show positive degree-days but no melting. This occurs because refreezing is allowed within the land model: snow melts during the day when the temperature is above zero but does not always have time to run off before nightfall when it refreezes. Second, note the spread of PDD values for a given melt rate (e.g., for 1000 kg m−2 yr−1 melt rate, there is a PDD range of 100–300; Fig. 9c). A value of 100 PDD corresponds to lower elevations (more red; 300 m) at the ice cap edge, and higher values of 300 PDD have points corresponding to higher elevations (more purple; ~500 m) closer to the center. Thus, use of PDD alone to represent melting does not capture terms in the surface energy budget that depend on surface elevation and tend to lead to less melting higher on the ice cap.
Finally, the cloud radiative forcing is shown in Fig. 10c and reveals that the initial snow cover extent substantially impacts the cloud radiative forcing, mostly through temperature and surface albedo changes. Changing the meteorology causes the magnitude of the net cloud forcing to change while leaving the functional dependence on elevation similar (cf. blue and purple solid lines, e.g., in Fig. 10c). Where initial snow cover is prescribed, the SW CRF is less negative than over the tundra regions without an initial snow cover (Fig. 10a), which is consistent with observations showing the CRF becoming more negative when winter snow melts away (Shupe and Intrieri 2004). The cloud properties are very similar for all simulations, and the SW CRF changes are as expected: increasing surface albedo reduces the reflection of shortwave radiation by clouds.
We have studied the glacial inception problem using a high-resolution, cloud-resolving atmospheric model, focusing on the importance of topography and clouds. Given present-day insolation and average meteorology, our model predicts the Penny Ice Cap is close to equilibrium (would retreat from 800- to 900-m elevation with an equilibrium line altitude of ~1400 m), indicating that our model is appropriate for the study of the possible ice cap expansion in different climates. We find that a combination of Milankovitch forcing, realistic topography, and moderately cold and wet meteorology allows mountain glacier growth. This is in contrast to some previous studies (e.g., Otieno and Bromwich 2009), suggesting that additional large cooling of 4°C or some additional feedbacks—for example, vegetation—are needed. However, it is possible that such feedbacks would be needed in later stages of inception when snow cover spreads outside of Baffin Island. We find that our choice of boundary conditions heavily constrains the heat and moisture fluxes to the region, and taking advantage of interannual variability allows us to use modern climate as a proxy for the slightly cooler climate experienced during glacial inception.
Because our model does not include active ice flow, we attempt to represent that effect by calculating the glacier mass balance integrated from the top downward. While the model does not show a buildup of perennial snow below the current level of the Penny Ice Cap, the accumulation over the ice cap increases under the appropriate inception conditions mentioned above, implying that simulations with 115-kya insolation and “cold” boundary conditions would allow the ice cap to expand downslope. Our cold boundary conditions are represented by the meteorology of years 1986–88 and have a cooling of only 0.5°C and precipitation increase of ~7%, leading to the only simulation with a positive mass balance over the ice cap of 0.2 m yr−1. We find this sufficient for glacier area growth of ~50%–100%—even without any height–mass balance feedback. This is a greater sensitivity of mass balance to climate than found in previous studies that required a much stronger cooling to get substantial glacier area expansion.
The topography in state-of-the-art GCMs is characterized by ~500-m maximal elevations in Baffin Island as opposed to the actual 2000-m topography, which is well represented in our study. We find that GCM-like representation of topography is insufficient and leads to biases, despite allowing for some progress on the glacial inception problem (Jochum et al. 2012). In particular, the lower topography leads to increased ablation and decreased precipitation over the area crucial to glacial inception. When we run our model with the GCM-like topography, it does not produce any net accumulation zone for average present-day meteorology. The lowered GCM-like topography leads not only to warmer conditions and thus more summer melt but also to reduced accumulation, as orographic precipitation—a critical process in this region—is eliminated. Smoothed topography also causes rainfall rather than snow during the summer, further contributing to theinability of GCMs to produce inception.
Moving beyond the initial glacier expansion that marks the first step of inception, we expect increases in snow and ice cover to lead to a regional snow-albedo feedback that amplifies changes in mass balance. We explored this feedback by testing the sensitivity of the mass balance to changes in initial snow cover, setting the land cover to snow/ice everywhere above an elevation of 800, 400, or 0 m and also considering both average and cold boundary conditions. We find that the initial snow cover impacts temperature and melting but changes precipitation very little. The mass balance deduced from these runs using cold meteorology indicates that once covered by snow, the area would support an equilibrium glacier extent above 300-m elevation, although an actual ice flow model would be needed to verify this conclusion and calculate the relevant response time scales of the ice caps.
Previous GCM work by Jochum et al. (2012) found that shortwave cloud radiative forcing acts as a negative feedback on Milankovitch forcing: when the insolation is reduced, low clouds let more SW radiation penetrate and weaken the tendency toward inception. We confirmed this conclusion using a model that does not rely on parameterized convection and clouds. We also find that clouds act as a negative feedback to topography changes: when the topography is lowered from its realistic values to a GCM-like smooth topography, clouds respond by letting less SW penetrate to the surface. Thus, although lowering the topography leads to surface warming that tends to prevent inception from happening, the cloud forcing makes the state more favorable to inception. This cloud response may allow inception to occur without realistic topography, representing a model bias and suggesting that inception in GCMs with smooth topography should be examined carefully.
We conclude that with realistic topography and 115-kya insolation values, glacier growth (the first step of inception) may be possible. Furthermore, the snow-albedo feedback can amplify the effects of a mere half-degree of cooling (relative to present-day average) and cause substantial snow growth downslope to ~300 m (the second phase of glacial inception). Thus, inception may occur during an unusually long sequence of moderately cold and wet years, and inception commences with downslope ice flow rather than snow accumulating spontaneously in unglaciated regions. Given that this required cooling is not extreme and the precipitation anomaly is only one standard deviation above average, it is plausible that such a sequence could happen because of normal decadal and millennial climate variability in the interglacial period preceding the inception. A succession of many cooler years could also be a result of forcing by other mechanisms such as changes in CO2, ocean circulation, or Milankovitch forcing. It is also possible for precipitation to increase during glacial inception because of increased storm activity (Jackson and Broccoli 2003), and we are interested in exploring the impact of insolation and ice sheets on large-scale circulation. Future work would need to employ an ice flow model in addition to the realistic topography used here to test these conclusions. The inclusion of an ice sheet model would also allow for the exploration of the height–mass balance feedback, the third step of glacial inception.
This work was funded by the Harvard Climate Change solutions fund, by the NSF climate dynamics program (Grant AGS-1303604), and by NSF P2C2 program (Grant OCE-1602864). Timothy Cronin was supported by a NOAA Climate and Global Change Postdoctoral Fellowship and by the Harvard University Center for the Environment. 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 acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the Community Earth System Model contributors for producing and making available their model output. For CMIP the U.S. Department of Energy’s Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. We would also like to thank Charles Jackson, Anthony Broccoli, and two anonymous reviewers for their insightful comments and guidance.
Insolation and Orbital Parameters in WRF
The incoming insolation varies over thousands of years according to changes in Earth’s orbital parameters (Berger and Loutre 1991): obliquity β, eccentricity ε, and precession ω. Within the WRF radiation module, we are able to change the obliquity explicitly, but to account for eccentricity and precession, WRF uses an eccentricity factor (Paltridge and Platt 1976), which is multiplied by the solar constant. This factor is a Fourier expansion of the square of the ratio of present-day mean distance from the sun to the actual sun–Earth distance :
where approximates the position in the orbit and t is the Julian date .
To modify Eq. (A1), we use the ratio of the semimajor axis (same as above) to the distance from the sun :
where again ε is the eccentricity and θ is the position in the elliptic orbit. With θ and ε, we fit a Fourier expansion to Eq. (A2) like Paltridge and Platt (1976) did. To apply precession, we simply shift the curve according to the difference in present-day precession angle and past precession angle. We validate this method with present-day orbital parameters from Table A1 and then apply this method to yield 115- and 128-kya insolation.