2. TWO MODELS OF CLIMATIC PERSISTENCE

We consider two commonly invoked statistical models of climatic persistence. The first is a first-order autoregressive process (≡ AR(1); e.g. Box and others, Reference Box, Jenkins and Reinsel2008). Let b _t be the annual-mean mass balance in year t. An AR(1) process can be written as

(1) $${b_t} = r{b_{t - \Delta t}} + {a_0}\Delta t\hskip1pt{\mu _t},$$

where Δt = 1a, r and  a ₀ are constants and μ _t is a random number drawn from a normally-distributed stochastic process. Thus b _t depends in part on its value in the previous year (and so has persistence), and in part on a random stochastic impulse occurring in the current year. r is the autocorrelation coefficient at lag 1 a. Equation (1) is the discrete form of the continuous equation: db/dt + b/τ _c = a ₀ μ(t), and so r can be related to a climatic persistence timescale (or memory),  τ _c:  r = 1 − Δt/τ _c. This model of persistence was introduced into climate science in a series of papers by Klaus Hasselmann and colleagues (e.g. Hasselmann, Reference Hasselmann1976; Frankignoul and Hasselmann, Reference Frankignoul and Hasselmann1977), which demonstrated that persistence in sea-surface temperatures was the result of the ocean-mixed layer acting to integrate rapid and essentially random, uncorrelated fluctuations of air-sea heat fluxes.

Note that one can also consider a higher-order autoregressive process, AR(p), which depends on the previous p states (e.g. Box and others, Reference Box, Jenkins and Reinsel2008), and which can be interpreted as the discrete form of a p ^th-order ordinary differential equation. Reichert and others (Reference Reichert, Bengtsson and Oerlemans2002) fit an AR(3) process to climate-model output for Nigardsbreen Glacier in Norway, for instance. The AR(3) coefficients they fit imply a persistence timescale of 1.5 a, and also an underlying damped oscillation with a 3 a period and a 1 a phase-memory. It is important to make a physical interpretation of the higher-order AR(p) process, and to establish the statistical necessity of the higher-order model via standard methods (e.g. Box and others, Reference Box, Jenkins and Reinsel2008). We only consider AR(1) here.

For an AR(1) process, the power spectral density as a function of frequency and lag-1 a autocorrelation (≡ P _b(f, r)), is given by:

(2) $${P_{\rm b}}(\,f,r) = \displaystyle{{{P_0}} \over {1 + {r^2} - 2r\cos (2\pi \Delta tf)}},\quad {\rm on}\quad 0 \le f \le \displaystyle{1 \over {2\Delta t}}.$$

where P ₀ is a constant, defined in the Appendix. P _b(f, r) increases towards lower frequencies, and by analogy with spectra for visible light, is often termed red noise.

The autocorrelation function (ACF) is the formal statistical description of the structure of the persistence in a time series (e.g. Box and others, Reference Box, Jenkins and Reinsel2008), and it characterizes the similarity between the data as a function of the time lag between them. Let  nΔt represent a lag of n years. The ACF is defined by ρ(nΔt) = 〈b _t b _t−nΔt 〉. From the Wiener-Kinchin theorem, the ACF is the Fourier transform of the power spectrum: $\rho (n\Delta t) = {{\int_0^{{f_0}} {P_{\rm b}}(f)\cos (2\pi n\Delta tf){\rm d}f} / {\int_0^{{f_0}} {P_{\rm b}}(f){\rm d}f}}$ .

For an AR(1) process the ACF is:

(3) $$\rho (n\Delta t) = {\left( {1 - \displaystyle{{\Delta t} \over {{\tau _{\rm c}}}}} \right)^n}{\rm \simeq} \exp \left( { - \displaystyle{{n\Delta t} \over {{\tau _{\rm c}}}}} \right),$$

and thus ρ asymptotes exponentially towards zero with an e-folding timescale, τ _c.

The second model for persistence characterizes the spectrum of a time series as a power law. This form is suggested by several analyses of long instrumental and paleoclimate proxy records for temperature (e.g. Pelletier, Reference Pelletier1997, Reference Pelletier1998; Fraedrich and Blender, Reference Fraedrich and Blender2003; Fraedrich and others, Reference Fraedrich, Luksch and Blender2004; Huybers and Curry, Reference Huybers and Curry2006), and also for precipitation and other hydrological variables (Hurst, Reference Hurst1951; Ault and others, Reference Ault2013). Let  ν be the slope of the power spectrum when plotted on logarithmic axes. So in this case, the spectrum P _b(f, ν), is defined by:

(4) $${P_{\rm b}}(\,f,\nu ) = {P_0}{({\,f_0}/f)^\nu}, $$

where P ₀ is the spectral density at frequency  f ₀ = 1/(2Δt), the maximum resolvable frequency from an annual-mean time series; and ν is a positive constant. Thus, spectral power continually increases towards low frequencies.

