Glacier mass balance provides a direct indicator of a glacier’s relationship with local climate, but internally generated variability in atmospheric circulation adds a significant degree of noise to mass-balance time series, making it difficult to correctly identify and interpret trends. This study applies “dynamical adjustment” to seasonal mass-balance records to identify and remove the component of variance in these time series that is associated with large-scale circulation fluctuations (dynamical adjustment refers here to a statistical method and not a glacier’s dynamical response to climate). Mass-balance records are investigated for three glaciers: Wolverine and Gulkana in Alaska and South Cascade in Washington. North Pacific sea level pressure and sea surface temperature fields perform comparably as predictors, each explaining 50%–60% of variance in winter balance and 25%–35% in summer balance for South Cascade and Wolverine Glaciers. Gulkana Glacier, located farther inland, is less closely linked to North Pacific climate variability, with the predictors explaining roughly 30% of variance in winter and summer balance. To investigate the degree to which this variability affects trends, adjusted mass-balance time series are compared to those in the raw data, with common results for all three glaciers; winter balance trends are not significant initially and do not gain robust significance after adjustment despite the large amount of circulation-related variability. However, the raw summer balance data have statistically significant negative trends that remain after dynamical adjustment. This indicates that these trends of increasing ablation in recent decades are not due to circulation anomalies and are consistent with anthropogenic warming.
a. Glaciers and climate variability
Variations in climate occur in response to both external forcing and internally generated variability. External forcings are generally defined as mechanisms outside of the climate system that change the underlying radiative balance of the planet; such forcings can be natural (e.g., changes in volcanic or solar activity) or anthropogenic (e.g., greenhouse gas and aerosol emissions and changes in land use) in origin (e.g., IPCC 2013). Internal variability is also fundamental to the climate system and occurs even in the absence of external forcing; these variations arise as chaotic fluctuations in oceanic and atmospheric circulation and are integrated by components of the climate system operating on a range of time scales (e.g., IPCC 2013). Though such variability is essentially stochastic in nature (e.g., Hasselmann 1976), there are preferred modes of variability that emerge on interannual to decadal time scales. A large body of research exists on these persistent patterns and their origins in the coupled atmosphere–ocean system (e.g., Deser et al. 2010, and references therein) as well as on the effects of time-varying atmospheric dynamics on hemisphere-scale temperature trends (Wallace et al. 1995). Model developments and increases in computational power over the past decade have allowed for improved understanding of how this internal variability differs from external forcing. For example, Deser et al. (2012) integrated an ensemble of identical global climate models with identical external forcing scenarios, but with minor perturbations in the initial climate state between ensemble members. The resulting divergence among the model integrations can be due only to internal variability. Deser et al. (2012) demonstrated that internal variability can dominate over external forcing in regional climate trends even on multidecadal time scales (see also Wallace et al. 2015).
Glaciers interact with their local climate via their mass balance: the accumulation and ablation (i.e., melt) of snow and ice each year. Worldwide glacier retreat is a highly visible and widely cited result of a changing climate. However, glaciers respond not only to external forcings but also to interannual variability in temperature and precipitation (e.g., Oerlemans 2001; Roe and Baker 2014) associated with internal climate variability. Large year-to-year fluctuations can make it difficult to identify and attribute trends in the global archive of mass-balance records. However, when combined with other climate data, the variability itself can yield information about the large-scale climate processes that drive glacier changes. For example, Bitz and Battisti (1999) examined correlations between the mass-balance records of several western North American glaciers and indices of prominent modes of climate variability, as well as meteorological data from local to synoptic scales. Their analyses provide insight into the dynamical links between North Pacific climate and the targeted glaciers as well as the relevant differences between the glaciers’ climatic settings. In this study, we build on these analyses of the signatures of large-scale climate variability in glacier mass-balance records, with the added objective of improving the identification of trends in mass balance.
b. Dynamical adjustment
For any climate record that spans only a few decades, identifying the effect of anthropogenic forcing can be problematic; trend estimates may lack statistical significance owing to the background noise of interannual variability. Furthermore, the trends themselves may be biased by limited temporal or spatial sampling of low-frequency, internally generated fluctuations (e.g., Casola et al. 2009; Wallace et al. 2012). “Dynamical adjustment” is a method that seeks to extract the component of the variance in a climate time series that is attributable to large-scale atmospheric circulation anomalies (i.e., dynamically induced variability) rather than to external forcing. With this component removed, the adjusted time series has less variance and may exhibit a different trend. For example, Wallace et al. (2012) dynamically adjusted observed surface air temperatures poleward of 40°N and concluded that 0.7°C of the ~1.7°C warming trend in wintertime temperatures from 1965 to 2000 was dynamically induced. With adjusted trends unbiased by internal effects, such analyses can clarify the anthropogenic signal in a target climate time series.
A number of studies in recent decades have developed and used this method. Although not all were referred to as dynamical adjustment, they nonetheless established the objective of evaluating the role of dynamically induced variability. These studies analyzed surface-temperature records using a variety of approaches, including regression of the temperature fields themselves (Wallace et al. 1995), with established climate indices (e.g., Hurrell 1996; Bitz and Battisti 1999; Thompson et al. 2000) and with sea level pressure (SLP) fields (Thompson et al. 2009). More recent studies have adjusted additional climate variables such as snowpack, air temperature, and hurricane activity (Smoliak et al. 2010; Wallace et al. 2012; Smoliak et al. 2015) and are based on the method of partial least squares (PLS) regression (see Abdi 2010). PLS regression decomposes a set of predictor variables into components of variability that are optimized to explain the variability in a variable of interest (the predictand). In the context of dynamical adjustment, the predictors are grid points of a time-varying field [e.g., sea level pressure and sea surface temperature (SST)] that is assumed to have a dynamical link to the predictand (e.g., regional temperature, snowpack). We present the PLS regression algorithm in section 3.
Note that in our context, the term “dynamical adjustment” should not be confused with the dynamics of a glacier’s geometrical adjustment to climate changes. In this study, dynamical adjustment will refer only to the methodology described above. We apply a PLS-based dynamical adjustment to glacier mass-balance records using SLP and SST, independently, as predictors. SLP has been used previously for dynamical adjustments (Smoliak et al. 2010; Wallace et al. 2012; Smoliak et al. 2015), and its variability is a widely used indicator of circulation changes on a range of time scales. SLP patterns control the direction and magnitude of near-surface winds and thus are strongly linked to variability in surface temperature and precipitation (e.g., Wallace and Hobbs 2006). This relationship makes SLP a clear candidate for dynamically adjusting glacier mass-balance records.
SST is related to glacier change through a number of pathways. It is established that the onshore flow of warm, moist marine air masses during storms links SSTs to snowpack in coastal mountains (Casola et al. 2009) and, by extension in this study, to the accumulation/winter balance of glaciers in western North America. Furthermore, SST and SLP variability are closely related via dynamical coupling between the ocean and atmosphere. For example, Deser and Phillips (2009) showed that SST variability is related to recent decadal trends in North Pacific atmospheric circulation. At the same time, SST also responds to higher frequency atmospheric pressure variations, as demonstrated by Johnstone and Mantua (2014), who found that the leading mode of monthly variability in North Pacific SST resembles a lagged response to the 11-month running mean of SLP variability. The coupling of atmospheric pressure and ocean temperature on time scales from months to decades motivates our investigation of both SLP and SST datasets as predictors for glacier mass-balance variability.
2. Study area and datasets
a. Target glaciers
We focus our study on the mass-balance records of three U.S. Geological Survey (USGS) benchmark glaciers: Wolverine Glacier in Alaska’s Kenai Range, Gulkana Glacier in the Alaska Range, and South Cascade Glacier in Washington State’s Cascade Range (Fig. 1). These glaciers have the longest continuous mass-balance records in North America, with winter balance Bw, summer balance Bs, and annual balance Ba data available from 1959 to 2011 for South Cascade Glacier (World Glacier Monitoring Service 2012, 2013) and from 1966 to 2015 for Wolverine and Gulkana Glaciers (Fig. 1) (O’Neel et al. 2016). Monitoring is ongoing, but more recent measurements for South Cascade were not released at the time of our analysis.
In addition to having long-term, high-quality mass-balance records, the glaciers exist in distinct climate settings. Accordingly, their balance records show some characteristic differences: Wolverine Glacier receives ample moisture from the Gulf of Alaska (approximately 50 km away) and has a large mean-winter balance rate [Bw = 2.2 meter water equivalent per year (mwe yr−1)] and also large variability (standard deviation: σ = 0.9 mwe yr−1); Gulkana Glacier is ~300 km inland and blocked from much of this precipitation by high coastal mountains (e.g., Rasmussen and Conway 2004) and thus experiences a continental climate with less precipitation and less variability (Bw = 1.3 mwe yr−1, σ = 0.3 mwe yr−1). Meanwhile at lower latitudes, South Cascade Glacier is 250 km inland from the Pacific but is only partially blocked from onshore moisture flow by the Olympic Mountains and still resides in a maritime climate as evidenced by its high winter accumulation and variability (Bw = 2.8 mwe yr−1, σ = 0.6 mwe yr−1).
The long-term mean mass-balance values are −0.4, −0.5, and −0.6 mwe yr−1 for Wolverine, Gulkana, and South Cascade, respectively, indicating that all three glaciers have been out of equilibrium with the average climate over the study period. However, as stated previously, Wolverine and South Cascade show significant variability in annual balance (σ = 1.2 and 1.0 mwe yr−1), such that years of positive annual balance are not uncommon. For Gulkana, the mean annual balance is comparable in magnitude to the standard deviation (σ = 0.6 mwe yr−1), and indeed there are only 7 years of positive balance in the 50-yr record, with the most recent in 2003.
b. Mass-balance data
The glacier-averaged balance values reported by the USGS are calculated from point measurements of accumulation and ablation, which are converted to meters-water-equivalent based on density measurements and then extrapolated over the whole glacier area using empirical altitude-dependent relations. Wolverine and Gulkana each have three index sites at which point mass balance is measured, though measurements at additional sites have been made intermittently to constrain results. Additionally, digital elevation models (DEMs) from aerial photogrammetry are used to update the area–altitude distribution used for extrapolation from point balance to glacierwide balance as the glacier’s geometry changes, linearly interpolating between years with available DEMs (Van Beusekom et al. 2010; O’Neel et al. 2014). South Cascade currently uses six fixed index sites (Bidlake et al. 2010). Snow depth, which is often more easily measured than ablation, is frequently measured in additional locations on the glaciers in order to mitigate the effects of local irregularities. Maximum accumulation and ablation values for each index site are often derived quantities because measurements rarely fall precisely at the transition between accumulation and ablation seasons, and furthermore, these dates may not be synchronous across the glacier. For all three glaciers, temperature and precipitation measurements from nearby meteorological stations are used to estimate the dates of these transitions and to model the accumulation/ablation that occurred between the measurement date and the seasonal transition.
Complete descriptions of the models and conventions used to generate the mass-balance datasets are available in USGS technical reports—for example, Bidlake et al. (2010) for South Cascade Glacier and Van Beusekom et al. (2010) for Wolverine and Gulkana Glaciers. For our purposes, we use Bw, Bs, and Ba in their available forms as estimates of glacierwide change from season to season. We discuss the potential effects of observational error in section 5.
c. Sea level pressure
The SLP predictor is derived from 2.5° × 2.5° monthly mean SLP data from the NCEP–NCAR reanalysis (Kalnay et al. 1996). To match the respective seasons of the mass-balance time series, averages of winter (October–March) and summer (April–September) SLP are used. The domain is restricted to 20°–70°N and 150°–250°E. Our results are not qualitatively sensitive to the domain, provided that most of the North Pacific is included. However, the fraction of mass-balance variability explained by the predictor declines somewhat with significantly larger or smaller domains.
d. Sea surface temperature
For the SST predictor, 1° × 1° monthly values come from the Met Office Hadley Centre Sea Ice and Sea Surface Temperature dataset (HadISST1; Rayner et al. 2003). The global mean is subtracted from the SST field at each time, and we use the resulting anomaly fields as the SST predictor. This removes the global, long-term warming signal associated with external radiative forcing (e.g., IPCC 2013) and in principle retains trends unique to the study domain, though this residual is small compared to the global mean trend. The domain is the same as that for SLP, except it is further restricted by removing all grid points that experience sea ice cover, even if only seasonally. Seasonal averages are taken as with SLP.
3. PLS regression method
Using a time-varying field (i.e., multiple predictor time series) to explain variability in the predictand preserves spatial information, but since the individual predictors are spatially indexed observations of the same climate variable, they may be highly correlated. While this could be problematic for some regression techniques, the PLS method eliminates this problem by decomposing the predictor set into orthogonal patterns of variability.
This idea is common to a number of statistical approaches for analyzing variability in climate data. For example, principal component analysis (PCA) decomposes a single dataset into a set of patterns that best capture its own variability, and maximum covariance analysis (MCA; see, e.g., Bretherton et al. 1992; von Storch and Zwiers 1999) identifies patterns that maximize covariance between two different variables. Like MCA, PLS regression identifies common patterns between different climate variables, but here the constraint is one way: the predictor is decomposed into patterns that explain the maximum amount of variance in the predictand. This makes PLS especially suitable for dynamical adjustment, where the primary goal is to identify the large-scale origins of variability in the predictand.
Dynamical adjustment using PLS regression proceeds as follows. Let be an n × m matrix with n observations in time at m spatial points (i.e., m predictors), which has been standardized to zero mean and unit variance in time. In the case of a gridded dataset such as SST or SLP, the grid is reshaped into a 1 × m vector for each observation time (where m is the product of the original grid’s length and width), yielding the n × m matrix for the whole observational period. Let Y be the predictand time series of n observations, also standardized to zero mean and unit variance. In our case, Y is an n × 1 vector, though PLS regression can be generalized to the case that the predictand is a matrix. First, a correlation map W is created using detrended time series:
where the prime denotes that the predictand and each predictor grid point have been detrended in time. This is done so that the correlation map is not biased by trends, which could more easily be correlated without a dynamical link. The quantity W is thus an m × 1 vector that can be reshaped back into the predictor’s physical (grid) dimensions to provide a map of the detrended correlations between predictor and predictand at each observation point (see Fig. 2a). The vector W is weighted by the cosine of latitude (denoted Wc) to equally distribute influence by area and finally normalized to unit magnitude. Brackets 〈A〉 hereinafter denote the following normalization for a generic vector A:
Next, the predictor is projected onto 〈Wc〉:
where t is an index in time (n × 1) that expresses the variations of the spatial correlation patterns. Here we note the similarity with PCA, in that t is analogous to the principal component of an empirical orthogonal function of , where here we simply enforce that the function is the correlation map Wc. The index t then determines regression coefficients that are used to remove variance from both the predictor and the predictand. Let P be the vector of regression coefficients for the set of predictor variables :
Let β be the regression coefficient for the predictand Y:
Then, the index is subtracted from the predictor and predictand weighted by P and β, which yields the dynamically adjusted variables:
The fractions of variance in and Y explained by t are given, respectively, by the following:
Note that, in general, the respective fractions of variance explained in and Y are not equal since PLS finds patterns that maximize variance in Y explained by .
At this stage, the PLS regression has removed from the component of its own variability that explains the most variance in Y. And from Y, it has removed the component of variability that can be attributed to variability in . In our context of PLS-based dynamical adjustment, the pattern regressed out of is interpreted as a mode of large-scale climate variability that drives changes in glacier mass balance. However, there may in fact be multiple independent modes whereby the predictor fields drive the predictand, and thus more associated covariability may remain between and Yadj. Much of the utility of PLS regression comes from the capability to reiterate the regression with and Yadj as predictor and predictand. Further iterations yield a series of indices t1, t2, …, tn. These indices are mutually orthogonal, and thus the total variance explained is additive.
Successive iterations typically explain progressively less variance and will eventually fail to yield significant or physically meaningful adjustments. Previous studies (Smoliak et al. 2010, 2015; Wallace et al. 2012) have used cross validation to determine the number of PLS modes that should be considered in analysis. We follow in this vein, using a procedure outlined in Abdi (2010), which evaluates the predictive skill of the PLS modes when applied to an observation left out of the regression. This method indicated that between one and four modes had predictive power depending on the variables used. However, this metric—especially with relatively short records—is not guaranteed to inform physical significance. Thus, for all dynamical adjustments presented here, we consider the two leading modes for the sake of consistency. This involves some risk of retaining insignificant modes (or omitting significant modes), but since the modes beyond the leading pattern explain a small amount of variance, they are unlikely to project strongly onto mass-balance trends and significantly bias our conclusions.
a. Raw data
We begin by characterizing the noise and trends in the raw mass-balance records using standard statistical metrics. We use a straightforward version of the t statistic to evaluate trend significance (see, e.g., Casola et al. 2009; Roe 2011). For a given mass-balance record, the t value is given by
where ΔB is the magnitude of total change in mass balance as estimated by a least squares linear fit to the data over the observing period, σ is the standard deviation of the detrended residual, and ν is the number of degrees of freedom. The critical t value for a given threshold of confidence and number of degrees of freedom can be found in standard statistics tables, and we can then solve for a critical ΔB/σ (the signal-to-noise ratio) needed to declare a trend significant.
A variable with persistence (i.e., serial correlation) on the order of the sampling interval Δt has fewer independent degrees of freedom than observations. To account for the possibility of persistence, we estimate ν from the autocorrelation function of the time series. One straightforward metric based on Bartlett (1946) states that the 95% confidence bound for a lag-1 autocorrelation rΔt indistinguishable from zero is , where N is the number of observations. After linear detrending, all of the mass-balance records fall below this threshold of ~0.27, indicating that the test does not provide evidence of persistence. Another method, based on Leith (1973), assumes a red noise process with a decorrelation (e-folding) time τ estimated from rΔt. In this case,
Following this metric, ν/N approaches 1 as rΔt approaches e−2 (~0.13). Again, rΔt for all records fall below this threshold. On this basis we conclude that the observations lack persistence, and thus each time series is consistent with N independent degrees of freedom. We discuss the sensitivity of our conclusions to the choice of ν in the next section.
The critical t value, and thus the critical signal-to-noise ratio, also depends on whether a one-tailed or two-tailed test is used—that is, whether trends breaching the chosen confidence level can come from one or both sides of the t distribution. A one-tailed test can enhance detectability provided the sign of the trend can reasonably be assumed a priori. However, considering that internal climate variability can have a considerable impact on decadal climate trends (e.g., Deser et al. 2012; Wallace et al. 2012), and lacking an a priori assumption about the sign of this trend contribution, we argue that a two-tailed significance test is most rigorous for this study. This allows for either positive or negative trends to be detected but requires a higher signal-to-noise ratio for significance at a given confidence level.
For the three target glaciers, we see a common pattern. All three have negative trends in summer balance, significant at the 5% level (Table 1). However, none of the glaciers have statistically significant trends in winter balance. Although several studies have documented declining snowpack in Washington’s Cascade Mountains (e.g., Mote et al. 2005; Casola et al. 2009; Stoelinga et al. 2010), winter balance depends on a different hypsometric signature on the landscape as well as other local effects specific to the glacier’s setting within the mountain range (e.g., avalanching, wind deposition, microclimates) and thus may not exhibit the same trends as regional snowpack.
Finally, driven by the summertime contribution, the annual mass balance has a negative trend for all three glaciers but is significant only for Wolverine and Gulkana.
b. Adjustment of winter balance
We now apply dynamical adjustment to investigate what fraction of each mass-balance record can be attributed to variability in large-scale atmospheric circulation. We first present results for winter mass balance, as they demonstrate the clearest dynamical link between predictors and predictand. For South Cascade Glacier, the leading SLP predictor pattern explains 53% of variance in winter balance. The spatial pattern of the leading SLP mode (Fig. 2a) shows that variability in winter accumulation is strongly correlated with the strength of the Aleutian low, a persistent wintertime pattern of low surface pressure and cyclonic circulation in the North Pacific. Variability in this feature reflects shifts in the latitude and intensity of the winter storm track and thus variability of precipitation over South Cascade Glacier (e.g., Bitz and Battisti 1999).
A second iteration of PLS regression yields a second pattern explaining an additional 10% of the variance in the original record (Fig. 2b). In general, spatial patterns for successive modes are more difficult to interpret physically since the remaining correlations after the preceding mode(s) have been removed are weaker and less spatially coherent. Additionally, the constraint that modes are mutually orthogonal can obscure physical interpretations, given that dynamical processes are not necessarily linearly independent (e.g., Hannachi et al. 2007). The two orthogonal indices t1, t2 generated from these patterns are shown in Fig. 2c, and their combined effect is seen in Fig. 2d; the adjusted time series has 63% less variance than the raw winter balance record.
Using SST as a predictor returned similar results for South Cascade winter balance; the leading patterns explained 47% and 12% of variance, respectively. The spatial correlation patterns are distinct for SST, with the leading pattern resembling the spatial signature of the Pacific decadal oscillation (e.g., Mantua et al. 1997) (see Fig. 3a). However, the indices (Fig. 3c) produced from SST are strongly correlated with those of SLP (r = 0.84 for leading mode), indicating that the dynamical adjustment extracts a common element of variance in South Cascade Glacier’s mass balance driven by the coupled ocean–atmosphere system. Similar to results using SLP, the adjustment with SST removes 59% of the original variance in Bw (Fig. 3d).
The dynamical adjustments for Wolverine Glacier accounted for similar, albeit slightly weaker, components of variance in winter balance. The first two SLP patterns explained 46% and 9% of variance, and the SST patterns explained 34% and 9%, respectively. Notably, the leading correlation patterns associated with SLP and SST have similar structure but opposite polarity to their respective counterparts for South Cascade Glacier (Fig. 4). As a result, the leading indices for Wolverine and South Cascade Glaciers are strongly negatively correlated (r ≈ −0.9 for SLP and −0.8 for SST). For the years that the mass-balance observations overlap for all three glaciers (1966–2011), the records themselves are moderately anticorrelated (r = −0.42), a relationship that has been attributed to circulation anomalies that shift winter storms toward southeast Alaska and away from Washington State, or vice versa (e.g., Walters and Meier 1989; Bitz and Battisti 1999; Rasmussen and Conway 2004). Additionally, Bitz and Battisti (1999) showed that winter balance has different temperature sensitivities between these settings, with warmer winters more favorable for accumulation in Alaska and the opposite for Washington. This may then be a coupled temperature and precipitation signal; it is difficult to say which dominates from the patterns alone, although the end effect is reflected in winter accumulation. The adjusted records for Wolverine and South Cascade Glaciers are effectively uncorrelated (r = −0.05 using SLP), reinforcing that the adjustment is indeed identifying and removing the dynamical contribution to mass-balance variability.
Winter mass balance for Gulkana Glacier showed a much weaker connection to the predictors; SLP and SST patterns explained 23% + 7% and 22% + 6%, respectively. Recalling that Gulkana’s winter balance record has comparatively little variability to begin with (Fig. 1), the year-to-year changes in accumulation that are driven by internal North Pacific climate variability are quite small, consistent with its more continental setting.
In summary, by applying dynamical adjustment we can account for almost two-thirds of the total variance in the maritime glaciers’ winter accumulation, indicating a close connection between internal variability in North Pacific circulation and winter mass balance; however, the case of Gulkana Glacier, with less than one-third of variance explained, suggests that this relation is muted for continental glaciers, which are buffered by coastal mountains and long overland fetches.
c. Adjustment of summer balance
For the maritime glaciers (South Cascade and Wolverine), warm-season (April–September) SLP and SST explain substantially less variance in the summer balance records. South Cascade Glacier is more closely linked with SLP than SST, with two predictor patterns explaining 28% + 7% of variance for SLP, compared with only 15% + 6% for SST. (The second modes may not in fact be significant but are reported here for the sake of consistency.) For Wolverine Glacier, both predictors explain roughly the same total variance, with 21% + 8% for SLP and 24% + 4% for SST. The variance explained by summer SLP for Gulkana Glacier (22% + 9%) is similar to that of the maritime glaciers, but it is notable in that it makes Gulkana the only glacier with comparable amounts of dynamically induced variability in summer and winter mass balance. This is broadly consistent with the rule of thumb that continental glaciers, being further removed from the ocean’s moisture-laden storms and moderating effect on temperature, are comparatively more sensitive to melt-season climate (e.g., Medwedeff and Roe 2016). It is perhaps not surprising that summer SSTs are a weak predictor for Gulkana Glacier’s summer balance, explaining only 15% of total variance.
By far the strongest summer SST relationship is the leading mode for Wolverine Glacier, explaining 24% of variance. The correlation map (Fig. 5a) shows what may be an intuitive result for this glacier situated closest to the Pacific: nearshore SST and summer balance are anticorrelated (i.e., years of anomalously warm ocean surface are associated with more summer melt). Correlation maps for the leading summer SLP patterns (associated with ~20% of variance for each glacier) all show moderately positive correlation with SLP to the south of each glacier (Figs. 5b–d). However, the links between summer climate anomalies and summer balance are weaker than for winter, consistent with generally less vigorous circulation in summer.
d. Adjusted trends
Our analyses show that a substantial amount (~30%–60%) of variance in mass balance is explained by variability in SLP and SST. We can now return to one of our main motivating questions: What is the role of natural variability in the observed trends in mass-balance records? It is important to note that, although the correlation map W [Eq. (1)] is generated from detrended data, the component of variance removed from the predictand may itself contain a linear trend. This is because the index t and the regression coefficients β and P contain data that have been standardized but not detrended [see Eqs. (3)–(5)]. Thus, low-frequency variance that might be projecting onto the linear trend can in principle be removed from the mass-balance records—indeed, this is one motivation for pursuing dynamical adjustment. Though all of the dynamically adjusted mass-balance records necessarily have reduced overall variance, removing this low-frequency variance may yield residual time series [Yadj in Eq. (6)] that have enhanced or reduced trends, affecting the interpretations of observed trends and their causes.
Since both the noise and the trend are affected by the dynamical adjustment, we revisit the t statistic [Eq. (8)] as a way to compare the raw and adjusted mass-balance records. Recall that the significance of a trend in a given number of independent observations depends on the signal-to-noise ratio ΔB/σ. This ratio is a useful metric for comparing the relative significance of trends, and for significance at the 5% level, our records requires that ΔB/σ > 0.98 for N = 53 years (South Cascade) and ΔB/σ > 1.00 for N = 50 years (Wolverine and Gulkana). Table 2 shows ΔB/σ for all raw records and adjustments, and we note that trends in the South Cascade and Wolverine winter records are not significant in raw or adjusted form. Gulkana’s adjusted winter balance trend barely breaches the threshold for significance at 95% when using SLP as a predictor but is just short of it using SST. However, all three of the summer records have significant trends in raw form and improved signal-to-noise ratios after adjustment. In other words, the trends in summer mass balance are not associated with low-frequency circulation variability and are more consistent with an external forcing.
a. Methods for testing trend significance
Several factors should be borne in mind when interpreting these results. First of all, recall that the threshold for trend significance is dependent on the number of degrees of freedom ν estimated from the time series’ autocorrelation [Eq. (9)]. Given that the autocorrelation functions are rather noisy for these short time series, it is worth considering significance estimates for the case ν < N to account for the possibility that sampling effects have obscured weak persistence. While this does increase the critical signal-to-noise ratio for significance, we find that our conclusions are robust to different choices of ν. Consider a decorrelation time twice the maximum allowed for our initial assumption that ν = N, which then cuts ν in half. If ν = 25 (26), the critical ΔB/σ for 95% confidence required by the t test is 1.49 (1.45) for a time series of 50 (53) years. Even in this case, the interpretation would remain that all summer trends are significant at 95%, both before and after dynamical adjustment (see Table 2).
As an alternative to the t test, we also applied the method of phase randomization to test for trend significance (see, e.g., Theiler et al. 1992; Vafeiadis et al. 2008). In this framework, a random phase rotation is applied to each Fourier frequency component of the detrended time series. After transformation back into the time domain, this yields a surrogate time series where the spectrum and autocorrelation structure have been preserved. A large ensemble (104–105) of these surrogate datasets can then provide a distribution of trends that arise merely from undersampling the spectrum of the original time series. The observed trends for raw and adjusted summer balance fell well beyond the central 95% of the distribution. However, the method of phase randomization has some caveats for short time series—namely, that the probability distributions of surrogate data are not well preserved. Thus, while consistent with the other trend tests, this method may not stand on its own in this application.
b. Choice of predictor variable
Additionally, we find that the trend changes after adjustments vary with the choice of predictor (however, the significance at 95% changes only for Gulkana Glacier’s winter balance, where both trends are near the threshold). While SLP and SST are dynamically linked and yield highly correlated indices (see section 4b), their trend adjustments are not consistent. For South Cascade and Wolverine Glaciers, SLP and SST adjust the mass-balance trends in opposite directions, in both the summer and winter cases (Fig. 6). This is consistent with the overall result that circulation variability has not contributed substantially to the trends over the period of observation.
In addition to the choice between SLP and SST, we also assessed the sensitivity to the choice of data source for the same climate variable. The Met Office HadSLP2 dataset is an alternative SLP dataset, providing monthly means on a 5° × 5° grid (Allan and Ansell 2006). Results using this predictor are qualitatively equivalent to those using the NCEP–NCAR reanalysis; the correlation patterns and respective amounts of variance explained are very similar, and the trend adjustments are also minimal.
c. Observational error
Both predictor and predictand datasets contain some degree of observational and sampling error. The datasets for the Alaskan glaciers have been reanalyzed, as described by Van Beusekom et al. (2010) and in the supplementary material of O’Neel et al. (2014). Work by the USGS addresses both systemic and isolated errors in reported balance data, and ongoing efforts include calibrating field measurements with geodetic methods using available DEMs.
What are the effects of remaining measurement error on dynamical adjustment? Since PLS regression operates only on anomalies, it is insensitive to systemic biases in mass-balance data (unlike, e.g., estimates of cumulative balance). Random errors, however, may confound correlation maps and trend estimates. But while a large error in one year’s reported mass balance certainly could bias an initial linear trend estimate, it is unlikely to be erroneously removed by dynamical adjustment, and thus the changes to trends (or lack thereof) should be relatively insensitive to this form of uncertainty. Furthermore, the cases in which the statistics are shored up by a physical interpretation, as in the relationship between South Cascade, Wolverine, and the Aleutian low (section 4b), suggest that observational or sampling errors in the mass balance or predictor datasets are not necessarily a barrier for identifying relevant dynamical patterns.
d. Conventional mass balance vs reference surface mass balance
While we have focused on partitioning the effects of internal climate variability and external forcing on mass-balance changes, it is also important to consider the effect that a changing glacier geometry can have on mass balance measurements. Terminus changes constitute a negative feedback by changing the ablation area, while thickness changes can counteract this by lowering or raising surface elevation across the glacier. These processes have motivated analyses of reference surface mass balances, which are calibrated to a constant glacier hypsometry (e.g., Elsberg et al. 2001; Huss et al. 2012). Reference surface values were available for Wolverine and Gulkana Glaciers through 2009 (Van Beusekom et al. 2010), and we found that dynamically adjusting these data produced qualitatively similar results to those using conventional balance. Trend estimates are slightly different for reference surface balances, but summer trends retained significance, suggesting that they are not principally a geometrical effect.
e. Length of time series
Finally, the length of the predictand time series has some interesting implications. We focus on the USGS benchmark glaciers in part because they offer long uninterrupted records and thus offer more information to be compared with the other climate datasets. However, the trend biases caused by the low-frequency components of natural variability can be more dramatic in shorter records. This is a familiar challenge associated with limited sampling, but it means that the skill of dynamical adjustment in terms of removing dynamically induced trends is sensitive to record length. As an illustrative example, we consider separate adjustments for two halves of the South Cascade winter record (Fig. 7). In the raw data, the periods from 1959 to 1984 and from 1985 to 2011 have markedly different trends, a feature that is lost when a linear trend is estimated for the whole record. The separate adjustments yield more dramatic trend changes than with the entire record. For example, winter balance from 1985 to 2011 has a strong positive trend (ΔB/σ = 1.35; significant at 90% but not at 95%), but this is nearly completely removed by dynamical adjustment, indicating that it is associated with circulation variability. Of course, the half-and-half split is arbitrary, and we do not attempt here to define an optimal time series length for adjustment. However, given the preponderance of one- to two-decade mass-balance records for glaciers around the world (Medwedeff and Roe 2016), these may be important considerations for future applications of dynamical adjustment.
6. Summary and conclusions
We have applied PLS-based dynamical adjustment to seasonal glacier mass balance and have found this to be a useful method for analyzing dynamically induced variability, explaining from about one-third of summer variance to well over half of winter variance in these records. Building on its previous applications to regional averages or fields (e.g., snowpack, surface air temperature, and hurricane activity), we find that the method of dynamical adjustment is also capable of identifying relationships between large-scale circulation and more localized processes (i.e., accumulation and ablation on individual glaciers).
For winter balance, dynamical adjustment shows that a large portion of accumulation variability is explained by North Pacific variability and is captured in comparable proportions by both SLP and SST reanalysis fields. The effect is strongest for the maritime glaciers (South Cascade and Wolverine), explaining over half of the variance. The adjustments also demonstrated that the primary source of dynamically induced variability for these glaciers is variability in the Aleutian low. For summer balance, dynamical adjustment using SLP as a predictor explained approximately one-third of variance for all glaciers but substantially less when using SST. The variability left unexplained by dynamical adjustment is likely the aggregate of many factors, including the nuances of regional climate, mountain meteorology, and the noise and sampling effects inherent in the climate and mass-balance datasets. Further analysis using meteorological data and/or regional climate models could provide insight into such processes, as well as the limits of the relationships that can be resolved with dynamical adjustment.
Dynamical adjustment also reveals whether trends in time series are related to circulation variability. We found that, despite the large degree of circulation-related variability in seasonal mass balance, the dynamical adjustment did not change the trends substantially. However, by accounting for the possibility of a dynamically induced bias, we have a refined view of mass-balance change for these glaciers. The negative trends in summer balance, which are significant in the raw data, persist with greater signal-to-noise ratios when dynamically induced variability is removed. This supports an interpretation that the trends are externally forced and associated with anthropogenic greenhouse warming. However, no statistically significant trends exist in adjusted winter balance (with the exception of Gulkana’s winter balance, but its significance is not robust to the choice of predictor). This strengthens the interpretation that anthropogenic trends in accumulation over the glaciers have yet to emerge.
The absence of winter trends is common to many winter mass-balance records (Medwedeff and Roe 2016) and may simply reflect the different temperature sensitivities of accumulation and ablation. In theory, warming winters would eventually cause a decline in accumulation as the average freezing level climbs, though in dry, cold settings the initial effect may instead be increased precipitation due to a greater atmospheric moisture capacity. Indeed, winter warming and elevated freezing levels have been observed in northwestern North America, but the actual effect on winter balance varies depending on the hypsometric profile of the glacier and the baseline climatology (Arendt et al. 2009). It follows that the time of emergence—and even the sign—of winter balance trends under warming may vary considerably by region and by glacier. Since such changes would inevitably exist amid natural variability, dynamical adjustment could be a useful method for diagnosing any new trends that emerge.
The general question of the relative importance of external forcing and internal variability in observed trends is an important one. Several studies have addressed it for glacier mass balance in different ways (e.g., Huss et al. 2010; Marzeion et al. 2014). Huss et al. (2010) concluded from running-mean correlations that up to half of the twenty-first-century declines in mass balance in the Swiss Alps could be accounted for by the Atlantic multidecadal oscillation. Dynamical adjustment would provide an alternative analysis method that could address the same question and would supply the associated correlation patterns that could be evaluated for dynamical mechanisms. Marzeion et al. (2014) used ensembles of global climate models to estimate the magnitude of natural mass-balance variability relative to simulated trends, aggregated into different regions. Dynamical adjustment would provide a complementary analysis here, too.
Finally, it is also important to understand how much of a model’s projection of a future climate trajectory may be attributable to internal variability (e.g., Deser et al. 2012). By identifying the variability associated with circulation patterns in the modern climate, dynamical adjustment could again be used to constrain the externally forced portion of the climate projection. Dynamical adjustment is thus a versatile tool, whose wider application to a variety of glacier mass-balance settings and problems may prove useful.
We gratefully acknowledge the USGS Climate and Land Use Change program for maintaining the benchmark glacier-monitoring program. We also thank Shad O’Neel for helpful discussion regarding mass-balance measurements and uncertainty. Finally, we thank editor Peter Huybers and three anonymous reviewers for comments that greatly improved the manuscript. J.E.C. was supported by the National Science Foundation Graduate Research Fellowship Program (DGE-1256082).