Introduction
Glacier retreat during the last century is one of the most prominent icons of anthropogenic climate change. At the global scale, glacier loss contributes substantially to sea level rise, but terminus retreat also has important impacts at local scales (e.g. Moore and others, Reference Moore2009). Changes to seasonal streamflow, stream temperature and geohazards can depend directly on the response of a few key glaciers. And despite the ubiquity of glacier retreat around the world (Leclercq and others, Reference Leclercq2014) and confident attribution to climate change (Roe and others, Reference Roe, Baker and Herla2017), differences exist in the scale and pace of observed retreat, even among glaciers in the same mountain range. Understanding the basis for these differences in transient glacier response is important for understanding local observations and impacts, and is a prerequisite for making reliable projections.
In this study, we investigate the differences in glacier evolution that arise from different response times among glaciers. The response time of the terminus position to climate variations is a fundamental aspect of glacier dynamics, and is important to understand in the context of a rapidly changing climate. Past studies have laid out the theoretical basis of response times, based on fundamental glacier characteristics (e.g. Jóhannesson and others, Reference Jóhannesson, Raymond and Waddington1989; Harrison and others, Reference Harrison, Elsberg, Echelmeyer and Krimmel2001; Oerlemans, Reference Oerlemans2001; Roe and Baker, Reference Roe and Baker2014). Response times can also be assessed empirically from terminus and climate observations (Harper, Reference Harper1993; Oerlemans, Reference Oerlemans2007) or numerical model output (e.g. Leysinger Vieli and Gudmundsson, Reference Leysinger Vieli and Gudmundsson2004; Zekollari and others, Reference Zekollari, Fürst and Huybrechts2014, Reference Zekollari, Huss and Farinotti2020).
Essentially all of these approaches give response times from 10 to 100 years for most mountain glaciers, with variations depending on glacier geometry and climate setting. The variety of glaciers within a given region can thus be expected to yield a variety of response times. This is an important distribution to assess, especially with growing interest in regional- to global-scale glacier simulations (e.g. Hock and others, Reference Hock2019). A number of previous studies have done so, using simplified scaling arguments (Haeberli and Hoelzle, Reference Haeberli and Hoelzle1995; Hoelzle and others, Reference Hoelzle2007), as well as more sophisticated models (Zekollari and others, Reference Zekollari, Huss and Farinotti2020). These studies both targeted the European Alps, and also New Zealand in the case of Hoelzle and others (Reference Hoelzle2007).
In this study, we focus on the glaciers of the Cascade mountain range of Washington State, USA. A number of previous studies have discussed the transient behavior of glaciers in the Cascades. Several have analyzed the climatic drivers and glacier characteristics that led a number of glaciers in the region to advance in the 1950s–1970s (Hubley, Reference Hubley1956; Harper, Reference Harper1993; Pelto and Hedlund, Reference Pelto and Hedlund2001; Rasmussen and Conway, Reference Rasmussen and Conway2001; Stevens and others, Reference Stevens, Conway, Kennard, Rasmussen and Koutnik2018). For example, Rasmussen and Conway (Reference Rasmussen and Conway2001) interpreted differences between South Cascade glacier and Blue Glacier (in the nearby Olympic Mountains) partly as a result of their different adjustment timescales. Additionally, Pelto and Hedlund (Reference Pelto and Hedlund2001) estimated response times for 21 glaciers in the Cascades, and drew conceptual links between these estimates and terminus behavior.
Our goal is to build upon these studies, assessing the full distribution of response times across the Cascades. We estimate response times with a simple scaling method, using geometric attributes provided in the Randolph Glacier Inventory version 6 (RGI Consortium, 2017, hereafter RGIv6), constrained by local mass-balance and ice-thickness observations. We then explore the dynamical implications of these estimates using theoretical tools developed in the last few years (Roe and Baker, Reference Roe and Baker2014; Christian and others, Reference Christian, Koutnik and Roe2018). The goal of this study is not to produce a detailed simulation of these glaciers, but rather to assess, based on robust physical principles, how transient responses vary within this population of glaciers. We analyze elements of glacier change where the response times found in the Cascades imply a consequential range of behaviors for individual glaciers. These include the degree of disequilibrium with current climate, the response to climate variability and changes in glacier runoff.
Methods
Transient models of glacier response
In this study, we use a response time τ proposed by Jóhannesson and others (Reference Jóhannesson, Raymond and Waddington1989):
where H is a characteristic glacier thickness and b t is the (negative) annual mass-balance rate near the glacier terminus. Equation (1) is a straightforward reservoir timescale: H scales the volume change associated with a given length change (assuming constant width), and b t scales the volume output rate where the reservoir is changing (i.e. the terminus). Furthermore, H and b t implicitly capture the essential dynamics of a glacier in steady state: it has evolved toward the thickness and terminus position necessary for ice flow to balance the pattern of accumulation and melt imposed by the landscape and climate.
Simplified models have been developed on the assumption that these same dynamics – implicitly represented in the glacier geometry – govern climate-driven departures from steady state (e.g. Jóhannesson and others, Reference Jóhannesson, Raymond and Waddington1989; Harrison and others, Reference Harrison, Elsberg, Echelmeyer and Krimmel2001; Oerlemans, Reference Oerlemans2001; Lüthi and others, Reference Lüthi, Bauder and Funk2010; Roe and Baker, Reference Roe and Baker2014). We use a linear model for glacier length that assumes a response time given by Eqn (1) (Roe and Baker, Reference Roe and Baker2014, hereinafter RB14). For simple geometries, the model accurately emulates the transient terminus behavior of flowline glacier models (Christian and others, Reference Christian, Koutnik and Roe2018, RB14), and provides useful analytical solutions. The RB14 model gives length anomalies (L′(t)) described by a third-order ordinary differential equation:
Here, τ is defined as above; note that it is not an e-folding timescale. That is, length perturbations do not decay exponentially, which would imply the terminus begins reacting immediately to a climate change. Instead, the RB14 captures an initial lag associated with thickness changes, as is seen in numerical models (RB14). β = A tot/(wH) where A tot is the total glacier area and w is the width near the terminus; and $\epsilon = 1/\sqrt {3}$
. b′(t) are surface-mass-balance anomalies assumed uniform over a fixed surface, and thus reflect only climate changes. The glacier-averaged balance (often denoted B) can be calculated from the RB14 model by incorporating the evolving geometry (Roe and others, Reference Roe, Christian and Marzeion2021), but our focus here is on length anomalies.
The equilibrium length sensitivity to mass-balance perturbations is
This is a simple statement of mass conservation, expressing the change in ablation area (w ΔL eq) that must accommodate a perturbation Δb applied over the entire glacier (Jóhannesson and others, Reference Jóhannesson, Raymond and Waddington1989).
An important consequence of multidecadal response times in an era of anthropogenic warming is that glacier-terminus positions are out of equilibrium with the current climate (e.g. Lüthi and others, Reference Lüthi, Bauder and Funk2010; Zekollari and others, Reference Zekollari, Fürst and Huybrechts2014; Christian and others, Reference Christian, Koutnik and Roe2018; Marzeion and others, Reference Marzeion, Kaser, Maussion and Champollion2018; Zekollari and others, Reference Zekollari, Huss and Farinotti2020). The RB14 model has an analytical solution for the response to a linear mass-balance trend (first derived in Roe and others, Reference Roe, Baker and Herla2017). Christian and others (Reference Christian, Koutnik and Roe2018) used these solutions to show that the degree of disequilibrium depends fundamentally on τ, and can be described by some simple metrics, illustrated in Figure 1. If a linear mass-balance trend $b'( t) = \dot {b}t$
begins at t = 0, the transient terminus response (L′) lags the instantaneous equilibrium response (L′eq, dashed line in Fig. 1a). Hereinafter, disequilibrium will refer to the length difference, L′ − L′eq (Fig. 1b). For a continuous trend, the RB14 solution predicts that disequilibrium asymptotes to a constant:
While this expression is a linear approximation, it shows three leading controls on disequilibrium: a factor of the response time (3ετ), the sensitivity (τβ) and the slope of the trend in mass-balance forcing ($\dot {b}$
, which has units of m a−2). Disequilibrium, as defined here, is equivalent to the additional committed retreat if the climate were to stabilize at a given time.
Schematic of idealized glacier retreat for a glacier with τ = 25 years. Adapted from Christian and others (Reference Christian, Koutnik and Roe2018), Fig. 3. (a) Transient (solid) and equilibrium (dashed) length responses to continuous negative trend in mass balance starting at t = 0. (b) Disequilibrium is defined as the difference between transient and equilibrium length responses, and corresponds to the additional committed retreat at any given time. (c) Fractional equilibration is defined as the ratio of transient to equilibrium responses.
Let fractional equilibration refer to the evolving ratio of L′ to L′eq (Fig. 1c; see Christian and others, Reference Christian, Koutnik and Roe2018). Length sensitivity cancels, and so fractional equilibration depends only on τ and the time since the onset of forcing:
This metric can be used to compare the state of glaciers independent of their length sensitivity, or to scale observed retreats (L′) to the equilibrium response (L′eq), provided that τ and t can be estimated.
We apply the RB14 model in two different ways in this study. First, after estimating individual response times for glaciers in the Cascades (described next), we apply Eqn (5) so that we can assess their fractional equilibration simply in terms of τ. The fractional metric is easy to compare between different glaciers because it is a normalized quantity, and it also eliminates the need to estimate additional parameters (e.g. β, which is related to the length sensitivity) for all glaciers. For the second set of model analyses, we use the model (Eqn (2)) to analyze additional consequences of these response times. We do this using two synthetic glaciers broadly representative of those in the Cascades, whose parameters are given in the Appendix.
Note that Eqns (4) and (5) come from the particular solution of the RB14 model to a linear trend. They are useful approximations for the component of glacier change due to long-term forcing, but of course other climate variations have also played a role for real glaciers. For example, many studies indicate that the rate of anthropogenic forcing accelerated in the in mid-20th century (e.g. Haustein and others, Reference Haustein2019), making Eqn (5) likely to overestimate the level of equilibration. We investigate the effects of changes in the rate of forcing, as well as climatic noise, in the second set of model analyses.
Estimating response times from inventory data
The Cascade Mountains of Washington State, USA, have a temperate, maritime climate and the largest glacierized area in the contiguous USA (e.g. Fountain and others, Reference Fountain, Glenn and Basagic IV2017). Most of the range's largest glaciers are found on stratovolcanoes: Mount Baker, Glacier Peak, Mount Rainier and Mount Adams (the smaller Mount Saint Helens, which erupted in 1980, hosts much smaller glaciers). The four major volcanoes are the highest peaks in the Cascades, and their glaciers are characterized by large elevation spans and steep slopes. The non-volcanic peaks of the Cascades host smaller valley and cirque glaciers.
RGIv6 reports 1709 glacier outlines for the Cascade range in Washington State. We consider only glaciers with a reported area >0.1 km2 and elevation span >250 m, reducing the sample to 383 glaciers. RGIv6 reports basic geometric parameters, but in general, observations of H and b t exist for only a small fraction of glaciers and are not compiled in the inventory. Haeberli and Hoelzle (Reference Haeberli and Hoelzle1995) and Hoelzle and others (Reference Hoelzle2007) used a simple algorithm to estimate H and b t from inventory data, which we adapt here to the Cascades.
Thickness
Following Haeberli and Hoelzle (Reference Haeberli and Hoelzle1995), we estimate characteristic ice thickness from the glacier's average slope (α) and an assumed basal shear stress (S b). Maximum basal shear stresses are typically on the order of 105 Pa across a wide range of geometries, because of the nearly plastic rheology of ice (e.g. Nye, Reference Nye1952; Cuffey and Paterson, Reference Cuffey and Paterson2010). Thus, assuming a characteristic value of S b for all glaciers, ice thickness is approximated as
where f = 0.8 is a shape factor, ρ i = 900 kg m−3 is the approximate density of glacier ice and g = 9.81 m s−2 is gravitational acceleration. α is estimated by $\alpha = \arctan {( Z_{\rm max}-Z_{\rm min}) /L_{\rm max}}$
where altitudes (Z max, Z min) and length (L max) are taken from RGIv6.
We set S b = 1.5 × 105 Pa, which yields reasonable agreement with observational estimates of mean thickness. Table 1 lists published thickness estimates for the few glaciers with observations from radar sounding or borehole depths. Average slope (α) is shown as well. Note that the observations are concentrated on relatively steep glaciers (α ~ 17 − 25°), with the exception of the much flatter South Cascade glacier (α ~ 10°). Though South Cascade's RGIv6 area is smaller than most of the others listed in the table, its ice thickness of up to ~200 m (Hodge, Reference Hodge1979; Fountain, Reference Fountain1994) supports that slope is indeed a key predictor for H. As τ is proportional to H in this framework (Eqn (1)), this is also consistent with empirical studies demonstrating the impact of slope on response time (Leysinger Vieli and Gudmundsson, Reference Leysinger Vieli and Gudmundsson2004; Zekollari and others, Reference Zekollari, Huss and Farinotti2020).
Comparison of published thickness observations, numerical estimates from Huss and Farinotti (Reference Huss and Farinotti2012) (HF12, data published in Farinotti and others, Reference Farinotti2019) and simple shear-stress scaling (Eqn (6)). Glacier slopes are also shown for reference.
Thickness observations come from aDriedger and Kennard (Reference Driedger and Kennard1986a), bDriedger and Kennard (Reference Driedger and Kennard1986b), cHarper (Reference Harper1992), dHodge (Reference Hodge1979), eFountain (Reference Fountain1994).
Ranges for observed h max correspond to histogram bins reported in b. The first group are glaciers on Mount Rainier, and the next are on Mount Baker. Note that estimates correspond to different time periods, which may also lead to some discrepancies between methods.
We also compared the scaling estimates to the global-scale thickness estimates of Huss and Farinotti (Reference Huss and Farinotti2012) (data published in Farinotti and others, Reference Farinotti2019, model 1), who used an inverse method to estimate distributed thickness from surface hypsometry. Means and maxima of these thickness maps agree well with observations on Mount Baker and Mount Rainier, but substantially underestimate observed thickness for South Cascade glacier (Table 1). We assessed data from the Huss and Farinotti (Reference Huss and Farinotti2012) model on its own, to compare our method with a single self-consistent dataset rather than the ‘consensus’ (weighted mean) from Farinotti and others (Reference Farinotti2019). However, the different models in Farinotti and others (Reference Farinotti2019) were reported to give similar overall results for the Cascades. In any case, since the global estimates are not optimized for the Cascades, more-detailed glacier modeling would need to evaluate assumptions within the numerical approach against local observations. It is important to remember that thickness observations also contain sampling uncertainty and constitute snapshots during an evolving response.
There is some inherent ambiguity in what ‘characteristic’ thickness means. Jóhannesson and others (Reference Jóhannesson, Raymond and Waddington1989) used maximum thickness on a simple flowline, while the RB14 model has agreed well with numerical models when using the mean flowline thickness (RB14, Christian and others, Reference Christian, Koutnik and Roe2018). However, both Jóhannesson and others (Reference Jóhannesson, Raymond and Waddington1989) and RB14 noted that topographic variations complicate these rules of thumb. Any single thickness metric, whether derived from scaling, numerical inversions or direct observations, unavoidably loses some of the detail of the real system. For first-order estimates of τ we proceed with Eqn (6), as it focuses on a key geometric constraint (slope) and is consistent with the simplicity of our method. A comparison with the Huss and Farinotti (Reference Huss and Farinotti2012) thickness estimates is provided in the Supplementary materials.
Mass balance
We must also estimate terminus mass-balance rates (b t) to estimate response times (Eqn (1)). For any glacier, b t depends on the local mass-balance gradient, and how far the glacier extends through this gradient. Haeberli and Hoelzle (Reference Haeberli and Hoelzle1995) and Hoelzle and others (Reference Hoelzle2007) estimated b t by assuming surface mass balance is near zero at the glacier's mean elevation ($\bar {Z}$
), and then extrapolating a vertical mass balance gradient to the terminus (Z t):
An alternative is to extrapolate a horizontal gradient (db/dx) along the glacier length, which would in principle allow for non-elevation factors (e.g. topographic shading, avalanching, wind effects) to affect the mass-balance profile:
We can assume a typical gradient based on direct observations of mass balance. These exist for a handful of glaciers in the Washington Cascades, although estimates of mass-balance gradients vary widely in the published literature and data (Table 2). Observations suggest db/dz of 5–10 m w.e. a−1 km−1 on several glaciers, but much steeper on South Cascade glacier (Meier and Tangborn, Reference Meier and Tangborn1965; Meier and others, Reference Meier, Tangborn, Mayo and Post1971; Baker and others, Reference Baker2018). Also, estimates for Nisqually Glacier on Mount Rainier vary between early reports (Meier and others, Reference Meier, Tangborn, Mayo and Post1971) and recent monitoring (e.g. Riedel and Larrabee, Reference Riedel and Larrabee2015, Supplementary material, Fig. S2). Furthermore, mass-balance gradients are hard to constrain on very small glaciers (e.g. Sandalee and Noisy Creek glaciers, L ~ 0.8 and 1.5 km, respectively).
Published or derived mass-balance constraints for several glaciers with observations, and comparison with scaling estimates for b t based on either vertical or horizontal gradients
Sources and approximate observational periods are as follows:
aMeier and Tangborn (Reference Meier and Tangborn1965), 1957–74.
bMeier and others (Reference Meier, Tangborn, Mayo and Post1971), 1960s.
cBaker and others (Reference Baker2018), 1986–2019.
dRiedel and Larrabee (Reference Riedel and Larrabee2015), 2002–12.
eRiedel and Larrabee (Reference Riedel and Larrabee2018), 1993–2012.
fRiedel and Larrabee (Reference Riedel and Larrabee2016), 1993–2009.
gNational Park Service provisional data (personal communication from Jon Riedel, 2020), 1993–2017.
hTangborn and others (Reference Tangborn, Fountain and Sikonia1990), 1947–61.
iPelto and Hedlund (Reference Pelto and Hedlund2001), 1984–1998 (varying by glacier).
Note that sources c−g refer to mass-balance stake data. For these, we estimated db/dz with a linear fit (see Supplementary materials, Fig. S2). In the table, ranges are given where sources suggest different values. In the first column, italics indicate glaciers on volcanoes.
Estimating b t for each glacier based on a single characteristic gradient, and assuming the equilibrium line is at the glacier's midpoint, is thus a highly simplified approach. However, it does capture first-order differences between glaciers based on their vertical or horizontal span. Table 2 compares observations of b t, which have been reported for an additional subset of glaciers, with our scaling estimates for b t, using either Eqn (7) or (8). We find that using Eqn (8) with db/dx = 2.7 m w.e. a−1 km−1 captures the heterogeneity in observed b t slightly better than applying a vertically defined gradient. Importantly, it captures the strongly negative rates (−5 to −10 m w.e. a−1) observed on larger glaciers with a range of geometries (e.g. compare South Cascade and Nisqually glaciers). However, a vertical gradient yields a similar overall distribution for the whole sample (Supplementary material, Fig. S3). We emphasize that in either case, errors may be substantial for individual glaciers.
Before moving to the results, it is also worth noting a more general limitation for these response time estimates and application of the model. Like any linearized model, the RB14 parameters are meant to reflect some initial equilibrium configuration, and solutions will become less accurate for large changes. In particular, τ corresponds to an equilibrium geometry in theory, yet most glacier observations have occurred well after the onset of anthropogenic warming (~1880; e.g. IPCC, 2013). Dates for the RGIv6 geometries range from 1959 to 1985 in this region, and most mass-balance observations correspond to only the last few decades. For how long is an estimate of τ valid? The answer likely varies by glacier. Christian and others (Reference Christian, Koutnik and Roe2018) showed that Eqns (1) and (2) compared well with nonlinear numerical models over the course of a 200-year, 2°C warming for simple geometries. On the other hand, RB14 showed that the linear assumptions break down when the terminus traverses significant geometric variations (e.g. slope breaks). Finally, this framework ignores the feedback between surface thinning and mass balance, which tends to lengthen response times (Harrison and others, Reference Harrison, Elsberg, Echelmeyer and Krimmel2001). We expect this to be a small effect for the relatively thin, steep glaciers common to the Cascades, but the true response time might be underestimated for flatter, thicker glaciers such as South Cascade (as we will see, there are only a few of these geometries). It must be kept in mind that τ will always be an inexact metric for real glaciers undergoing large changes. Nevertheless, if focus is directed toward the basic physical tendencies that it captures, τ can be a useful metric (in conjunction with the model) for comparing the responses of glaciers with different geometries and scales.
Results
Response-time estimates
Figure 2a shows the estimated distribution of characteristic thicknesses for the 383 glaciers in our sample. The distribution peaks around 40 m, which is unsurprising for a collection of relatively small glaciers, and is consistent with previous applications in other regions (Hoelzle and others, Reference Hoelzle2007). The maximum estimate for H occurs on South Cascade Glacier (123 m), and only two other glaciers have estimates exceeding 100 m (Honeycomb and Columbia glaciers). Recall that the method appears to capture a characteristic value closer to the mean than maximum thickness (Table 1).
Estimated glacier parameters and response times. (a) Distribution of characteristic thickness using shear-stress scaling. (b) Estimated b t, using db/dx = 2.7 m w.e. a−1 km−1. (c) Estimated response times according to Eqn (1). (d) As for (c), but weighted by glacier area reported by RGIv6. (e) Map of each glacier in our sample. Dot size corresponds to area and color indicates τ. Clusters of glaciers on the volcanoes are indicated, as is South Cascade Glacier. Note many large glaciers are clustered on the volcanoes, and dots may overlap. The basemap shows topography from the Shuttle Radar Topography Mission (Jarvis and others, Reference Jarvis, Reuter, Nelson and Guevara2008).
Figure 2b shows the resulting distribution of b t for the sample. The small scale of most glaciers gives modest values for b t (i.e. their terminus is not far from their equilibrium line), while the volcanoes host a few glaciers with b t near or exceeding 10 m w.e. a−1 in magnitude. An important caveat for some of these cases is that several glaciers on Mount Rainier have extensive debris cover, which can lower melt rates (e.g. Moore and others, Reference Moore, Nelson and Groth2019). While observations do support strongly negative b t, our estimates may be too negative for debris-covered glaciers.
Combining estimates of H and b t yields estimates for τ, whose distribution is shown in Figure 2c. This distribution represents the central estimate for the region-wide scaling, and we again emphasize that uncertainties are large for individual glaciers. We discuss the impact of errors in τ in a later section. Note that we have presented mass-balance rates in water-equivalent units, as is standard for reported observations. To calculate response times, however, we convert b t to ice-equivalent units for consistency with H. Figure 2e shows the locations of individual glaciers, where the marker size indicates area and color indicates τ. τ falls between 10 and 60 years for >$90\percnt$
of glaciers. However, when the distribution is weighted by reported glacier area (Fig. 2d), it skews towardshorter response times. The skewness of the area-weighted distribution is due in large part to the volcanoes, which host clusters of large, fast-responding glaciers. Their steep slopes lead to thinner ice, and their low termini have high ablation rates, which together lead to short response times. These tendencies on H and b t are supported by observations (Tables 1 and 2), and short response times on the volcanoes have also been proposed based on 20th-century fluctuations (e.g. Harper, Reference Harper1993; Pelto and Hedlund, Reference Pelto and Hedlund2001). Our estimates of τ, while uncertain for individual glaciers, set us up to investigate what the range of response times implies for the state of glaciers throughout the Washington Cascades.
Current glacier disequilibrium
We now use our estimates of τ and the RB14 model to assess the adjustment of glaciers in the Cascades to anthropogenic warming, which began in the late 19th century (e.g. IPCC, 2013). We assume that a linear trend in mass balance, associated with warming surface air temperatures, captures the majority of the anthropogenic signal. Assuming a linear trend allows us to use the analytical solution for ‘fractional equilibration’ (Eqn (5)) from the RB14 model to analyze glacier response as a function of τ. Recall that fractional equilibration is defined with respect to an equilibrium response that is continually evolving with the climate trend (Fig. 1c). We would need additional constraints on glacier sensitivity or length history to assess absolute length disequilibrium (Eqn (4); Fig. 1b) for all glaciers. The advantage of fractional equilibration is that it depends only on τ and the duration of the climate trend, so it allows us to compare glaciers with different geometries and sensitivities. Additionally, though we approximate the forcing as a linear trend, fractional equilibration does not depend on the magnitude of the trend, or the relationship between temperature and mass balance. The following results are not a comprehensive reconstruction of glacier retreat, but a comparison of how different response times affect the component of glacier evolution associated with long-term warming.
Figure 3a shows the evolution of fractional equilibration for all 383 glaciers, assuming a linear trend beginning in 1880. Figure 3b shows the current distribution of fractional equilibration. The most important and broadly applicable result is that a fairly typical range of individual glacier response times implies a wide range in the current fractional equilibration of those glaciers. For example, τ = 10 years yields a fractional equilibration of 88% after 140 years of forcing, while τ = 40 years gives only 51%. In other words, the retreat of fast-responding glaciers, such as those on the volcanoes, has nearly kept pace with warming thus far. But those with multidecadal response times might have retreated only half of their full response, implying additional committed retreat even with no further warming.
Ratio of transient to equilibrium length response evolves according to τ. (a) Fractional equilibration of each glacier (Eqn (5)) assuming a climate trend began in 1880. (b) The distribution of L′/L′eq today, 140 years into the forcing. A typical range of response times means a large range in current equilibration. (c) Fractional equilibration for Nisqually (blue) and South Cascade glaciers (red). The shading corresponds to the range of τ consistent with available observations for each glacier (5–14 years for Nisqually; 20–41 years for South Cascade), while the darker line corresponds to the scaling estimate (5 and 25 years, respectively). (d) Fractional equilibration can be used to estimate the current equilibrium length based on observed retreat. Despite similar amounts of observed retreat over the last century (dots), Nisqually and South Cascade glaciers have very different amounts of additional committed retreat. Dashed lines show the range of committed retreat corresponding to uncertainty in τ.
Bearing in mind the limitations of the linear model discussed earlier, we can analyze two well-observed glaciers to illustrate this disparity in current disequilibrium. Both Nisqually and South Cascade glaciers have more than a century of length observations, and thickness and mass-balance measurements provide more-direct constraints on their response times. We consider a plausible range for τ corresponding to the observational estimates of mean and maximum thickness (Table 1), and published values for b t (Table 2). This gives a range of 5–14 years for Nisqually, and 20–41 years for South Cascade. The scaling method, using the region-wide parameters, estimates 5 years for Nisqually and 25 years for South Cascade. Figure 3c shows estimated fractional equilibration for each, which are distinct despite uncertainty in τ. We can use current fractional equilibration to estimate their total committed retreat, based on their retreat histories. Leclercq and others (Reference Leclercq2014) report a 2.2 km retreat over 1885–2001 for Nisqually and 1.8 km over 1900–2007 for South Cascade glacier. Taking these as L′, we can use Eqn (5) to solve for L′eq based on τ and t. For the range of τ described above, this implies 180–230 m of additional retreat committed for Nisqually in 2001, and 730–2300 m for South Cascade in 2007 (Fig. 3d). For South Cascade glacier, the upper bound would mean much of the remaining glacier is lost. This is consistent with Rasmussen and Conway (Reference Rasmussen and Conway2001), who estimated an ‘equilibrium topography’ based on the pattern of mass balance, which is a small fraction of current area. While Eqn (5) would break down as a glacier approaches complete loss, the point of this example is to demonstrate the current disequilibrium associated with multi-decadal response times vs shorter decadal response times.
Implications of short vs long response times
The simple scaling arguments based on mass balance and geometry suggest that response times in the Cascades range from less than a decade to several decades, which corresponds to a notable spread in current equilibration (Fig. 3a, b). We now illustrate some further contrasting behavior implied by this range of response times. Since our overall focus is on the population of glaciers in the Cascades, we want to highlight general consequences of short vs long response times, which would be applicable to multiple glaciers in the inventory. Therefore, instead of simulating individual glaciers (which would also require further observational constraints), we introduce two generic, synthetic glaciers to illustrate these points. We chose parameters giving τ = 12 years and 48 years for the two glaciers, in order to roughly bracket the distribution of response times in the Cascades. Additional parameters are characteristic of small maritime glaciers (described in the Appendix), but the exact values are not critical for these comparisons; the main point is to illustrate the effect of different response times. We will refer to these as the ‘fast’ and ‘slow’ synthetic glaciers.
Climate variability
First of all, the calculation of fractional equilibrium via Eqn (5) assumes a linear climate trend. However, glaciers also integrate year-to-year climate anomalies (e.g. Oerlemans, Reference Oerlemans2001; Roe and Baker, Reference Roe and Baker2014), adding variability to the overall retreat due to anthropogenic warming (Roe and others, Reference Roe, Baker and Herla2017). The glacier state at any given time reflects both the disequilibrium associated with anthropogenic forcing since ~1880 (hereafter, the forced disequilibrium) as well as with these natural fluctuations, and it is important to understand their relative magnitude.
For mass-balance anomalies b′ consistent with white noise (i.e. equal power at all frequencies and no persistence), RB14 derived the standard deviation of length anomalies σL as
where σb is the standard deviation of b′; $\psi ( \tau ) = \sqrt {[ ( 1-\kappa ) ( 1 + 4\kappa ^2 + \kappa ^4) ] /( 1 + \kappa ) ^5}$
; and κ = 1 − Δt/ετ. Here, Δt = 1 year and β, τ and ɛ are defined as above. ψ(τ) simplifies to a more-manageable $\sqrt {3 \Delta t / ( 16 \varepsilon \tau ) }$
for τ ≫ Δt. The salient aspect of ψ(τ) is that it is a decreasing function of τ, here expressing the damping caused by glacier memory. The factors in Eqn (9) show that length variability (σL) depends on glacier sensitivity (βτ), the glacier's damping effect (ψ(τ)) and the magnitude of the imposed climate variability (σb).
We illustrate these fluctuations with the idealized glaciers in Figure 4. We apply a noisy trend in mass balance of $\dot {b} = 1$
m a−1 century−1, beginning in 1880 (Fig. 4a). We add white-noise anomalies with σb = 1 m a−1, consistent with the 59-year record on South Cascade Glacier (Baker and others, Reference Baker2018). The trend thus has a signal-to-noise ratio (SNR) Δb/σb ~ 1 after one century of forcing. Recall that b′ (Fig. 4a) reflects the climate forcing alone, which is the input to the model. The glacier-averaged balance (not shown) would include the effect of changing glacier geometry. This levels off as the glacier retreats, but remains negative as long as the climate trend continues because the glacier remains in disequilibrium.
(a) Synthetic mass-balance forcing, which is a linear trend starting in 1880 (dashed line) plus white-noise anomalies (thin red line). The trend has an SNRb of ~1 after one century. Bold red shows the 30-year running mean. (b) Length responses to the anomalies in (a) for glaciers with τ = 12 years (orange) and 48 years (blue). Shading correspond to ± 1σL bounds around the response to the trend. Dashed lines show equilibrium response (without variability). (c) The ratio of disequilibrium to σL in the limit t ≫ τ. The relationship is shown for different choices of SNRb (solid vs dashed lines). Stars correspond to the glacier states in (b) in 2020. (d) April–September temperature from Rohde and others (Reference Rohde2013) and October–March precipitation from Matsuura and Willmott (Reference Matsuura and Willmott2018) for the northwest Cascades (~121° W, 48.5° N), with 30-year running means (bold). (e) Length responses of the idealized glaciers to a mass-balance trend with a 30-year break (inset). Dashed lines show the equilibrium response.
The glacier responses, simulated with the RB14 model, are shown with bold lines in Figure 4b. The shading shows the ± 1σL bounds given by Eqn (9). Climate variability makes the ‘instantaneous’ equilibrium L eq harder to define, but we illustrate it here based on the forced trend (dashed lines). For real glaciers, it can be conceptualized as the long-term mean terminus position, were the external forcing trend to stop.
It is visually clear that the disequilibrium compared to σL is different for the fast and slow glaciers, but we can generalize this using the RB14 model. The forced disequilibrium approaches $3\epsilon \tau ^2 \beta \dot {b}$
for t ≫ τ (Eqn (4)), which is a decent approximation for most response times after ~140 years of anthropogenic forcing. Dividing Eqn (4) by Eqn (9) gives the ratio of the long-term forced disequilibrium to the standard deviation of natural terminus variability. The sensitivity βτ cancels, leaving
The first ratio on the right-hand side is an increasing function of τ, and $\dot {b}/\sigma _{\rm b}$
is the ratio of the forced change in b to the natural variability (i.e. the SNR per unit time, 1 century−1 here). Figure 4c shows Eqn (10) for fixed $\dot {b}/\sigma _{\rm b}$
. The values vary significantly over the range of typical response times. Using the synthetic glaciers as examples, forced disequilibrium equals ~1σL and 9σL by 2020 (orange and blue stars, respectively). $\dot {b}/\sigma _{\rm b}$
may vary somewhat by glacier, but mass-balance anomalies are generally quite regionally coherent (Pelto, Reference Pelto2006; Huybers and Roe, Reference Huybers and Roe2009). We show a range from 0.5 to 2 century−1 for reference (dotted lines).
The point is that the range of response times in the Cascades includes glaciers where long-term, forced disequilibrium may be obscured by natural fluctuations, as well as glaciers where disequilibrium stands well beyond the noise. This helps put the differences in fractional equilibration (Fig. 3) in a more practical light. One may ask: if climate change paused today, would we notice the additional committed retreat in the coming decades? The answer is a clear yes for multi-decadal response times, but a careful accounting of climate variability would be needed for short response times.
Our calculation of fractional equilibration also assumes that the external forcing is described by a linear trend, but global-mean temperature records show a slowdown in warming roughly from the 1940s to 1970s, and faster warming thereafter (e.g. IPCC, 2013). Local climate is noisier, but data from the Pacific Northwest show similar changes in the rate of warming, and ample precipitation (Fig. 4c). Regardless of the natural vs anthropogenic contributions, the mid-20th century clearly included a period relatively favorable for glacier mass balance in the Cascades (Rasmussen, Reference Rasmussen2009), and several glaciers briefly advanced (e.g. Hubley, Reference Hubley1956; Harper, Reference Harper1993; Pelto and Hedlund, Reference Pelto and Hedlund2001).
It is important to note that variations in decadal trends occur even in uncorrelated white noise (see the running means in Fig. 4a) which can cause glaciers with short memories to advance (Fig. 4b). However, distinct phases of forcing are often invoked heuristically in the literature, and so we consider what this assumption implies for different glaciers. To illustrate the effects on disequilibrium, we apply a forcing broken into two linear trends separated by a constant climate from 1940 to 1970 (Fig. 4d). A 30-year break in forcing is enough time for the glacier with τ = 12 years to nearly equilibrate, while the glacier with τ = 48 years remains far out of equilibrium, and its retreat rate is essentially unchanged. While we are not proposing that a step-wise trend is the optimal model, Figure 4d simply illustrates that τ limits how much glaciers could have adjusted to multi-decadal changes in the rate of forcing.
Uncertainties in current equilibration
Interannual variability and multidecadal changes in the rate of forcing, both of which are relevant to glaciers in the Cascades, imply potential errors in estimates of equilibration from Eqn (5), which is the solution for a linear trend. Uncertainty in τ itself is obviously another source of uncertainty, but is associated with the glacier rather than the forcing. In Figure 5, we use the idealized glaciers to compare these different sources of uncertainty.
Errors in estimated fractional equilibration (L′/L′eq) from three different sources. (a) Climate variability drives natural glacier fluctuations, which temporarily drive a glacier toward or away from its long-term equilibrium length. Shaded bounds correspond to L′ ± 1σL/L′eq for glaciers with τ = 12 years (orange) and 48 years (blue). (b) L′/L′eq for the case with a forcing break from 1940 to 1970. The difference compared to a linear forcing (dotted lines) is significant during the break, but minimal after a few decades of resumed forcing. (c) The spread in L′/L′eq for a range of $\pm 50\percnt$
in τ. Uncertainty in the response time means persistent uncertainty in fractional equilibration.
First, if we assume L′eq is described by a linear trend in the long term but allow for a spread in L′ due to interannual variability, the spread in fractional equilibration is (L′ ± σL)/L′eq. Figure 5a shows this spread (shaded bounds), as well as L′/L′eq for the particular realization of noise in Figs 4a, b. σL swamps L′ ± σL/L′eq early on, when L′ and L′eq are small. As expected from Figure 4, the uncertainty in current fractional equilibration is greater for short response times, because forced disequilibrium can be overwhelmed by natural terminus variations.
Next, Figure 5b compares L′/L′eq for the trend with a break in forcing (solid lines) vs a linear trend (dashed lines). Fractional equilibration increases more rapidly during the break (L′eq stops changing), but tends back toward the linear case after a few decades. While the glacier with τ = 12 years nearly equilibrates during the break, it regains its (small) level of forced disequilibrium quickly. The result for both glaciers is that, provided a break in forcing was several decades ago, the effects on current fractional equilibration are small. L′/L′eq is simply dominated by the total forcing and response since 1880.
Finally, uncertainty in τ leads to a persistent uncertainty in disequilibrium and fractional equilibration (Christian and others, Reference Christian, Koutnik and Roe2018). Available observations (Tables 1 and 2) suggest that for individual glaciers, a substantial uncertainty in τ should be considered. Even for those with direct observations, the ambiguity in defining a ‘characteristic’ thickness implies a conservative range from $\bar {h}$
to h max, which is often a factor of 2 or more (Table 1). Errors in b t may be similar without direct observations (especially for small b t; Table 2). Thus, Figure 5c shows a very broad spread in τ of $\pm 50\percnt$
(i.e. a factor of 3 between upper and lower bounds). In contrast to errors from neglecting short-term fluctuations, uncertainty in τ is more consequential for long response times.
Figure 5 shows three qualitatively different uncertainties associated with our approach to estimating the current fractional equilibration of glaciers in the Cascades. These synthetic examples suggest that simplifying the external forcing to a linear trend is a reasonable approach for estimating the fractional equilibration to the total anthropogenic forcing thus far. The caveat, of course, is that forced disequilibrium simply never emerges far from the noise for glaciers with very short response times (Fig. 4c). For glaciers with longer response times, the more serious issue for estimating disequilibrium with current climate is uncertainty in individual response times.
Runoff changes
Finally, we consider how fast vs slow glacier responses affect their contributions to downstream hydrology, which are important for some drainages in the Cascades, especially in the late summer (e.g. Riedel and Larrabee, Reference Riedel and Larrabee2016). A common expectation is that climate warming causes glacier runoff to increase initially, and later decline as glacier area diminishes, leading to the concept of ‘peak runoff’ in glacierized watersheds (e.g. Jansson and others, Reference Jansson, Hock and Schneider2003). This phenomenon depends fundamentally on transient glacier dynamics. Carnahan and others (Reference Carnahan, Amundson and Hood2019) showed that, over a wide range of parameters, the peak in simulated glacier-melt runoff occurred ~1τ after the onset of a climate trend, using a metric for τ very similar to Eqn (1) (see Harrison and others, Reference Harrison, Elsberg, Echelmeyer and Krimmel2001).
Mountain hydrology also depends on snowpack, precipitation and other local factors, which can affect the timing of peaks in total streamflow. As a local example, Frans and others (Reference Frans, Istanbulluoglu, Lettenmaier, Fountain and Riedel2018) simulated a complex picture of recent and future changes for several glacierized basins in the Pacific Northwest, with runoff trajectories varying substantially by basin. We do not address all of these factors here, but focus on the glacier-melt contribution to runoff. This ignores trends in precipitation (which are much less robust than changes in temperature), and other processes such as refreezing or evapotranspiration. Direct glacier melt is the runoff component most dependent on glacier dynamics, via transient changes in melt area as the glacier advances or retreats. And, in the Cascades, it can play an important role in late summer, when precipitation and snowmelt make minimal contributions (e.g. Riedel and Larrabee, Reference Riedel and Larrabee2016; Frans and others, Reference Frans, Istanbulluoglu, Lettenmaier, Fountain and Riedel2018).
By simply calculating the total melt associated with the glacier's geometry in the RB14 model, we can illustrate differences between the fast and slow glaciers that would be relevant for assessing hydrological contributions from different glaciers in the Cascades. We take the total glacier melt flux M g (m3 a−1) to be equivalent to the local melt rate m (m a−1), integrated over the glacier's surface area. If we assume a constant width w, this is:
where x is the coordinate from the glacier head (x = 0) to terminus (x = L). We assume a linear melt gradient based on a standard temperature-melt factor, lapse rate and glacier slope (see the Appendix), and we also assume that the melt anomaly driven by warming is spatially uniform across the glacier. These assumptions simplify the integral to the mean value between melt at the glacier head and terminus, at any time t. The melt flux is then
Because long-term mass-balance trends and glacier retreat are dominated by warming rather than precipitation trends (e.g. Medwedeff and Roe, Reference Medwedeff and Roe2017), we assume that specific melt anomalies are equivalent to the mass-balance anomalies (b′) that drive glacier retreat in the model. Note that this ignores changes in the proportion of rain in the precipitation on the glacier.
Figure 6a shows melt-flux changes following a trend in mass balance (and therefore specific melt), with L(t) given by the RB14 model for the idealized glaciers. M g peaks at t ~ τ for both glaciers (peaks are at 15 and 49 years), consistent with Carnahan and others (Reference Carnahan, Amundson and Hood2019). The role of the response time is conceptually straightforward: for long τ, a slower initial retreat means specific melt anomalies can increase more before much melt area is lost, and the peak in M g is thus later than for short τ. Note that M g drops steeply after this peak, as area losses become large.
Melt changes for idealized glaciers forced by a linear warming trend. (a) Top panel shows local melt anomaly, and bottom shows integrated melt (M g) as the glacier responds. Solid lines show the actual transient melt, and dotted lines show the melt that would occur over the instantaneous equilibrium geometry. The timing and magnitude of peak melt depends on τ. (b) As for (a), but with a break in forcing from 1940 to 1970. A long response time means that a second peak is still affected by earlier warming.
This basic mechanism causing M g to peak and decline is straightforward at the onset of an idealized trend, but variability in the rate of climate forcing introduces additional transient changes. Figure 6b shows M g in response to the two-step trend considered earlier, which, again, is an idealized representation of the multidecadal rates of warming in the Pacific Northwest (Figs 4d, e). In this scenario, there are some additional differences between the fast- and slow-responding glaciers. Because the slow glacier (τ = 48 years) remains far from equilibrium during the 30-year break in warming (Fig. 4e), its melt flux during the second stage of warming is compounded by its continued adjustment to earlier warming. In this case, there is no second peak in M g, because the loss of melt area is dominant and nearly uninterrupted. In contrast, the fast-responding glacier yields two similar peaks to both phases of warming, since it can keep up with multi-decadal variations in forcing.
Finally, the committed retreat associated with glacier-length disequilibrium (e.g. Fig. 3d) means a committed loss of ablation area and thus of melt contributions to runoff. The dashed lines in Figure 6 show the melt flux that would occur integrated over the equilibrium glacier length, illustrating the total committed change as the forcing evolves.
It should be noted that multiple runoff metrics exist (e.g. O'Neel, Reference O'Neel, Hood, Arendt and Sass2014), and the timing and magnitude of peaks depend on the metric chosen. For instance, total runoff from a fixed basin area peaks later than that from the evolving glacier surface (e.g. Carnahan and others, Reference Carnahan, Amundson and Hood2019; Rounce and others, Reference Rounce, Hock and Shean2020). The most appropriate metric depends on the question at hand, and obviously other variables are important for overall basin hydrology. However, the glacier melt contribution (M g) is useful to analyze on its own as well, since it directly tracks the evolving area of glacier ice that provides runoff once seasonal snowpack has waned. It is also at these times when the committed retreat of slow-responding glaciers is most consequential for committed hydrological changes.
Discussion
As noted previously, a number of studies have analyzed the transient terminus behavior of glaciers in the Cascades (e.g. Hubley, Reference Hubley1956; Harper, Reference Harper1993; Pelto and Hedlund, Reference Pelto and Hedlund2001; Rasmussen and Conway, Reference Rasmussen and Conway2001; Roe and O'Neal, Reference Roe and O'Neal2009; Stevens and others, Reference Stevens, Conway, Kennard, Rasmussen and Koutnik2018). In particular, Pelto and Hedlund (Reference Pelto and Hedlund2001) estimated response times for 21 glaciers in the Cascades using multiple metrics including Eqn (1), although for unmeasured glaciers they assumed H was either 50 or 100 m. They linked these estimates to differences in terminus responses, including the failure of some glaciers to advance in the 1950s–1970s. Additionally, Rasmussen and Conway (Reference Rasmussen and Conway2001) and Harrison and others (Reference Harrison, Elsberg, Echelmeyer and Krimmel2001) pointed out the apparently long response time of South Cascade glacier, and other studies have inferred short response times for glaciers on Mount Baker (Harper, Reference Harper1993; Roe and O'Neal, Reference Roe and O'Neal2009).
Our conclusions are qualitatively consistent with these studies, namely that different behaviors, linked to the response time, should be expected throughout the range. Our method offers an advance by taking advantage of data from RGIv6 to estimate τ for more glaciers, and also applying the RB14 model to explore additional implications of these response times. While individual response times may be uncertain, the model provides a formalism for understanding implications of the range of response times found in the Cascades.
One such implication is that the dynamic state of one glacier may be substantially different from that of its neighbors, even among glaciers with response times in the typical range of one to a few decades. For example, a pause in retreat (or advance) for one glacier may indicate that it is near equilibrium with local climate, but other glaciers with longer response times may stay far from equilibrium during the same period (Fig. 4). We have illustrated this with a simple linear model that uses τ explicitly, but the point holds for any model meant to capture transient glacier dynamics. This should be considered when calibrating models that simulate many glaciers. If such models are initialized well after the onset of industrial-era warming, differences in initial disequilibrium would be important for projecting accurate responses to additional forcing. This also applies to modeling a peak in glacier runoff, which depends on the lag between warming and retreat. These issues should be considered as hydrological studies incorporate increasingly complex ice dynamics. For example, distributed ice-flow models require a spin-up period to capture the current transient state and disequilibrium of glaciers, but assumptions built in to model initialization vary in the literature (e.g. Naz and others, Reference Naz, Frans, Clarke, Burns and Lettenmaier2014; Clarke and others, Reference Clarke, Jarosch, Anslow, Radić and Menounos2015; Frans and others, Reference Frans, Istanbulluoglu, Lettenmaier, Fountain and Riedel2018). Estimates of response times, as we have used here, could help ensure that spin-ups are commensurate with characteristics of the glaciers being simulated.
A final point relates to estimating parameters for glaciers with minimal direct observations. Our method simplified this to H and b t, but these encapsulate key parameters for any glacier model. There has been substantial progress in estimating ice thicknesses around the world (e.g. Farinotti and others, Reference Farinotti2019), but it is important to remember that mass-balance gradients are similarly important for a glacier's transient response and overall length sensitivity. This is implicit in Eqn (1) (via b t), but also affects the response of numerical models that don't rely on a stipulated response time. Yet at the same time, modeled equilibrium length and thickness are not highly sensitive to the mass-balance gradient, provided that the overall mass budget is roughly correct. This point has been made before (Oerlemans, Reference Oerlemans2001), but is worth reemphasizing, as it highlights a challenge for tuning and initializing models: mass-balance gradients – and thus important aspects of the future response – are not necessarily constrained by matching the glacier-wide mass budget or an observed geometry. Fortunately, existing observations indicate that a glacier's mass-balance gradient is fairly constant from year to year, despite variability in glacier-wide balance (e.g. Oerlemans, Reference Oerlemans2001). Progress in constraining gradients on high-priority glaciers, and the heterogeneity within a region, might thus be possible from new observations even without long-term records. Accounting for specific sources of heterogeneity in mass-balance gradients, such as variable lapse rates (Minder and others, Reference Minder, Mote and Lundquist2010) and local topography (Florentine and others, Reference Florentine, Harper, Fagre, Moore and Peitzsch2018), could also help improve extrapolations to glaciers that lack direct observations. Recent advances in calculating geodetic changes in seasonal timescales (e.g. Bushan and others, Reference Bhushan, Shean, Alexandrov and Henderson2021; Hugonnet and others, Reference Hugonnet2021) may also provide new tools to constrain the patterns of mass balance, and thus important glacier-response characteristics at regional scales.
Conclusions
The strong geometric controls on a glacier's response time (τ = −H/b t) mean that first-order estimates can be made from basic parameters available in inventories like RGIv6. We have adapted established scaling methods to 383 glaciers in the Washington Cascades. Uncertainties may be large for individual glaciers, but our main goal was to perform an overall survey of the region. When calibrated to the few available thickness and mass-balance observations from the Cascades, the scaling yields a distribution of decadal-to-multidecadal response times (Fig. 2c), which is generally comparable to results for other mountain ranges (Haeberli and Hoelzle, Reference Haeberli and Hoelzle1995; Hoelzle and others, Reference Hoelzle2007; Zekollari and others, Reference Zekollari, Huss and Farinotti2020). A result specific to the Cascades is that the largest glaciers, which reside on the major stratovolcanoes, tend to have short response times (~10 years).
Differentiating between fast and slow response times is important in the context of century-scale anthropogenic warming trends. We analyzed three issues where a typical range of response times implies markedly different interpretations among individual glaciers.
First, we analyzed the current fractional equilibration of glaciers, which is the ratio of their transient response (i.e. observed retreat) to the full equilibrium response to a climate trend (Eqn (5)). After a warming trend of ~140 years, there is a wide spread in this fractional equilibration among individual glaciers the Cascades (Fig. 3b).
Second, a glacier's response to natural climate variability is important context for this disequilibrium, and also depends on the response time. For glaciers with long τ, the disequilibrium with current climate far outweighs fluctuations due to interannual climate variability. For short-τ glaciers, noise-driven fluctuations can be similar in magnitude to long-term disequilibrium. Interannual climate variability would thus be more important for interpreting the current state, or making near-term projections, of these fast-responding glaciers. In the Cascades, this is an important consideration for the volcano glaciers.
Finally, a range of response times has implications for changes in the glacier-melt contribution to runoff. We used a simplified model to account for the leading effects of increased melt and decreasing glacier area following a warming trend. After a warming trend exceeds a few multiples of τ (which is the case for the vast majority of Cascades glaciers), the loss of melt area from cumulative retreat provides a strong negative tendency on total melt that may overwhelm transient gains from increased local melt rates. This glacier-dynamics perspective is consistent with empirical assessments, which have linked declining glacier area to declining late-summer melt contributions in the Cascades’ Skagit watershed (Riedel and Larrabee, Reference Riedel and Larrabee2016), and the nearby upper-Columbia River basin (Moore and others, Reference Moore, Pelto, Menounos and Hutchinson2020). Frans and others (Reference Frans, Istanbulluoglu, Lettenmaier, Fountain and Riedel2018) projected declining late-summer glacier runoff in the coming decades in several Pacific Northwest drainages, although they also noted the possibility of transient gains in high-elevation basins from increased precipitation, which we did not consider here. The concept of ‘peak runoff’ in response to warming is physically robust, but quickly becomes nuanced when climate variability and a population of glaciers are considered. Our main point is that transient ice dynamics play a key role, and can lead to differences in the melt contributions of individual glaciers.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2021.133.
Acknowledgments
J.E.C. was supported by a National Science Foundation Graduate Research Fellowship (DGE-1256082), a U. Washington Innovation Award (PI: D. Shean) and an internship at the USGS Washington Water Science Center in Tacoma, WA through the NSF Graduate Research Internship Program (GRIP). Coordination of GRIP at USGS is through the Youth and Education in Science programs within the Office of Science Quality and Integrity. Conclusions are those of the authors. E.C. was supported by fellowships from the U. Texas Provost Office and Jackson School of Geosciences. We thank Andrew Fountain, Shad O'Neel and David Shean for helpful insights, and Jon Riedel and Michael Larrabee for providing stake data for glaciers monitored by the National Park Service. We are grateful to many individuals from USGS, NPS and Nichols College for decades of mass-balance monitoring in the Cascades. Finally, we thank two anonymous reviewers for detailed and insightful comments, and Scientific Editor Hester Jiskoot for helpful and effective handling of the manuscript. Code to reproduce these analyses is available as a public repository at https://github.com/johnerich/cascades-responsetimes.
Appendix: Synthetic glacier parameters
Table 3 lists parameters used for the synthetic glaciers discussed in the main text. To calculate length anomalies, the RB14 model (Eqn (2)) only requires β and τ directly (in addition to mass-balance forcing). Here, we present the geometric and climate parameters that gave us β and τ for these cases (see RB14 for a more general description). We chose values roughly representative of small maritime glaciers, but note that the same response times can be achieved with different climatic and geometric parameters. The main point here was to compare synthetic glaciers with different response times. Note, however, that the runoff calculations use a number of these parameters explicitly (e.g. L 0, w), so they have a direct role beyond their effects on β and τ.
Geometric and climate parameters for the synthetic glaciers
In RB14, β = A tot/(wH), where A tot is the total glacier area, w is the width near the terminus and H is the characteristic thickness. We assume a constant width w, which simplifies the expression to β = L 0/H, where L 0 is the steady-state length. A melt factor μ relates melt-season temperature to melt rate, and a standard atmospheric lapse rate Γ describes the vertical temperature gradient. The vertical mass-balance gradient is thus db/dz ≡ μΓ. Given melt-season temperature at the top of the glacier (T Zmax) and solid precipitation P, the mass balance at the terminus is b t = P − μ(T Zmax + ΓL 0 tan θ), where θ is the average slope. Combined with H, this defines the response time τ = −H/b t.
Open access