For the ACF there is an exact analytical solution in the form of an ordinary hypergeometric function.

(5) $$\rho (n\Delta t{) =\; _2}{F_1}\left[ {\displaystyle{{1 - \nu} \over 2},\left\{ {\displaystyle{1 \over 2},\displaystyle{{3 - \nu} \over 2}} \right\}, - {{(2\pi n\Delta t{\,f_0})}^2}} \right],$$

which looks intimidating but the key point is that, at large lags, it decays towards zero as (nΔt) ^ν−1 (e.g. Percival and others, Reference Percival, Overland and Mofjeld2001; Box and others, Reference Box, Jenkins and Reinsel2008). It thus decays much more slowly than an AR(1) process (Eqn 3), and so exhibits greater persistence at long timescales. For this reason Eqn (4) is referred to as a ‘long memory’ process. This behavior will be shown in Figure 2.

Huybers and Curry (Reference Huybers and Curry2006) present spatial maps of ν for 100 year-long temperature observations, and find values between 0 and 0.5 over land, with the highest values confined to maritime climate settings. Values can exceed 1 over the ocean, particularly in the tropics. Ault and others (Reference Ault2013) fit Eqn (4) to spectra of observed precipitation in western North America and find a patchwork pattern with ν ranging from ~−0.2 to +0.5, but with marginal statistical significance. Persistence evaluated from paleoclimate proxy time series should be regarded with some caution, since many of the proxies themselves act as integrators of climate, which together with dating uncertainty and post-depositional processes can act to redden the proxy time series independent of the driving climate. But such datasets do indicate that internal variability can be characterized as a power-law spectrum (Pelletier, Reference Pelletier1997, Reference Pelletier1998; Fraedrich and Blender, Reference Fraedrich and Blender2003; Fraedrich and others, Reference Fraedrich, Luksch and Blender2004; Huybers and Curry, Reference Huybers and Curry2006). Furthermore there are theoretical grounds to think that power-law behavior should be expected across decadal to millennial frequencies: for quasi-diffusive heat uptake in the deep ocean, progressively more of the ocean becomes involved at low frequencies, increasing the effective inertia of the climate system (Hoffert and others, Reference Hoffert, Callegari and Hsieh1980; Pelletier, Reference Pelletier1997; Fraedrich and others, Reference Fraedrich, Luksch and Blender2004), which is also diagnosed in comprehensive climate models (e.g. Zhu and others, Reference Zhu, Fraedrich, Liu and Blender2010; MacMynowski and others, Reference MacMynowski, Shin and Caldeira2011). Power-law spectra can also be generated by a fractionally differenced process, which is equivalent to an aggregation of AR(1) processes with many different timescales (Beran, Reference Beran1994; Percival and others, Reference Percival, Overland and Mofjeld2001; Box and others, Reference Box, Jenkins and Reinsel2008).

2.1. Observed variability and persistence in glacier mass-balance records

We first evaluate the variability and persistence of three of the longest annual-mean mass-balance observations available. Our goal is not to present a comprehensive evaluation of glacier mass-balance records, which is done in Medwedeff and Roe (Reference Medwedeff and Roe2016), but the results shown here are typical. The longest continuous mass-balance record is that on Clariden glacier, Switzerland (since 1914), although the measurements are only at two snow stakes, and so not necessarily representative over the whole glacier; Storbreen (since 1949) is the longest record in Norway; and South Cascade glacier in Washington State (since 1959) is the longest in North America.

Figure 1a presents these records after they have been linearly detrended to remove any anthropogenic trend. The standard deviations in mass balance (≡ σ _b) and 95% bounds for each glacier are: Clariden, 0.74 {0.65, 0.86} m a⁻¹; Storbreen, 0.69 {0.58, 0.83} m a⁻¹; South Cascade, 1.04 {0.87, 1.3} m a⁻¹, where the confidence bounds have been calculated assuming the n-years sample variance follows a χ ² distribution with n − 1 degrees of freedom (e.g. VonStorch and Zwiers, Reference vonStorch and Zwiers1999). All three datasets are consistent with a normal distribution at better than 95% confidence based on a Jarque-Bera test (e.g. Steinskog and others, Reference Steinskog, Tjøtheim and Kvamstø2007).

Fig. 1. Three long records of observed annual-mean glacier mass balance. Clariden (upper stake), Switzerland (95 a, σ _b = 0.74 m a⁻¹); Storbreen, Norway (63 a, σ _b = 0.69 m a⁻¹); and South Cascade, USA (54 a, σ _b = 1.04 m a⁻¹). (a) Time series of the observations, with the linear trend and mean removed. (b) Autocorrelation function. The horizontal dashed lines are drawn at $2/\sqrt n $ , where n is the length of each record in years – these lines indicate the level above which the lag 1 a autocorrelation should lie to reject the null hypothesis of white noise at better than 95% confidence (Bartlett, Reference Bartlett1946; Box and others, Reference Box, Jenkins and Reinsel2008). (c) Power spectra (unwindowed periodogram) of mass-balance data. Dashed lines show the best-fit slopes (based on a least-squares linear regression). All slopes are positive, but none are significantly different from zero at better than 95% confidence.

