Modeling the evolution of the Juneau Icefield between 1971 and 2100 using the Parallel Ice Sheet Model (PISM)

Abstract

We study the evolution of the Juneau Icefield, one of the largest icefields in North America (>3700 km²), using the Parallel Ice Sheet Model (PISM). We test two climate datasets: 20 km Weather Research and Forecasting Model (WRF) output, and data from the Scenarios Network for Alaska Planning (SNAP), derived from spatial interpolation of observations. Good agreement between simulated and observed surface mass balance was achieved only after substantially adjusting WRF precipitation to account for unresolved orographic effects, while SNAP's climate pattern is incompatible with observations of surface mass balance. Using the WRF data forced with the RCP6.0 emission scenario, the model projects a decrease in ice volume by 58–68% and a 57–63% area loss by 2099 compared with 2010. If the modeled 2070–99 climate is held constant beyond 2099, the icefield is eliminated by 2200. With constant 1971–2010 climate, the icefield stabilizes at 86% of its present-day volume. Experiments started from an ice-free state indicate that steady-state volumes are largely independent of the initial ice volume when forced by identical scenarios of climate stabilization. Despite large projected volume losses, the complex high-mountain topography makes the Juneau Icefield less susceptible to climate warming than low-lying Alaskan icefields.

1. INTRODUCTION

Glaciers in Alaska and adjacent Yukon and British Columbia (~87 000 km²; Kienholz and others, 2015) are substantial contributors to sea level rise. During the period 2003–09 their rate of mass loss was $ \sim \!\!50\,{\rm Gt} {{\rm \;a}^{ - 1}}$ , ~1/5 of the mass loss of the global glaciers excluding ice sheets (Gardner and others, 2013). The largest mass losses are found in the maritime climates around the Gulf of Alaska (Arendt and others, 2002, 2013; Berthier and others, 2010). The mass loss of Alaska's glaciers is expected to continue; for example, Radić and others (2013) project losses between 18 and 45% by the end of the 21st century.

Large icefields mantling complex mountainous topographies with strong relief are common in Alaska (e.g. Harding Icefield, Juneau Icefield, Stikine Icefield). Accurate projections of their mass changes are important for assessing the future contribution of Alaska's glaciers to sea level rise. There are two main feedbacks competing in an icefield's response to a warming climate (e.g. Harrison and others, 2001; Huss and Farinotti, 2012). As outlet glaciers retreat, low-lying, high-ablating parts of the icefield are lost, leading to less negative icefield-wide specific mass balances. This stabilizing (negative) feedback competes with a destabilizing (positive) feedback due to glacier-wide thinning commonly referred to as the Bodvardsson effect (Bodvarsson, 1955) or climate-elevation feedback. As the glacier thins, the surface is exposed to higher temperatures leading to increased melt and further thinning. While for relatively thin mountain glaciers in steep terrain the stabilizing feedback usually dominates, ice caps can fall victim to the destabilizing feedback. For different initial ice thicknesses these feedback mechanisms can then lead to different stable states under the same climate forcing. Juneau Icefield combines the characteristics of an ice cap and mountain glaciers leading us to question which of the two feedbacks will dominate under future climate conditions.

So far, due to the global scope of existing studies, projections of Alaska's icefields have relied on simple scaling methods rather than on physically-based flow modeling to account for glacier front variations (Marzeion and others, 2012; Radić and others, 2013). In this study we simulate the recent and future evolution of the Juneau Icefield, the second largest icefield in Alaska, using the open-source Parallel Ice Sheet Model (PISM; Khroulev and the PISM Authors, 2014) forced by output from the regional Weather Research and Forecasting (WRF) Model. PISM is a hybrid shallow ice, shallow shelf model, which has been extensively applied at a wide range of spatial scales. Applications include the entire Northern hemisphere (Ziemen and others, 2014), the Last Glacial Cordilleran Ice Sheet (Seguinot and others, 2014), the Antarctic and Greenland ice sheets (e.g. Martin and others, 2011; Aschwanden and others, 2013), their sub-regions (e.g. Habermann and others, 2013), the New Zealand Southern Alps (Golledge and others, 2012) and individual glaciers (van Pelt and others, 2013). The goal of our work is to project the 21st century evolution of the icefield. We explore the sensitivity of the modeled icefield to model parameters and climatological variables in order to assess the robustness of our results. Finally, we perform a series of model experiments to investigate the stability of the icefield under hypothetical climate stabilization scenarios, and assess the ability to regrow from an ice-free state.

2. PHYSICAL SETTING

The Juneau Icefield is located in the northern Coast Mountains at the border between southeast Alaska and British Columbia (Fig. 1). The icefield presently covers an area of ~3700 km² and ranges in elevation from sea level at the western and southern margins to 2300 m a.s.l. (Kienholz and others, 2015). Oriented northwest-southeast, roughly parallel to the coast, the icefield is split into a maritime western part that receives 3–4 m a^–1 of precipitation (Pelto and others, 2013) and a much drier eastern part. Melt occurs across the entire icefield in summer (Ramage and others, 2000), and rainfall and melt are common in the lower parts in winter. All glaciers of the icefield are presently land- or lake-terminating.

Fig. 1. Modeling domain including the Juneau Icefield and surrounding glaciers (in total ~4150 km² based on present-day ice cover; Kienholz and others, 2015). Colors distinguish individual outlet glaciers. The glaciers south of the Lynn Canal (shown in gray) are excluded in our study.

The largest outlet glacier is Taku Glacier with an area of 735 km² and a length of 60 km (Kienholz and others, 2015). Taku Glacier has experienced tidewater glacier cycles (Meier and Post, 1987) in past centuries that yielded an advance/retreat pattern largely asynchronous with regional glacier change. Following a rapid tidewater retreat for about a century, in 1850 the glacier began a steady advance that continues today, in sharp contrast to the general retreat and thinning of the other glaciers of the icefield (Motyka and Begét, 1996; Truffer and others, 2009). The glacier has built its own dam by pushing a moraine into the Taku Inlet, protecting itself from marine calving (Motyka and Begét, 1996; Kuriger and others, 2006). Therefore the glacier is presently land-terminating. Seismic profiles have revealed local ice thicknesses exceeding 1.4 km and basal troughs reaching down to 600 m below sea level (Nolan and others, 1995).

Llewelyn Glacier is the second largest (450 km², 37.5 km long), and drains substantial parts of the drier eastern side of the icefield. Mendenhall Glacier (109 km², 25 km long) close to Juneau is a well-studied glacier that retreated 3 km during the 20th century (Motyka and Echelmeyer, 2003). The glacier was historically land-terminating but now ends in a proglacial lake. The calving flux into the lake accounts for only 3–6% of the total net mass loss, but plays an important role in the ongoing retreat. Maximum ice thickness near the terminus is ~70–120 m and reaches ~500–600 m further up-glacier (Motyka and Echelmeyer, 2003; Boyce and others, 2007; Molnia, 2007).

