1. Introduction
Mountain glacier mass loss rates in the Alaska region are among the highest on Earth (Gardner and others, Reference Gardner2013; Zemp and others, Reference Zemp2019). 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, Reference Zemp2013; O'Neel and others, Reference O'Neel2019). 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, Reference Zemp2013). 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, Reference Zemp2013). 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, Reference Cogley2011). Such synthesis allows long-term trends to be constrained by geodetic analyses while interannual variations are determined from glaciological observations (Cox and March, Reference Cox and March2004).
Of the ~27 000 glaciers in the Alaska region (including neighboring Canada; Keinholz and others, Reference Kienholz2015), 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, Reference Criscitiello, Kelly and Tremblay2010). This dichotomous response to the same regional climate forcing provided an impetus for glaciological measurements beginning in the mid-20th century (LaChapelle, Reference LaChapelle1954). 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−1 (Roth and others, Reference Roth2018).
Taku Glacier (725 km2) is a tidewater glacier that began advancing in the late 19th century (Motyka and Echelmeyer, Reference Motyka and Echelmeyer2003). It actively calved into Taku Inlet until 1950 (Motyka and Beget, Reference Motyka and Beget1996), when the advance produced a subaerial shoal that prevented direct marine influence (Motyka and Echelmeyer, Reference Motyka and Echelmeyer2003). 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, Reference Nolan, Motkya, Echelmeyer and Trabant1995).
Land-terminating Lemon Creek Glacier is much smaller, with a present-day area of 9.7 km2 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, Reference Marcus, Chambers, Miller and Lang1995; Larsen and others, Reference Larsen2015). 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, Reference LaChapelle1954; Wilson, Reference Wilson1959; Miller, Reference Miller1975). 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, Reference Pelto and Miller1990; Miller and Pelto, Reference Miller and Pelto1999; Pelto and others, Reference Pelto, Kavanaugh and McNeil2013).
Existing glacier-wide annual mass balance estimates for Taku and Lemon Creek Glaciers rely exclusively on JIRP data (Pelto and others, Reference Pelto, Kavanaugh and McNeil2013). 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, Reference Pelto and Miller1990; Pelto and others, Reference Pelto, Kavanaugh and McNeil2013). 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−1; Lemon Creek Glacier = ±0.20 m w.e. a−1) were estimated solely from accumulation data, and did not incorporate errors due to spatial extrapolation or area (Pelto and Miller, Reference Pelto and Miller1990; Miller and Pelto, Reference Miller and Pelto1999). While recent mass balance studies (e.g. Zemp and others, Reference Zemp2013; Beedle and others, Reference Beedle, Menounos and Wheate2014) have incorporated increasingly robust uncertainty assessments, no attempt has been made to improve the error estimates for the JIRP time series (Pelto and others, Reference Pelto, Kavanaugh and McNeil2013). 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, Reference Cogley2011), 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, Reference Marcus, Chambers, Miller and Lang1995; Berthier and others, Reference Berthier, Larsen, Durkin, Willis and Pritchard2018) have provided coarse validation of the JIRP time series (Pelto and others, Reference Pelto, Kavanaugh and McNeil2013), 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, Reference Larsen2007, Reference Larsen2015), seasonal variability in acquisitions dates (e.g. Arendt, Reference Arendt2006; Melkonian and others, Reference Melkonian, Willis and Pritchard2014), inconsistent glacier boundaries and non-uniform application of elevation data co-registration (e.g. Larsen and others, Reference Larsen, Motyka, Arendt, Echelmeyer and Geissler2007; Berthier and others, Reference Berthier, Larsen, Durkin, Willis and Pritchard2018). For example, differences in Taku Glacier's area estimates range by 7% (~50 km2; e.g. Pelto and others, Reference Pelto, Kavanaugh and McNeil2013; 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, Reference Westoby, Brasington, Glasser, Hambrey and Reynolds2012) 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, Reference Neigh, Masek and Nickeson2013; Shean and others, Reference Shean2016). 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, Reference Larsen, Motyka, Arendt, Echelmeyer and Geissler2007). 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, Reference Rignot, Echelmeyer and Krabill2001). 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, Reference Berthier, Larsen, Durkin, Willis and Pritchard2018), 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, Reference Nuth and Kääb2011) 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, Reference Klug2018). Margin identification along Taku Glacier's nine ice divides relied on flow-field divergence (Burgess and others, Reference Burgess, Forster and Larsen2013; Kienholz and others, Reference Kienholz2015), 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, Reference Elsberg, Harrison, Echelmeyer and Krimmel2001). For years between DEMs, we linearly interpolated changes in the hypsometry to prevent step changes in time (Harrison and others, Reference Harrison, Cox, Hock, March and Pettit2009; Huss and others, Reference Huss, Hock, Bauder and Funk2012).
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, Reference Shean2016; Sass and others, Reference Sass, Loso, Geck, Thoms and Mcgrath2017) 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, Reference Larsen2015). Glacier-wide average surface elevation changes were converted to mass changes assuming a constant density of 850 kg m–3 (Huss and others, Reference Huss2013).
4.1.3 Area and geodetic mass balance uncertainties
We used an inverse power-law model (Pfeffer and others, Reference Pfeffer2014) 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, Reference Kienholz2015). 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 km2), 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 km2).
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, Reference O'Neel2019). 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, Reference Shean2016). 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−3 (Huss and others, Reference Huss2013). These various uncertainty terms are then summed in quadrature to yield the total uncertainty (σ G) for any geodetic measurement:
where $\sigma _\rho$
represents density uncertainty as a unitless conversion factor of 0.85 (i.e. 850/1000 kg m−3), ${\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−3; Beedle and others, Reference Beedle, Menounos and Wheate2014; O'Neel and others, Reference O'Neel2019).
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, Reference Motyka, Truffer, Kuriger and Bucki2006). 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−1, 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, Reference O'Neel2019). 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, Reference O'Neel2019), 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, Reference Kuriger, Truffer, Motyka and Bucki2006). 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, Reference Kienholz, Hock, Truffer, Bieniek and Lader2017). 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, Reference Van Beusekom, O'Neel, March, Sass and Cox2010; O'Neel and others, Reference O'Neel2019) to adjust mid-season measurements to a consistent time system. A positive degree-day (PDD) model (e.g. Hock, Reference Hock2003) estimates ablation on a site-by-site basis as:
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:
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, Reference Van Beusekom, O'Neel, March, Sass and Cox2010) 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:
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, Reference Menne, Durre, Vose, Gleason and Houston2012). 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−1. The optimal lapse rate was selected using R 2, 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, Reference Mayo, Meier and Tangborn1972; Cogley and others, Reference Cogley2011; O'Neel and others, Reference O'Neel2019). 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, Reference O'Neel2019) 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, Reference Zemp2013; O'Neel and others, Reference O'Neel2019). 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, Reference Zemp2013).
Two time-windowing approaches for geodetic calibration have been commonly applied in previous glacier mass balance studies (O'Neel and others, Reference O'Neel2019). The first is a ‘global’ approach (Van Beusekom and others, Reference Van Beusekom, O'Neel, March, Sass and Cox2010; O'Neel and others, Reference O'Neel, Hood, Arendt and Sass2014) 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, Reference Zemp2013; Andreassen and others, Reference Andreassen, Elvehøy, Kjøllmoen and Engeset2016) 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, Reference O'Neel2019). 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, Reference O'Neel2019), 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, Reference Machguth, Purves, Oerlemans, Hoelzle and Paul2008). 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, Reference Pelto and Miller1990; Pelto and others, Reference Pelto, Kavanaugh and McNeil2013), due to low spatial variability in both snow depth (McGrath and others, Reference McGrath2015) and mean snow densities throughout the column (LaChapelle, Reference LaChapelle1954; Pelto and others, Reference Pelto, Kavanaugh and McNeil2013). 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, Reference Mernild2013). For these reasons, we applied ±0.2 m w.e. as a nominal glaciological uncertainty (Lliboutry, Reference Lliboutry1974; Beedle and others, Reference Beedle, Menounos and Wheate2014; Sass and others, Reference Sass, Loso, Geck, Thoms and Mcgrath2017). 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, Reference Sold2016) and the site-index method (an area-weighting approach long used by USGS; March and Trabant, Reference March and Trabant1996; Van Beusekom and others, Reference Van Beusekom, O'Neel, March, Sass and Cox2010).
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, Reference Larsen, Motyka, Arendt, Echelmeyer and Geissler2007; Criscitiello and others, Reference Criscitiello, Kelly and Tremblay2010; Pelto and others, Reference Pelto, Kavanaugh and McNeil2013; Berthier and others, Reference Berthier, Larsen, Durkin, Willis and Pritchard2018), 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 km2 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 km2 a−1) to the 1979–2013 interval (−0.05 km2 a−1) and further accelerated by a factor of seven (−0.14 km2 a−1) 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−1 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−1 at Taku Glacier and −1.05 m w.e. a−1 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 2 (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−1 for Taku Glacier and ±0.38 m w.e. a−1 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−1 (2-year interval) to ±0.05 m w.e. a−1 (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−1) (Lemon Creek Glacier: r = 0.90; p < 0.001; MAE = 0.33 m w.e. a−1). 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, Reference Zemp2013; Beedle and others, Reference Beedle, Menounos and Wheate2014; Andreassen and others, Reference Andreassen, Elvehøy, Kjøllmoen and Engeset2016). 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−1 at Taku Glacier and −0.60 ± 0.06 m w.e. a−1 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−1, respectively.
The maximum mass loss rate for both glaciers occurred during 2013–18 (Taku Glacier = −0.84 m w.e. a−1; Lemon Creek Glacier = −1.88 m w.e. a−1). 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, Reference Pelto, Kavanaugh and McNeil2013) 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−1 (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, Reference Motyka, O'Neel, Connor and Echelmeyer2002; Criscitiello and others, Reference Criscitiello, Kelly and Tremblay2010; Bieniek and others, Reference Bieniek2012; Bieniek and others, Reference Bieniek, Walsh, Thoman and Bhatt2014; O'Neel and others, Reference O'Neel2019). 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, Reference O'Neel2019). 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, Reference Criscitiello, Kelly and Tremblay2010) 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, Reference Bieniek, Walsh, Thoman and Bhatt2014).
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, Reference Fountain and Vecchia1999; Huss and others, Reference Huss, Hock, Bauder and Funk2010; Pelto and others, Reference Pelto2018), while divergent mass balance trends among nearby glaciers have also been reported (e.g. Larsen and others, Reference Larsen2015; Berthier and others, Reference Berthier, Larsen, Durkin, Willis and Pritchard2018). Taku and Lemon Creek Glaciers are located within a single climatic division (Bieniek and others, Reference Bieniek2012), 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−1 decrease to the average mass balance of Taku Glacier and a +0.51 m w.e. a−1 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, Reference Jiskoot, Curran, Tessler and Shenton2009). 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−1 decrease of the average mass balance of Taku Glacier and a +0.76 m w.e. a−1 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, Reference McGrath, Sass, O'Neel, Arendt and Kienholz2017), 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−1 increase to the average mass balance of Taku Glacier and a −0.36 m w.e. a−1 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, Reference McGrath2015) and modeled accumulation patterns across the Juneau Icefield (Roth and others, Reference Roth2018).
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, Reference Larsen2015; Berthier and others, Reference Berthier, Larsen, Durkin, Willis and Pritchard2018). 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, Reference Post1975; Brinkerhoff and others, Reference Brinkerhoff, Truffer and Aschwanden2017) described by Post and Motyka (Reference Post and Motyka1995). Previous advances and retreats of the glacier have been forced by both climate (synchronous) and dynamics (asynchronous) (Motyka and Beget, Reference Motyka and Beget1996). 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, Reference Larsen2015).
Taku's most recent retreat initiated in about 1750, coincident with other Juneau Icefield and neighboring Glacier Bay glaciers (Lawrence, Reference Lawrence1950; Post and Motyka, Reference Post and Motyka1995). 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, Reference Nolan, Motkya, Echelmeyer and Trabant1995; Post and Motyka, Reference Post and Motyka1995), 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, Reference Motyka and Beget1996); (2) a pinning point near the 1890 terminus (Fig. 1; Nolan and others, Reference Nolan, Motkya, Echelmeyer and Trabant1995) and (3) high sedimentation rates (Brinkerhoff and others, Reference Brinkerhoff, Truffer and Aschwanden2017).
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, Reference Post and Motyka1995). 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, Reference Larsen, Motyka, Arendt, Echelmeyer and Geissler2007; Larsen and others, Reference Larsen2015; Berthier and others, Reference Berthier, Larsen, Durkin, Willis and Pritchard2018).
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−1; Table S3) drove rapid expansion (+0.48 km2 a−1; Table 2). During 1989–2013, near-equilibrium mean annual mass balance (−0.04 ± 0.11 m w.e. a−1) slowed the expansion to +0.13 km2 a−1. Since 2013, the average surface mass balance has become negative (−0.84 ± 0.16 m w.e. a−1). The associated rise in ELA reduced the average (1946–2011) AAR from 0.85 (Pelto and others, Reference Pelto, Kavanaugh and McNeil2013) 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−1 (Table 5) and retreat averaged 59 m at Taku and 85 m at Hole-in-the-Wall (Fig. 3), with a −1.9 km2 total reduction in glacier area. We do not know the precise date of retreat initiation, but a 2 August 2015 image (Larsen, Reference Larsen2015) 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, Reference Truffer, Motyka, Hekkers, Howat and King2009). 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, Reference McAfee, Walsh and Rupp2014; Zeiman and others, Reference Ziemen2016), 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, Reference Cogley2011; Sutherland and others, Reference Sutherland2019), subject to further positive feedbacks of dynamics if retreat enters deep water (e.g. Post, Reference Post1975; Meier and Post, Reference Meier and Post1987; Post and others, Reference Post, O'Neel, Motyka and Streveler2011; Brinkerhoff and others, Reference Brinkerhoff, Truffer and Aschwanden2017). However, shoal preservation will likely favor a lake-calving retreat, similar to Mendenhall (Boyce and others, Reference Boyce, Motyka and Truffer2007) or Yakutat Glaciers (Trussel and others, Reference Trüssel, Motyka, Truffer and Larsen2013), 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, Reference Larsen2015). 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, Reference Nolan, Motkya, Echelmeyer and Trabant1995), the current and projected warming trend (Bieniek and others, Reference Bieniek, Walsh, Thoman and Bhatt2014; McAfee and others, Reference McAfee, Walsh and Rupp2014) 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, Reference Brinkerhoff, Truffer and Aschwanden2017). 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, Reference Royer and Grosch2006) 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, Reference Nolan, Motkya, Echelmeyer and Trabant1995), 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, Reference Zemp2019), similar to the retreat at Columbia Glacier (~1000 km2 pre-retreat) (Meier and Post, Reference Meier and Post1987; Post and others, Reference Post, O'Neel, Motyka and Streveler2011; Larsen and others, Reference Larsen2015).
7. Conclusions
This study revealed significantly different annual mass balances, with a mean of +0.25 ± 0.28 m w.e. a−1 at tidewater Taku Glacier and −0.60 ± 0.15 m w.e. a−1 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−1) 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.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2020.22
Acknowledgements
Hundreds of people have contributed to the collection of these datasets during the past seven decades. Specific acknowledgement needs to be made of William O. Field, Richard Hubley, Edward LaChapelle, Austin Post and Maynard Miller for initiating long-term data collection on such a large spatial scale. Logistical support by the Juneau Icefield Research Program has been invaluable for data collection both historically and currently. Additional gratitude to Jason Amundson and Martin Truffer must be conveyed for supplying data from lower Taku Glacier in 2014 and 2015. Comments from Caitlyn Florentine, Roman Motyka, and an anonymous reviewer immensely improved this manuscript. Finally, unwavering support from Kuma McNeil was critical to this project between 2012 and 2019. We thank the Geologic Society of America and the Arctic Institute of North America for partially funding this research. Primary funding for the work came from the US Geological Survey under the guise of many programs throughout the years. Previous support came through the USGS National Institutes of Water Resources. Today, support is granted by the USGS Land Resources Mission Area, Research and Development Program. Any use of trade, firm or product names is for descriptive purposes only and does not imply endorsement by the US Government.
Author contributions
Christopher McNeil (CM), Shad O'Neel (S0), Michael Loso (ML), Mauri Pelto (MP), Louis Sass (LS), Emily H. Baker (EB), Seth Campbell (SC). Concept and design CM, SO, LS, ML, MP. Programming CM, LS. Formal analyses CM. Data curation EB, CM. Writing CM, SO, ML, LS, MP, EB, SC. Edit and review: CM, SO, ML, LS, MP, EB, SC. Visualization CM. Funding acquisition and project administration/supervision SO, ML, CM, SC
Open access