Most mass-balance records are not as long as those we selected, of course, and it is important to emphasize that, for shorter records, the uncertainty estimates on σ _b are quite large. For a 20 a mass-balance record, for instance, the 95% bounds are {0.76, 1.46} times the standard deviation of the data sample. In other words, for many records, there is approximately a factor of two uncertainty in σ _b (Medwedeff and Roe, Reference Medwedeff and Roe2016).

Figure 1b shows the ACFs for the three records. The null hypothesis of no persistence can only be rejected at 95% confidence if the lag one-year autocorrelation exceeds $2/\sqrt n $ (Bartlett, Reference Bartlett1946; Box and others, Reference Box, Jenkins and Reinsel2008). As is clear from Figure 1, on the basis of this test one concludes that all three records cannot be distinguished from a null hypothesis that the fluctuations are uncorrelated (i.e. they are white noise). This is confirmed using standard algorithms for selecting and fitting AR(p) models (e.g. Box and others, Reference Box, Jenkins and Reinsel2008), which identify AR(0) as the best description for all three records. Burke and Roe (Reference Burke and Roe2014) fit AR(p) models to long instrumental records and glacier mass-balance records in Europe, and also concluded persistence could not be definitively established.

An alternative measure of persistence is to evaluate whether the power spectra of the records have non-zero slopes. Figure 1c shows the power spectra for the three records (unwindowed periodograms), together with the best-fit slopes, based on a least-squares linear regression. The results for the slopes, ν, are: Clariden, 0.25 ± 0.32 (2σ); Storbreen, 0.23 ± 0.61 (2σ); South Cascade, 0.09 ± 0.64 (2σ). For all the records the slopes are positive, but the 2σ ranges all straddle zero and so none are statistically different from ν = 0 (Clariden gets close with a p value of 0.12). Therefore, by this measure also, one concludes the observations cannot be distinguished from white noise.

Medwedeff and Roe (Reference Medwedeff and Roe2016) evaluate the full, global mass-balance dataset available from the World Glacier Monitoring Service (WGMS, Reference Zemp2013 and earlier issues) and find these results apply more broadly: observed annual-mean and season-mean mass-balance records are consistent with white noise (plus a linear trend).

We are therefore left with the situation that, while paleoclimate proxies and long instrumental records both suggest that there is a degree of persistence in climate time series (which varies as a function of location and climate variable), glacier mass-balance records are too short to detect such persistence statistically, or to accurately determine its type (i.e. AR(1) or power law) or degree (e.g. Percival and others, Reference Percival, Overland and Mofjeld2001). Nor can persistence in mass balance be easily inferred from other climate variables. Persistence is different for different variables. Glacier mass balance reflects a complex amalgam of many variables, and also depends in part on other stochastic processes such as avalanching and drifting. Nonetheless, as has been preliminarily explored by Reichert and others (Reference Reichert, Bengtsson and Oerlemans2002), climatic persistence has the potential to significantly enhance glacier length fluctuations. In this paper we take a what-if approach: what if climatic persistence, of a degree that is hard to detect, exists in glacier mass balance? What would be the impact of such persistence on the amplitude of glacier length fluctuations? The implications for the interpretation of past glacier records are broached in Section 4.

2.2. Synthetic climate time series with persistence

To explore the impacts of possible undetected persistence in mass balance on the interpretation of glacier-length records, we first generate synthetic mass-balance records with differing degrees of persistence. We then compare the theoretical glacier responses to these synthetic mass-balance histories with the responses that would be predicted under the erroneous assumption that the mass-balance fluctuations represent white noise.

We generate five synthetic 10 ka mass-balance time series using the method outlined by Percival and others (Reference Percival, Overland and Mofjeld2001). For a given choice of r or ν, the time series is generated from the Fourier transform of the spectrum amplitude (the square root of either Eqn (2) or (4)), multiplied by a random phase at each frequency. We use the same set of random phases for each time series, to facilitate comparison.

Our control mass-balance time series is white noise, with standard deviation σ _b = 1 m a⁻¹. For power-law persistence, we take values for ν of 0.25 and 0.4, since these sit within the range suggested for land-based climate variables. The lag one-year autocorrelation coefficients for such power-law persistence are 0.17 and 0.28 (Eqn 5), and these are the two values we use for r for AR(1)-type persistence. The memory, τ _c, associated with these r values is 1.2 and 1.4 a, respectively. Based on the ~95% confidence intervals (Bartlett, Reference Bartlett1946), for a given r it takes n ~ (2/r)² years to detect persistence. Thus, for our choices for ν of 0.25 and 0.4 it would require ~150 and 50 a of data, respectively, before persistence could be established. Such durations are comparable with some of the longest available instrumental and glacier mass-balance records. In other words we are emulating a realistic situation: we are introducing persistence with a magnitude that is hard to detect in observations.