3. PREVIOUS WORK

The Juneau Icefield has been the subject of a variety of glaciological studies. Most of them focused on individual outlet glaciers (e.g. Heusser and Marcus, 1964; Pelto and others, 2013) or included the icefield as part of a much larger domain (e.g. Lawrence, 1950; Larsen and others, 2007; Berthier and others, 2010). Only Ramage and others (2000) and Melkonian and others (2014) explicitly dealt with the entire Juneau Icefield as such. Here we review relevant previous studies on mass balance and ice velocity, since these data are used for model calibration.

The most prominent study is the long-term monitoring of surface mass balance with focus on Taku and Lemon Creek glaciers (9.5 km²) by the Juneau Icefield Research program, a student training program in operation since 1946 (e.g. Pelto and others, 2013). Surface mass balance is measured in snow pits in 17 (Taku) and 5 (Lemon Creek) fixed locations between 950 and 1800 m a.s.l. in July and August and adjusted to represent annual balances at the end of the melt season. Ablation stake measurements over a few years revealed net surface mass losses of 12–14 m w.e. a^–1 close to the terminus (Pelto and Miller, 1990; Motyka and Echelmeyer, 2003) consistent with losses of 10–14 m w.e. a^–1 found close to the terminus of neighboring Mendenhall Glacier between 1998 and 2005 (Boyce and others, 2007). Glaciological measurements on Mendenhall Glacier indicate an area-averaged mass balance of $ - 0.95 \pm 0.11\,{\rm m} {\rm \;w}. {\rm e}. {{\rm \;a}^{ - 1}}$ from 1948 to 1995 (Motyka and Echelmeyer, 2003).

Several studies have reported geodetic mass balances over various periods based on sequential DEM comparisons (Table 1). Larsen and others (2007) computed glacier volume changes in southeast Alaska by differencing the Shuttle Radar Topography Mission (SRTM) DEM (acquired in February 2000) with DEMs derived from maps referring to 1948 for Alaska (61% of the Juneau Icefield), and 1982 (17%) or 1987 (22%) for Canada (yielding ~1962 as an area-averaged date). Berthier and others (2010) subtracted the same early-date DEM from a DEM for 2006 derived from SPOT5 and ASTER images, but most of Taku Glacier and large parts of Llewelyn Glacier were excluded due to data gaps. Melkonian and others (2014) estimated mass change between 2000 and 2009/2013 by applying a weighted linear regression on a pixel-by-pixel basis to elevations from ASTER DEMs stacked with the SRTM DEM.

Table 1. Previously reported specific mass balance rates $\dot B$ for different periods for the entire Juneau Icefield and Taku Glacier. The Larsen and others (2007) value for Taku is taken from Melkonian and others (2014)

While there is consensus that the icefield has lost mass over all investigated periods, a less negative rate is found for the period 2000–2009/2013 compared with the earlier periods (Table 1). The difference between the periods can be seen all across the icefield. It must therefore either be caused by different climatic conditions or by systematic differences between the methods. For Taku Glacier the geodetically derived balance rates are generally positive, which is consistent with its current advance, but rates vary between studies. One factor that contributes to the differences between Pelto and others (2013) and the other studies are different outlines used for Taku Glacier (Table 1). The glaciological method yields a positive balance for Taku Glacier for the period 1948–2000 and slightly negative mass change rates for more recent periods, although the lack of glaciological measurements in the ablation zone could be a significant source of error (Pelto and others, 2013).

Melkonian and others (2014) combine ASTER and Japanese ALOS pixel tracking with ERS InSAR-based measurements from various dates between 1995 and 2012 to obtain velocity fields for the icefield. They find maximum speeds of ~1.5 m d^–1 along the main trunk of Taku Glacier, with maximum speeds at other outlet glaciers typically ranging from 0.4 to 0.8 m d^–1, consistent with maximum speeds of ~0.5–0.6 m d^–1 measured by Motyka and Echelmeyer (2003) 1–3 km up-glacier from the terminus of Mendenhall Glacier using GPS between 1997 and 2000.

Using GPS measurements McGee and others (2007) did not detect significant interannual velocity variability in the accumulation area of Taku between 1993 and 2006, but found a slight acceleration along the main trunk at four profiles of 6–7% between 1998/2000 and 2007, an increase too small to be detected by the pixel-tracking results in Melkonian and others (2014). Speeds closer to the terminus of Taku Glacier exhibit pronounced seasonal variability. Truffer and others (2009) found summer speeds ~0.2–0.3 m d^–1 higher than the annual average at its terminal kilometer between 2003 and 2005, consistent with observations along the same profiles between 2009 and 2010 made by Melkonian and others (2014).

4. DATA

PISM requires bed elevation and initial ice thickness fields. In addition, we force the model with monthly-mean near-surface (2 m) air temperature and precipitation fields. For model calibration we use the available mass balance data and velocity fields reported in the literature. Additionally we use equilibrium-line altitudes (ELAs) that we derive from late-summer snow lines for individual glaciers of the icefield from satellite data (Section 4.3).

4.1. Ice thickness and bed elevation

Current ice thickness (referring to the year 2000) was derived by M. Huss (personal communication, 2014) using the method of Huss and Farinotti (2012). This approach uses a highly parametrized mass balance to determine the volumetric balance flux along the glacier, which in turn is used to compute spatially distributed ice thicknesses using the isothermal shallow ice approximation applied to the SRTM DEM. Here the method of Huss and Farinotti (2012) is updated to allow ice thicknesses larger than zero along ice divides between individual glaciers. The glacier outlines are taken from Kienholz and others (2015). The resulting ice thickness distribution is shown in Figure 2. The average thickness is 270 m and the maximum thickness 1283 m. Total volume is 1120 km³. Given the simplicity of the model, the maximum thickness is in good agreement with the more than 1400 m obtained by Nolan and others (1995), and calculated thicknesses compare well with observations on Taku Glacier (Huss and Farinotti, 2012).

Fig. 2. Spatial distribution of ice thickness computed with an updated version of the method by Huss and Farinotti (2012).

We derive the basal topography by subtracting the ice thickness from the SRTM surface topography, and keep it constant. Ice thickness and bed topography grids have the nominal resolution of 30 m of the SRTM DEM.

4.2. Climate data

We use the regional atmospheric model of WRF (Skamarock and others, 2008) to dynamically downscale near-surface air temperature and precipitation (Figs. 3, 4b) from one of the Coupled Model Intercomparison Project Phase 5 (CMIP5) simulations by the Community Climate System Model 4 (CCSM4; Gent and others, 2011). Coarse-scale CMIP5-CCSM4 simulations from the ensemble member Mother of all Runs (MOAR) for the historical period of 1971–2005 and projections for 2006–99 are downscaled with WRF to a 20 km resolution grid covering Alaska and adjacent Canada following the approach by Zhang and others (2007). The projections are forced with the greenhouse gas emission scenario of Representative Concentration Pathway 6.0 (RCP6.0), which is considered a ‘middle-of-the-road’ emission scenario. WRF predicts an increase in near-surface air temperature by 3.6 K averaged over the icefield in summer and by 2.5 K in winter. Changes in precipitation are minor (Table 3).

