Modeling the impacts of climate change on mass balance and discharge of Eklutna Glacier, Alaska, 1985–2019

Abstract

Alaska's largest city, Anchorage, depends on Eklutna Glacier meltwater for drinking water and hydropower generation; however, the 29 km² glacier is rapidly retreating. We used a temperature-index model forced with local weather station data to reconstruct the glacier's mass balance for the period 1985–2019 and quantify the impacts of glacier change on discharge. Model calibration involved a novel combination of in situ, geodetic mass-balance measurements and observed snowlines from satellite imagery. A resulting ensemble of 250 best-fitting model parameters was used to model mass balance and discharge. Eklutna Glacier experienced a significant negative trend (−0.31 m w.e. decade⁻¹) in annual mean surface mass balance (mean: −0.62 ± 0.06 m w.e.). The day of the year when 95% of annual melt occurs was five days later in 2011–19 than in 1985–93, demonstrating a prolongation of melt season (May–September). Modeled mean specific discharge increased at 0.14 m decade⁻¹, indicating peak water, the year when annual discharge reaches a maximum due to glacier retreat, has not been reached. Four of the five highest discharge years occurred since 2000. Increases in discharge quantity and melt season length require water resource managers consider future decreased discharge as the glacier continues to shrink.

1. Introduction

The mountain cryosphere provides clean drinking water and hydropower resources to an estimated 1.6 billion people living downstream of mountainous areas (Immerzeel and others, 2020). The thinning and retreating of glaciers worldwide can have an immediate socio-economic implication in addition to the longer-term glacier contribution to sea level rise (Hock and others, 2019). This is true even in Alaska, where the largest city, Anchorage, is critically dependent upon the meltwater of Eklutna Glacier in the western Chugach Mountains for both drinking water (~80% of the city's supply) and hydropower generation (10–15% of the city's supply; Moran and Galloway, 2006). The regional area-average glacier mass balance for all the western Chugach Mountains from 1962 to 2006 was −0.64 ± 0.07 m w.e. a⁻¹ (Berthier and others, 2010), and the geodetic mass balance of Eklutna Glacier calculated for the time periods 1957–2010 and 2010–15 was −0.52 ± 0.46 and −0.74 ± 0.10 m w.e. a⁻¹, respectively (Sass and others, 2017a). These results indicate an acceleration in mass loss. Reconstructing the annual variations of Eklutna Glacier's historic mass balance and resulting glacier discharge can help anticipate and mitigate future impacts on water resources.

Here, we calibrate a temperature index model with observations from 2011–15 and then reconstruct multi-decadal mass-balance variations of Eklutna Glacier to quantify the impacts of glacier change on discharge patterns. A novel approach combining in situ mass-balance measurements and observed snowlines from satellite imagery in conjunction with a geodetic mass balance is applied to calibrate the model parameters and identify an ensemble of the 250 best-performing model parameter combinations. This ensemble is then used to reconstruct the glacier's mass balance and discharge over the 35 mass-balance (hydrological) years 1985–2019 (e.g. hydrological year 1985 spans from 1 October 1984 to 30 September 1985).

2. Study site

Ranging in elevation from 600 to 1700 m a.s.l., Eklutna Glacier (61.21°N, 148.98°W) is a ~29 km² glacier located in the Western Chugach Mountains of south-central Alaska (Fig. 1). The glacier's tongue is fed by two branches: a larger main (east) branch (~16 km²) and a west branch (~13 km²). Ground-penetrating radar performed in 2010 indicated a mean ice thickness of 139 m with a maximum thickness of 430 m in the upper basin of the main branch (Sass, 2011). The glacier is in a transitional maritime climate. An automated weather station (AWS) in coastal Girdwood (76 m a.s.l.), 25 km south of the glacier, shows a mean annual precipitation of 1907 mm from 1984 to 2014 (NOAA Cooperative Station #500243, www.wrcc.dri.edu). Measurements from an AWS operated close to the equilibrium line altitude (ELA) on the glacier from 2009 to 2015 indicate a mean melt season (May–September) air temperature of 2.2°C.

Fig. 1. Map of Eklutna Glacier, its watershed and observation sites. The glacier area (year 2019) is shown in light grey, while the earlier 1985 extent consistent with legend is blue. Black dots represent stakes with the number of measurements over 2011–2015 calibration period. Elevation (2010) contours on glacier surface are at 100 m intervals. Black border depicts the model domain watershed. Tick marks represent easting/northing (UTM 6N, WGS84). Upper left inset map shows the location relative to Anchorage, Girdwood, Wolverine Glacier, Eklutna Glacier, and the surrounding Eklutna watershed boundary.