The inclusion of persistence enhances the variance of a time series. However if the presence of persistence went undetected, the measured standard deviation would be taken to be the standard deviation of a white-noise process. The parameters of the assumed white noise used to predict glacier length fluctuations would then be based on this standard deviation. We therefore apply a correction, normalizing each time series so the standard deviation is equal to that of the control series (σ _b = 1 m a⁻¹). This ensures that the variance we use to drive the glacier models is consistent with what would be determined from observations – the only thing that is different is the persistence present in each time series. This correction is conservative in the sense that it reduces the resulting glacier variance. Without the correction, adding persistence would increase the standard deviation of mass balance (and glacier response) by 1.5, 4.2, 15 and 29%, for r = 0.17, 0.28, ν = 0.25, 0.4 respectively.

Figure 2a presents a 100 a segment of each of these synthetic climate time series. The series are virtually indistinguishable by eye, but the ACFs (Fig. 2b) reveal the differences. The AR(1) processes asymptote exponentially to zero (Eqn 3), whereas the ACFs of the power-law processes decline more slowly (Eqn 5) and thus represent a climate with a small-but-persistent chance that a series of successive years have the same sign anomaly (i.e. a long memory). The spectra are shown in Figure 2c and show that the synthetic time series accurately match the theoretical power spectra.

Fig. 2. Synthetic mass-balance records, with differing degrees of persistence. (a) a 100 a time slice from the 10 ka realizations used to drive the glacier models. Time series for white noise, ν = 0.25, 0.40 and r = 0.17, 0.28 are shown, and have been normalized so that σ _b = 1.0 m a⁻¹ in each case. For clarity, the time series have been successively offset by 2000 mm on the y-axis. (b) Autocorrelation function for the five time series, indicating the persistence present in each record. (c) Power spectra (windowed periodogram) of the five time series (thick curves), and the theoretical spectra for each process (thin curves, see Eqns (2) and (4)). The discrepancies at the lowest frequencies are artifacts of the windowing.

In Section 3, we drive two different glacier models with these synthetic mass-balance time series. We also derive algebraic formulae that capture the impact of climatic persistence on glacier length.

3. IMPACT OF CLIMATIC PERSISTENCE ON GLACIER LENGTH

The first glacier model we use is the numerical flowline model from Roe and O'Neal (Reference Roe and O'Neal2009), with shallow-ice dynamics and basal sliding. It uses the same set of parameters as Roe and Baker (Reference Roe and Baker2014), characteristic of the small Alpine glaciers around Mt. Baker in Washington State, U.S.A. Readers are referred to Roe and Baker (Reference Roe and Baker2014) for the equations and more details.

The second model is the three-stage model, also of Roe and Baker (Reference Roe and Baker2014), a linear model for length fluctuations around a prescribed mean state, which was shown to closely emulate the behavior of the numerical model under a wide range of parameter choices. Its relatively simple form makes it an extremely versatile tool. It provides analytic solutions for important metrics of the response, such as a glacier's power spectrum, ACF, variance and excursion probabilities in response to climate variability, and a glacier's sensitivity to climate trends and step-changes. The three-stage equation is

(6) $${\left( {\displaystyle{{\rm d} \over {{\rm d}t}} + \displaystyle{1 \over {{\epsilon} {\tau _{\rm g}}}}} \right)^3}L^{\prime}(t) = \displaystyle{\beta \over {{{\epsilon} ^3}\tau _{\rm g}^2}} b(t),$$

where L′ is the perturbation in length away from a prescribed mean state; τ _g = −H/b _term is the glacier response timescale originally proposed by Jòhannesson and others (Reference Jòhannesson, Raymond and Waddington1989), where H is the mean thickness, and b _term is the (negative) net mass balance at the glacier terminus, both in the mean state; and ε is a constant ( $ = 1/\sqrt 3 {\rm \simeq} 0.58$ ). The glacier geometry also determines β = A _tot /(wH), where A _tot is the total glacier area and w is the characteristic valley width at the terminus. As in Roe and Baker (Reference Roe and Baker2014), for our standard control glacier, we take τ _g = 6.74a, and β = 178 (unitless).

Equation (6) reflects the three stages of glacier adjustment identified and described in Roe and Baker (Reference Roe and Baker2014): changes in mass balance first cause changes in thickness, which then causes changes in ice flux at the terminus, which finally causes changes in length. These three stages are robust to the spatial pattern of mass-balance changes, and are intrinsic to the high-frequency response of a glacier. Roe and Baker (Reference Roe and Baker2014) further showed that Eqn (6) could also be written as a discrete equation in time steps of  Δt = 1a. The discrete equation more closely matches the nature of actual mass balance, coming as it does in discrete winter and summer seasons. Roe and Baker (Reference Roe and Baker2014) showed that, for white-noise forcing, the length variance ( $ = \sigma _{\rm L}^2 $ ) for the discrete three-stage model is:

