Glacier mass and area changes on the Kenai Peninsula, Alaska, 1986–2016

Abstract

Glacier mass loss in Alaska has implications for global sea level rise, fresh water input into the Gulf of Alaska and terrestrial fresh water resources. We map all glaciers (>4000 km²) on the Kenai Peninsula, south central Alaska, for the years 1986, 1995, 2005 and 2016, using satellite images. Changes in surface elevation and volume are determined by differencing a digital elevation model (DEM) derived from Advanced Spaceborne Thermal Emission and Reflection Radiometer stereo images in 2005 from the Interferometric Synthetic Aperture Radar DEM of 2014. The glacier area shrunk by 543 ± 123 km² (12 ± 3%) between 1986 and 2016. The region-wide mass-balance rate between 2005 and 2014 was −0.94 ± 0.12 m w.e. a⁻¹ (−3.84 ± 0.50 Gt a⁻¹), which is almost twice as negative than found for earlier periods in previous studies indicating an acceleration in glacier mass loss in this region. Area-averaged mass changes were most negative for lake-terminating glaciers (−1.37 ± 0.13 m w.e. a⁻¹), followed by land-terminating glaciers (−1.02 ± 0.13 m w.e. a⁻¹) and tidewater glaciers (−0.45 ± 0.14 m w.e. a⁻¹). Unambiguous attribution of the observed acceleration in mass loss over the last decades is hampered by the scarcity of observational data, especially at high elevation, and by large interannual variability.

1. Introduction

Alaskan glaciers account for ~12% of the total global glacierized area excluding the Greenland and Antarctica ice sheets (Pfeffer and others, 2014; RGI Consortium, 2017) and are an important contributor to global sea level rise (Dyurgerov and Meier, 1997; Arendt and others, 2002, 2013; Meier and Dyurgerov, 2002; Berthier and others, 2010; Arendt, 2011; Gardner and others, 2013; Larsen and others, 2015; Zemp and others, 2019). Gardner and others (2013) found a mass change of −50 ± 17 Gt a⁻¹ for the period 2003–2009 for all glaciers in Alaska and adjacent Yukon Territory. Subsequent estimates based on gravimetric methods (Wouters and others, 2019), glaciological and geodetic measurements (Zemp and others, 2019) and a combination of these methods (Box and others, 2018) vary between −48 ± 9 and −73 ± 17 Gt a⁻¹ for the period 2006–2016. Johnson and others (2013) found highly variable mass change rates for the glaciers in the Glacier Bay region (~6400 km²) during the period 1995 and 2011 with no clear trend, while Das and others (2014) found a decrease in mass balance of all glaciers in the Wrangell Mountains (~5000 km²) from −0.07 ± 0.19 m water equivalent (w.e.) a⁻¹ in 1957–2000 to −0.24 ± 0.16 m w.e. a⁻¹ in 2000–2007 although uncertainty ranges overlap. Wastlhuber and others (2017) found accelerated elevation change for the glaciers (~700 km²) in the Susitna River basin (−0.14 ± 0.07 m a⁻¹ in 1951–2005 compared to −1.20 ± 0.75 m a⁻¹ in 2005–2010). Also, all five benchmark glaciers in Alaska (Gulkana, Wolverine, Lemon Creek, South Cascade and Sperry Glacier) have lost mass with average mass-balance rates ranging from −0.58 to −0.30 m w.e. a⁻¹ since mid-20^th century (O'Neel and others, 2014; O'Neel and others, 2019).

Several studies report on glacier changes on the Kenai Peninsula in south-central Alaska. An early study showed that glaciers in this region have experienced widespread recession since the Little Ice Age Maximum (late 1700s through late 1800s) (Wiles and Calkin, 1992). Most studies on Kenai glaciers have been carried out on parts of the whole glacierized region. Aðalgeirsdóttir and others (1998) studied the volume change of Harding Icefield (subregion II in Fig. 1) through comparing the United Stations Geological Survey (USGS) topography maps in the 1950s and airborne Lidar elevations along center-lines in 1994/1996, and found the Harding Icefield had lost ~34 km³ of its total volume, corresponding to an area averaged elevation change of −21 m over the ~43 years. Another study using Landsat images shows that the Harding Icefield and Grewingk-Yalik Glacier Complex (subregions I and II, Fig. 1) lost ~78 km² (~3.6%) of their glacier area between 1986 and 2002 (Hall and others, 2005). A comparison between USGS maps (1950) and SRTM DEMs also indicates that glaciers in the same region (subregions I and II) were thinning at a rate of 0.61 ± 0.12 m a⁻¹ for the period 1950–1999, and the thinning rate for Harding Icefield had further accelerated by a factor of 1.5 during the mid-1990s–1999 relative to 1950 to the mid-1990s (VanLooy and others, 2006).

Fig. 1. Map of studied glaciers on the Kenai Peninsula. Black dots show the ICES at footprint on stable terrain (i.e. terrain without glacier, water and vegetation), while red stars refer to weather stations (Table S1). The dashed polygons refer to four subregions analyzed separately in this study. The pie charts show the regional fraction covered by land-terminating, lake-terminating and tidewater glaciers.

This study aims to quantify recent area and mass changes of all glaciers on the Kenai Peninsula aggregated in four subregions defined by the major icefields of the peninsula. First, we compile four new glacier inventories based on multiple satellite images for the years 1986, 1995, 2005 and 2016, which are more comprehensive and detailed than those used in previous studies (Hall and others, 2005; VanLooy and others, 2006; Le Bris and others, 2011; Kienholz and others, 2015). Second, from these inventories we assess the glacier area changes. Finally, we determine the glacier volume and mass changes between 2005 and 2014 using DEMs derived from satellite imagery.