Fig. 3. Mean summer (April–September), annual, and winter (October–March) near-surface air temperatures from the WRF model over the Juneau Icefield and surrounding glaciers (Fig. 1).

Fig. 4. Winter (October–March) precipitation fields averaged over the period 1971–2000. (a) 2 km resolution SNAP data based on PRISM. (b) 20 km resolution WRF data. (c) Adjusted WRF data; red lines are isolines of the precipitation correction factor field by which the WRF data are multiplied prior to forcing PISM (Section 6.2). All the data are interpolated to the 300 m ice sheet model grid. Black outline shows the present-day ice extent of the Juneau Icefield.

Table 3. Winter (October–March), summer (April–September) and annual mean near-surface air temperature and precipitation averaged over the present-day icefield and four periods. The data are based on WRF with the precipitation adjusted as detailed in Section 6.2. The first column refers to the constant climate split-off (S) and regrowth (R) experiments. ΔV and ΔA refer to the volume and area changes relative to the 2010 state of reference run REF. t₉₀ is the number of model years when 90% of the final volume change has occurred, t_s is the number of model years when a steady state has been reached according to Eqn (4) and is given relative to 2010 for S-experiments and to the start of the run for the R experiments.

In WRF, the processes that occur at the sub-grid scale are described by various physical parametrizations that are selectable by the user depending on their particular needs. The selection of the most suitable and best-performing parametrizations schemes for a particular region is critical for achieving the best possible performance from WRF. To address this challenge, a series of sensitivity tests to optimize the WRF configuration for application over Alaska and surrounding areas have been carried out (Liu and others, 2013; Zhang and others, 2013).

We bilinearly interpolate the WRF downscaled near-surface temperature and precipitation data from 20 km to the 300 m resolution PISM grid. Air temperatures are corrected for the elevation difference between the WRF topography and the evolving surface elevations of the icefield using a constant lapse rate of $ - 5 {\rm K} {\rm \;k}{{\rm m}^{ - 1}}$ (Abe-Ouchi and others, 2007; Ziemen and others, 2014). In the following we refer to this dataset as the WRF data.