(7) $$\sigma _{\rm L}^2 = {\beta ^2}\tau _{\rm g}^2 {\sigma _{\rm b}^2} \displaystyle{{(1 - \kappa )(1 + 4{\kappa ^2} + {\kappa ^4})} \over {2\Delta t{{(1 + \kappa )}^5}}},$$

where κ = (1 − (Δt/ετ _g)). In the limit as ετ _g ≫ Δt this simplifies to:

(8) $$\sigma _{\rm L}^2 = \displaystyle{{3\Delta t{\tau _{\rm g}}} \over {16{\epsilon}}} {\beta ^2}{\sigma _{\rm b}^2}. $$

For σ _b = 1 m a⁻¹, Eqn (7) predicts σ _L = 284 m, which agrees closely with the output from numerical flowline model, σ _L = 285 m. The time series of length fluctuations from the numerical model and the three-stage model are essentially identical (grey and black lines in Fig. 3), confirming that the three-stage model emulates ice-flow dynamics very effectively.

Fig. 3. Glacier length fluctuations driven by climates with differing persistence. (a) 2 ka time slices of glacier-length anomalies for the five climate time series. The thick curves are the output from the numerical flowline model, the thin curves are calculated from the linear three-stage model of Roe and Baker (Reference Roe and Baker2014). For clarity, the time series have been successively offset by 1000 m. The greater variance for persistent climates is evident by eye. (b) The autocorrelation function of the glacier response. Note the expanded x-axis scale compared with Figure 2, and that the white noise and AR1 curves plot nearly on top of each other. (c) Power spectra of glacier response (windowed periodogram). Thick curves are from the numerical flowline model, thin curves are the theoretical curves from the three-stage model of Roe and Baker (Reference Roe and Baker2014).

Let R be the ratio of the glacier variance with-and-without persistence in mass-balance forcing. In the Appendix we derive expressions from the three-stage model. For persistence based on an AR(1) process, we find:

(9) $${R_r} \equiv \displaystyle{{{{\left. {{\sigma _{\rm L}^2}} \right \vert} _r}} \over {{{\left. {{\sigma _{\rm L}^2}} \right \vert} _{r = 0}}}} = \displaystyle{{(1 - {r^2})} \over {{{(1 - r)}^2}}} \cdot \displaystyle{{{\epsilon} {\tau _{\rm g}}(3{{\epsilon} ^2}\tau _{\rm g}^2 + 9{\epsilon} {\tau _{\rm g}}{\tau _{\rm c}} + 8\tau _{\rm c}^2 )} \over {3({\epsilon} {\tau _{\rm g}} + {\tau _{\rm c}}{)^3}}}.$$

In the limit of  ετ _g ≫ τ _c, this simplifies to:

(10) $${R_r} = \displaystyle{{(1 - {r^2})} \over {{{(1 - r)}^2}}}.$$

For persistence based on a power-law process, we find:

(11) $${R_\nu} \equiv \displaystyle{{{{\left. {\sigma _{\rm L}^2} \right \vert} _\nu}} \over {{{\left. {\sigma _{\rm L}^2} \right \vert} _{\nu = 0}}}} = \displaystyle{{{\pi ^\nu} {{({\epsilon} {\tau _{\rm g}}/\Delta t)}^\nu} (1 - {\nu ^2})(\nu + 3)\sec (\nu \pi /2)} \over 3}.$$

These expressions can be used to calculate the standard deviation of glacier length:

(12) $${\sigma _{\rm L}}{ \vert _{r,\nu}} = \sqrt {{R_{r,\nu}}} \cdot {\sigma _{\rm L}}{ \vert _{wn}},$$

where σ _L| _wn is the response to white-noise forcing Eqn (7) or (8). Equations (9) and (11) allow for efficient computation across a wide range of parameters.

For the mass-balance time series generated by AR(1) processes with r = 0.17 and 0.28, the numerical flowline model responds with σ _L = 336 and 377 m, respectively. These are increases of 18 and 32% over σ _L for white-noise forcing, and Eqn (9) gives good agreement, predicting 17 and 31%, respectively. The length time series for numerical and three-stage models are shown in Figure 3a, and the increased length variability for increasing r can be seen by eye, particularly if you focus on adjacent maxima and minima. Furthermore, we see that the three-stage model emulates the numerical flowline model very well. The ACFs (Fig. 3b) look very similar to the white-noise case, but the spectra show the enhanced power at low frequencies (Fig. 3c). At frequencies <1/100 a⁻¹, the glacier dynamics no longer damp the response, and the glacier exists in near-equilibrium with the climate forcing (e.g. Roe and Baker, Reference Roe and Baker2014).