Eklutna Glacier has been subject to numerous research projects (Brabets, 1993; Larquier, 2011; Sass, 2011; Sass and others, 2017a). A field study conducted by the United States Geological Survey (USGS) for the mass-balance years 1984/85 to 1987/88 found that the annual glacier-wide mass-balance rate varied from −0.34 to 0.42 m w.e. a⁻¹ (mean: −0.1 m w.e. a⁻¹) and the accumulation-area ratio (AAR) ranged from 0.58 to 0.75 (mean: 0.63) (Brabets, 1993). For the hydrological years 1986–88, annual specific discharge (i.e. discharge per unit area of watershed area) at a gauging station (USGS 15277800 WF Eklutna) 3 km downstream from the glacier terminus ranged from 1.39 to 1.47 m (mean: 1.43 m). Alaska Pacific University has monitored both the glacier's mass balance and discharge since 2008 (Larquier, 2011; Sass, 2011). Sass and others (2017b) found winter mass balance varied from 1.4 to 2.5 m w.e. a⁻¹ (mean: 1.7 m w.e. a⁻¹) and summer mass balance from −1.4 to −2.1 m w.e. a⁻¹ (mean: −1.7 m w.e. a⁻¹) over 2008–15.

3. Data sources and use

To model Eklutna Glacier's mass balance and discharge, a digital elevation model (DEM) of the glacier and DEM-derived slope, aspect and topographic shading, as well as climate datasets, were required. The model was forced with daily mean air temperature and precipitation observations from AWS. For calibration and validation, we used in situ surface mass-balance observations (point balances), a geodetic mass balance and snowline observations digitized from optical satellite imagery (Sass and others, 2017a).

3.1 DEM, terrain layers and glacier boundaries

A LIDAR mission flown September 2010 with an Optech Gemini Airborne scanning system created a point-cloud with a nominal post spacing of 1.9 m and a vertical accuracy of 0.3 m (Sass and others, 2017a). The USGS processed the point-cloud resulting in a 2.5 m gridded DEM (Sass and others, 2017a). The DEM was used to derive surface slope, aspect, a hillshade layer, daily means of potential direct solar radiation and to delineate a watershed boundary (ArcGIS, vers. 10.6, ESRI, 2019). The hillshade layer was used to support the digitization of the 2010 glacier boundary. All layers were resampled (cubic convolution method) to a 25 m resolution for use as model input for computational efficiency while allowing spatial variations in mass balance to be captured. Glacier boundaries were manually digitized from Landsat satellite imagery for years 1985, 1993, 1999, 2006 and 2013 (Table S1). Specific years were chosen to establish a ~7-year interval over the simulation period allowing us to account for the retreat of the glacier in the simulations. Between 1985 and 2013, the glacier lost ~7% of its initial area.

3.2 Climate and weather data

Our model relies on two sources for temperature and precipitation: data from an AWS operated on the Eklutna glacier only during melt seasons (May–September) over 2011–19 and data from the nearest long-term year-round weather stations during 1984–2019. Transfer functions allowed the extension of the on-glacier temperature series to the entire period 1984–2019 for year-round model forcing.

The Eklutna Glacier AWS near the ELA (~1390 m a.s.l.) recorded hourly air temperature at 2 m above surface (Campbell Scientific, model CS215 temperature sensor) and liquid precipitation (Campbell Scientific, model TE525 tipping bucket precipitation sensor). We refer to these data as the ‘Eklutna data’.