We also tested a high-resolution gridded monthly climate dataset (Fig. 4a) generated by the Scenarios Network for Alaska Planning (SNAP), University of Alaska (retrieved March 2014 from http://www.snap.uaf.edu), which has been created specifically for Alaska. SNAP downscales the Climate Research Unit global reconstruction (Harris and others, 2014) and CMIP5 climate model projections to a higher resolution using the PRISM (Parameter-elevation Relationships on Independent Slopes Model) climate mapping system (Daly and others, 2002). PRISM applies a climate/elevation regression for each grid cell. Stations included in the regression are assigned weights based primarily on the physiographic similarity of the station to the grid cell depending on factors such as slope, aspect and coastal proximity.

The SNAP data have a high nominal spatial resolution of 2 km, but the very few weather stations available in the study area are concentrated at lower elevations, leaving PRISM relatively uninformed at higher elevations. Hence PRISM precipitation fields largely depend on altitude (Fig. 4a), resulting in a precipitation maximum on the highest parts of the icefield. This conflicts with our expectations for orographic effects of large mountain chains (e.g. Roe, 2003) which suggest a precipitation maximum on western and southern flanks of the icefield. A precipitation maximum on the flanks is also necessary for our mass balance modeling to match the existing mass balance measurements (Section 6.2).

The WRF data, in contrast, have a resolution of 20 km, but WRF directly models the atmosphere dynamics and places the precipitation maximum in the southwestern part of the icefield. The WRF and SNAP precipitation data differ both in total amounts and spatial distribution. WRF's total winter (October–March) precipitation is 2.2 m averaged over the entire icefield for 1971–2000, 0.40 m less than the SNAP data for the same period. Comparing each grid cell's winter precipitation, the RMS difference between the SNAP and WRF datasets is 0.74 m, indicating substantial differences in spatial pattern. The icefield average near-surface air temperature is higher in the SNAP dataset by 5.7 K in summer (April–September) and 5.8 K in winter (October–March). The RMS is 5.8 K in summer and 6.1 K in winter. The differences are smaller in the low-lying areas (<50 m a.s.l.), where weather stations contribute to PRISM (and thus the SNAP dataset), with summer temperatures 2.9 K, and winter temperatures 0.9 K higher in the SNAP dataset.

4.3. Equilibrium line altitudes

We map snowlines for individual glaciers of the icefield from nine Landsat 5 Thematic Mapper (TM) satellite scenes as well as three Landsat 7 Enhanced Thematic Mapper (ETM+) and four Landsat 8 Operational Land Imager (OLI) scenes, covering the period 1996–2014. We aim at annual (end-of-summer) snowlines but due to frequent cloud cover we analyze all cloud-free images from August to September. In case several snowlines are available for the same glacier and year, we retain the highest. We determined each snowline's characteristic altitude by calculating the median of the SRTM altitudes within 30 m of the snowline. Since we assume superimposed ice to be negligible, the snow line altitudes are equivalent to ELAs (Cogley and others, 2011), thus allowing us to compare them with modeled ELAs for model calibration. The method of deriving the snowlines is detailed in Appendix A.

5. MODEL DESCRIPTION AND APPLICATION

5.1. The Parallel Ice Sheet Model

The Parallel Ice Sheet Model (PISM) includes a hybrid stress balance model (Bueler and Brown, 2009; Aschwanden and others, 2012) combining the Shallow Ice Approximation (Hutter, 1983) for vertical deformation and the Shallow Shelf Approximation (Morland, 1987) for longitudinal stretching. The ice rheology is described using Glen's flow law

(1) $$D = E A \vert \tau { \vert ^{n - 1}}\tau, $$

where D is the strain rate tensor, E is the flow enhancement factor, a tuning parameter used to account for unresolved aspects of the ice flow, A is the rate factor, τ is the deviatoric stress tensor and n = 3 is the flow law exponent. The whole Juneau Icefield is assumed to be temperate; consequently we run PISM in an isothermal configuration and use a constant rate factor of $A = 4.5293 \times {10^{ - 24}} P{a^{ - 3}} {s^{ - 1}}$ , the result of the Paterson and Budd (1982) rate factor equation for ice at the melting point. We refrain from lumping A and E into one single parameter in Eqn (1) to facilitate the interpretation of the corresponding sensitivity experiments.

We parametrize basal sliding using a power law with a quadratic dependence of the basal velocity ${\vec u_{\rm b}}$ on the basal stress ${\vec \tau _{\rm b}}$ ,

(2) $${\vec u_{\rm b}} = - {u_{\rm c}}\displaystyle{{ \vert {{\vec \tau} _{\rm b}} \vert {{\vec \tau} _{\rm b}}} \over { \vert {\tau _{\rm c}}{ \vert ^2}}},$$

where τ _c (units of stress) and u _c (units of velocity) are two tuning parameters representing only one degree of freedom. Using two parameters instead of one simplifies the interpretation of the sliding law. We set ${u_{\rm c}} = 500\,{\rm m} {{\rm a}^{ - 1}}$ , while τ _c, a pseudo-yield-stress, is part of the calibration procedure (Section 6.1).

Surface mass balance is computed as the difference between accumulation and ablation. Precipitation is split into solid and liquid fractions using a linear transition between − 10 and +7°C monthly mean surface air temperature based on observations analyzed by (Marsiat, 1994). The large range is explained by the use of monthly mean data rather than instantaneous temperatures. The solid fraction of precipitation is treated as snow accumulation and converted to ice at the end of each mass balance year (end of September). The liquid fraction is discarded as runoff. Ablation is computed from monthly near-surface air temperature using a positive degree day (PDD) scheme (Hock, 2003). Temperatures within each month are described by a normal distribution and the PDD integral is evaluated analytically using the error function (Calov and Greve, 2005). Following Ziemen and others (2014) we implement this approach by providing PISM with monthly maps of temperature mean and standard deviation computed from the hourly WRF output.

The model is set up on Alaska Albers equal area coordinates with a horizontal resolution of 300 m. Interpolation from the SRTM grid to the PISM grid is done conservatively. The modeling domain covers the Juneau Icefield and surrounding unconnected glaciers (Fig. 1). For this domain, the adaptive time stepping yields time steps of ~8 h. Each run took ~36 h on 64 CPUs including spinup. All simulations used PISM version 0.6. Due to the limited spatial and temporal scale of the experiments we did not model isostatic rebound.

5.2. Model spinup

We perform a multi-stage spinup to avoid spurious mass changes resulting from coupling shocks at the start of the experiments. These would interfere with the tuning and validation. As a starting point for the spinup, we compute an initial ice thickness field by projecting the year 2000 ice thickness field (Section 4.1) back in time by 60 a to 1940 using the thickness change rates from Larsen and others (2007). As climate forcing we use the 1971–2000 climate data, because this is the first 30 a interval for which we have downscaled the climate data available. We choose 30 years to average out short-term variability.

During the spinup, we nudge the model to the 1940 ice thickness using Newtonian relaxation (e.g. Jeuken and others, 1996). After each time step the deviation of modeled ice thickness from the targeted ice thickness ΔH of each grid cell is scaled by a relaxation factor γ, and that value is subtracted from the following time step's surface mass balance. For example, for γ = 0.1 a⁻¹ and a thickness deviation of ΔH = 1 m, the surface mass balance rate computed with the PDD scheme is reduced by $\gamma \Delta H = 0.1 \,{\rm m}{{\rm \;a}^{ - 1}}$ .

We run the model for 90 years repeating the 1971–2000 climate data three times. For the first 60 years we employ a relaxation constant $\gamma = 0.1 \,{{\rm a}^{ - 1}}$ . This allows the icefield to evolve into a state that is consistent with the model physics, while keeping the surface elevations very close to the reconstructed 1940 state. For the remaining 30 years we decrease the constant to $0.01 \,{{\rm a}^{ - 1}}$ to arrive at the 1971 initial thickness field. The reduction of γ in the last cycle allows the icefield to evolve almost freely to avoid a coupling shock when starting the experiments. The short spin-up is possible because the ice flow does not depend on ice temperature in our setup.

6. CALIBRATION

The model is calibrated by performing hindcasts for the period 1971–2010 and comparing the hindcasts with four sets of observations: surface speeds (Melkonian and others, 2014), the mean 1946–86 surface mass balance gradient (Pelto and Miller, 1990), our satellite-derived transient ELAs for individual outlet glaciers and multi-year geodetic mass change (Table 1). Since the WRF data are produced by a free running climate model, we cannot expect individual model years or even decadal means to line up with corresponding observations, and therefore only compare long-term averages. The model run using the calibrated model parameters is referred to as reference run (REF) in the following.

6.1. Ice flow

To calibrate the two ice dynamics parameters, E (Eqn (1)) and τ _c (Eqn (2)), we run PISM's surface relaxation scheme with a relaxation constant of $0.1 \,{\rm m}{{\rm \;a}^{ - 1}}$ . This isolates the effect of the flow parameters from the effects of the surface mass balance and the surface evolution. At the same time, the surface can slightly adjust to the flow, thereby smoothing out local inconsistencies between prescribed surface and flow physics. We find E = 1 and τ _c = 0.5 MPa provide reasonable agreement between the modeled and observed ice velocities (RMS difference of $69 \,{\rm m}{{\rm \;a}^{ - 1}}$ ; Fig. 5).

Fig. 5. (a) Measured surface speed (Melkonian and others, 2014), (b) modeled surface speed and (c) modeled basal sliding speed.

6.2. Surface mass balance

In a first step we calibrate the surface mass balance model by adjusting the PDD factors for ice and snow to match the observed mass balance profile over Taku Glacier (Fig. 6).

Fig. 6. Surface mass balance vs surface elevation for Taku Glacier. The black curve indicates the average profile derived from measurements between 1948 and 1986 (Pelto and Miller, 1990). Green dots mark modeled surface mass balances of all Taku Glacier grid cells averaged over the period 1971–2000.

Using the 2 km resolution SNAP dataset (Section 4.2; Fig. 4a), it was impossible to reproduce the observed mass balance profile, while keeping PDD factors within reasonable ranges. While the SNAP data show high winter precipitation, especially at high altitudes (Fig. 4a), it also shows relatively high air temperatures over the icefield (Supplementary Figure S1), yielding a too negative mass balance at all elevations, even for unrealistically low PDD factors (Supplementary Figure S2a). We investigated possible systematic temperature biases. To increase the surface mass balance, we decreased the temperature by 2.5 K uniformly in the whole domain. With this temperature field, a reasonable fit to the Pelto and Miller (1990) mass balance profile could only be obtained for a snow PDD factor of $0.5{\kern 1pt} \,{\rm mm}\,{{\rm K}^{ - 1}}{\kern 1pt} {{\rm d}^{ - 1}}$ (Supplementary Figure S2b) and an ice PDD factor of $14{\kern 1pt} \,{\rm mm}\,{{\rm K}^{ - 1}}{\kern 1pt} {{\rm d}^{ - 1}}$ . However, the PDD factor for snow is unrealistically low and indicates that the high altitudes are still too warm. We therefore decreased the air temperatures even further by a total of 5 K and found best agreement with the observed mass balance profile using ${\rm PD}{{\rm D}_{{\rm snow}}} = 2\,{\rm mm}\,{{\rm K}^{ - 1}}{\kern 1pt} {{\rm d}^{ - 1}}$ and ${\rm PD}{{\rm D}_{{\rm ice}}} = 14\,{\rm mm}\,{{\rm K}^{ - 1}}{\kern 1pt} {{\rm d}^{ - 1}}$ (Supplementary Figure S2b). These factors are too far apart from each other to be plausible (Hock, 2003). They indicate a too small temperature difference between accumulation and ablation areas. Overall, this provides further circumstantial evidence that the spatial distribution of air temperature and precipitation over the icefield is not realistic in the SNAP data.

In contrast, we are able to obtain a good fit to the observed mass balance profile by using the dynamically downscaled WRF data (Section 4.2). The fit to the Pelto and Miller (1990) mass balance profile is best for snow and ice factors of 4 and 10 mm w.e. a^–1 respectively, consistent with the range of values typically found in the literature (Hock, 2003). Observed ELAs range between ~800 m a.s.l. in the western part and ~1600 m a.s.l. in the eastern part (Supplementary Figure S3). Compared with the observations, modeled ELAs are considerably higher in the western part and lower in the eastern part (Fig. 7b). We attribute this discrepancy to the WRF model not resolving the pronounced precipitation gradient across the icefield in the west-east direction, which is caused by orographic effects on a spatial scale too small to be resolved by the climate model (Fig. 4b).

Fig. 7. Modeled annual ELAs for the period 1990–2019 vs observed ELAs for the period 1996–2014 for individual glaciers of the Juneau Icefield. The time window for the model data is chosen as 30 years, centered on the observation period. Dark filled circles show the median of all available annual modeled and observed values for each glacier. Large circles mark glaciers >100 km². Small dots connected by lines mark the ELAs for each individual year (blue: modeled; green: observed). Also shown are 1 : 1-line (light gray) and linear fit of the medians (dark line). (a) Model run REF with precipitation adjustment to account for systematic bias (Section 6.2). (b) Model run NOGRAD using the WRF data without precipitation adjustment.

We account for this apparent systematic bias by adjusting the precipitation fields so that precipitation near the coast is increased, while it is decreased on the eastern side of the icefield. Similar to Aðalgeirsdóttir and others (2003), we fit a plane, $z = a x + b y + c$ with a = 8.41 × 10⁻³, b = 9.36 × 10⁻³ and c = −58 616 m through the observed ELAs, where x and y are the coordinates in Alaska Albers projection. From this field z, we derive a precipitation correction factor ζ for each grid cell:

(3) $$\zeta = \alpha (a x + b y + \beta ) + 1.$$

This parametrization maintains the direction of the gradient, but adjusts the slope and the offset to obtain a field with values decreasing from ~1.5 near the western margin to 0.5 near the eastern margin of the icefield (red lines in Fig. 4c). The offset $\beta = 1200 \,{\rm m} - c = - 59816\,{\rm m}$ was chosen to roughly conserve the total amount of winter precipitation over the present-day icefield, and to provide a smooth gradient across the whole ice-covered area. The slope factor α is determined from tuning. The optimized value is $\alpha = - 1.5 \times {10^{ - 3}} {{\rm m}^{ - 1}}$ . The line for ζ = 1, where precipitation stays unaltered, is independent of α. To maintain a minimum precipitation in the northeast, all values of ζ below 0.1 are set to 0.1. Multiplication with this field results in a new adjusted precipitation field (Fig. 4c) used as input for PISM. Forcing the model with this precipitation field yields ELAs that correlate well with observed ELAs (Fig. 7a) while maintaining a good fit to the observed mass balance profile on Taku Glacier (Fig. 6).

The modeled ELAs are on average ~100 m above the observed ones possibly because the observed ELAs represent snapshots in late summer, but not necessarily the maximum level during the season. They therefore provide a lower estimate of the annual ELAs. Lowering the modeled ELAs by reducing the degree-day factors (sensitivity experiment SMELT↓, Section 7.2) leads to a better fit (Supplementary Figure S4), but also increases the icefield's total mass balance for the period 1971–2010, which is inconsistent with observations (Table 1). With the slightly elevated ELAs, the modeled icefield-wide mass balance rate for the period 1971–2010 is consistent with the observed glacier-wide balances found in previous studies. The spatial pattern of the results during the reference period largely agrees with the observations, (Supplementary Figure S5) although the model shows unrealistic growth on Gilkey Glacier. Due to the large discrepancies we put less emphasis on matching the exact geodetic mass balance observations, but aim to reproduce the large-scale spatial patterns of change with negative balance rates for the calibration period within the observed range of −0.62 to $ - 0.13{\kern 1pt} \,{\rm m\ w}.{\kern 1pt} {\rm e}.{{\rm \;a}^{ - 1}}$ (Table 1) for the entire icefield, while maintaining a positive mass balance for Taku Glacier (Supplementary Figure S5).

7. RESULTS

7.1. Icefield evolution 1971–2099

The ice volume evolution of the reference run REF for the period 1971–2099 is shown in Figure 8 (black curve). During the calibration period 1971–2010, ice volume increases slightly by 2% between 1971 and 1996, which is followed by a pronounced decrease, leading to an overall volume loss for 1971–2010 of $62{\kern 1pt} \,{\rm k}{{\rm m}^3}$ (5%). The corresponding specific balance rate of $ - 0.33{\kern 1pt} \,{\rm m}\,{\rm \; w}.{\kern 1pt} {\rm e}.{{\rm \;a}^{ - 1}}$ falls within the range of reported multi-decadal estimates (Table 1). The decrease in ice volume following 1996 is consistent with simultaneous pronounced summer warming (Fig. 3). Mass gains are modeled on Taku and Gilkey Glacier (Supplementary Figure S5).

Fig. 8. Modeled evolution of annual ice volumes over the period 1971–2100 for the reference run and the sensitivity experiments (Table 2). The legend is sorted by final ice volume. Experiments in which the same parameter is varied share the same color.

From 2010 to 2040 volume loss rates are almost linear and all outlet glaciers except for Gilkey and Taku are shrinking. Volume loss rates accelerate and by 2070, most outlet glaciers have negative glacier-wide balances. By the end of 2099 the entire icefield has lost 63% of its volume and 62% of its area compared with 2010. Most of the outlet glaciers have retreated substantially and the northern parts of the icefield have disappeared (Fig. 9d). In contrast Taku and Gilkey Glacier have thinned substantially, but show only a slight retreat. The surface mass balance averaged over the period 2070–99 is negative everywhere except for the high-elevation areas of Taku and Llewelyn glaciers. Despite their positive surface mass balances, these regions are losing mass by dynamic thinning.

Fig. 9. Modeled ice thicknesses in year (a) 2011, (b) 2041, (c) 2071 and (d) 2099.

7.2. Uncertainty and sensitivity analysis

Uncertainties in projections arise from inaccuracies in input data, inadequacies in model physics and the calibration procedure. Error propagation cannot be applied in this study because the validity of inherent assumptions, such as normality or independence of parameters, cannot be tested. Therefore, we evaluate uncertainties by assessing the sensitivity of projected volume changes to the choice of model parameters: the flow enhancement factor E (Eqn (1)), the yield stress τ _c (Eqn (2)), the PDD factors for snow and ice, the slope factor α (Eqn (3)), the initial ice thickness and the grid resolution. For most parameters or variables we analyze the effects of a positive and negative change relative to the reference run REF and leave all other parameters unaltered (Table 2). For the ice thickness experiments, the initial ice thickness field (Section 5.2) is multiplied by a constant factor prior to spin-up, and basal topography is adjusted accordingly, leaving the surface elevations unchanged.

Table 2. Sensitivity experiments including varied parameters, and key results. f _α is a factor by which α is multiplied to vary the strength of the precipitation adjustment (Eqn (3)), PDD_s and PDD_i are the positive degree-day factors for snow and ice, respectively, E and τ _c are flow parameters (Eqns (1) and (2)), f _thk is the factor by which the initial ice thickness field is multiplied. Displayed results are icefield-wide specific surface mass balance rate $\dot B$ , volume change ΔV and area change ΔA for the calibration period 1971–2010 and the projection period 2011–99, and 2099 volume relative to the volume of Ref in 2099 V _rel. Parameters varied in each experiment are in bold. Modeled $\dot B$ consistent with the observed range (–0.62 to –0.13 m w.e. a^–1, Table 1) are in italic.

Sensitivity experiments are carried out by running each experiment first through the spin-up and then from 1971 to 2099. Consequently, each experiment has a different ice volume at the beginning of the prognostic simulation. This procedure avoids shocks in the model that would result from sudden parameter changes.

Initial ice volumes in 1971 are within ~10% of the REF volume and most experiments show similar patterns of volume change throughout the simulation period, 1971–2099. Experiments that start at an elevated ice volume in 1971 also stay at an elevated volume throughout the experiment (and vice versa). In most cases there are no substantial deviations in the pattern of changes from those of the reference run. For those experiments where there are significant deviations, we list them in the following and show characteristic examples in Figure 10. Section 8.3 includes a discussion of the results.

Fig. 10. Ice thickness differences in 2099 for two sensitivity experiments. (a) NOGRAD−REF. (b) RES ↓ −REF. Note the different color scales.

Experiments that strongly deviate from the volume loss rate in REF are those where we changed the ice thickness. Experiment THICK starts with 18% greater ice volume, but due to a substantially larger volume loss rate, the ice volume by 2099 is similar to that of REF. The increased ice thickness causes an acceleration in ice flow and therefore increased transport of ice to lower, high-ablating elevations. In addition the decreased basal topography allows for more glacier area at lower elevations. Despite a considerably lower initial ice volume (18%) than REF, the experiment THIN also drops to a similar 2099 ice volume as REF. The year 2099 volumes of THIN and THICK are within 2% of REF.

Varying the enhancement factor by a factor of 1.5 (SLOW, FAST) shows only modest impact on the evolution of the ice volume. The differences in the 2100 ice thickness resulting from variations in the flow parameters are spatially very homogeneous, and do not show any substantial differences between outlet glaciers and the interior of the icefield (not shown).

Varying the PDD factor for ice by ±1 mm K^–1 d^–1 (IMELT↓, IMELT↑) has little effect on the ice volume evolution. In contrast, varying the snow PDD factor by the same amount strongly affects the ice volume evolution (SMELT↓, SMELT↑).

Compared with REF, running the model without the precipitation adjustment (NOGRAD) results in a 55% higher ice volume in 2100 and yields considerably greater ice thicknesses and ice extent in the north-eastern side of the icefield, in particular on Llewelyn Glacier, and thinner ice in the south. As a result the icefield is shifted to the northeast (dipole pattern in Fig. 10a). In contrast, a higher gradient (GRAD↑) in the precipitation field causes a westward migration (not shown).

Decreasing the model grid resolution from 300 to 600 m (RES↓) leads to an increase in ice thickness, especially on narrow outlet glaciers (Fig. 10b). As outlet glaciers in narrow valleys are less well resolved, the ice slows down leading to dynamic thickening at higher surface elevations, which in turn increases the icefield-wide mass balance. The unrealistically strong growth of Gilkey Glacier may also result from more inflow from the interior of the icefield. The 600 m resolution grid does not properly resolve the mountain ranges separating the glacier from the interior of the icefield and more ice can spill over into the valley. This is consistent with a band of no or negative thickness change around Gilkey Glacier (Fig. 10b). Similar effects may explain the greater ice thicknesses on Mendenhall Glacier.

Eight of the sensitivity experiments yield a mean mass balance rate over the calibration period 1971–2010 that is within the geodetically derived rates of −0.62 to $ - 0.13\,{\kern 1pt} {\rm m}\,{\rm \;w}.{\kern 1pt} {\rm e}.{{\rm \;a}^{ - 1}}$ (Table 2). These experiments also show reasonable agreement between modeled and observed ELAs, Taku Glacier's mass balance profile (Supplementary Figure S4) and observed patterns of mass change during the calibration period (Supplementary Figure S5). Hence, we consider these projections equally plausible as the reference run. Volume changes from 2010 to 2099 for these runs range between 58 and 68%, and area changes range between 57 and 63%.

7.3. Constant-climate split-off experiments

All projections above are based on the RCP6.0 scenario and show continuous volume losses throughout the 21st century, suggesting a likely disappearance of the icefield, if the temperature and precipitation trends were to continue. Here we investigate whether or not the icefield is able to survive and reach a new steady state, if air temperature and precipitation stabilize earlier in this century. We perform four constant-climate split-off experiments. In the first experiment PISM is forced repeatedly with the climate of the calibration period (1971–2010). Three additional experiments start in 1971 with the transient climate scenario but the model is forced repeatedly with the 2011–40, 2041–70 and 2071–99 time series beyond the years 2040, 2070 and 2099 respectively. We refer to these four experiments as S2010, S2040, S2070 and S2099 in the following. Simulations run for 1500 years or until a steady state is reached, whichever occurs first. We define that a steady state has been reached at time t, if the criterion

(4) $$\displaystyle{{\sum\nolimits_{i,j} { \vert {H_{i,j}}(t) - {H_{i,j}}(t - \Delta t) \vert}} \over {\Delta t\sum\nolimits_{i,j} {1/2({H_{i,j}}(t) + {H_{i,j}}(t - \Delta t))}}} \lt \displaystyle{{1\%} \over {100 {\rm a}}},$$

is fulfilled, where H _i,j is ice thickness at grid point (i,j) and t and t − Δt are the end times of two consecutive climate forcing cycles. This criterion is different than requiring the volume change rate to be lower than 1% per century. The use of the absolute value of the local differences prevents a steady state from being inferred when opposing mass trends in different areas of the icefield cancel.

Key results and the statistics on temperature and precipitation for the four stabilization scenarios are given in Table 3 and Figures 11 and 12. The experiments are forced with segments of the transient climate, and therefore experience cycles of climate warming that are followed by step-like cooling at the beginning of the following cycle. This leads to cyclic growth and retreat of the icefield during the split-off experiments (Fig. 11).

Fig. 11. Modeled annual ice volumes for the constant climate split-off (S) and regrowth (R) experiments. Experiments with the same climate share shades of the same color. S-experiments start with present-day ice volumes, while R-experiments start from an ice-free state at nominal year 2000. Except for REF, the legend is sorted by final ice volume. Runs end either when a steady state is reached (Table 3) or in year 3500, whichever comes first.

Fig. 12. Modeled ice thickness (m) at the end of constant-climate split-off experiments (a) S2010, (b) S2040, (c) S2070 and (d) of the regrowth experiment R2070.

The experiment forced with the climate from the calibration period (S2010) shows a stabilization of the icefield volume at slightly above 1000 km³, corresponding to a loss of 14% of the ice volume in 2010. Substantial area losses and glacier retreats are projected, in particular, in the northern part of the icefield but also in the east (Fig. 12a). Gilkey and Taku Glacier are projected to advance.

Experiments S2040 and S2070 show substantial volume and area losses, but the icefield survives in the higher reaches of the mountains, reaching a steady state after 1110 and 1050 years respectively. Most of the volume changes occur much earlier. While it takes S2040 460 years to reach 90% of the volume change, S2070 passes that mark after 280 years. The earlier stabilization in S2070 compared with S2040 occurs because the warmer climate causes an expansion of the ablation areas, which results in a more rapid retreat to higher elevations than in the cooler scenario. Assuming the 2071–99 climate remains constant (S2099), the icefield's volume rapidly declines to 4% of its 2010 volume by 2220, and continues to decrease thereafter, but at a greatly reduced rate, with some ice still remaining after 1500 years.

7.4. Regrowth experiments

To further investigate the multi-stability of the icefield, and to shed light on the potential of the Juneau Icefield for glaciation, we also perform regrowth experiments (R2010, R2040, R2070). Experiments start from zero ice thickness and are run with repeated climate segments of the corresponding split-off experiments. A regrowth experiment with the 2071–99 climate is not performed since this climate leads to the destruction of the icefield in the split-off experiment S2099. To simplify interpretation we set all experiment start times to year 2000. In all the three experiments, the icefield nucleates in the high-elevation areas and subsequently fills lower-elevation valleys. R2010 and R2040 reach states that are very similar to S2010 and S2040 respectively, and would converge further for longer experiment duration. In contrast, R2070 grows to a volume 24% smaller than S2070 since it is not able to re-occupy the deep troughs of Taku and Llewelyn valleys (Figs. 11, 12c, d).

8. Discussion

8.1. Icefield evolution

Our reference run, forced by the RCP6.0 emission scenario, indicates substantial volume losses by 2100 (64% relative to 2010). Volume losses range from 58 to 68% for the eight sensitivity experiments that are within observational uncertainty in the calibration period. This exceeds the range of 18–45% projected for all Alaskan glaciers with 14 general circulation models for RCP4.5 and RCP8.5 scenarios by Radić and others (2013).

Under the 2071–2100 climate (S2099) the icefield almost completely disappears. All other climate split-off experiments indicate substantial volume losses but the icefield is able to reach a new steady state by retreating to higher elevations.

Assuming the present climate conditions, the icefield stabilizes at a volume only 14% lower than the 2010 volume, which is in sharp contrast to the low-elevation Yakutat Icefield (342 km², ~200 km to the northwest). When forced by projections from the same climate model and emission scenario, Yakutat Icefield collapses completely by 2070 and, under the present climate conditions, by 2110 (Truessel and others, 2015). The different behavior is due to distinctly different topographies. The Juneau Icefield is able to reach a new steady state due to significantly higher maximum elevations where the glaciers can retreat (2000 vs 650 m a.s.l.). Thus, in Yakutat Icefield, the destabilizing mass balance – elevation feedback dominates, while in Juneau Icefield the stabilizing feedback prevails.

Overall, our experiments show that the icefield is not simply bi-stable as climate warms, but exhibits a variety of possible retreat states. These retreat states are largely controlled by the climate and are independent of the initial geometry. In two out of three scenarios, the icefield state obtained by starting from the present-day geometry is practically identical to the state obtained by starting from a completely removed icefield. Only in one of the three climate forcings (R2070), some parts of the icefield cannot be regrown due to the elevation-mass balance feedback.

The survival of parts of the icefield under warming scenarios and the regrowth from an ice-free state are possible because of the complex topography of the Juneau Icefield. There are a multitude of narrow valleys and several regions where the glaciers can retreat up into high-elevation mountains, which allows them to maintain areas with positive surface mass balance even under the applied warming scenarios. However, even under the fairly moderate S2040 climate scenario, the icefield loses two-thirds of its volume and more than half of its area in the long-term. This shows that, independent of the exact climate trajectory, we have to expect large glacier mass changes in this region.

8.2. Model performance and input data

Good agreement with surface mass balance observations was not achieved using unaltered SNAP and WRF climate data directly, although both datasets were specifically designed for Alaska. The 2 km SNAP data product is based on extrapolation of weather station data, but largely relies on elevation dependency due to sparsity of weather stations in the region. Hence it cannot resolve the intricate precipitation patterns in this region. WRF simulates the dynamics of the atmosphere, and thus more adequately captures spatial patterns other than elevation dependence. However, the grid spacing is too coarse to fully resolve the slope precipitation on the steep flanks of the Juneau Icefield. It therefore cannot capture the strong gradient in precipitation from the coast to the interior.

Only by substantially increasing the WRF precipitation in the west and decreasing it on the eastern parts of the icefield were we able to reproduce the observed surface mass balance pattern. Our study emphasizes the need for high-elevation weather station data to better inform gridded climate products in the study region, in addition to climate modeling at a spatial resolution sufficient to capture orographic effects in southeastern Alaska's complex topography. A moisture tracking scheme (e.g. Smith and Barstad, 2004), forced with regional climate model output, might prove to be a good compromise between the computational cost of running regional climate models at km-scale resolution and the need for representing the physics of the cloud processes. Observations of the cross-icefield surface mass balance gradient over several years are necessary to provide better tuning and validation targets.

With the adjusted WRF precipitation PISM's surface mass balance scheme is able to reproduce the observed surface mass balance profile on Taku Glacier and the spatial variability in ELAs across the icefield, although modeled ELAs are on average ~100 m above the observed (Figure 7 Supplementary Figure S3). The modeled mass gain for Taku Glacier and mass loss for the remainder of the icefield is consistent with observations (Supplementary Figure S5). However, available observations, even for identical periods, differ greatly from each other (Table 1), emphasizing the need for further studies on the present state of the icefield to better constrain glacier models.

Although it mimics the basic pattern of orographic precipitation, our method of fitting a plane is simplistic and only serves as a first-order approximation. For example, the modeled rapid deglaciation in the area of Meade Glacier (Fig. 9) may be overestimated as a result of exaggerated reduction in precipitation in this area, although satellite images show the glacier in an unhealthy state with the terminal section (2.3 km) disintegrating in a newly formed proglacial lake between 1986 and 2014 and the loss of tributaries.

The present-day ice thickness was derived by an updated version of the method of Huss and Farinotti (2012) using the shallow ice approximation to estimate the ice thickness necessary to keep the ice flux in equilibrium with a highly parametrized surface mass balance. Ice flow grows approximately with the fifth power of the ice thickness (Eqn (5.110) of Greve and Blatter, 2009); therefore, an over- or underestimation of the surface mass balance, and thus the ice flux, by a factor of 2 only leads to a thickness error of about 15%. Our experiments involving varying the initial ice thickness by $ \pm 20\% $ , and adjusting the basal topography accordingly, show that these systematic changes yield final ice losses within or close to the range of losses of the other sensitivity tests. Furthermore, increasing the ice thickness, and thus the ice losses (THICK), leads to increased volume losses in the calibration period and in the projection period, and vice versa for THIN. The uncertainty due to ice thickness on the ice losses is therefore within the uncertainty covered by calibrating the model against geodetic mass balance estimates.

While ice thickness estimates could be further improved by assimilating more precise mass balance estimates and velocity observations in mass conservation approaches, these methods provide minimal constraints on the ice thickness at the ice divides. This kind of information needs to be collected with electromagnetic, or in areas with especially thick, warm ice, with seismic transects (Nolan and others, 1995). Better information on the ice thickness at the ice divides is crucial for improving future projections of the evolution of the locations of ice divides, and thus of the distribution of meltwater runoff into the different drainage basins. In the case of the Juneau Icefield this can make the difference between meltwater flowing directly into the Gulf of Alaska or following the Yukon River into Bering Strait.

In our model runs, Gilkey Glacier advances, while observations indicate retreat. This may be attributed to overestimation of the influx of ice from the high plateaus of the icefield in combination with the elevation-mass balance feedback. The surface topography of the accumulation area of Gilkey Glacier is very complex and not fully resolved in the model. Since Gilkey Glacier grows during the calibration period and in S2010, while it should be retreating in both cases, we expect our experiments to over-estimate the volume of Gilkey Glacier.

8.3. Model sensitivity

We find that the 21st century volume and area reduction is not very sensitive to the variations applied to the flow parameters or to doubling the grid cell spacing. The transport of ice from the accumulation to the ablation area is dominated by shear flow through steep narrow valleys. The enhancement factor is roughly a constant multiplier of the velocity field, thus a doubling in speed leads to a doubling in ice flux. Increased drawdown is, in part, compensated by thickness changes that negatively feed back on the ice velocity. Due to this negative feedback, the modeled volume changes only weakly depend on the details of the ice dynamics, as long as they are kept consistent throughout spinup, calibration and projection. Supports the findings of the simplified-geometry experiments by Leysinger Vieli and Gudmundsson (2004) that, in the absence of significant basal sliding, the response of alpine glaciers to climatic changes can be projected by lower-order approximations to the Stokes Equations. For the Juneau Icefield, the changes in projected surface air temperature lead to an expansion of the ablation area across most of the icefield (Fig. 13) and decrease the importance of ice transport.

Fig. 13. Modeled surface mass balance over the ice-covered area (m w.e. a^–1) (a) averaged over the period 1971–2000, and (b) averaged over the period 2070–99.

While varying the PDD factor for ice by ±1 mm K^–1 d^–1 (IMELT↓, IMELT↑) has little effect on the ice volume evolution, the modeled volume is highly sensitive to changing the PDD factor for snow. This is probably due to the large fraction, at least initially, of area above the ELA, typical of ice caps and icefields, to which the snow factor is applied. In addition the effect on the ice factor is expected to be lower since variation of the PDD factor by an equal amount has a smaller effect on a larger PDD factor.

The experiments that show strongest deviations from REF in the final ice volume (SMELT↓, SMELT↑, NOGRAD, NOSLIDING, SLIDING↑) also show strong deviations of the same sign in the ice volume evolution during the calibration period. The effect of these parameters can therefore be constrained by tuning the model against mass balance observations.

9. CONCLUSIONS

We used PISM to model the evolution of the Juneau Icefield over the 21st century. The eight model runs that agreed well with observations for 1971 to 2010 project volume losses of 58–68% and area losses of 57–63% between 2010 and 2099, when forced with CCSM4 climate downscaled by WRF and driven by the RCP6.0 emission scenario. This scenario projects a near-surface air temperature increase of 3.6 K averaged over the present-day icefield in summer and of 2.5 K in winter. The modeled loss exceeds the Alaska-average projected losses of Radić and others (2013) under RCP4.5 as well as RCP8.5. The outlook still is more positive than the complete destruction of neighboring Yakutat icefield projected by Truessel and others (2015) using the same climate model and emission scenario.

Sensitivity experiments suggest that within reasonable ranges of parameters, model results are generally more sensitive to the choice of the surface mass balance parameters than the flow parameters, indicating a dominance of surface mass balance processes in the icefield's evolution under the applied warming scenario. Uncertainties in icefield evolution may be most effectively reduced by better constraining the PDD factor for snow. Hence, surface mass balance observations are essential to constraining the model.

We were not able to calibrate PISM with the SNAP data, which is based on spatial interpolation of observations. Calibration of the model using WRF data was successful, but only after correcting precipitation using a scaling that substantially increases precipitation on the western side of the icefield and reduces it on the eastern side. Our results suggest that both the SNAP and WRF data do not adequately represent the spatial pattern of precipitation over the icefield. This emphasizes the need for more measurements in this data-scarce region. In addition, climate models should be run at a spatial resolution adequate to resolve the strong precipitation gradient across the mountain range.

Additional model experiments indicate that the icefield almost completely disappears if the climate towards the end of the century is assumed to continue. If the climate stabilizes earlier, the icefield loses substantial volume but eventually reaches a steady state. The steady state icefield geometries are largely independent of the initial icefield geometry. Very similar icefield geometries can be obtained by starting the model from the present-day geometry and from a completely ice-free state. The bed topography of Juneau Icefield is sufficiently high to allow the stabilizing mass balance feedback due to glacier retreat to counteract the destabilizing mass balance elevation feedback due to thinning.

In summary, the simulations show that the Juneau Icefield will lose substantial amounts of mass during the 21st century. The most relevant factors for improving projections are spatially distributed mass balance measurements for constraining the model and improved climate projections that resolve the local temperature and precipitation patterns. With these datasets, more sophisticated methods that exploit the full spatial information, could be employed for tuning the model.