For the mass-balance time series generated by power-law processes with ν = 0.25 and 0.4, the numerical flowline model responds with σ _L = 399 and 489 m, respectively. These are increases of 40 and 72% over σ _L for white-noise forcing, and Eqn (11) gives good agreement, predicting 43 and 79%, respectively. Again the increase in variance is clear by eye in the time series in Figure 3a, with the three-stage model again accurately reproducing the behavior of the numerical flowline model. In contrast with the AR(1) case, the ACFs for the power-law climate (Fig. 3b), clearly show the extra persistence at large lags. These ACFs and the larger increase in variance are consistent with the spectra (Fig. 3c) that show power continuing to increase towards lower frequencies.

Reichert and others (Reference Reichert, Bengtsson and Oerlemans2002) fit an AR(3) process to a climate-model simulation of mass balance for Nigardsbreen, Norway, and generated an increase in σ _L of 36%. However our results are hard to directly compare because of the different AR(p) models, and the different glacier geometries.

We have deliberately implemented a degree of persistence small enough to be hard to detect in instrumental records, and our analyses show it causes an increase in σ _L of between 17 and 79%, for our chosen glacier parameters.

3.1. The impact of persistence on excursion probabilities

A practical measure of how much persistence matters for understanding past glacier length fluctuations is the impact on the probabilities of the glacier undergoing fluctuations of a given size. Roe and Baker (Reference Roe and Baker2014) showed that the three-stage model could be used to accurately estimate excursion probabilities. We extend that analysis here by incorporating persistence. Under the assumptions that (1) the excursion amplitudes of L are normally distributed, and (2) this distribution is independent of the distribution of $\dot L\ (={\rm d}L/{\rm d}t)$ . The expected return time of a given glacier advance, L ₀, beyond the equilibrium position is given by:

(13) $$R({L_0}) = 2\pi \displaystyle{{{\sigma _{\rm L}}} \over {{\sigma _{\dot {\rm L}}}}}\exp \left[ {\displaystyle{1 \over 2}{{\left( {\displaystyle{{{L_0}} \over {{\sigma _{\rm L}}}}} \right)}^2}} \right].$$

where ${\sigma _{\dot {\rm L}}}$ is the standard deviation of the time rate of change of L. The curves in Figure 3 give everything needed to calculate R(L ₀). Figure 4a shows that the expected return time of an advance past zero (i.e. equilibrium) is 42 a for white-noise forcing, increasing slightly with persistence to 58 a for ν = 0.4. However there is a dramatic reduction in the return time of large advances because Eqn (13) has an exponential sensitivity to σ _L. For instance, for our control glacier, a 1 km advance is expected to occur on average once every 20 ka for white noise, but once every 400 a for ν = 0.4 (Fig. 4a).

Fig. 4. (a) Return times of a given advance, as a function of climatic persistence, calculated from Eqn (13) for the control glacier parameters. Note the logarithmic y-axis. (b) Excursion probabilities as a function of climatic persistence. It shows the probability of exceeding a given total excursion (i.e. maximum minus minimum extent) in any 1 ka period of constant climate. Curves are calculated from Eqn (14) for the control glacier parameters. The presence of even a small degree of persistence can dramatically increase the likelihood of seeing a given excursion.

Another relevant metric is the probability of a given total excursion (i.e. maximum − minimum extent) in a given period of time, say 1 ka. It is useful to evaluate, for example, whether a given moraine might have been deposited simply because of internal climate variability, or whether the moraine must have required a climate change. Such calculations can also be used to correct estimates of past climate change from moraines records, recognizing that the moraines represent extrema in the lengths of fluctuating glaciers, and that the terminus position associated with the average climate during any given period likely lies somewhat up-valley (Anderson and others, Reference Anderson, Roe and Anderson2014; Rowan and others, Reference Rowan2014).

Using standard results from extreme-value statistics (e.g. Vanmarcke, Reference Vanmarcke1983), Roe (Reference Roe2011) derived the probability of a glacier exceeding a maximum total excursion, ΔL, in any given interval of time, ${\rm {\cal T}}$ , assuming the extreme-excursion return times following a Poisson distribution:

(14) $$\eqalign{& p({L_{max}} - {L_{min}} \gt \Delta L)({\rm {\cal T}}) = \cr & \quad \int_0^\infty {\Delta L \over R ({L_1})\sigma _{\rm L}^2} {e^{ - {{\cal T}}/R({L_1})}}\left( {1 - {e^{ - {\rm {\cal T}}/R({L_1}-\Delta L)}}} \right){\rm d}{L_1}.} $$

Figure 4b shows the probabilities of exceeding a given total excursion in any 1 ka period, calculated from the return times shown in Figure 4a. Including persistence has a dramatic impact. For instance, a 2 km excursion has only a 1% chance of happening for white noise, but a 98% chance of happening for ν = 0.4. These results represent a severe complication in evaluating whether a particular moraine on a landscape required a climate change or not, given that the degree of persistence we have used is hard to detect in mass-balance records: that is, the statistical inferences are acutely sensitive to some poorly constrained assumptions.