2. Study area

The Kenai Peninsula is located in south-central Alaska, between the Cook Inlet and the Gulf of Alaska (Fig. 1). The Kenai Mountains are an effective barrier for wet airflow (Le Bris and others, 2011), resulting in frequent cloud cover, and high precipitation on the windward south-eastern side of the mountains, while the leeward side lies in the rain shadow. Hence, the glaciers experience a predominantly maritime climate along the south and east side of the peninsula, while a more continental climate prevails in the northern and western parts. The annual total precipitation at a weather station located in the northeast (station A3, Fig. 1) is ~1800 mm, while stations A1 and A2 on the west coast of the Peninsula receive annual totals of ~600 and ~500 mm, respectively, averaged over the period 1986–2016 (https://gis.ncdc.noaa.gov/maps/ncei/summaries/monthly, Table S1).

There are ~1460 glaciers on the Kenai Peninsula, spanning from sea level to ~2000 m a.s.l. and covering 4165 km², which corresponds to ~5% of the total glacierized area in Alaska according to the Randolph Glacier Inventory (RGI 6.0; RGI Consortium, 2017). Of all glaciers, only 11 are tidewater and 14 are lake-terminating, but they cover 948 km² (22.7% of total area) and 968 km² (23.2%), respectively. The Sargent and Harding icefields are the largest ice masses. Especially the Harding icefield has been the focus of several glaciological studies (Aðalgeirsdóttir and others, 1998; Hall and others, 2005; VanLooy and others, 2006).

3. Data and methods

3.1 Glacier outlines extraction and attributes calculation

We used data from Landsat complemented by Sentinel-2 and ASTER images (Table S2), where Landsat data were not available or relevant parts of the Landsat scenes were too cloudy or snowy to extract glacier outlines with sufficient accuracy to compute area changes. A total of 12 orthorectified satellite scenes acquired at the end of the ablation season (04 August to 26 September) were used to delineate glacier outlines for the years 1986, 1995, 2005 and 2016 (Fig. S1). Only images with minimum cloud cover and seasonal snow were chosen. All images (Table S2) were retrieved from the USGS (https://glovis.usgs.gov).

The glacier outlines for clean ice were delineated using a semi-automatic procedure based on the band ratio segmentation method (Paul, 2000; Albert, 2002; Guo and others, 2015; Li and others, 2017). A red/short-wave infrared (R/SWIR) band ratio with a threshold of 2–2.5 (TM3/TM5 of Landsat TM and ETM+, TM4/TM6 of Landsat OLI imagery) and an additional threshold on the blue band (band 1 of TM/ETM+ and band 2 of OLI) digital number were adopted and thresholds were selected interactively (Paul and Kääb, 2005; Raup and others, 2007). We employed sectionalized glacier outline delineation by using moving windows with different optimal band ratio thresholds on the same satellite image to minimize the influence of various snow/ice ablation levels on the optimal thresholds (Guo and others, 2015). We also applied a median filter (3 by 3 kernel) to reduce noise in shadowed regions and remove isolated pixels outside the glaciers, although it might reduce the size of small glaciers to some extent (Paul, 2002; Andreassen and others, 2008). Only glaciers with an area ≥0.01 km² were mapped. Following previous studies (Hall and others, 1992; Racoviteanu and others, 2008a; Pan and others, 2012; Burns and Nolin, 2014), debris-covered glacier margins were delineated manually from visual inspection of the images. Clean ice outlines were visually inspected using all available contrast-enhanced false-color composite images as well as Google Earth^TM images to correct errors caused by cast shadow, light clouds and seasonal snow.

Following the top-down method (Guo and others, 2015), which identifies the ice divides by aspect differences between the two sides of mountain ridges, we derived ice divides by using the IFSAR DEM, within a 1500 m buffer distance from all glaciers. Results were scrutinized with the help of the shaded relief data and the Landsat scenes, while gross errors were corrected manually where necessary. The ice divides were used to split the glacier complex into individual glaciers.

We extracted a series of glacier-specific attributes including Global Land Ice Measurements from Space (GLIMS)-IDs, coordinates, surface area, debris-covered area, mean slope and aspect, maximum, mean and minimum elevation using the IFSAR DEM and glacier outlines, and added glacier names (if existing) and glacier terminate types (land-terminating, lake-terminating and tidewater) from RGI 6.0. Mean aspect was derived from the arctangent of the mean values of the sine and cosine of the aspect of all DEM cells within a glacier.

3.2 DEM generation and mass change derivation

DEMs derived from ASTER L1A stereo images in 2005 (Fig. 1) and the IFSAR DEM in 2014, were used to compute surface elevation and glacier volume changes over the period 2005–2014. The ASTER data were acquired between 9 August and 14 August, while the IFSAR data were acquired between 29 August and 12 September 2014 (http://viewer.nationalmap.gov/basic/ access date: December 2018, Table S3).

3.2.1 IFSAR DEM and validation

As part of the Alaska Statewide Digital Mapping Initiative, which started in 2010, the IFSAR DEM was acquired from airborne radar operating in the X- and P- bands, with a native resolution of 5 m mosaics. In this study, we used the elevation data product, IFSAR Digital Surface Model (DSM), which is based only on the X-band, since the X-band signal generally penetrates the glacier surface less deep than the P-band (Bert and others, 2011). X-band radar (9.7 GHz) has been shown to penetrate snow and ice (e.g. Rignot and others, 2001; Gardelle and others, 2012; Gusmeroli and others, 2013) depending on the grain size and surface wetness (Mätzler and Schanda, 1984), as well as frequency (Davis and Poznyak, 1993; McNabb and others, 2019). X-band penetration depths are not well known, but in cold and dry snow penetration depths of several meters have been found in Antarctica (Davis and Poznyak, 1993; Surdyk, 2002) and at high elevations in the European Alps (Millan and others, 2015; Dehecq and others, 2016). Some geodetic mass-balance studies have therefore applied multi-meter corrections to account for radar penetration (e.g. Melkonian and others, 2014; Round and others, 2017; Lambrecht and others, 2018). However, over wet glacier surfaces radar penetration is limited to a few centimeters below the surface (Mätzler and Schanda, 1984).

To estimate surface wetness, we analyze air temperature records from three weather stations at high elevation (990–1420 m a.s.l.) on or close to the ice cap during and preceding the period of IFSAR acquisition (Fig. S2). Daily mean and minimum air temperatures were generally well above 0°C leading us to conclude that meltwater was present at or immediately below the surface during this period. Hence, we assume radar penetration to be negligible in this study.

Elevation data derived from Ice, Cloud and Land Elevation Satellite Geoscience Laser Altimeter System (ICESat/GlAS) laser altimeter were used to validate the IFSAR DEM. We used the data from GLAH14 (release 634), available from the National Snow and Ice Data Center (NSIDC). The data has a 10.6 ± 4.5 m horizontal and 0.34 m vertical accuracy (Magruder and others, 2007), and is considered to be the most consistent globally available elevation data (Zwally and others, 2014). To minimize the influence of land surface changes, we only compared the IFSAR and GLAH14 data over stable terrain, i.e. we excluded all glacier, lake, ocean and vegetation pixels as extracted from the land cover map of Alaska (USGS, 2011, https://gapanalysis.usgs.gov/gaplandcover/data/download/).

The Ellipsoid of the ICESat data was first transformed from TOPEX/Poseidon Ellipsoid to WGS84 Ellipsoid and then to NAD83, while the vertical datum was transferred from WGS84 to NAVD88 (same as IFSAR DEM) through an online tool (https://vdatum.noaa.gov/vdatumweb/). For each ICESat footprint, the corresponding IFSAR elevation was extracted by calculating the area-weighted average of the elevations of all pixels fully or partially enclosed by a circle with 35 m radius around the footprint's center (Brun and others, 2017). ICESat/GLAS data with slope larger than 5° were omitted (VanLooy and others, 2006; Le Bris and Paul, 2015), and outliers were removed based on the ‘Interquartile range’(IQR), where only data between Q ₁  − 1.5 ⋅ IQR and Q ₃  + 1.5 ⋅ IQR were used (Q ₁ is the middle value in the first half of the rank-ordered dataset; Q ₃ is the middle value in the second half of the rank-ordered dataset), and the rest discarded as outliers (McGill and others, 1978). The mean and Std dev. of the elevation differences between ICESat data and IFSAR DEM is 0.53 ± 1.38 m, indicating reasonable agreement between the two data sets. For comparison differences between the IFSAR and ICESat/GLAS data for slopes larger than 5° are shown in Table S4.

3.2.2 Derivation of ASTER DEMs and bias corrections

The ASTER DEMs were generated from ASTER L1A stereo images (3N and 3B bands) in early August 2005, using the DEM Extraction Wizard of ENVI 5.1 with the geographic reference system NAD83 UTM subregion 5N or 6N (NASA LP DAAC, 2015; Table S3, data downloaded from https://lpdaac.usgs.gov). Since ground control points were not available, more than 200 tie points for each scene, covering different altitudes and surfaces including ice and snow, were automatically collected and then manually examined (removing those influenced by clouds or water) for all stereo-pairs. The maximum Y parallax errors of the tie points were constrained to less than one pixel (±15 m). The DEMs were generated with 30 m resolution. Artifacts, which mainly resulted from clouds and other noises, were removed using a median filter. All extracted ASTER DEMs cover the glacierized area and surrounding terrain within ~10 km from the ice margin.

The resulting DEMs may contain planimetric and altimetric biases (Moore and others, 2002; Molnia, 2007; Berthier and others, 2010). The DEM biases caused by misregistration between DEMs were corrected by the co-registration method suggested by Nuth and Kääb (2011), which depends on the elevation differences, slope and aspect over stable terrain. The IFSAR DEM was used as master DEM after resampling to 30 m spatial resolution using pixel assembly, and all DEMs (ASTER DEMs and IFSAR DEM) were referenced to the same datum and projection. Only the DEM differences of pixels that do not contain glaciers, water bodies and vegetation (determined by the land cover map of Alaska, USGS (2011)) were used to derive the adjustment coefficients. Then the IQR (see section 3.2.1) was adopted to remove outliers which were typically caused by clouds, data gaps and DEM edges. All iterations of the co-registration adjustment were stopped when the improvement in Std dev. of elevation difference on stable terrain was <2%. Then the along/cross track biases of all ASTER DEMs related to satellite acquisition geometry were corrected using higher order polynomial fittings (5^th to 8^th order; Nuth and Kääb, 2011).

Figure 2 illustrates the co-registration of the ASTER DEM (image AST-4, Table S3) to IFSAR DEM. Table 1 and Figure 3 show the co-registration vectors and correction results of each DEM pair.

Fig. 2. Results of the co-registration of the ASTER DEM (AST-4) to the IFSAR DEM. (a, b) Elevation differences normalized by the slope tangent vs aspect (a) before co-registration and (b) after co-registration, including error bars (blue) and cosine-fit (red line). The fitted equation and coefficient of determination r ² are given in the lower right corner. (c) Residual elevation difference with error bars in the along-track direction. (d) Normal probability density curve of elevation differences. The bin size is 0.0175 (5°) for panels (a) and (b), and 500 and 3 m for panels (c) and (d), respectively.

Fig. 3. Frequency distribution of elevation changes dh over stable terrain and the glacier area during the study period. Bin size is 3 m. Mean and Std dev. is given for various stages of co-registration (Fig. 2). Panels a, b, c and d refer to the subregions I, II, III and IV (Fig. 1), respectively.

Table 1. Statistics of elevation differences between ASTER DEMs (AST-1, AST-2a, b, AST-3, AST-4) and IFSAR DEM over stable terrain after co-registration and after additional along/cross track correction. The grid resolution is 30 m. N is the number of pixels, $\overline {\Delta h}$ is the mean elevation difference (m) and σ _non is the Std dev. of elevation difference (m), dx, dy, dz are the three components of the full co-registration adjustment vector (in meters) between the ASTER DEMs and IFSAR DEM.

3.2.3 Elevation changes

The co-registered ASTER DEMs were subtracted from the IFSAR DEM to calculate elevation difference (dh) over the glacierized areas of the four subregions I–IV (Fig. 1). After scrutinizing the elevation difference maps and difference histograms, we regard pixels with dh dt ⁻¹ values >5 m a⁻¹, < −15 m a⁻¹ for elevations below 1000 m a.s.l. (−10 m a⁻¹ for subregion I) and < −5 m a⁻¹ for elevations higher or equal 1000 m a.s.l. as outliers. This procedure proved to be the most suitable approach in this region compared to other methods, such as empirical value (Berthier and others, 2010; Bolch and others, 2011), IQR (O'Neel and others, 2014; Fischer and others, 2015; Kienholz and others, 2015; Pieczonka and Bolch, 2015).

Data gaps were filled with the mean value of the elevation change in each 50 m elevation band (Fig. 4), which has been proven as an effective method to interpolate voids in DEM difference maps and can tolerate a rather high percentage of data voids, up to ~60% (McNabb and others, 2019).

Fig. 4. Elevation change rates (blue) for all glacierized cells of ASTER DEMs for which data were available (median ± 1 Std dev.) and area-altitude distribution (50 m bin size) of total glacier area (black), and the area with data gaps (pink). The panels a, b, c and d represent the four subregions I, II, III and IV, respectively (Fig. 1).

3.2.4 Volume and mass change

Glacier volume changes for each subregion and period were computed by multiplying the mean surface elevation changes by the glacier area in 2005 (Tables 2 and 3). A volume to mass conversion factor of 850 ± 60 kg m⁻³ (Huss, 2013) was used.

Table 2. Glacier area in 1986, 1995, 2005, 2016 for four regions (Fig. 1) and the area change 1986−2016

Table 3. Glacier-wide mean rates of elevation change (dh dt ⁻¹), volume change (dV dt ⁻¹), specific mass change (dM dt ⁻¹) in each subregion (Fig. 1) for 2005–2014. Area refers to the glacier area in 2005. The region-wide average is the area-weighted average of the four subregions

^a Based on original acquisition dates of IFSAR (29 August−12 September 2014) and ASTER DEM (09 August–14 August 2005).

^b Includes correction for seasonal difference in acquisition data, so that the ASTER data refers to 6 September (Section 3.2.4).

3.2.5 Correcting for difference in DEM acquisition dates

To account for the temporal difference in DEM acquisition (9 August and 14 August 2005 for ASTER and between 29 August and 12 September 2014 for IFSAR), we approximated the mass change that occurred between 11 August and 6 September 2005, and adjusted the area-averaged specific mass change rates obtained above accordingly for each subregion. Following Van Beusekom and others (2010), we used a simple degree-day mass balance model forced by daily temperature and precipitation data from the Wolverine weather station (990 m a.s.l.) to compute idealized mass-balance profiles for the elevation range of the glacier-covered area (0–2000 m a.s.l.). Since the model is not well-constrained, we applied Monte Carlo methods to randomly vary model parameters within reasonable limits and generated 50 profiles. Model parameters and applied ranges are given in Supplementary Table S5. The mass-balance profiles, which provide a mass-change value for each elevation, were assumed constant across the corresponding glacier regions, and used to compute a mass change for each glacierized gridcell. Area-averaged mass changes were computed for each subregion. Over the 27 days, the total mass change in the period 11 August to 6 September 2005 ranged from −0.86 ± 0.17 m w.e. (subregion III) to −0.73 ± 0.17 m w.e. (subregion II) with a region-wide average of −0.78 ± 0.10 m w.e..

3.3 Uncertainty estimation

3.3.1 Glacier outlines

Glacier area errors normally include technical errors (misregistration), interpretation errors and methodological constraints (Raup and others, 2007; Paul and Andreassen, 2009). In this study, the band ratio method was used to derive almost all glacier outlines from Landsat data and no distinct horizontal shift between images was found when visually checking over stable landforms (e.g. mountain ridges, peaks or lateral moraines), so the technical error was considered negligible and ignored (Bolch and others, 2010; Guo and others, 2013). Interpretation errors are mostly associated with the definition of a glacier, while methodological errors largely depend on the resolution of the source imagery, and the experience of compilers. Following Guo and others (2015), we incorporated all error sources into one term, assuming ±10 and ±30 m accuracies for clean-ice and debris-covered glacier outline delineation respectively, and ±5 m for ice divides.

The buffer method was selected to calculate the uncertainty of glacier area generated from Landsat images, which included both the length of glacier outlines and positional accuracies (e.g. Rivera and others, 2005; Granshaw and Fountain, 2006; Kienholz and others, 2015; Tielidze and Wheate, 2018). The area error of each individual glacier E _a, was defined as,

(1)$$E_a = \mathop \sum \limits_{i = 1}^n L_iE_{\,pi}$$

where, i is the type of glacier outlines (clean-ice, debris-covered, ice divides, n = 3) and L _i is the length of glacier outline and E _Pi is the position error. The errors of the position of the ice divides were omitted for considering the whole glacier area of each subregion, and the errors of boundary between clean and debris-covered ice were omitted in the single and regional glacier area error estimation.

The uncertainty of area change E _ΔA over a period t ₁–t ₂, was calculated by

(2)$$E_{\Delta A} = \sqrt {E_{A1}^2 + E_{A2}^2 } $$

where E _A1 and E _A2 are the errors of the glacier area in t ₁ and t ₂, respectively.

3.3.2 Elevation change

Following Zemp and others (2013) and Huber and others (2020), the uncertainty of area-averaged elevation change (σ) was estimated from three components: the uncertainty related to spatial autocorrelation in elevation differences (σ _autocorr), the uncertainty related to the residual elevation errors after co-registration (σ _coreg) and the uncertainty due to data voids (σ _void), although this approach may underestimate the total error (Berthier and others, 2012, 2016).

(3)$$\sigma = \sqrt {\sigma _{autocorr}^2 + \sigma _{coreg}^2 + \sigma _{void}^2 } $$

The Std dev. of the elevation differences between the ASTER DEMs and IFSAR DEM over nonglacierized pixels (stable terrain; σ _non) can be used as a first estimate of the σ _autocorr, if the spatial correlation of the elevation differences is accounted for (Rolstad and others, 2009). Following Gardelle and others (2013), we calculated the standard error (SE) of the mean elevation change:

(4)$${\rm SE} = \displaystyle{{\sigma _{non}} \over {\sqrt {N_{eff}} }}$$

where N _eff is the effective number of the pixels after de-correlation defined by,

(5)$$N_{eff} = \displaystyle{{N_{tot}\cdot R} \over {2d}}$$

where N _tot is the total number of pixels on the nonglacierized area (stable terrain); R is the pixel size (30 m), d is the de-correlation distance assumed to be 600 m as suggested by Bolch and others (2011) and also used in previous geodetic mass-balance studies (Brun and others, 2017; Menounos and others, 2019; Pelto and others, 2019). σ _autocorr was then calculated using SE and the mean elevation differences of the nonglacier region $\lpar \overline {\Delta h} \rpar$:

(6)$$\sigma _{autocorr} = \sqrt {{\rm S}{\rm E}^2 + {\overline {\Delta h} }^2} $$

σ _coreg was calculated using the triangulation method by Nuth and Kääb (2011). Based on three elevation data sets (ASTER DEM, IFSAR DEM and ICESat), we did co-registration by the correction vectors between each of them and then used the residual between the triangulation of vertical vectors to estimate σ _coreg (Table S6).

The DEM difference maps include data voids encompassing 7, 19, 41 and 1% of the four subregions I, II, III and IV, respectively. Based on findings by McNabb and others (2019), who derived differences to true elevation changes as a function of void percentage, we used the constant values of 0.05 and 0.10 m a⁻¹ to estimate σ _void for subregions II and III, respectively, and assumed zero σ _void for subregions I and IV due to nearly complete coverage.

3.3.3 Errors in volume and specific mass changes

Following standard error propagation, we finally obtained the error of the volume change E _Vol and specific mass balance E _Mass (m w.e. a⁻¹), assuming that error terms are independent of each other by

(7)$$E_{Vol} = \sqrt {{\lpar {dh\cdot E_A} \rpar }^2 + {\lpar {A\cdot \sigma } \rpar }^2} $$

where A is the glacier area in 2005, E _A is the corresponding error.

(8)$$E_{Mass} = \sqrt {{\lpar dh\cdot E_\rho \rpar }^2 + {\lpar \sigma \cdot \rho \rpar }^2} $$

where ρ is the assumed volume to mass conversion factor (850 kg m⁻³), $E_\rho$ is the corresponding error (±60 kg m⁻³); dh is area-averaged surface elevation change, while σ is the corresponding elevation change error (Eqn (3)). The mean annual uncertainty is then given by dividing E _Mass by the number of years.

4. Results and discussion

4.1 Glacier inventory in 2016

In the 2016 Landsat inventory (L2016), the 1357 glaciers larger than 0.01 km² cover an area of 3857 ± 87 with 35 ± 12 km² (~1%) of the surface debris covered. Half of the glaciers are smaller than 0.1 km², but they contribute only to <1% to the total area. Only 76 glaciers (~6%) are larger than 5 km², however, they make up about 85% of the total area (Fig. 5a).

Fig. 5. Characteristics of glacier distribution of Kenai Peninsula in 2016. (a) Distribution of glaciers in different area classes; (b) scatterplot of glacier size vs maximum and minimum elevation (Max. elev., Min. elev.). The solid lines (blue and pink) give mean values for distinct size classes while the dark line shows the median elevation per area class for the same area classes as shown in (a); (c) Mean elevation vs aspect per glacier; the solid line refers to the mean per cardinal sector; (d) Mean slope vs glacier size per glacier. The solid line shows the trend (2^nd polynomial fit). X-Axis in (b) and (d) is logarithmically transformed.

Figure 5b indicates that the glaciers' minimum elevation decreases while their maximum elevation increases as the glacier size increase, indicating that glaciers reaching sea-level require a large elevation range, similar to the findings on southern Baffin Island (Paul and Svoboda, 2009). The median elevation of all glaciers fluctuates slightly around ~1100 m a.s.l., regardless area class which may indicate that the equilibrium line of the Kenai Peninsula glaciers is at a similar elevation. There is no correlation between mean glacier elevation and aspect (Fig. 5c). As expected, the mean glacier slope decreases with glacier size, but the scatter is large for smaller glaciers (Fig. 5d) due to local topography (Haeberli and Hoelzle, 1995; Paul and others, 2011).

Figure 6 illustrates the area distribution per elevation, slope and aspect class for each subregion as well as the entire region. Nearly all ice-covered area (>99%) is situated between 150 m and 1750 m a.s.l. (Fig. 6a). The area-elevation distribution of subregion II peaks at a higher elevation than the remaining regions. As shown in Figure 6b, glacier area with slopes <30° occupies 90% of the total area. Figure 6c indicates that there is a tendency for glacierized pixels in all subregions to face towards the northwest or east.

Fig. 6. Area distribution of topography variables for glacierized area in each subregion (I, II, III and IV) and the entire region in 2016 including (a) hypsography with bin size 50 m, and (b) slope with bin size of 1°; (c) aspect, in which numbers represent the glacier area with unit km² in 22.5° aspect bins.

4.2 Area change

Total glacier area has shrunk by 543 ± 123 km² (12 ± 3%) between 1986 and 2016. Glaciers in different regions retreated at different speeds (Table 2) and the highest area loss rate occurred in subregion IV (total area loss of 20 ± 4%).

Larger glaciers tend to have smaller relative area losses than smaller glaciers, but the scatter especially for the smaller glaciers is large (Figs 7a, 8b), consistent with previous studies (Maisch and others, 1999; Serandrei Barbero and others, 1999; Kääb and others, 2002; Paul, 2002; DeBeer and Sharp, 2009; Bolch and others, 2010). The average area loss for glaciers <0.5 km² is nearly 56%, ranging from 0 to 100%, while it is 8% for glaciers larger than 10 km². Approximately 270 glaciers (all <0.5 km²) disappeared during 1986–2016. Larger retreat rates of smaller glaciers can be attributed to a stronger impact of local topographic characteristics on the processes controlling the glacier mass balance (Florentine and others, 2018). For example, a moderate rise in the equilibrium line may turn the entire glacier to ablation zone (Kaser and Osmaston, 2002; Racoviteanu and others, 2008b). A given specific mass change may also result in larger relative area changes for smaller and thinner glaciers than larger glacier with steeper flanks (Kääb and others, 2002). Generally, the area loss rate decreases with increasing elevation (Fig. 7b) except for elevations above 1700 m a.s.l. where the sample size is very small and errors most likely larger due to steep terrain, cast shadows and seasonal snow.

Fig. 7. (a) Relative (blue) and absolute (red) glacier area change between 1986 and 2016 (a) for each glacier (n = 1660) as a function of glacier area, and (b) as a function of elevation (based on all glacier pixels falling within sequential 50 m bins). Solid lines in (a) show the mean values per area class.

Fig. 8. (a) Glacier area-elevation distribution for four different years, and (b–d) glacier area-change rates for different periods between 1986 and 2016 vs (b) area size class, (c) mean slope and (d) aspect. The bin size in panel a is 50 m. Numbers in panel d refer to area change rates in % (10 a)⁻¹.

For glaciers smaller than 0.1 km² (nearly 50% of the total number of glaciers) with a more notable reduction in area than others, the area loss during 2005–2016 is markedly lower than during the earlier periods of 1986–1995 and 1995–2005 (Fig. 8b). One possible explanation is that small glaciers may retreat into more sheltered locations that facilitate their preservation (DeBeer and Sharp, 2009).

We investigated the role of the slope and aspect in spatial variations in area change. Results show that slope played an important role in relative glacier area reduction. Glaciers located on steeper terrains experienced larger relative area loss than less steep glaciers (Fig. 8c) consistent with previous studies (Salerno and others, 2008; Racoviteanu and others, 2015; Tielidze and Wheate, 2018). Southeast-south-southwest facing glaciers experienced generally more pronounced shrinkage than northwest-north-northeast facing glaciers. This is similar to the pattern found in the Himalaya by Ahmad and Rais (1999), which was attributed to more solar radiation on the south-facing than north-facing slopes.

4.3 Glacier elevation and mass changes

Figure 9 shows the spatial distribution of surface elevation changes for the four subregions during 2005–2014. All regions experienced most pronounced thinning in the lower parts especially below 1000 m a.s.l., and less pronounced thickening or thinning in the glaciers' upper reaches.

Fig. 9. Spatially distributed surface elevation changes rates for the period 2005–2014. Dark gray areas represent the cloud error mask of the ASTER DEM, while white area with blue dashes mark the areas with no data. The red line is the 1000 m contour line. The panels a, b, c and d represent the four subregions I, II, III and IV, respectively (Fig. 1).

In subregion I (Grewing-Yalik glacier complex), most thickening occurred above 1000 m a.s.l., but the lower tongue of Portlock Glacier and some smaller glaciers in the southeast show some thickening at lower elevations (Fig. 9a). In subregion II (Harding Icefield), especially the upper reaches of some glaciers in the southwest and southeast experienced thickening. It remains unclear whether this thickening is caused by net snow accumulation, changes in ice dynamics or by errors due to the noise of ASTER DEM (Fig. 9b). Large parts of McCarty Glacier (<1000 m a.s.l.) show a positive surface elevation change, similar to results found by Aðalgeirsdóttir and others (1998) for the period 1950s–1990s.

The two lake-terminating glaciers (Falling and Excelsior glacier), experienced extreme thinning close to their termini (nearly 15 m a⁻¹). The lower tongue of Ellsworth Glacier shows glacier thickening during the studied period (Fig. 9c). However, based on recent satellite imagery, no evidence of a surge could be detected that could explain the local thickening. More studies, such as the surface velocity and its recent dynamics, or the distribution of its thickness, are needed to reveal its causes. Errors in the DEMs may also explain this anomaly.

The regional elevation and mass changes during 2005–2014 for the four subregions are given in Table 3. Results emphasize the need to correct for differences in acquisition dates of any two DEMs to be differenced in regions with considerable mass change over short periods. Without accounting for the 27-days elevation difference between the ASTER DEMs and IFSAR DEM acquisition date, the area-average region-wide mass-balance rate would have been underestimated by 0.09 m w.e. a⁻¹ (~10%).

On average the glaciers on the Kenai Peninsula experienced specific mass changes of −0.94 ± 0.12 m w.e. a⁻¹. Mass-balance rates were most negative for lake-terminating glaciers (−1.37 ± 0.13 m w.e. a⁻¹, ~980 km²), followed by land-terminating glaciers (−1.02 ± 0.13 m w.e. a⁻¹, ~2166 km²) and tidewater glaciers (−0.45 ± 0.14 m w.e. a⁻¹, ~941 km²). This pattern cannot be explained with differences in mean elevation (Table 4). Many lake- or marine-terminating glaciers experienced pronounced thinning especially on their tongues (Wosneseski, Petrof, Skilak, Bear, Pedersen, Falling, Excelsior, Spencer and Trail glaciers; Fig. 9), but not all (for example McCarty glacier). This trend is consistent with airborne laser altimetry derived results across Alaska (Larsen and others, 2015), which indicated more rapid thinning of many lake-terminating glaciers near their termini compared to land-terminating glaciers.

Table 4. Area (% of total), mean elevation and specific mass balance (dM dt ⁻¹ in m w.e. a⁻¹) during the period 2005–2014 for three glacier types (land-terminating, lake-terminating and tidewater glaciers) in each subregion (Fig. 1). Area refer to the inventory from 2005 and elevations to the IFSAR DEM in 2014

Mean specific mass change rates varied considerably among the four subregions (Table 3). Subregion I (Grewing-Yalik glacier complex) experienced the most negative mean mass change (−1.63 ± 0.28 m w.e. a⁻¹), and subregion III (Sargent Icefield) the least negative rate (−0.69 ± 0.21 m w.e. a⁻¹). These differences are consistent with the relative distribution of glacier types. Subregion I has roughly twice as much lake-terminating glacier area (in %) than the other subregions while tidewater glaciers are lacking. In contrast subregion III has a considerably higher percentage of tidewater glacier area compared to the other regions (Fig. 1, Table 4).

Wolverine glacier (~17 km², Fig. S3), located in subregion IV and part of the USGS benchmark glacier monitoring project, had an average glacier-wide mass-balance rate of −0.55 m w.e. a⁻¹ during the period 2005–2014 (O'Neel and others, 2019), is similar to the averaged mass-balance rate of subregion IV (−0.73 ± 0.13 m w.e. a⁻¹). However, we find a less negative mass-balance rate for Wolverine glacier for the same period (−0.20 ± 0.13 m w.e. a⁻¹) and attribute this discrepancy to artifacts in the DEMs which are introduced by the steep local topography (Koblet and others, 2010; Le Bris and Paul, 2015).

4.4 Comparison with other studies

Our new Landsat derived glacier inventory of the Kenai Peninsula for 2005 includes 1421 glaciers covering an area of 4087 ± 84 km², which agrees well with RGI 6.0 (1457 glaciers, ~4175 km², Kienholz and others, 2015) compiled in 2005–2007. Hall and others (2005) found an area reduction of 3.6% (~78 km²) for the Harding Icefield and Grewingk-Yalik Glacier Complex (excluding the glaciers in the surroundings that are not connected to these ice fields) over the time interval 1986–2002. Their results are similar to our results for 1986–2005 for the same ice masses (area loss of 91 ± 22 km², 4.2%). O'Neel and others (2019) found that the Wolverine glacier experienced a 1.5 km² area reduction between 1969 and 2018, which is similar to the area loss of 1.4 km² found in this study for the period 1986–2016.

Aðalgeirsdóttir and others (1998) calculated a surface thinning rate for the Harding Icefield (subregion II) of −0.49 ± 0.12 m a⁻¹ for the period 1950s–1996, which is close to the rate of −0.47 ± 0.01 m a⁻¹ found by VanLooy and others (2006) for roughly the same period (1950 to mid-1990s). The latter study found an accelerated rate of surface elevation change (−0.72 ± 0.13 m a⁻¹) in the period mid-1990s to 1999. For the glaciers in subregions I and II combined, they report a surface elevation change rate of −0.61 ± 0.12 m a⁻¹ for the entire period 1950–1999. Correspondingly, we find accelerated elevation change rates of −1.97 ± 0.29 and −1.26 ± 0.25 m a⁻¹ for subregions I and II (the area-average value is −1.39 ± 0.21 m a⁻¹), respectively, for the period 2005–2014.

For the entire Kenai Range (subregions I-IV, Fig. 1), our results (−0.94 ± 0.12 m w.e. a⁻¹) are considerably more negative than the mass-balance rates of −0.06 ± 0.40 m w.e. a⁻¹ from 1994/1996 to 1999/2001 reported by Arendt and others, 2009, and −0.45 ± 0.11 m w.e. a⁻¹ between 1962 and 2006 found by Berthier and others, 2010, suggesting that glaciers mass loss in the Kenai Peninsula has strongly accelerated since at least the mid-2000s. These findings are in agreement with observed acceleration of glacier mass loss in Western Canada between the periods 2009–2018 and 2000–2009 (Menounos and others, 2019).

Our results are also more negative than regional-scale estimates in other regions in Alaska during similar periods, for example, on the Juneau Icefield and the Stikine Icefield (−0.68 ± 0.15 and −0.83 ± 0.12 m w.e. a⁻¹, respectively) over the period 2000–2016 (Berthier and others, 2018) and also averaged over all glaciers in Alaska and adjacent Yukon (−0.69 ± 0.18 m w.e. a⁻¹) for the period 2006–2015 (Hock and others, 2019).

4.5 Regional climate variability and trends

To investigate the causes of the accelerated glacier mass losses we analyzed air temperature and precipitation records on the Kenai Peninsula (Fig. 10, Table S7). Only four weather stations (A1, A2, A3, A4; Fig. 1) have longer than 30 years of largely uninterrupted observations through at least 2016, while data from the other five stations close to the ice masses (B1–B5, Fig. 1) are only available since 2005. All records show large interannual variation for both air temperature and precipitation (Fig. 10). We applied a moving t-test (Afifi and Azen, 1972; Fu and others, 1999) to the summer and winter temperature data to explore possible abrupt changes in temperature. Results indicate a significant change (at significance level α = 0.05) in 2005 for all stations, except for summer temperature at A1 and A4, where the change occurred in 2002 and 1997, respectively (Fig. 11).

Fig. 10. Time series of (a) summer (May to September) temperature, (b) winter (October to April) temperature and (c) winter (October to April) precipitation for eight weather stations in the Kenai Peninsula (Fig. 1) in the period 1986–2016. Dotted lines are the mean values before and after 2005 when an abrupt change was detected with statistical significance α = 0.05 (Fig. 11). Numbers in parentheses indicate the elevation of each weather station (m a.s.l.). Precipitation data are not available for B3. Data from NOAA, except A4 (Wolverine weather station), from the USGS.

Fig. 11. The moving t-test curve of the air temperatures of four weather stations (A1, A2, A3 and A4) on the Kenai Peninsula. Dashed horizontal lines indicate the 95% confidence level (α = 0.05). Results indicate an abrupt change in 2005 (or 2002 in panel a, 1997 in b). Panels a–d refer to summer temperatures, while panels e–h refer to winter temperatures.

We calculate temperature trends for the period 2005–2016 for the five stations (B1–B5) that lack prior data and over the period 1986–2016 for the four stations (A1–A4) with longer records. Except for station A1 (1986–2016) and station B1 (2005–2016) linear warming trends in summer are not significant (Table S7). Three stations (B3, B4 and B5) show pronounced positive linear trends in winter temperature over the period 2005–2016, however, trends for the stations (A1–A4) going back to 1986 are not significant. In fact, mean summer and winter temperatures of the periods 1986–2005 and 2005–2016 are similar (Fig. 10). Precipitation variations are complex with both positive and negative trends, but only two stations (A1 and A4) have trends (negative) over the period 1986–2016 that are significant (p < 0.05, Table S7). Although decreasing precipitation is consistent with accelerated mass loss, the precipitation records are too short and scarce, especially at high elevation, to elucidate the role of precipitation in the observed glacier changes.

Overall, the large interannual variability and mostly short air temperature and precipitation records with largely insignificant changes of most of the investigated records hamper unambiguous attribution of the observed accelerated specific mass losses of the glaciers on the Kenai Peninsula during 2005–2014 compared to earlier periods covered in previous studies. In addition, changes in frontal ablation (calving and submarine melt) of tidewater glacier due to changes in ocean properties close to the calving fronts may also have contributed to accelerated mass loss, but their role remains unknown.

5. Conclusion

Four new glacier inventories of the Kenai peninsula glaciers in 1986, 1995, 2005 and 2016 were compiled from Landsat images indicating substantial area loss between 1986 and 2016 (543 ± 123 km², ~12%). Despite substantial scatter, relative area losses were generally considerably higher for smaller than larger glaciers, consistent with previous studies elsewhere. Geodetic mass-balance estimates derived from the IFSAR DEM in 2014 and DEMs generated from ASTER images in 2005 reveal substantial thinning and mass loss between 2005 and 2014 (−0.94 ± 0.12 m w.e. a⁻¹). Mass-balance rates vary strongly among the four subregions ranging from −0.69 ± 0.21 m w.e. a⁻¹ (Sargent Icefield) to −1.63 ± 0.28 m w.e. a⁻¹ (Grewing-Yalik glacier complex). These rates are considerably more negative than those found in previous studies for various periods between the early 1960s and late 1999s indicating strong acceleration of mass loss in this region. Thinning is generally most pronounced at lower elevations and both thinning and thickening is observed at higher elevations. The area-averaged specific mass-balance rate of the peninsula's lake-terminating glaciers is almost three times as negative than the rate of the tidewater glaciers and about twice as negative than the rate of the land-terminating glaciers.

Although the acquisition day of the year of the ASTER and IFSAR DEMs differed only by ~27 days, the calculated average region-wide mass-balance rate would have been underestimated by 0.09 m w.e. a⁻¹ (~10%) if it was not corrected for, highlighting the importance of correcting for seasonal differences in acquisition dates in this case. Trends in available temperature and precipitation data on the Kenai Peninsula are mostly insignificant. Records are scarce and often short hampering unambiguous attribution of the acceleration in mass loss over the last decades. Further analysis is needed to determine the exact drivers of this acceleration. Process-based glacier modeling informed by the available data may help to attribute the observed acceleration in glacier mass loss in this region to their causes.