The nearest long-term and year-round meteorological observations are from Girdwood, Alaska ~25 km south of Eklutna Glacier. We use the data from two stations as neither have a complete record of precipitation and air temperature data over the simulation period. An AWS located at the base of Alyeska Ski Resort in Girdwood recorded daily mean air temperature (76 m a.s.l., NOAA Cooperative Station #500243, www.wrcc.dri.edu) for the period October 1984–2016. A National Resource Conservation Service (NRCS) SNOTEL site also in Girdwood on Mt. Alyeska (470 m a.s.l., SNOTEL #1103, www.wcc.nrcs.usda.gov) has daily mean air temperature from 2010 to 2019. This SNOTEL site also has daily precipitation for the period October 1984–2019. We refer to the air temperature and precipitation data from the two sources as the ‘Girdwood data’.

To extend the Eklutna data temperature series (2012–19) to the entire period 1984–2019, we compared the Eklutna data to the Girdwood data to build transfer functions based on the common melt seasons of each of the Girdwood stations. An air temperature sensor failure for Eklutna data prevented the inclusion of the 2011 data observations. Since neither the NOAA nor the SNOTEL site data cover the entire period 1984–2019, we first regressed the Eklutna data against the NOAA Cooperative Station data (n = 479 days) and applied this transfer function to the NOAA Cooperative Station data 1984–2015 (Fig. 2a). Next, we regressed the Eklutna temperature data against the SNOTEL data (Fig. 2b) and applied this transfer function to the SNOTEL data for the melt seasons 2016–19. The SNOTEL daily precipitation data were used as model input without any adjustments and assumed to refer to the location (elevation) of the Eklutna AWS. Instead, biases were accounted for through a constant precipitation correction factor that was derived by model calibration (Section 2.1).

Fig. 2. Daily mean air temperatures (°C) during the melt season (May–September) at the Eklutna AWS, T _AWS, versus synthetic Eklutna temperatures, T _Synthetic, derived from (a) the NOAA Cooperative Station (76 m a.s.l, 2012–2015) and (b) the SNOTEL site (470 m a.s.l., 2016–2019) in Girdwood. The synthetic data refer to temperatures after application of the transfer functions derived from regressing the measured Eklutna AWS temperatures with the NOAA and SNOTEL temperatures for the overlapping years. Transfer functions for each dataset are provided. The dashed line depicts the 1:1 line.

3.5 In situ point mass balances

Alaska Pacific University has maintained a glacier mass-balance monitoring program on both branches of Eklutna Glacier since 2008. Index sites and locations have evolved over the course of the program, with four to seven sites per year generally located near the glacier centerline and spanning the accumulation and ablation zones of each branch (Fig. 1). Point mass balances were derived at each site from measured ablation stake height, snow depth and snow pit densities (Sass and others, 2017a, 2017b). The sites were visited at least twice each year. In spring, snow depth was measured and snow density derived from snow pit measurements to calculate winter balances. Stakes were also installed in spring and measured again in fall to observe the summer balance. Seasonal mass balances thus refer to the floating (summer balance) and combined time system (winter balance, Cogley and others, 2011). For model calibration, we used observations from 2011–15, resulting in a total of 50 summer (n = 25) and winter point balances (n = 25). The point mass-balance data for the period 2011–2015 are given in Sass and others (2017a, 2017b).

3.6 River discharge

River discharge was monitored on the West Fork Eklutna River ~3 km downstream from the glacier terminus. The corresponding subwatershed is 46% glacierized. A non-vented submersible pressure transducer (Onset Hobo, U20) mounted to a boulder in the river recorded the stage (water level) at hourly intervals and was verified to gage datum by standard wire-weight gage on the bridge (Rantz and others, 1982). Water pressures were adjusted to reflect atmospheric pressure variations (Onset Hobo, U20L) measured concurrently at the gage site. Discharge was measured by mid-section velocity measurements using a mechanical current meter (Price AA) and sounding weight (34 kg) or wading per USGS standards (Rantz and others, 1982). We constructed a stage–discharge relationship from 23 observations over the years 2015–19 (Fig. 3); the relationship between stage and discharge is well established for discharge <25 m³ s⁻¹. An hourly discharge time series for the melt seasons were aggregated as daily averages to compare with model output. Historic discharge data from the USGS for 1985–1988 were also used for validation (Brabets, 1993; Site #15277800, https://waterdata.usgs.gov/nwis).

Fig. 3. Stage-discharge rating curve constructed using 23 measurements (diamonds) at the watershed's gauging station (Fig. 1) of the West Fork Eklutna River for 2015–19. Dashed lines indicate ±1 sd.

3.7 Snowline position

Transient snowlines were digitized from satellite imagery over the June to September period between 1985 and 2015 (Table S1). Cloud-free satellite scenes with clearly visible snowlines were selected from Landsat 4–8 (30–60 m spatial resolution) and SPOT (5 m) as available from the USGS Earth Explorer website (earthexplorer.usgs.gov). Within a GIS (ArcGIS Ver. 10.6, ESRI, 2019), snowlines were manually digitized for each satellite melt season period scene (n = 61 snowlines; Table S1; June: 23 scenes, July: 26 scenes, August: 12 scenes). For model calibration, we created a centerline profile along the main (east) branch (Fig. 1). We then determined the centerline distance from the high point of the centerline at the glacier head to the intersection of each snowline with the centerline. This distance is used to compare modeled results with observations (henceforth referred to as snowline position). This centerline method avoids potential border effects along glacier margins from shadows, avalanche input and wind deposition. Because of difficulty discriminating snow from firn surfaces, we delineated snowlines only below the firn line (Kienholz and others, 2017). Field observations and late season satellite imagery allowed identification of the approximate location of the firn line. We assume a ±60 m horizontal digitizing accuracy, reflecting the coarsest satellite resolution.

4. Methods

4.1 Mass-balance and discharge model

We used the open access Distributed Enhanced Temperature Index Model (DETIM, Hock, 1999; http://regine.github.io/meltmodel) to recreate a 35-year (1985–2019) record of surface mass balance and discharge from Eklutna Glacier. At each gridcell of the 25 m resolution DEM and at a daily time step, DETIM simulates snow accumulation, melt through a temperature index method and resulting discharge via a linear reservoir approach (Baker and others, 1982). Accumulation is computed from precipitation using a threshold temperature to discriminate between rain and snow. Precipitation is distributed across the glacier using a linear precipitation gradient relative to the site and elevation of the on-glacier Eklutna weather station. Melt, M (mm d⁻¹), is calculated from daily average air temperature, T (°C), as

(1)$$M = \left\{{\matrix{ \!\!{( {\,f_{\rm m} + r_{{\rm ice/snow}}\;R} ) \;T\;\colon \;\;\;\;T > 0} \cr {0\hskip4.5pt\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\colon\;\; \;\;T \le 0\;\;} \cr } } \right.\;\;\;$$

where f _m is a melt factor (mm d⁻¹°C⁻¹), r _ice/snow (mm m² W⁻¹ d⁻¹°C⁻¹) represents a radiation factor for snow/firn and ice surfaces, and R is the potential clear-sky direct solar radiation (W m⁻²).

We compute mean daily discharge from the entire watershed defined by the gauging station (i.e. including the glacier and non-glacierized areas; Fig. 1) using a linear reservoir approach (Baker and others, 1982). The method assumes that, at any time, t, the reservoir volume, V(t), is proportional to the reservoir discharge, Q(t)

(2)$$V( t ) = kQ( t ) , \;$$

where k is a reservoir-specific storage constant with the unit of time. Solving Eqn (2) for Q and accounting for volume changes through input to the reservoir from melt and rainwater (I), hourly discharge from the reservoir is given by:

(3)$$Q( {t_2} ) = Q( {t_1} ) {\rm e}^{{-}1/k} + I( {t_2} ) -I( {t_2} ) {\rm e}^{{-}1/k}.$$

Following Hock and Noetzli (1997), we use three linear reservoirs defined by surface type (firn, snow, ice) to account for the markedly different hydraulic properties and associated water transit velocities of these zones. The firn reservoir reflects the glacier area above a user-defined firn line which is kept constant. The snow reservoir includes any area (on or outside the glacier) with snow present at the surface but excluding the glacier firn area; the ice reservoir includes the bare ice and snow-free non-glacierized area. Ice and snow-free areas outside the glacier are treated as one single reservoir since flow velocities from ice surfaces and the generally rocky and vegetation-free non-glacierized terrain are assumed to be of similar magnitude, and we aimed to keep the number of model parameters to a minimum. On the glacier, the firn reservoir has the greatest k-value, in the order of hundreds of hours, while values for snow typically are tens of hours, and ice <30 h (Hock and others, 2005). K-values are obtained by calibration (Section 4.2).

For each time step, melt and rain over all catchment grid cells defining the firn, snow and ice reservoirs are summed up, and Eqn (3) is applied separately for each reservoir. Finally, each reservoir's total discharge is summed to obtain the total discharge at the location of the gauging station.

4.2 Model calibration

The mass-balance model was calibrated by optimizing six model parameters over the period 2011–15 including the temperature lapse rate (ϒ), melt factor (f _m), radiation factors for ice (r _ice) and snow (r _snow), precipitation correction factor (p _cor) and a precipitation gradient (p _grad). A grid search for the best-performing mass-balance parameter combinations was applied by running the model with all parameter combinations inside a prescribed parameter space of defined increments established by preliminary testing. This resulted in a total of 257 040 model runs (Fig. 4). Each model run covered the period 1 October 2010 to 30 September 2015 (i.e. the five mass-balance years 2011–15). This period is closely aligned with the period of the geodetic mass balance reported by Sass and others (2017a, 2017b) (16 September 2010–24 August 2015).

Fig. 4. Parameter values of the 250 best-performing parameter combinations superimposed on the search parameter space. ϒ is the temperature lapse rate, f _m is the melt factor, r _snow and r _ice are the radiation factors of ice and snow (Eqn (1)), p _grad is the precipitation gradient, and p _cor is the precipitation correction factor. Numbers on the left and right sides indicate the range of the parameter space and black dots mark the parameter values tested. Grey lines connect the parameter combinations of each of the 250 best-performing parameter sets. The values above each point reflect the number of successful combinations through a tested parameter (only shown if >0).

To select the best-performing parameter combinations, we first selected all parameter sets that calculated a geodetic mass balance rate that matched (within uncertainties) the 2010–2015 geodetic mass balance (−0.74 ± 0.10 m w.e. a⁻¹) found by Sass and others (2017a, 2017b). This reduced the number of potential parameter combinations to 25 062 from our original total of 257 040. Next, we performed a multi-criteria calibration by comparing the results of each remaining DETIM parameter set with both summer and winter point balance measurements and snowline positions. For each stake location, the point balance was extracted from the continuous model runs for the exact same period the measurement refers to, thus allowing direct comparability of the 50 point balances. We calibrated the model to the combined dataset of summer and winter point balances to capture the effects that melt parameters have on winter balance, for example, due to winter melt events, and precipitation parameters on summer balance, for example, due to summer snowfall events.

A total of 20 snowline positions obtained from satellite imagery were available during the calibration period (Table S1). For both variables, we quantified the agreement between modeled and observed values by calculating standard z-scores as:

(4)$$Z_i = \displaystyle{{{\rm RMSE}_i-\overline {{\rm RMSE}} } \over {{\rm sd}}},$$

where RMSE_i is the root mean square error from the ith parameter set, $\overline {\rm RMSE}$ is the mean RMSE averaged across all remaining 25 062 parameter sets, and sd is the std dev of RMSE. Hence, the unitless z-score represents each parameter set's departure from the mean performance in std dev. for each variable, i.e. the point balance at index sites (in m w.e.) and distance along the glacier centerline (in m) for the snowline positions. A z-score of zero indicates an ‘average’ error; a negative z-score indicates a parameter set with less error than average, and vice versa.The z-scores of better performing parameter sets (i.e., negative z-scores) were then normalized to range between zero and one, both to allow equal weighting between variables of different units and so that larger positive values represented better fitting parameter sets than those closer to zero.

We found that a total of 7051 parameter sets (~28% of those reproducing the observed geodetic balance by Sass and others (2017a, 2017b)) performed better than average for both variables. The mean z-score (±1 sd) for these 7051 parameters sets of stakes and snowline location were 0.47 ± 0.25 and 0.23 ± 0.16, respectively. For each variable, the normalized z-score maximum value of one indicates the best possible agreement between modeled and observed. Since the observations have errors and the model is overparameterized, it is not possible to determine a single best model run. Different parameter combinations can perform equally well in reproducing the observations. For example, overestimation of melt due to an overestimated degree-day factor can be compensated by an underestimated precipitation correction factor. Therefore, following previous studies (Trüssel and others, 2015; Kienholz and others, 2017), we use an ensemble of the best-performing parameter sets for further analysis and present results in terms of ensemble mean and std dev. We chose the best 250 parameter sets, thus considerably more than the ensemble of 15 and 40 best parameter sets used by Trüssel and others (2015) and Kienholz and others (2017), respectively.

The 250 parameter sets were chosen so that the same minimum z-score threshold was exceeded for both variables. Figure 5 demonstrates how some parameter sets perform well for point balances (high z-score) and poorly for snowlines (low z-score), and vice versa. Thus, we only retain parameter sets that yield high z-scores for both point balances and snowlines. A threshold of 0.5 yielded 250 parameter sets (Fig. 5).

Fig. 5. Normalized z-scores for ablation stakes versus normalized z-scores of snowline positions for 7051 parameter sets. Values >0.5 for both variables (marked in red) corresponded to the 250 best-performing parameter sets. Grey lines mark the threshold of 0.5.

Figure 4 depicts the frequency of the top 250 best-performing parameter sets for the tested parameter space of each parameter. The temperature lapse rate among the best parameter sets ranged from −0.6 to −0.2°C (100 m)⁻¹ with a mean of −0.3°C (100 m)⁻¹ and a mode of −0.2°C (100 m)⁻¹. All melt factors ranged from 3.75 to 6.00 mm d⁻¹°C⁻¹ with the most frequent values being 5.75 and 6.00 mm d⁻¹°C⁻¹. The radiation factor for ice was distributed across three values between 0.0242 and 0.414 m² W⁻¹ mm d⁻¹°C⁻¹ while the radiation factor for snow was most frequent in the range of 0.005 and 0.0098 m² W⁻¹ mm d⁻¹°C⁻¹. The precipitation gradient of 15% (100 m)⁻¹ and a precipitation correction factor of 0 and 5% dominated among the best parameter sets.

In a final step, we calibrated the three storage constants k (Eqn (2)) using the 250 best-performing sets of mass-balance parameters. The coefficient does not affect the annual amounts of discharged water but modifies the seasonality of discharge by increased (higher k-values) or decreased (lower k-values) flow speed of water through the firn, snow and ice reservoirs. The k-values were determined by varying them within established ranges (Hock and others, 2005). Tested values for k _firn, k _snow and k _ice ranged from 240 to 400 h (using 20 h intervals), 30–200 h (10 h intervals) and 5–25 h (1 h interval), respectively. The k-values were calibrated over 1 September 2010 to 30 August 2015 (3402 daily discharge values). A Nash–Sutcliffe model efficiency coefficient (R ², Nash and Sutcliffe, 1970) was calculated between modeled and observed daily discharge to assess model performance. R ² values are typically used to assess the efficiency of hydrological model results (Krause and others, 2005). Values can be negative, thus differing from the coefficient of determination r ². These best-performing k-values (k _firn = 300 h, k _snow = 70 h and k _ice = 15 h) were then used for all model runs.

4.3 Model validation

We cross-validated the 250 best-performing model parameter sets using (i) snowline positions, and (ii) discharge observations over periods excluded from calibration. Modeled snowline locations were compared to observations on 41 melt season days between 1985 and 2010 (Figs 6, 7). There is a tendency for the model to over-predict snow cover extent early in the season and to under-predict at the end of the season. This discrepancy may at least partially be caused by the darkening of the snow surface over the melt season, for example, through particulate and biota accumulation which is not accounted for in the model (Skiles and others, 2018). Additionally, daily mean modeled discharge was compared to observed discharge over the periods 1985–88 (Brabets, 1993) and 2016–19 (Fig. 7). Nash-Sutcliffe efficiency values R ² for the hydrological years ranged from 0.66 to 0.85 (mean of best-performing 250 parameter sets) with a value of 0.77 averaged over all 8 years.

Fig. 6. Modeled and observed snowline locations for 41 days during the melt seasons 1985–2010 for the best-performing parameter combination (ϒ = −0.2 °C (100 m)⁻¹, M _f = 5.5 mm°C⁻¹ d^−1, r _ice = 0.0414 m² W⁻¹ mm d⁻¹(°C)⁻¹, r _snow = 0.0098 m² W⁻¹ mm d⁻¹ (°C)⁻¹, p _cor = 15% and p _grad = 25% (100 m)⁻¹). Modeled snow-covered glacier area is shown in grey, observed snowlines are depicted by blue lines, and centerline is white line. Lower right plot depicts the 41 observed (x-axis) versus modeled (y-axis) snowline positions as measured along the centerline profile (in units of 1 000 m) including the 1:1 line (grey).

Fig. 7. Daily mean discharge during the melt seasons (May–September) for 1985–88 and 2016–19 calculated from the 250 best-performing parameter sets. The grey shading represents the range of modeled discharges with the black line representing observed discharge. Nash–Sutcliffe efficiency values (R ²) and the ratio of the modeled and observed discharge expressed in percent are provided. The blue bars reflect unaltered daily precipitation (SNOTEL #1103, www.wcc.nrcs.usda.gov). Ticks mark the first day of each month.

5. Results

5.1 Mass balance

Using the 250 best-performing mass-balance parameter sets and forcing the model with the synthetic meteorological time series derived from the observations from the Girdwood stations, we hindcast surface mass balance for the mass-balance years from 1985 to 2019 (Fig. 8). A fixed date system was applied for computing winter and summer balances with 12 May marking the end of winter and 17 September the end of summer. These dates were determined from the average annual maxima and minima of modeled cumulative mass balances (daily resolution) derived from preliminary model runs.

Fig. 8. Modeled winter, summer and annual surface mass balance (m w.e.) for the mass-balance years 1985–2019. Results refer to mean values from the 250 best-performing parameter sets. Vertical black lines show ±1 sd (only winter and summer balances).

The mean annual surface mass balance averaged over the 250 parameter sets was −0.62 ± 0.06 m w.e. (±1 sd). The maximum annual balance (0.83 m w.e.) occurred in 1988 and the minimum (−2.3 m w.e.) in 2004 (Fig. 8). Throughout the time series, results suggest the annual surface mass balances were mostly negative. From 1985 to 2019, there is a statistically significant negative trend in annual mass balance (−0.31 m w.e. decade⁻¹, p = 0.02; Fig. 9), due mostly to a significant negative trend in summer balance (−0.26 m w.e.  decade⁻¹, p < 0.01). No trend in winter balance was found (−0.05 m w.e. decade⁻¹, p = 0.45).

Fig. 9. Modeled winter (blue), summer (red) and annual (black) mass balance (m w.e.) for the mass-balance years 1985–2019. Dots show the mean of the 250 best-performing parameter sets (±1 sd) and lines show the linear trends.

Modeled glacier melt increased steadily over the near-decadal periods from 1985 to 2019 (Fig. 10). Changes in mean annual cumulative melt (hydrological year) over four consecutive 8–9-year periods (1985–93, 1994–2001, 2002–10 and 2011–19) demonstrated a 17% increase from the earliest to the latest period. In addition, the day of the year when 5 and 95% of melt have occurred is 7 days earlier in spring and 5 days later in fall in the latest period than in the earliest period, demonstrating a prolongation of the melt season. Like the near-decadal trends in annual melt, we found considerable changes over time in transient and annual AAR, computed as the ratio of snow-covered and total glacier area (Fig. 11). Annual AARs (defined by the annual AAR minima for the year) decreased from ~65% during the earliest period (1985–93) to ~45% during the latest period (2011–19). Transient AARs indicate that the glacier's winter snow cover extent is depleted faster in each consecutive period.

Fig. 10. Modeled cumulative glacier melt between 25 April and 30 September averaged over four consecutive periods from 1985 to 2019. Cumulative melt for each year is relative to the start of the mass-balance year, i.e. 1 October of the previous calendar year. Lines give ensemble means and shading indicates ±1 sd for the 250 best-performing parameter sets. Black dots indicate when 5 and 95% of annual cumulative melt is reached. Ticks mark the first day of each month.

Fig. 11. Modeled transient accumulation-area ratio, AAR (%) averaged over four consecutive periods between 1985 and 2019. Lines show the ensemble mean of the 250 best-performing parameter sets and shading indicates ±1 sd for the 250 best-performing parameter sets. Black dots depict the date and value of the annual AAR at the time of its minimum. Ticks mark the first day of each month.

5.2 Discharge

The modeled daily mean discharge over the melt season (May–September) increased by 17% between the 1985–1993 and the 2011–2019 interval (Figs 12, 13). The mean specific discharge averaged over the melt season period from 1985 to 2019 was 2.1 ± 0.11 m (±1 sd) with the greatest modeled specific discharge most recently (2.7 ± 0.14 m in 2019) and the lowest in 1985 (1.5 ± 0.07 m). Four of the five highest discharge years for the model period occurred after 2000 (2004, 2005, 2016 and 2019). Despite strong day-to-day variability due to differences in amount and timing of daily precipitation and melt, we find a significant trend of increasing modeled specific discharge within the melt season (0.14 m per decade, p < 0.01). Modeled melt season mean discharge represented 88 ± 3% of annual discharge over the 1985–2019 time period. The strongest increase in discharge occurred in the main melt season and fall. Summer precipitation showed no trend (p = 0.83) while temperature had an increasing trend (0.23°C per decade, p < 0.01; Fig. 12). This reinforces the notion that increasing melt dominated the discharge trend.

Fig. 12. Mean specific discharge (a), mean air temperature (b), and total precipitation (c) during the melt season (May–September) from 1985 to 2019. Specific discharge refers to the mean of the 250 best-performing parameter sets (±1 sd). Air temperature and precipitation data reflect the data used to force the model (Section 1.5), i.e. air temperature depicts the record for the AWS site on the Eklutna Glacier extended to the entire period based on transfer functions with nearby weather stations, and precipitation refers to the nearby weather station record prior to applying the calibrated precipitation correction factor (Section 3.2). Lines show the linear trends. Vertical dotted lines mark the four averaging periods used in Figures 10, 11, 13.

Fig. 13. Modeled daily mean discharge (m³ s⁻¹) for the 250 best-performing parameter sets during the melt seasons (May–September) of the period 1985–2019. Four intervals are depicted: 1985–93, 1994–2001, 2002–10, 2011–19. Mean melt season precipitation for each interval is given in the legend. Shading indicates ±1 sd for the 250 best-performing parameter sets. Peaks in late fall reflect large precipitation driven flood events on 21 September 1995 and 3 October 2003. Ticks mark the first day of each month.

The entire time series was used to calculate standard hydrological metrics following Fleming and Clarke (2005), including annual median daily discharge, total annual discharge, annual maximum daily discharge and the centroid of annual hydrograph (half the year's total flow volume). Both the median daily discharge evaluated over the melt season and total annual discharge have positive significant correlations (median discharge = 0.319 × year – 57.9; p = 0.02: total annual flow = 11.8 × year – 112.77; p < 0.01). The annual maximum daily discharge and the centroid of the annual hydrograph both showed a positive trend, though not significantly correlated (max daily discharge = 0.109 × year – 191.4; p = 0.44; centroid day of year (DOY) = 0.046 × year + 114.92; p = 0.54).

6. Discussion

6.1 Model calibration

Results from model calibration indicate the value of using multi-criteria validation that includes the use of a geodetic mass-balance constraint, point balances and snowline positions. Even after a geodetic constraint, the point balances alone were not well-constrained, as the comparison between modeled and predicted point balances found 95% of individual parameter set model runs had an r ² > 0.90 (n = 25 062). This suggests that selecting the best model parameters based solely on point balances is not sufficient. Point balances have frequently been used for melt model calibration (Hock, 1999; Schuler and others, 2005), while the use of snowlines is becoming more common (Evans and others, 2008; Ziemen and others, 2016; Kienholz and others, 2017; Barandun and others, 2018). Using only snowline observations would result in selecting different model parameter sets; however, this is expected to lead to decreased agreement between modeled and measured point mass balances as evident from Figure 5. Barandun and others (2018) found the inclusion of snowlines within mass-balance modeling significantly narrowed uncertainty. Calculating a normalized z-score for each variable (point balances and snowlines) allowed equal weighting between the two variables of different units. Setting the threshold value of 0.5 allowed selecting an ensemble of 250 parameters where agreement was good for both point balances and snowlines (Fig. 4). The inclusion of both variables in addition to the geodetic balance proved necessary to provide a more robust calibration.

Model parameter values vary widely within the set of 250 best-performing parameters; however, the range in the modeled variables of the ensemble is relatively narrow for most variables (Figs 9–12). This provides confidence that the modeled trends are robust. Obviously, individual model parameter values combine to yield similar responses.

6.2 Mass balance

Eklutna Glacier mass-balance results were compared to the USGS benchmark glaciers Wolverine and Gulkana, which have been studied since the 1960s (O'Neel and others, 2019). Wolverine Glacier (15.6 km²) is located within a maritime climate (Bieniek and others, 2012) on the Kenai Mountains, ~90 km south of Eklutna Glacier. Gulkana Glacier (16 km²) is in a continental climate in the eastern Alaska Range, ~290 km northeast of Eklutna. The ensemble mean of the modeled annual mass balance for Eklutna Glacier is strongly correlated with the observed mass balance of Wolverine Glacier from 1985 to 2019 (r ² = 0.73, p < 0.01, n = 35 years) and moderately correlated to the observed mass balance of Gulkana Glacier (r ² = 0.33; p < 0.01, n = 35 years; Figs 14, 15). Some of the scatter can be attributed to the different methodologies in determining the balances including different time systems. Nevertheless, the significant correlations indicate spatially coherent climatic drivers of annual mass change. The lower correlation with Gulkana is expected given its distant location and different climatic setting. Eklutna and Wolverine glaciers are 90 km apart and each are >35 km from the Gulf of Alaska. Eklutna Glacier's annual mass balances are on average slightly higher than Wolverine. This is possibly due to Eklutna's northerly versus Wolverine's southerly orientation, and the generally cooler, somewhat transitional climate of Eklutna compared to the more maritime climate of Wolverine Glacier. While both are near the ocean, Eklutna is on the lee side of the Chugach Mountains and Wolverine on the coastal side of the Kenai Mountains.

Fig. 14. Modeled annual mass balance (m w.e.) for Eklutna and observed balances for Wolverine and Gulkana glaciers (O'Neel and others, 2019) over the period 1985–2019. Balances for Eklutna 2008–15 based on the glaciological method are also shown (Sass and others, 2017a) correlating well with the modeled annual balances (r ² = 0.89, p < 0.01, bias = −0.05 m w.e.).

Fig. 15. Modeled annual mass balance (m w.e.) for Eklutna Glacier versus reported balances from (a) Wolverine Glacier and (b) Gulkana Glacier (O'Neel and others, 2019) over mass-balance years 1985–2019. The grey dashed line depicts the 1:1 line.

Disagreements between modeled Eklutna and Wolverine Glacier mass-balance results in 1990 and 2009 coincide with Redoubt volcanic eruptions. Redoubt Volcano (~220 km away) erupted throughout the winter of 1989/1990 and spring of 2009 and deposited substantial amounts of volcanic ash on glacier surfaces across parts of Alaska (Waythomas and Nye, 2002; Schaefer and Wallace, 2012). The largest deposit occurred in 1992 from Spurr Volcano with an estimated 500–1 000 g m⁻² deposition on Eklutna Glacier and no deposition on Wolverine Glacier (McGimsey and others, 2001). Such ash deposits may influence surface albedo across multiple years. DETIM as implemented here did not consider these altered surface conditions.

6.3 Discharge

The modeled increase in annual discharge over the simulation period is consistent with the negative trend in summer mass balances coincident with increasing temperatures (Fig. 12). Since summer precipitation shows no significant trend, we attribute the increase in melt season discharge to the increasingly negative glacier mass balances, and thus to the additional water released from the glacier. Changes in seasonal snow storage or precipitation over the non-glacierized catchment area may also affect the discharge trend; however, model results do not support an increase in snow accumulation outside the glacier that could have increased melt season discharge.

Results also show a lengthening of melt season (Figs 10, 13). A positive feedback may be occurring where more melt causes larger bare ice areas, exposing lower albedo ice and leading to more melt. More exposed glacier ice surface also leads to a faster throughflow of water (Hock and others, 2005). Hence, one can expect increases in flood peaks, especially when coincident with heavy precipitation events.

Our model results are consistent with the observations and also with other glacier discharge studies in Alaska (O'Neel and others, 2014; Beamer and others, 2016). Beamer and others (2016) found glacier volume loss within the Gulf of Alaska watershed contributed 760 km³ a⁻¹ (8%) to mean annual discharge over the period 1980–2014. O'Neel and others (2014) found Wolverine Glacier's mass loss over the 1966–2011 period caused a 23% increase in discharge. Annual glacier discharge typically first increases as mass balances become more negative but then tends to decrease as the glacier becomes smaller (Huss and Hock, 2018).

The positive trend in annual discharge indicates that peak water, i.e. the year when annual discharge reaches a maximum due to glacier retreat, has not been reached within the Eklutna watershed. This is consistent with the model and observational studies elsewhere in Alaska (see Hock and others (2019) for summary). Valentin and others (2018) project that in Eastern Alaska's Copper River watershed (60 800 km²), glacier discharge may peak in 2070 using the moderate emission scenario RCP 4.5. Global-scale projections indicate that peak water in Alaska will be reached later in the 21st century than elsewhere in the world (Huss and Hock, 2018).

The anticipated peak water has implications for Eklutna watershed water resources. The negative mass balance currently generates water for power generation and fresh water supply above precipitation inputs; however, this surplus will decrease as glacier volume further diminishes. The earlier spring melt will provide additional water in the near term, but will not be sustainable after peak flow, as melt depletes the glacier ice reservoir. A recent study of headwater glaciers of the Columbia River found that meltwater contribution there has already reached peak flow (Moore and others, 2020). A future study could apply the model developed here to forecast when peak flow is reached in the Eklutna watershed under various climate change scenarios.

7. Conclusion

We calibrated the mass-balance and discharge model DETIM using a combination of geodetic and point mass-balance observations as well as snowline data. Results indicate that a multi-criteria optimization including diverse types of observations is necessary to constrain model parameters. The calibrated model allowed us to reconstruct the mass balance and glacier discharge history of Eklutna Glacier over the 35-year period 1985–2019.

Despite a wide range of model parameter values within the set of 250 best-performing model parameters, the range of the ensemble's modeled mass balance and discharge was relatively small compared to the modeled changes over time. Significant positive trends in modeled discharge are consistent with positive trends in air temperature, modeled melt season water production and prolongation of the melt season, as well as modeled negative trends in summer mass balance and AAR. The lack of significant trends in winter balance and summer precipitation suggests that the increased discharge from the highly glacierized catchment is driven by the loss of glacier mass. The modeled increase in melt season discharge throughout the simulation period indicates that peak water has not been reached. The negative trend in annual mass balance and the associated increase in discharge show that warming summers ‘mined’ the glacier’s stored water. This will have implications for water managers who seek to maximize water resources for hydropower needs as glaciers retreat and thin under a warmer climate (Hock and Jansson, 2005). Continued monitoring of both Eklutna Glacier mass balance and discharge as model input and validation will better inform predictive models of future glacier mass and discharge as well as peak water timing.

While we found fairly good agreement between observed and predicted discharge, improvements could include the effect of variable albedos to capture the volcanic tephra deposited on Eklutna Glacier during several eruptions during the simulation period (McGimsey and others, 2001; Schaefer and Wallace, 2012). Such efforts will require documenting the extent and amounts of the tephra and their impact on melt.