We note that the Poisson distribution governs events that occur independently, and that this is only an approximation when applied to total glacier excursions (i.e. both maxima and minima) and when climatic persistence is included. In other analyses (not shown) we found the Poisson distribution held well for 1 ka periods for ν < 0.5, but care should be taken applying the formula for shorter periods. Its validity depends on the glacier response time and the degree of climatic persistence. The three-stage model is simple enough that the validity of the assumption can always be checked using Monte Carlo methods.

3.2. Spanning a broad range of parameters

So far, we have used one standard set of glacier parameters and varied only the persistence. Figures 5a and b shows how σ _L varies as a function of glacier timescale, τ _g, and r or ν, using Eqns (7), (9), and  (11); and recognizing that τ _g for our control glacier (~7 a) is quite small relative to values estimated for many other glaciers (e.g. Leclercq and Oerlemans, Reference Leclercq and Oerlemans2012).

Fig. 5. Contours of the standard deviation of glacier response σ _L (m) as a function of glacier response time, τ _g, and climatic persistence (either r or ν). (a) AR(1) persistence, r. (b) Power-law persistence, ν. For these calculations the value of β and σ _b are held constant.

Increasing τ _g increases σ _L (i.e. Eqn 8). A longer glacier response time implies a weaker restoring tendency for a given L′ (e.g. Eqn 6), and so the glacier will undergo larger and longer excursions in response to stochastic forcing. The impact of persistence grows with τ _g, more so for power-law persistence (i.e. Eqn 11) than for AR(1) persistence (i.e. Eqn 10). The reason is that a glacier with a larger τ _g can be thought of as a low-pass filter with a lower-frequency cut-off than a glacier with a smaller τ _g. The differences between the white-noise and power-law spectra increase as frequency decreases, and those lower frequencies contribute more to the variance of the large-τ _g glacier than the small-τ _g glacier. So for our control glacier with τ _g ~ 7 a, going from ν = 0 to ν = 0.4 has a less than two-fold impact on σ _L; however for τ _g ~ 50 a the effect is nearly three-fold (Fig. 5b). Finally, we note that the strong sensitivity at larger values of ν is because of the dependence of σ _L on secant (νπ/2) (i.e. Eqn 11), which reflects the fact that the variance can become unbounded for large enough ν (Eqn A13).

We have also held β and σ _b constant. From Eqn (7) we see that σ _L is linearly proportional to each of them, so our results can easily be applied to other glaciers by reading off the contours in Figure 5. It should also be noted that τ _g and β are related via the glacier thickness, H (see the discussions in Section 3 after Eqn (6)).

4. SUMMARY AND DISCUSSION

Observational records of glacier mass balance are short, and this makes it hard to establish the degree of climatic persistence that pertains to them. Typically, they are short enough that the null hypothesis of no persistence (i.e. white noise) cannot be rejected. And yet there are other observational datasets and also theoretical grounds for being confident that a degree of climatic persistence does, in fact, exist. We have explored the impact of climatic persistence on glacier length fluctuations. For two different statistical models we have found that, even for a degree of persistence that is hard to detect in instrumental records, the magnitude of glacier length fluctuations can be significantly enhanced over that in the case of no persistence. The effect becomes even more important when one considers the impact on excursion probabilities. For extreme excursions, persistence can turn an ‘exceptionally unlikely’ event (<1%) into a ‘virtually certain’ event (>99%, using IPCC (Reference Stocker2013) definitions and a slight extrapolation from Fig. 4b).

Our results compound with the other challenges of identifying past climate change from the glacier length record. If one accepts that internal climate variability will produce stochastic fluctuations of glacier length then, ipso facto, identifying climate change in the glacial record is an exercise in signal-to-noise detection. In turn, such an exercise depends on our physical models for what the glacier ‘noise’ is. For most glaciers, observed mass-balance variability is not constrained to within a factor of two (Section 2.1). The alternative – modeling mass balance from meteorological data – relies on uncertain parameters, and in any event can only ever be calibrated and evaluated against the direct observations. Because the moraine record is intrinsically fragmentary, we must estimate glacier response to internal climate variability using models. In principle, analytical tools like the three-stage model of Roe and Baker (Reference Roe and Baker2014) are efficient vehicles for building intuition about the relative importance of forced and unforced climate variability. Many parameters are accurately observed, but some important parameters like glacier thickness, H, remain poorly constrained. Methodology is improving, but Huss and Farinotti (Reference Huss and Farinotti2012) report a standard error of ~30%. More comprehensive numerical models are not obviously more skillful, depending as they do, on uncertain rheologies and crude parameterizations of subglacial hydrology. In the present study, we have added to this litany of vexation by noting that persistence for glacier mass balance is poorly constrained, while demonstrating it can exert a large influence on the magnitude of glacier length fluctuations. Of course, there is no question but that the Pleistocene Ice Age constituted a climate change that left its mark on the moraine record, a conclusion that survives all modeling uncertainty (e.g. Anderson and others, Reference Anderson, Roe and Anderson2014; Rowan and others, Reference Rowan2014). However the interpretation of Holocene moraines, for which climate changes may be more marginal, is less clear. Events such as the putative Little Ice Age may not stand out relative to the noise of internal climate variability (e.g. Matthews and Briffa, Reference Matthews and Briffa2005).

