Explaining mass balance and retreat dichotomies at Taku and Lemon Creek Glaciers, Alaska

Abstract

We reanalyzed mass balance records at Taku and Lemon Creek Glaciers to better understand the relative roles of hypsometry, local climate and dynamics as mass balance drivers. Over the 1946–2018 period, the cumulative mass balances diverged. Tidewater Taku Glacier advanced and gained mass at an average rate of +0.25 ± 0.28 m w.e. a^–1, contrasting with retreat and mass loss of −0.60 ± 0.15 m w.e. a⁻¹ at land-terminating Lemon Creek Glacier. The uniform influence of regional climate is demonstrated by strong correlations among annual mass balance and climate data. Regional warming trends forced similar statistically significant decreases in surface mass balance after 1989: −0.83 m w.e. a^–1 at Taku Glacier and −0.81 m w.e. a^–1 at Lemon Creek Glacier. Divergence in cumulative mass balance arises from differences in glacier hypsometry and local climate. Since 2013 negative mass balance and glacier-wide thinning prevailed at Taku Glacier. These changes initiated terminus retreat, which could increase dramatically if calving begins. The future mass balance trajectory of Taku Glacier hinges on dynamics, likely ending the historic dichotomy between these glaciers.

1. Introduction

Mountain glacier mass loss rates in the Alaska region are among the highest on Earth (Gardner and others, 2013; Zemp and others, 2019). Although methods for remotely measuring regional glacier mass change at seasonal, annual and decadal timescales are emerging rapidly, glaciological records reveal processes unresolved by remotely-sensed solutions, and validate modeled estimates of glacier change (Zemp and others, 2013; O'Neel and others, 2019). Glaciological records form the only annual observations of glacier mass change prior to the satellite era, yet comparisons among these records, and to remote sensing or model outputs, are limited by methodological inconsistencies (Zemp and others, 2013). This motivates the consistent, repeatable approach for basin-scale annual glacier surface mass balance that we developed and used to interpret the divergent histories of Taku and Lemon Creek Glaciers.

Given changes in data availability and methods during the past half century, mass balance reanalysis now forms a critical step required prior to interpreting field-based (or glaciological) data (Zemp and others, 2013). Reanalysis identifies and corrects bias identified by comparison of glaciological and geodetic measurements, often incorporates previously unavailable data or methods, and presents a thorough error analysis (Cogley and others, 2011). Such synthesis allows long-term trends to be constrained by geodetic analyses while interannual variations are determined from glaciological observations (Cox and March, 2004).

Of the ~27 000 glaciers in the Alaska region (including neighboring Canada; Keinholz and others, 2015), four have continuous mass balance records exceeding 50 years in length. Taku and Lemon Creek Glaciers present the longest records in North America, measured annually since the mid-1940s and early 1950s. Despite the close proximity (~30 km) of these glaciers (Fig. 1), their mass balance histories have diverged. During the past 70 years, Taku Glacier gained mass and advanced, while Lemon Creek Glacier lost mass and retreated (Criscitiello and others, 2010). This dichotomous response to the same regional climate forcing provided an impetus for glaciological measurements beginning in the mid-20th century (LaChapelle, 1954). Today, Taku Glacier has one of the best documented tidewater glacier advances globally, while Lemon Creek Glacier is the sole benchmark glacier for the World Glacier Monitoring Service and US Geological Survey (USGS) in southeast Alaska. The length and proximity of these records, combined with their differences in response, glacier type and hypsometry constitute a unique glacier mass balance dataset – unmatched in the Alaska region.

Fig. 1. (a) Glaciers (dark shading) of the Alaska Region, with the arrow indicating the location of the Juneau Icefield. (b) Map of the Juneau Icefield highlighting Taku (purple) and Lemon Creek (orange) Glaciers. Symbols in (b) represent weather stations, with elevations provided in the legend. (c) Lemon Creek Glacier with the 1948 area outlined in black and 2018 in orange. JIRP mass balance sites are shown with black circle crosses. (d) Taku Glacier with same annotations as (c). The scale for (d) is an order of magnitude larger than for (c). The inset (e) provides a detailed view of both termini of Taku Glacier with the terminus positions for 1750 (green), 1890 (blue), 1948 (black) and 2018 (purple). The same colors are used for each glacier in all figures.

In this study, we present a reanalysis of surface mass balance time series for Taku and Lemon Creek Glaciers over the 1946–2018 interval. The reanalysis required four primary steps: (1) homogenize geodetic measurements of the area, glacier hypsometry and mass balance; (2) reprocess all recoverable field data; (3) perform geodetic calibrations to remove systematic errors from the glaciological time series; and (4) quantify uncertainties on annual and decadal timescales. We use the reanalysis results to identify the end of the 130-year Taku Glacier advance, and to better understand mass balance and retreat dichotomies in the context of glacier hypsometry, local climate and dynamics (fast flow, frontal ablation and submarine melt).

2. Study area

Taku and Lemon Creek Glaciers are located ~30 km apart on the southern flank of southeast Alaska's Juneau Icefield (Fig. 1). The regional climate is predominantly maritime, but transitions to a more continental climate on the eastern side of the Juneau Icefield. This transition results in a precipitation gradient ranging between 4400 and 2500 mm a⁻¹ (Roth and others, 2018).

Taku Glacier (725 km²) is a tidewater glacier that began advancing in the late 19th century (Motyka and Echelmeyer, 2003). It actively calved into Taku Inlet until 1950 (Motyka and Beget, 1996), when the advance produced a subaerial shoal that prevented direct marine influence (Motyka and Echelmeyer, 2003). The glacier descends from 2130 m to sea level with a median elevation of 1374 m and an average slope of 5°. The glacier has a general southerly aspect with four main tributary branches (Demorest, Matthes, Northwest and Southwest) that converge to form the main branch. The main branch splits and terminates as two piedmont lobes in the Taku River valley (Hole-in-the-Wall and Taku; Fig. 1). Sparse ice thickness measurements indicate the bed remains below sea level ~40 km upstream of the terminus with ice thickness approaching 1500 m (Nolan and others, 1995).

Land-terminating Lemon Creek Glacier is much smaller, with a present-day area of 9.7 km² and a northerly aspect. Surface elevations range from 660 to 1500 m, with a median elevation of 1117 m and a mean slope of 10° (Fig. 1). Lemon Creek Glacier thinned and retreated since studies began (Marcus and others, 1995; Larsen and others, 2015). Thinning has reached such extensive levels that most tributaries have disconnected from the main trunk during the past two decades.

3. Previous studies

The combination of logistical ease and scientific relevance has resulted in a rich research record at the Juneau Icefield. In particular, efforts initiated over 70 years ago by the Juneau Icefield Research Program (JIRP) provide long-term context. During the 1950s and 60s, JIRP made extensive accumulation and ablation observations during mid-September at Taku and Lemon Creek Glaciers (LaChapelle, 1954; Wilson, 1959; Miller, 1975). Since the late 1960s, the bulk of glaciological data collected consist of midsummer accumulation (from the previous winter) and short-duration snow ablation measurements above 1000 m elevation (Fig. 1; Fig. S2). This is due in part to the logistical challenges of accessing the lower part of both glaciers. Additional measurements, made by other research organizations, are discussed in the Supplementary Material (Section S-1).

Beginning in 1946, observations of Accumulation Area Ratio (AAR) and Equilibrium Line Altitude (ELA) occurred at the approximate time of annual mass minimum and provided an essential component in assessing the glacier-wide annual mass balance. Historically, aerial photographs were obtained to estimate the end-of-summer snowline elevation. Spaceborne imagery replaced aerial imagery starting in 1984 (Pelto and Miller, 1990; Miller and Pelto, 1999; Pelto and others, 2013).

Existing glacier-wide annual mass balance estimates for Taku and Lemon Creek Glaciers rely exclusively on JIRP data (Pelto and others, 2013). The JIRP time series consist of glacier-wide annual mass balances, constructed using a mass balance profile, constrained by midsummer measurements, historical data and the annually observed ELA. The annual mass balance profile slope in the accumulation zone is derived annually from midsummer accumulation measurements, while the slope in the ablation zone is held constant, constrained by 1950s and 60s ablation data (Pelto and Miller, 1990; Pelto and others, 2013). Each year, ELA observations connect the accumulation and ablation zone mass balance profile sections. JIRP mass balance profiles were integrated over a fixed glacier hypsometry (1948 for Taku Glacier and 1989 for Lemon Creek Glacier) to determine the glacier-wide annual mass balance in a floating date stratigraphic time system. Previous JIRP time series uncertainty estimates (Taku Glacier = ±0.14 m w.e. a⁻¹; Lemon Creek Glacier = ±0.20 m w.e. a⁻¹) were estimated solely from accumulation data, and did not incorporate errors due to spatial extrapolation or area (Pelto and Miller, 1990; Miller and Pelto, 1999). While recent mass balance studies (e.g. Zemp and others, 2013; Beedle and others, 2014) have incorporated increasingly robust uncertainty assessments, no attempt has been made to improve the error estimates for the JIRP time series (Pelto and others, 2013). Although the timing and distribution of JIRP data collection are sub-ideal, an issue compounded by analysis methods that deviate from standard techniques (e.g. Cogley and others, 2011), a high level of consistency characterizes the entire period of record yielding an analytically consistent dataset. However, a lack of direct incorporation of geodetic mass balance measurements leaves systematic errors in these mass balance records unconstrained.

Standalone geodetic analyses (e.g. Marcus and others, 1995; Berthier and others, 2018) have provided coarse validation of the JIRP time series (Pelto and others, 2013), but to date they have not been integrated with the glaciological data due to analytical inconsistencies in geodetic methods. Geodetic techniques and time intervals are inconsistent during the 1946–2018 period (Table S1). Additionally, these previous geodetic studies exhibit considerable range in associated uncertainties (±0.3 to ±15 m; Larsen and others, 2007, 2015), seasonal variability in acquisitions dates (e.g. Arendt, 2006; Melkonian and others, 2014), inconsistent glacier boundaries and non-uniform application of elevation data co-registration (e.g. Larsen and others, 2007; Berthier and others, 2018). For example, differences in Taku Glacier's area estimates range by 7% (~50 km²; e.g. Pelto and others, 2013; Randolph Glacier Inventory (RGI) Version 3), primarily due to poorly constrained locations of the glacier's nine ice divides.

4. Methods

All data used and produced in this study are publicly available (U.S. Geological Survey Benchmark Glacier Program, 2020). Our glacier-wide mass balance calculations require data of several types: geodetic data, including area and elevation change; meteorological data, including temperature and precipitation; and glaciological data, including accumulation and ablation measurements. Additional data used in this study were downloaded from https://earthexplorer.usgs.gov/ and http://ifsar.gina.alaska.edu/.

4.1 Geodetic methods

The geodetic component of this reanalysis required homogenization of geodetic measurements of the area, hypsometry and mass balance during the 1946–2018 period. To achieve this, we: (1) collected, produced and co-registered Digital Elevation Models (DEMs) and orthorectified imagery (orthoimages); (2) delineated glacier boundaries to produce time-variable glacier hypsometries; and (3) differenced the DEMs to estimate geodetic mass balance.

4.1.1 Geodetic data and co-registration

We used 22 DEMs and orthoimages, 18 of which we produced for this study (Table 1). These were produced using Structure-from-Motion photogrammetry methods (SfM; Westoby and others, 2012) on scanned stereo aerial photographs acquired before 1990. Publicly available DEMs and orthoimages in 2000 and 2013 were derived from Synthetic Aperture Radar (SAR). After 2013, DEMs and orthoimages were produced using Digital Globe imagery and the Ames Stereo Pipeline (Neigh and others, 2013; Shean and others, 2016). Low-resolution imagery, poor accumulation zone contrast and/or incomplete glacier coverage prevented DEM production for several orthoimages. These images were only used to quantify glacier area.

Table 1. Geodetic data used to reanalyze Taku (TG) and Lemon Creek (LCG) Glaciers

The reference DEM is marked with an *.

Use of the 1948 USGS National Elevation Dataset (NED; topographic map-derived) was complicated due to the terminus area of Taku Glacier being updated with more recent imagery (Larsen and others, 2007). To remedy this, a DEM was produced of Taku Glacier's terminus from the original 1948 aerial stereo images. This DEM was inserted into the reprocessed region in the NED to provide a time-consistent initial DEM.

The C-band SAR used to produce the 11 February 2000 Shuttle Radar Topography Mission (SRTM) DEM contains 4 ± 2 m cold snow penetration over temperate glaciers (Rignot and others, 2001). Rather than attempt to correct the SRTM DEM for acquisition date bias, we considered the DEM representative of the 1999 end-of-season glacier surface (i.e. assuming 4 ± 2 m of snow had already accumulated on the glacier by midwinter; Fig. S4). Our assumption is supported by recent comparisons of the SRTM DEM, ASTER DEMs and laser altimetry data collected at Taku Glacier between August 1999 and May 2000 (Berthier and others, 2018), which resolved similar amounts of snow penetration. We used a 27 August 1999 Landsat 5 image for an orthoimage roughly coincident with the SRTM DEM.

DEMs and orthoimages were co-registered using a universal method (Nuth and Kääb, 2011) to minimize systematic misalignment errors. Co-registration was applied to manually selected regions of stable, snow-free ground, covering a broad range of aspect, elevation and area limited to slopes <40° within DEM pairs. Co-registration shifts were also applied to coinciding orthoimages. Each DEM was co-registered to and differenced from the highest quality, or ‘reference’ DEM (Table 1).

4.1.2 Area, hypsometry and geodetic mass balance

Well-defined glacier margins were manually delineated from orthoimages. To minimize the issues arising from seasonal snow cover, we used orthoimages acquired between midsummer and the mass minimum date (Table 1). Tributaries present in 1948 were included throughout the analysis regardless of their connectivity in later years. This prevented step changes in the area and mass balance when tributaries disconnect (Klug and others, 2018). Margin identification along Taku Glacier's nine ice divides relied on flow-field divergence (Burgess and others, 2013; Kienholz and others, 2015), and each divide was assumed to remain stationary throughout the analysis period.

We combined glacier boundaries and co-registered DEMs to produce time-variable glacier hypsometries (Elsberg and others, 2001). For years between DEMs, we linearly interpolated changes in the hypsometry to prevent step changes in time (Harrison and others, 2009; Huss and others, 2012).

Each co-registered DEM was bilinearly resampled to the coarsest resolution DEM, then differenced from the reference DEM (Table 1) to calculate the average surface elevation change (Shean and others, 2016; Sass and others, 2017) over the glacier larger area between each DEM pair. DEMs with unresolved glacier areas >5% (due to low snow contrast; Table 1) were interpolated using an empirical relationship between elevation change in respect to the reference DEM's elevation (Fig. S4; Larsen and others, 2015). Glacier-wide average surface elevation changes were converted to mass changes assuming a constant density of 850 kg m^–3 (Huss and others, 2013).

4.1.3 Area and geodetic mass balance uncertainties

We used an inverse power-law model (Pfeffer and others, 2014) to estimate uncertainty in the glacier area. Uncertainties result from sensor resolution and landscape features (e.g. perennial snow or debris) that impair the exact determination of the glacier margin. These errors scale inversely to the glacier area (Kienholz and others, 2015). Consequently, the relationship between glacier area and error produces proportionally greater uncertainties as glacier area decreases. Additionally, our assumption that ice divides remain stationary on Taku Glacier may not be strictly true, which could increase systematic errors. Hence for Taku Glacier we used an uncertainty of 2% (14 km²), corresponding to the upper bound of the RGI uncertainty model for that glacier area. Lemon Creek Glacier does not have flow divides, so the best estimate from the RGI area uncertainty model was used, which coincidentally is also 2% (0.21 km²).

We assessed geodetic mass balance uncertainty from DEM accuracy and alignment, data gaps and density assumptions using methods identical to the USGS Benchmark Glacier Project (O'Neel and others, 2019). Errors from accuracy and alignment were quantified with the normalized median absolute deviation (NMAD) of post co-registration residual differences over stable, snow-free bedrock, between a given DEM and the reference DEM (Shean and others, 2016). For DEMs with data gaps exceeding 5% of the glacier area, we used the mean absolute error (MAE) of observed versus predicted values (Table 1). In these cases, we combined those errors as the area-weighted average of the NMAD and MAE. Finally, given a lack of firn density observations, we assumed a density and uncertainty of 850 ± 60 kg m⁻³ (Huss and others, 2013). These various uncertainty terms are then summed in quadrature to yield the total uncertainty (σ _G) for any geodetic measurement:

(1)$$\sigma _{\rm G} = \sqrt {\sigma _\rho ^2 \Delta {\overline {{\rm dh}} }^2 + \sigma _{{\rm dh}}^2 \rho ^2} \comma \;$$

where $\sigma _\rho$ represents density uncertainty as a unitless conversion factor of 0.85 (i.e. 850/1000 kg m⁻³), ${\rm \Delta }\overline {{\rm dh\;}}$ is the area-averaged surface elevation change, σ _dh the uncertainty of measured and interpolated area-averaged surface elevation change and ρ the glacier-wide density as a unitless conversion factor of 0.06 (i.e. 60/1000 kg m⁻³; Beedle and others, 2014; O'Neel and others, 2019).

Finally, we assessed geodetic mass balance bias at Taku Glacier's main terminus where basal erosion and sediment excavation (unmeasured by DEM differencing) occurred as the glacier advanced into its own sediment shoal during the past 70 years (Motyka and others, 2006). Ice thickness measurements made in 2003 suggest that the volume of basal erosion since 1948 is less than or equal to the subaerial volume estimate over the region of terminus advance. To estimate an uncertainty bound for this unmeasured volume change, we doubled elevation change values from the 1948–2018 DEM difference map in the terminus advance region (Fig. 2). This change increased the mean mass balance by +0.04 m w.e. a⁻¹, which we included in the upper geodetic uncertainty bounds for Taku Glacier.

Fig. 2. Illustration of sensitivity testing to estimate the sediment excavation bias in the Taku Glacier geodetic assessment. The dashed black line represents sea level and the grey area above sea level is the amount of surface elevation change measured by DEM differencing. The grey area below sea level represents the amount of potentially unmeasured change in the geodetic assessment.

4.2 Glaciological methods

Our glaciological approach follows recently developed methods for the USGS Benchmark Glaciers with two minor changes (O'Neel and others, 2019). We used an optimized mass balance model to accommodate the lack of late-summer ablation observations and constrained the mass balance profile with transient snowline (TSL) observations to fill data gaps in the ablation zone. Missing field data before 1995 at Taku Glacier and 1998 at Lemon Creek Glacier required using pre-existing estimates (similar to South Cascade Glacier; O'Neel and others, 2019), following geodetic calibration and uncertainty evaluation, to extend the records prior to the full reanalysis.

4.2.1 Glaciological measurements

A principal part of this reanalysis involved digitizing midsummer accumulation and ablation observations from JIRP field notebooks and compiling these with other existing glaciological measurements (Section 3; Figs S1, S2). Mid-summer JIRP snowpit data include snowpack thickness and density at fixed locations on the glacier surfaces. Measurements are made between late June and mid-August, with the number of snowpits varying annually. Ablation data consisted of short-term (between several days to weeks) and seasonal duration measurements on both glaciers. Additionally, both short duration and seasonal ice ablation were measured at Taku Glacier during 2004–05 (Kuriger and others, 2006). Finally, seasonal and sub-seasonal ablation measurements of both snow and ice were made on Taku Glacier (2014–15) and Lemon Creek Glacier (2016–18) for this study.

We addressed the lack of ablation zone observations using TSL observations (Kienholz and others, 2017). This involved using available georeferenced imagery (e.g. Landsat, Sentinel, Digital Globe) to extract snowline elevations between April and October each year. Point measurement of 0 m w.e. was assigned at the intersection of the TSL and glacier centerline before the TSL rose above the firn line. TSLs below the firn line were clearly identified at the transition between snow and bare glacier ice. We did not measure TSLs once it approached or passed the firn line to avoid the ambiguity of discerning seasonal snow from firn in low-resolution satellite imagery.

Ablation measurements were exclusively used for mass balance model calibration. Accumulation and TSL measurements were used to initialize the mass balance model each year.

4.2.2 Mass balance model parameter optimization

We used a mass balance model (Van Beusekom and others, 2010; O'Neel and others, 2019) to adjust mid-season measurements to a consistent time system. A positive degree-day (PDD) model (e.g. Hock, 2003) estimates ablation on a site-by-site basis as:

(2)$$a = kT^ + \comma \;$$

where a is the snow or ice ablation (m w.e.), k is an empirically-derived degree-day factor for snow (k _s) or ice (k _i) and T ⁺ is the sum of PDDs (degrees >0 °C). Accumulation was estimated as:

(3)$$c = P_{\lpar T \lt 1.7\rpar }m\comma \;$$

where c represents snow accumulation (m w.e.), P is measured precipitation when lapse rate-adjusted temperatures from the Juneau Airport were below 1.7 °C (Van Beusekom and others, 2010) and m is an empirically-calibrated scaling ratio between the weather station and each glacier. Finally, the change in mass balance was calculated between the measurement date and the glacier-wide mass extrema by summing the two model outputs as:

(4)$$\Delta M = a + c\comma \;$$

where ΔM is the mass balance at a specific site, determined from the daily sum of w.e. ablation, a and accumulation, c.

The model was forced with daily total precipitation and average temperature measured at the Juneau Airport at sea level (Fig. 1; Station ID: USW00025309; ftp.ncdc.noaa.gov; Menne and others, 2012). We optimized model parameters using all available observations from the Juneau Icefield, including short-term and discontinuous stations that are not robust enough to drive the model. First, we used temperature observations from seven sites located on the Juneau Icefield (1100–2100 m a.s.l; Fig. 1) to estimate lapse rates. To determine the optimal lapse rate used in our model, we compared PDDs measured at each weather station and those predicted by scaling the Juneau Airport data, exploring lapse rates from −4.0 to −6.5 °C km⁻¹. The optimal lapse rate was selected using R ², MAE and slope of the linear regression between observed and lapsed PDDs. Next, degree-day factors for snow (k _s), ice (k _i) and precipitation ratios (m) were derived using accumulation and ablation observations (Section 4.2.1) and the inverse of Eqn (2) and (3), with lapsed Juneau Airport data.

4.2.3 Reanalyzed annual mass balance

We used the optimized mass balance model to adjust all midsummer measurements to the mass minimum date in a common floating-date stratigraphic time system (Mayo and others, 1972; Cogley and others, 2011; O'Neel and others, 2019). Mass change after each point measurement date was determined using Eqn (4), applying k _s for point measurements >0 m w.e. (snowpits) and k _i for those ≤0 m w.e. (TSLs). The annual point mass balances were resolved by adding model-predicted mass changes to observed point measurements.

To resolve the glacier-wide annual mass balance, piece-wise linear annual mass balance profiles (O'Neel and others, 2019) were fit to the model-adjusted annual point mass balances and integrated over the glacier hypsometry. The time-variable hypsometry resulted in a conventional, stratigraphic floating-date, glacier-wide annual mass balance for each glacier, henceforth referred to as the reanalyzed time series.

4.2.4 Geodetic calibration and time series assembly

Geodetic calibration refers to the detection and removal of systematic errors in glaciological mass balances that cause a divergence between the cumulative glaciological time series and geodetic time series (Zemp and others, 2013; O'Neel and others, 2019). The glaciological time series is adjusted to the geodetic time series with a constant calibration coefficient, describing the annualized systematic error removed from the glaciological time series (Zemp and others, 2013).

Two time-windowing approaches for geodetic calibration have been commonly applied in previous glacier mass balance studies (O'Neel and others, 2019). The first is a ‘global’ approach (Van Beusekom and others, 2010; O'Neel and others, 2014) which resolves a single calibration coefficient for the entire record. This method's strength is that the calibration is weighted by the geodetic uncertainty of each DEM. However, this approach does not allow for calibrations to vary as either the glacier's geometry or measurement program evolves. This can lead to a poor post-calibration agreement between the glaciological and geodetic time series if the systematic errors vary throughout the record. Alternatively, a ‘sequential’ approach (Zemp and others, 2013; Andreassen and others, 2016) resolves multiple calibrations for discrete geodetic measurement intervals. This allows for time-variable calibration coefficients, as the measurements or glacier geometry changes. However, it does not account for variable geodetic uncertainties, and errors in the geodetic mass balance, which can alias the time series.

Balancing the strengths of those approaches, we used a ‘breakpoint’ geodetic calibration (O'Neel and others, 2019). Prior to calibration, we use the mass balance model to adjust the glaciological data to match the geodetic acquisition date. This choice results in minimal use of the model and removes any volume difference due to variable measurement timing. Breakpoints are inserted when calibration uncertainties exceed the difference between sequential calibrations. This method allows maximum use of geodetic data without the risk of short calibration periods accentuating or aliasing interannual variability.

Missing field data prevented a full glaciological reanalysis during the early part of the study interval, requiring a careful examination of the JIRP time series before a continuous time series could be constructed. Using the breakpoint approach, we geodetically calibrated the JIRP and reanalyzed estimates separately for both glaciers. We then compared the calibrated records for each glacier during the reanalysis period and determined that both approaches exhibited similar interannual variability. We finally extended the reanalysis time period using the calibrated JIRP time series, to produce continuous annual mass balance time series at Taku and Lemon Creek Glaciers during 1946–2018 and 1953–2018, respectively.

4.2.5 Glaciological uncertainties

We estimated random annual mass balance errors using a Monte Carlo simulation. This simulation mirrored sensitivity analyses for the USGS Benchmark Glaciers (O'Neel and others, 2019), but further incorporated glaciological measurement, model adjustment and extrapolation errors. We estimated mean uncertainty from the variance of the Monte Carlo results by randomly varying all components simultaneously.

For each component, a normal distribution of perturbation was defined with the component uncertainty (described below) as 1σ of the distribution curve (Machguth and others, 2008). This simulation was run 1000 times randomly perturbing every component, for each year, in the reanalysis time series. The iterations were geodetically calibrated using each of the approaches described in Section 4.2.4, in addition to our preferred breakpoint approach to produce a range of possible glacier-wide mass balance solutions annually for each glacier. Our final glacier-wide stochastic uncertainty was determined by the average std dev. of glacier-wide annual mass balance solutions resulting from the Monte Carlo simulation during the reanalyzed portion of the time series.

The initial uncertainty source in our reanalyzed time series resulted from annual point mass balance estimation, initially obtained via JIRP snowpits and TSL observations. Previous studies concluded a ±0.08 m w.e. uncertainty for JIRP snowpit measurements (Pelto and Miller, 1990; Pelto and others, 2013), due to low spatial variability in both snow depth (McGrath and others, 2015) and mean snow densities throughout the column (LaChapelle, 1954; Pelto and others, 2013). However, potential observer errors (e.g. previous year's summer surface misidentification) were never included and may increase errors in the accumulation zone. Additionally, we lacked ground control for TSL measurements, limiting quantitative uncertainty estimates (Mernild and others, 2013). For these reasons, we applied ±0.2 m w.e. as a nominal glaciological uncertainty (Lliboutry, 1974; Beedle and others, 2014; Sass and others, 2017). This represented a conservative error bound for all point measurements (both snowpit and TSL) in the Monte Carlo simulation.

Further uncertainty arose from late-summer adjustments of point measurements made using the mass balance model. The majority (96%) of adjustments involved the ablation portion of the model (Eqn 2), and error bounds were therefore represented with the MAE of the PDD calibration. Additionally, the simulation includes a 2% uncertainty in glacier area (Section 4.1.3), distributed randomly across all 100 m elevation bins of the AAD.

Spatial extrapolation of the mass balance profile over the glacier area was the most challenging uncertainty to quantify. Whether the profile captures mass balance (along and across any elevation bin) is generally unknown. The profiles may be only sparsely populated, especially at their tails, and fitting methods may not adequately capture true variability. To explore the sensitivity of extrapolating the mass balance profile over the glacier area, we varied both point balance values and the mass balance profile type. First, to simulate errors in lateral extrapolations, we used the MAE of the point mass balances (±0.38 m w.e. at Taku Glacier; ±0.36 m w.e. at Lemon Creek Glacier) compared to the mass balance profile to perturb the entire mass balance profile. We simulated longitudinal extrapolation errors for the lower and upper portions of the mass balance profile by varying the type of profile used. In addition to our preferred piece-wise mass balance profile (Section 4.2.3), we used two other recently cited mass balance profiles: a linear mass balance profile (Sold and others, 2016) and the site-index method (an area-weighting approach long used by USGS; March and Trabant, 1996; Van Beusekom and others, 2010).

5. Results

We reanalyzed area and annual surface mass balance time series for Taku and Lemon Creek Glaciers from 1946 to 2018. Taku Glacier gained mass during most of the 20th century before entering an equilibrium period and finally transitioning to a state of negative mass balance. Sustained mass loss characterized Lemon Creek Glacier, with accelerating area and mass loss rates throughout the study interval. We show that our results generally agree with previous studies at both glaciers (e.g. Larsen and others, 2007; Criscitiello and others, 2010; Pelto and others, 2013; Berthier and others, 2018), but have reduced and constrained uncertainty. A comparison of our geodetic efforts to pre-existing studies is provided in the Supplementary Material (Section S-4).

5.1 Area

Taku Glacier's area increased between 1948 and 2013 (Table 2), during which time both termini advanced. However, after 1979, the rate of area increase was 73% lower than in the early record. By 2018, a transition to area loss was evident. The glacier lost −1.9 km² since 2013, and both termini retreated (Taku averaged 59 m and Hole-in-the-Wall 85 m; Fig. 3). These recent changes at both termini mark the first retreat since the late 18th century.

Fig. 3. Taku Glacier terminus change between 4 September 2013 (yellow) and 1 October 2018 (purple). Panels (a–e) labeled counterclockwise, specific regions of change to the center panel (f) as indicated with arrows. (a) meltwater pooling between ice and terminal moraine at the western flank of the main terminus, (b) area where Norris River migrated toward the terminus initiating iceberg calving, (c) broad and continuous retreat of the central terminus, (d) outlet stream eroding sediment shoal along the eastern terminus, (e) meltwater pooling between the terminal moraine and the Hole-in-the-Wall Terminus. Orthorectified image from 1 October 2018.

Table 2. Area of Taku and Lemon Creek Glaciers

Conversely, Lemon Creek Glacier's area continuously decreased over the study interval (Table 2). Area loss doubled from the 1948–79 interval (−0.02 km² a⁻¹) to the 1979–2013 interval (−0.05 km² a⁻¹) and further accelerated by a factor of seven (−0.14 km² a⁻¹) between 2013 and 2018 (Table 2). Since 1999, most tributary branches of Lemon Creek Glacier have fragmented from the main branch of the glacier, leaving only 77% of the total glacier area directly flowing to the terminus (Fig. 1).

5.2 Mass balance

A lapse rate of −5.0 °C km⁻¹ to extrapolate Juneau Airport temperature measurements best predicted observed summer PDDs at both glaciers (Fig. 4; Table 3). Model calibration resolved different k _s, k _i, and m values for Taku and Lemon Creek Glaciers (Fig. 5; Table 4). The results suggest similarities in PDD forcing at both glaciers, but variable precipitation and ablation rates. Annual mass balance profiles were steeper on Lemon Creek Glacier, with ELAs averaging 115 m higher than at Taku Glacier (Fig. 6). Through the coincident reanalysis periods (1998–2018), the mean glacier-wide annual mass balance was −0.24 m w.e. a⁻¹ at Taku Glacier and −1.05 m w.e. a⁻¹ for Lemon Creek Glacier (Fig. 7).

Fig. 4. (a–f) Predicted versus observed PDDs at Juneau Icefield weather stations during the 1998–2018 period. Temperatures were reconstructed using constant lapse rates between −6.5 and −4.0 °C km^–1 (labeled) applied to Juneau Airport weather data. Reconstructed temperatures were then compared to those measured at various elevations and locations (symbols) on the Juneau Icefield. The lapse rate with the highest R ² (d) value is −5.0 °C km^–1 and was used in subsequent analyses.

Fig. 5. (a–d) Regressions of PDDs measured at the Juneau Airport with a −5.0 lapse rate and ablation measurements collected on Taku (purple) and Lemon Creek (orange) Glaciers during the 1946–2018 interval. (a–b) Snow ablation measurements, (c–d) ice ablation measurements.

Fig. 6. Mean annual mass balance profiles for (a) Taku Glacier (1995–2018) and (b) Lemon Creek Glacier (1998–2018). Grey dots represent model-adjusted annual point mass balances from snowpit, TSL observations and our mass balance model. The solid line shows the fitted piece-wise linear mass balance profile. The mass balance gradient in the ablation zone is provided in the bottom left corner of each panel, and the mass balance gradient in the accumulation zone in the upper right corner. Black vertical lines show the mean ELA. The 2018 hypsometry is plotted as a histogram with 100 m bins.

Fig. 7. Cumulative annual mass balance time series for Taku Glacier (purple) and Lemon Creek Glacier (orange) for the 1946/1953–2018 interval. Geodetic mass balances are shown for Taku Glacier with circles, and squares for Lemon Creek Glacier (Table 5). Monte Carlo simulation results are shown in gray, which also indicates the portion of the records that was completely reanalyzed.

Table 3. Temperature forcing comparison for variable lapse rates, including fit quality, fit slope and mean absolute error (MAE) of observed verse predicted positive degree days (PDDs)

Table 4. Mass balance model calibration, where n is the total number of data points to derive melt coefficients k _s (snow) and k _i (ice), and precipitation ratios m, and the mean absolute error (MAE) of observed versus predicted ablation or accumulation

The Monte Carlo output of ±0.45 m w.e. a⁻¹ for Taku Glacier and ±0.38 m w.e. a⁻¹ for Lemon Creek Glacier represent our best stochastic uncertainty estimate for any given year. However, this uncertainty is significantly reduced during the geodetic calibration process where average mass balances are evaluated over decadal-scale time intervals determined by the geodetic mass balance. The quadratic sum of geodetic uncertainties between any two measurement intervals range from ±0.66 m w.e. a⁻¹ (2-year interval) to ±0.05 m w.e. a⁻¹ (50-year interval) but vary depending on the measurement interval and glacier (Table 5). Lower uncertainty in the long-term mean mass balance at Lemon Creek Glacier results from higher precision geodetic estimates constraining the mass balance time series.

Table 5. Geodetic mass balances for Taku and Lemon Creek Glaciers between 1948 and 2018

The reference DEM is marked with an *.

Table 6 shows the geodetic calibrations applied to both the JIRP and reanalyzed time series. Both time series are significantly correlated at each glacier (Taku Glacier: r = 0.91; p < 0.001; MAE = 0.24 m w.e. a⁻¹) (Lemon Creek Glacier: r = 0.90; p < 0.001; MAE = 0.33 m w.e. a⁻¹). The MAEs are comparable to glacier-wide annual mass balance uncertainty estimates found in this study, as well as other published error assessments (Zemp and others, 2013; Beedle and others, 2014; Andreassen and others, 2016). The significant correlation between both time series substantiates the decision to splice the JIRP and reanalyzed time series to form a continuous record.

Table 6. Geodetic calibrations for reanalyzed and JIRP time series

Each calibration refers to the systematic error removed from the glaciological time series during the state period. Systematic errors are assumed to accumulate linearly over time, hence, are reported as an annualized quantity.

Over the entire concurrent study interval (1953–2018), we found the mean glacier-wide annual mass balance to be +0.25 ± 0.28 m w.e. a⁻¹ at Taku Glacier and −0.60 ± 0.06 m w.e. a⁻¹ at Lemon Creek Glacier (Fig. 7; Table S3). Despite the difference in the mean annual mass balance of the two glaciers, the anomaly time series are highly correlated (r = 0.91; p < 0.001; MAE = 0.26 m w.e.), demonstrating similar interannual variations at both glaciers. Using an ANOVA test between the 1953–88 and 1989–2018 periods, Taku and Lemon Creek Glaciers' mean annual mass balance decreased significantly (p < 0.001) by −0.83 and −0.81 m w.e. a⁻¹, respectively.

The maximum mass loss rate for both glaciers occurred during 2013–18 (Taku Glacier = −0.84 m w.e. a⁻¹; Lemon Creek Glacier = −1.88 m w.e. a⁻¹). As mass balance declined, the ELA rose, resulting in the complete removal of Lemon Creek Glacier's firnpack. At Taku Glacier, the average AAR decreased from 0.85 (Pelto and others, 2013) to 0.57 by the fall of 2018. Moreover, between 2013 and 2018, Taku Glacier thinned over its entire elevation range, at −1.3 m ice equivalent a⁻¹ (Table 5).

6. Discussion

For the past 70 years, Taku Glacier advanced devoid of marine influence, which allows us to isolate the surface mass balance signal during this portion of the tidewater glacier cycle. We contrast the influence of glacier hypsometry and local climate on the surface mass balances at these two glaciers in the context of a regional warming trend. Finally, we discuss the implications of continued warming and future dynamics on the historic mass balance and retreat dichotomies of Taku and Lemon Creek Glaciers.

6.1 Regional climate influence

The Juneau Airport weather station (Fig. 1) is the longest, continuous daily meteorological record available for the region and has been extensively used to represent regional climate forcing on Juneau Icefield glaciers (e.g. Motyka and others, 2002; Criscitiello and others, 2010; Bieniek and others, 2012; Bieniek and others, 2014; O'Neel and others, 2019). From 1953 to 2018, the mean annual and mean summer (April–September) temperatures have increased by 0.02 and 0.03 °C a^–1 (O'Neel and others, 2019). Congruent with the Juneau Airport temperature history, annual mass balance covaried and decreased at both glaciers. Correlations are strongest in summer, with mean summer temperatures inversely correlated with annual mass balances (Taku r = −0.67, p < 0.001; Lemon Creek r = −0.69, p < 0.001; Fig. 8). The correlations are weaker for annual average temperatures, but still significant (Taku, r = −0.49, p < 0.001; Lemon Creek, r = −0.55, p < 0.001). No significant correlations exist between total winter (October–March) precipitation and the annual mass balance (Taku, r = −0.14, p = 0.25; Lemon Creek, r = −0.20, p = 0.11). The correlation with increasing temperatures and the covarying annual mass balance variability suggest a similar climatic forcing at each glacier and provides a physical mechanism driving the coincident decreases in surface mass balance. While others have ascribed forcing mechanisms (e.g. PDO; Criscitiello and others, 2010) to annual variations in climate and mass balance for the Juneau Icefield, our goal is to explore the relative contributions of glacier geometry and local climate to the observed mass balance response, in light of regional warming trends (Bieniek and others, 2014).

Fig. 8. Comparison of mean summer temperatures measured at Juneau International Airport and glacier-wide annual surface mass balance for Taku Glacier (panel a) and Lemon Creek Glacier (panel b) during the 1953–2018 period. The color of each value indicates the year.

6.2 Hypsometry and local climate forcing

Covariance among annual mass balances within similar climatic regions are substantiated by previous work (e.g. Fountain and others, 1999; Huss and others, 2010; Pelto and others, 2018), while divergent mass balance trends among nearby glaciers have also been reported (e.g. Larsen and others, 2015; Berthier and others, 2018). Taku and Lemon Creek Glaciers are located within a single climatic division (Bieniek and others, 2012), yet dissimilar ELAs and mass balance profiles suggest different local climates. In addition to substantial contrasts in glacier hypsometry (Fig. 6), these glacier-specific characteristics likely dictate their respective surface mass balances given the same regional warming trend. To explore the relative control of glacier hypsometry and local climate on surface mass balance, we produced three scenarios, where hypsometry, ELA and mass balance profiles were exchanged between the glaciers. Each scenario was run over the coincident portions of the reanalyzed record (1998–2018) and then compared to the original reanalysis results (Section 5.2; Fig. 9). In the first scenario, we exchanged glacier hypsometries when integrating point balances to the glacier-wide scale. In the second, we shifted the ELA lower on Lemon Creek Glacier and higher on Taku Glacier each year by the mean difference in ELAs (115 m). For the third scenario, we exchanged the slope of the mass balance profile.

Fig. 9. Mass balance distributions resulting from imposed changes in ELA, hypsometry and annual mass balance profile for (a) Taku Glacier (purple) and (b) Lemon Creek Glacier (orange). ELA refers to shifting the ELA by the mean difference between the glaciers (Taku ELA shifted higher; Lemon Creek lower); Hypsometry refers to an exchange of hypsometries; Balance Profile refers to an exchange of annual mass balance profiles. (c) Differences between each scenario and the uncalibrated annual mass balance (Section 6.2).

Exchanging the hypsometries between the glaciers resulted in a −0.40 m w.e. a⁻¹ decrease to the average mass balance of Taku Glacier and a +0.51 m w.e. a⁻¹ increase to the mean mass balance of Lemon Creek Glacier. Hypsometric indices at both glaciers are equidimensional, meaning each glacier's median surface elevation is centrally located between the respective maximum and minimum elevations (Jiskoot and others, 2009). However, Taku Glacier spans a greater elevation range and has a higher median elevation than Lemon Creek Glacier (Fig. 6). Given the average ELAs, Taku Glacier has a 36% larger AAR than Lemon Creek Glacier, largely due to its greater maximum elevation.

Exchanging the ELAs between the glaciers resulted in a −0.59 m w.e. a⁻¹ decrease of the average mass balance of Taku Glacier and a +0.76 m w.e. a⁻¹ increase of the mean mass balance of Lemon Creek Glacier. In contrast to the expectation of higher inland ELAs due to decreased precipitation (McGrath and others, 2017), the Taku Glacier ELA averaged 115 m lower than at Lemon Creek Glacier (Fig. 6; Table S3). Weather station data from several sites on the icefield show consistent temperature lapse rates (Fig. 4; Table 3) and thus a similar temperature at respective elevations of the two glaciers. However, mass balance model calibrations imply higher snow ablation rates (larger melt coefficients) at Lemon Creek Glacier (Fig. 5; Table 4) despite its more favorable northerly aspect. The difference in coefficients is preserved regardless of the lapse rate tuning parameter in the mass balance model. Attaining similar melt coefficients at both glaciers requires dissimilar lapse rates between the glaciers, which is unsupported by the temperature data (Fig. 4). This suggests a significant difference in energy balances between the two glaciers, possibly due to some possible combination of shorter transport distances for off-glacier dust and debris at the smaller Lemon Creek Glacier, the greater influence of katabatic winds, or more frequent clear and cold nights at Taku Glacier's more continental setting.

Exchanging the annual mass balance profiles between the glaciers resulted in a +0.19 m w.e. a⁻¹ increase to the average mass balance of Taku Glacier and a −0.36 m w.e. a⁻¹ decrease to the mean mass balance of Lemon Creek Glacier. Taku Glacier's accumulation zone is located ~30 km further inland than Lemon Creek Glacier and experiences drier conditions than comparable elevations on Lemon Creek Glacier (Fig. 1; Table 4). Steeper annual mass balance profiles characterize Lemon Creek Glacier (Fig. 6), resulting from higher snow accumulation rates and steeper accumulation gradients in its more coastal location. This is consistent with ground-penetrating radar surveys that show an inverse correlation between glacier accumulation gradients and distance from the coast (McGrath and others, 2015) and modeled accumulation patterns across the Juneau Icefield (Roth and others, 2018).

At both glaciers, the exchange scenarios demonstrate that the ELA exerts the greatest control over mass balance, followed by the glacier hypsometry and finally the slope of the mass balance profile (Fig. 9). The lower ELA and higher peak in hypsometry at Taku Glacier promote positive mass balances by increasing the AAR, compared to those components at Lemon Creek Glacier. However, the slope of the mass balance profile counteracts the influence of the ELA and hypsometry at both glaciers. Applying the lower slope mass balance profile from Taku Glacier to Lemon Creek Glacier reduces accumulation and drives negative mass balances even further into deficit. Similarly, increasing the slope of the mass balance profile drives Taku Glacier's positive mass balance upwards (Fig. 9). This implies that while the ELA and hypsometry of Lemon Creek Glacier favor more negative mass balances, the higher accumulation rates contribute positive mass balance, buffering the net mass loss.

These results highlight the complexity of glacier response to climate, and in part, reveal forcing for the highly variable spatial pattern of glacier change (e.g. Larsen and others, 2015; Berthier and others, 2018). In addition to the importance of hypsometry, the influence of the ELA (likely driven by local climate energy balances) has significant implications on regional mass balance patterns. We have shown that each of those components has a deterministic role in surface mass balance. However, we further explore the implications of dynamics (fast flow, frontal ablation and submarine melt) in the following section, as the future response of Taku Glacier to regional climate warming will possibly depend on the initiation of calving.

6.3 Taku Glacier dynamics and mass balance

Through the past three millennia, Taku Glacier advanced and retreated both synchronously and asynchronously with other Juneau Icefield glaciers in tidewater glacier cycles (Post, 1975; Brinkerhoff and others, 2017) described by Post and Motyka (1995). Previous advances and retreats of the glacier have been forced by both climate (synchronous) and dynamics (asynchronous) (Motyka and Beget, 1996). In its advanced ‘stable’ position, Taku Glacier is one of only four large coastal glaciers in the Alaska region with lower reaches grounded below sea level. As such, it is currently susceptible to rapid calving retreat (Larsen and others, 2015).

Taku's most recent retreat initiated in about 1750, coincident with other Juneau Icefield and neighboring Glacier Bay glaciers (Lawrence, 1950; Post and Motyka, 1995). However, between 1793 and 1890, the glacier began re-advancing asynchronously with nearby glaciers and while the terminus was still in deep water (Nolan and others, 1995; Post and Motyka, 1995), which is unexpected in the simplest model of the tidewater glacier cycle. Several explanations of the premature conclusion of retreat have been suggested: (1) a positively biased surface mass balance resulting from hypsometry and ELA (Section 6.2; Motyka and Beget, 1996); (2) a pinning point near the 1890 terminus (Fig. 1; Nolan and others, 1995) and (3) high sedimentation rates (Brinkerhoff and others, 2017).

Since about 1950, sediment deposition from both Taku Glacier and the Taku River produced a subaerial shoal extensive enough to prevent calving and fast flow (Post and Motyka, 1995). As a result, surface mass balance took a dominant role over dynamic mass balance through the last 70 years. The intersection of the ELA with the glacier hypsometry promoted thickening and advance, despite all other Juneau Icefield glaciers having thinned and retreated (Larsen and others, 2007; Larsen and others, 2015; Berthier and others, 2018).

After more than 130 years of advance, Taku Glacier appears to be transitioning into retreat. From 1948 to 1989, positive surface mass balance (+0.56 ± 0.28 m w.e. a⁻¹; Table S3) drove rapid expansion (+0.48 km² a⁻¹; Table 2). During 1989–2013, near-equilibrium mean annual mass balance (−0.04 ± 0.11 m w.e. a⁻¹) slowed the expansion to +0.13 km² a⁻¹. Since 2013, the average surface mass balance has become negative (−0.84 ± 0.16 m w.e. a⁻¹). The associated rise in ELA reduced the average (1946–2011) AAR from 0.85 (Pelto and others, 2013) to 0.57 by the fall of 2018 (Section 5.2; Fig. 10). Between 2013 and 2018, the area-averaged thinning was −1.3 m ice equivalent a⁻¹ (Table 5) and retreat averaged 59 m at Taku and 85 m at Hole-in-the-Wall (Fig. 3), with a −1.9 km² total reduction in glacier area. We do not know the precise date of retreat initiation, but a 2 August 2015 image (Larsen, 2015) shows the Hole-in-the-Wall terminus had retreated from its terminal moraine in excess of seasonal length changes. However at this time, the main Taku terminus was still within the range of seasonal variability (10–15 m w.e.; Fig. 6; Truffer and others, 2009). By October 2018, both termini receeded well beyond possible the seasonal range with extensive meltwater ponds between the ice face and terminal moraines (Fig. 3).

Fig. 10. (a) Location of the ELA (Sentinal II imagery, 29 September 2018) in blue and the 1989 ELA in yellow (transferred from Landsat imagery, 16 September 1989). (b) Surface elevation changes resulting from DEM differencing (Section 4.1.2) using the 4 September 2013 and 1 October 2018 DEMs (Table 1). Taku Glacier is outlined in purple.

Despite a glacier hypsometry and ELA that favors positive mass balances, Taku Glacier's transition into retreat is being driven by progressive regional climate warming (Section 6.1). Regional climate is projected to continue warming through the 21st century (McAfee and others, 2014; Zeiman and others, 2016), implying that the negative surface mass balance forcing for retreat will also persist. These changes raise the obvious question, will Taku Glacier enter a calving retreat again? The answer depends on the stability of the subaerial shoal currently protecting the terminus from direct ocean interaction and associated frontal ablation. The transition to calving retreat largely depends on the sediment shoal's stability. A breach of the shoal will subject the terminus to calving and submarine ablation (Cogley and others, 2011; Sutherland and others, 2019), subject to further positive feedbacks of dynamics if retreat enters deep water (e.g. Post, 1975; Meier and Post, 1987; Post and others, 2011; Brinkerhoff and others, 2017). However, shoal preservation will likely favor a lake-calving retreat, similar to Mendenhall (Boyce and others, 2007) or Yakutat Glaciers (Trussel and others, 2013), with lower submarine ablation and calving rates relative to marine-terminating counterparts. In either scenario, a retreat of Taku Glacier will likely be faster than typical land-terminating glaciers (Larsen and others, 2015). Furthermore, this increased mass loss driven by dynamics would reverse the historical dichotomy of Taku Glacier as a positive mass balance outlier within the region, to one of extreme of negative mass balance. However, the complex interactions of climate and dynamics challenge the predictions of the pace and scale of future retreat.

Although Taku Glacier's last retreat did not span the entire ~40 km long fjord (Nolan and others, 1995), the current and projected warming trend (Bieniek and others, 2014; McAfee and others, 2014) increases the likelihood that this retreat will. Recent work on marine-terminating glacier retreat, using Taku Glacier as a test case, suggests upstream migration of the grounding line during successive tidewater glacier cycles (Brinkerhoff and others, 2017). During each advance of the glacier, the bed is further eroded producing a deeper fjord during each retreat, favoring terminus buoyancy and instability. Furthermore, regional ocean warming (Royer and Grosch CE, 2006) will likely increase the submarine ablation rates when compared to those during the late-19th century. Combined with the fact that the most over-deepened section is immediately upstream from the 1890 terminus position (Nolan and others, 1995), a calving retreat at Taku Glacier will likely be more extensive than the previous retreat during the 19th century. Consequently, such a retreat would increase the already disproportionate sea-level rise contribution from Alaska Glaciers (Zemp and others, 2019), similar to the retreat at Columbia Glacier (~1000 km² pre-retreat) (Meier and Post, 1987; Post and others, 2011; Larsen and others, 2015).

7. Conclusions

This study revealed significantly different annual mass balances, with a mean of +0.25 ± 0.28 m w.e. a⁻¹ at tidewater Taku Glacier and −0.60 ± 0.15 m w.e. a⁻¹ at land-terminating Lemon Creek Glacier during the 1953–2018 period. However, the time series exhibit strong interannual covariance (r = 0.91; p < 0.001; MAE = 0.26 m w.e. a⁻¹) and correlation with regional warming trends (Taku r = −0.67, p < 0.001; Lemon Creek r = −0.69, p < 0.001). Our findings indicate that coinciding decreases in surface mass balance at both glaciers are related to climate forcing, while the dichotomy in mass balance and glacier response (retreat vs advance) stems from distinct differences in local climate and glacier hypsometry.

Since 2013, a warming regional climate has overwhelmed the influence of local climate and hypsometry that allowed Taku Glacier to gain mass and advance. Persistent negative mass balances and associated thinning during the past 5 years have initiated retreat. Predicting the pace and scale of future retreat is complex, as it will be governed by both climate and dynamics. Approximately 40 km of the glacier is grounded below sea level with ice thicknesses approaching 1500 m, but the present-day terminus is protected by a subaerial shoal. As long as the shoal protects the terminus from warm ocean water, the retreat will likely resemble a lake-terminating glacier. If the terminus is exposed to the ocean, a full tidewater retreat may ensue. Regardless of the character of retreat, we expect Taku Glacier's mass balance to become dominated by dynamic feedbacks, ending the historic mass balance and retreat dichotomies with Lemon Creek Glacier and further increasing the already disproportionate role of the Alaska region glaciers in the global sea-level budget.