What to do? In this study we have focussed exclusively on the magnitude of the fluctuations of an individual glacier, and we have made the case that, with existing observations, the minimum uncertainty in the magnitudes of the fluctuations due to internal climate variability remains a factor of two or three (at least until observational records grow long enough to constrain σ _b better). However other information can also be brought to bear on interpreting glacier history. In particular, the spatial extent, or the temporal duration, of contemporaneous glacier length fluctuations can make for an independent test of a climate-change hypothesis. In other words, is the duration or spatial coherence of glacier length fluctuations surprising, given a null hypothesis that they are unconnected? The statistical machinery exists (e.g. Bretherton and others, Reference Bretherton, Widmann, Dymnikov, Wallace and Blad1999), to assess the spatial and temporal correlations in climate and glacier variability. Such tools have yet to be applied to interpreting the glacier record, and an important caveat is that the spatial coherence varies with timescale (e.g. Deser and others, Reference Deser, Alexander, Xie and Phillips2010; Meehl and others, Reference Meehl, Hu, Arblaster, Fasullo and Trenberth2013). However, some principles can be laid out. Glacier length fluctuations within individual mountain ranges will be highly coherent, since they experience essentially the same climate and so do not constitute independent records of climate change (e.g. Huybers and Roe, Reference Huybers and Roe2009). Observations and theory suggest that the degree of climatic persistence is likely greater for maritime settings than for continental settings (e.g. Pelletier, Reference Pelletier1997; Huybers and Curry, Reference Huybers and Curry2006). On larger spatial scales and longer timescales, interpretations of geologic observations hinge critically on the accuracy with which moraine dating can be made, and there is an ongoing and vibrant debate as to what is achievable (e.g. Davis and others, Reference Davis, Menounos and Osborn2009 and references therein; Balco, Reference Balco2009; Schaefer and others, Reference Schaefer2009; Kirkbride and Winkler, Reference Kirkbride and Winkler2012; Schimmelpfennig and others, Reference Schimmelpfennig2012; Solomina and others, Reference Solomina2015).

We have followed in the path of Reichert and others (Reference Reichert, Bengtsson and Oerlemans2002) by including persistence in the climate variability. Does such persistence by itself imply climate changes? In the case of AR(1) persistence, the answer is an emphatic no because the timescales are too short and reflect only, for example, the integrating effects of the ocean mixed layer. For power-law persistence the answer is a little more nuanced. Power-law variability is a constant climate in the sense that it requires neither changes in external forcing, nor changes to the dynamics of the climate system. However variability continues to increase out to lower frequencies and eventually merges into a portion of the spectrum that is most definitely a response to external forcing (e.g. Huybers and Curry, Reference Huybers and Curry2006). The fact the background spectrum of climate variability is a continuum means that what constitutes climate ‘change’ will always depend on somewhat arbitrary definitions.

A related question is whether climate variability is itself a function of mean state. If the climate change is large enough that there are significant shifts of major climatic zones, then there must, at some level, be changes in local variability, especially at the margins of those zones. However for the as yet relatively small climate changes of the last century (~1 K in the global mean), there are no theoretical grounds for expecting that changes in the variability (i.e. the second and higher-order moments of the statistical distribution) have had as big an impact as shifts in the mean (i.e. the first moment). This is largely borne out in observations and models (e.g. Simolo and others, Reference Simolo, Brunetti and Nanni2011; Donat and Alexander, Reference Donat and Alexander2012; Rhines and Huybers, Reference Rhines and Huybers2013). So while possible changes in variability are worth keeping in mind for the larger climate changes of past millennia, the modern instrumental record of variability provides a good foundation for building intuition about the impact of variability in past climates.

The amplitude and nature of internal climate variability is one of the most important challenges in climate science. It governs the predictability of future climate (e.g. Hawkins and Sutton, Reference Hawkins and Sutton2009; Deser and others Reference Deser, Phillips, Bourdette and Teng2012), and the interpretation and statistical significance of past changes, an example of which was our focus here. While we have concentrated on glaciers, the impact of climatic persistence on the response of all paleoclimate proxies is important to assess. There is currently a mismatch between the climatic persistence inferred from paleoclimate proxies and that generated by comprehensive climate models, which may reflect a problem with the proxy records, or the models, or both (Laepple and Huybers, Reference Laepple and Huybers2014). It remains to be seen what the resolution will be.