Glacier response to climate perturbations: an accurate linear geometric model

Abstract

In order to understand the fundamental parameters governing glacier advance and retreat, and also the spectral properties of fluctuations in glacier length in response to noisy weather, we examine outputs of a numerical flowline model solving the shallow-ice equations with sliding. The numerical results reveal a surprising simplicity: the time evolution and spectral shape of glacier excursions depend on a single parameter, a time constant determined by the geometrical properties of the glacier. Furthermore, the numerical results reveal that perturbations in mass balance over the glacier surface set in motion a sequence of events that can be roughly described as occurring in three overlapping stages: (1) changes in interior thickness drive (2) changes in terminus flux, which in turn drive (3) changes in glacier length. A simple, third-order linear differential equation, which extends previous models in the literature, successfully captures these important features of the glacier flow. This three-stage linear model is readily invertible to recover climate history. It provides clear physical insight and analytical expressions for some important metrics of glacier behavior, such as variance, sensitivity and excursion probabilities. Finally, it facilitates uncertainty analysis. The linear model can also be adapted for arbitrary catchment geometry, and is applied to Nigardsbreen, Norway.

1. Introduction

The basic physics of a glacier’s response to small perturbations in climate forcing is of longstanding and intrinsic interest in glaciology, and it is also of some applied importance for predicting future glacier change and for inferring the climate history that drove past glacier fluctuations. The exponential advances in computing power of the past few decades have allowed for efficient numerical integration of the governing equations. Nye (1965), Budd and Jensen (1975) and Oerle-mans (1986) are early examples. Most numerical models solve the shallow-ice equations (which neglect longitudinal stresses; e.g. Hutter, 1983) and incorporate a representation of basal sliding. For a one-dimensional flowline following the longitudinal profile of a glacier, a standard version of these equations can be written as

(1a)

(1b)

where h (x) is glacier thickness at position x; F (x) is the vertically integrated flux of ice, dz_s= dx is the surface slope, ḃ(X) is the local mass balance, ρ is the density and f _d and f _s are dynamical coefficients governing ice deformation and Weertman-style ice sliding, respectively (Budd and others, 1979; Oerlemans, 2001). The first equation represents local mass conservation, while the second represents the translation and deformation of ice associated with shear stresses. In combination, the equations have the form of a nonlinear diffusion equation in thickness. While the appropriateness of neglecting longitudinal stresses is sometimes debated, it has been shown that for parameters appropriate for typical alpine glaciers, the solutions to Eqns (1) closely match those from numerical models solving the full stress field (e.g. Leysinger Vieli and Gudmundsson, 2004).

A second strand of research has striven for simpler descriptions of glacier response to climate (e.g. Jóhannesson and others, 1989, hereafter JRW89; Oerlemans, 2001; Harrison and others, 2003; Lüthi, 2009; Roe and O’Neal, 2009). The reasons for doing this are several-fold. Most importantly, one gains physical insight into the glacier’s response. Simple models also provide understanding of how the numerical solutions depend on glacier geometry, and can be used to estimate the response of more glaciers and across a wider range of parameter uncertainty than is computationally tractable with numerical models. Finally, such models can be used to try to recover the climate history from glacier-length records (e.g. Klok and Oerlemans, 2004; Oerlemans, 2005; Lüthi, 2013).

We seek the simplest description that captures the time-dependent glacier response to climate change as well as the spectral characteristics of the glacier fluctuations that result from noisy weather. It is natural to begin with the classical ‘one-stage’ linear model (JRW89; Oerlemans, 2001; Roe and O’Neal, 2009), which is described in the next section. We will compare the results from this model with the output of the numerical flowline model (Eqns (1)) to identify the physical shortcomings of the one-stage model, and to develop a more satisfactory extension. For step changes to a new equilibrium climate and for long-term climate trends, the glacier response is determined by its geometry. On shorter timescales, glacier dynamics plays an important role. Adjustment proceeds in three overlapping stages: (1) changes in interior thickness drive, (2) changes in flux at the terminus, which in turn drive (3) changes in glacier length. We show that simple representations of each of these three stages can be combined into one third-order linear differential equation for glacier length that successfully emulates the behavior of the numerical model across a wide range of parameters.

2. One-Stage Models

Roe and O’Neal (2009) derived a model for how perturbations in glacier length, L ′, linearized about a mean length, , are driven by changes in melt-season temperature, T’, and annual-mean precipitation, P’. For a simplified geometry (constant thickness, terminus width and bed slope), and taking ablation proportional to melt-season temperature, a linearization of the mass-conservation equation yields

(2)

where

(3a)

(3b)

(3c)

H is the thickness, w is the terminus width, tan ϕ is basal slope; A _tot, A_{T> 0} and A _abl are the total area, the melt area (where the melt-season temperature exceeds 0°C) and the ablation area where there is net melting (i.e. below the equilibrium-line altitude (ELA)), respectively. Γ is atmospheric lapse rate and μ is the melt factor, i.e. the proportionality constant relating ablation to melt-season temperature.

Equations similar to Eqn (2) have been developed, or applied, in many other studies. JRW89 made a geometric argument that where is the net ablation rate at the glacier terminus. For a glacier of uniform width on a constant slope, it can be shown that ) (see Appendix), in which case the Roe and O’Neal (2009) and JRW89 timescales are identical. JRW89’s expression for τ is widely used to estimate the characteristic ‘memory’ or response time of glaciers, and has been evaluated against the output of numerical glacier models (e.g. JRW89; Oerlemans, 2001; Leclercq and Oerlemans, 2012). Such studies often report longer response times than suggested by JRW89. In the present study we demonstrate how to reconcile this apparent discrepancy: τ is indeed the correct timescale, but it must not be misinterpreted as an e-folding timescale. Raper and others (1996, 2000) incorporate glacier hypsometry, thickness and width, and Bahr and others (1998) and Harrison and others (2001) include feedbacks between glacier shape and mass balance. Other variants and extensions include those of Harrison and others (2003), Lüthi (2009) and Harrison (2013). Echoing Nye (e.g. 1960, 1963, 1965; which are summarized by Paterson, 1981), Oerlemans (2001) takes a different, semi-empirical approach and defines τ ~ /u, where u is a characteristic velocity. Whereas Nye took u to be the speed of kinematic surface waves at the terminus, Oerlemans relates u to mass flux and, via a scale analysis, to simple functions of glacier geometry, then finally calibrates it to the output of numerical models.

One appeal of a linear equation, such as Eqn (2), is that glacier length history can be easily inverted to recover climate history. Oerlemans and colleagues have used a version of it to make global temperature reconstructions (Klok and Oerlemans, 2004; Oerlemans, 2005; Leclercq and Oerlemans, 2012). Although the functional form is the same as earlier work, the specific model given by Eqns (2) and (3) offers a few advantages: it is formally derived from a mass-conservation equation; it distinguishes between climate forcing due to melt-season temperature and that due to precipitation; and the model parameters are cleanly and clearly related to glacier geometry. However, all glacier models of the form of Eqn (2) share a single dynamical stage: any changes in the mass balance on the glacier are immediately converted to a tendency on the terminus length. For this reason we refer to them as ‘one-stage’ models.

There are some simple solutions for Eqn (2). In the case of no climate anomalies (P ′, T ′ = 0), an initial length perturbation, L’(0), decays as L ′(t) exp (-t/τ). For the one-stage model then, τ can thus be identified as the ‘e-folding’, characteristic timescale over which the glacier ‘remembers’ previous states (although, as we show below, the concept of an e-folding timescale is in fact incorrect for glaciers). This memory is also expressed in the autocorrelation function (ACF), which is the correlation coefficient of the time series with itself lagged by a time interval, γ (ACF(γ) <L’(t)L’(t + γ)>/<L’(t)L’(t)>; e.g. Box and others, 2008):

(4)

Starting from equilibrium, if step-functions of ΔT and ΔP are applied at t ¼ 0, then the solution is

(5)

where

(6)

is the equilibrium length change. If trends in melt-season temperature, Ṫ, and precipitation, Ṗ, are applied at t = 0, then

(7)

For t ≫ τ this solution asymptotes to a straight line with a slope of τ(α and an intercept of τ on the time axis.

Equation (2) also discretizes naturally into time steps of Δt = 1 year:

(8)

Roe and O'Neal (2009) and Roe (2011) focused on the glacier response to stochastic climate variability: random year-to-year fluctuations in T′ and P′ that occur due to the vagaries of weather, even in a constant climate. Building on the earlier work of Oerlemans (2000) and Reichert and others (2002), they demonstrated that century-scale, kilometer-scale glacier fluctuations could be driven by such interannual climate variability. If the variability in T′ and P′ is Gaussian white noise (i.e. uncorrelated and normally distributed), with standard deviations of σ _T and σ _P, then the variance of glacier length is given by

(9)

Thus, glacier sensitivity to a step change in climate (Eqn (6)), glacier response to climate trends (Eqn (7)), and glacier variance driven by stochastic climate fluctuations (Eqn (9)) are all proportional to τ, making τ an important number to constrain. Moreover these three aspects of glacier behavior are inextricably interwoven: a high sensitivity to climate change goes hand-in-hand with a large natural variability. This is a general property of dynamical systems and also applies to the three-stage model developed below.

The power spectrum for the discrete length equation (Eqn (8)) is given by a standard formula (e.g. Box and others, 2008), which applies for frequencies 0 f 1/ (2Δt):

(10a)

(10b)

(10c)

H(f) is the complex frequency-response function, from which the phase of the response can be calculated, and P₀ is the amplitude of the power spectrum in the limit f → 0. Equations (10) can be combined to yield

(11)

It can be shown that σ_L ² = ʃ P(f) df, for the limit τ ≫ Δt.

2.1. Performance of a one-stage model

We first compare the performance of the one-stage model (Eqn (8)) with a numerical model of a flowline (Eqns (1)) on a constant slope using standard numerical techniques. Grid resolutions of between 10 and 100 m were used, though results presented are insensitive to this choice. The control-case glacier uses all the same parameters as Roe (2011), provided in Table 1, and which are representative of the glaciers around Mount Baker, Washington, USA. Additional parameters for the one-stage model are matched to the numerical model in equilibrium (Table 1), leading to α = –99: 5 m a⁻¹ °C⁻¹, β = 177 and τ = 6.73 years. We evaluate the responses of the two models for two climate-forcing scenarios: (1) step-function accumulation changes of ΔP = ±0: 5 m a⁻¹ (Fig. 1a) and (2) random, uncorrelated, normally distributed, interannual fluctuations in T ′ and P ′, with σ _T = 0: 8°C and σ _P = 1: 0 m a⁻¹ (consistent with observations; e.g. Roe and O’Neal, 2009). For the fluctuating climate we run each model for 10 000 years, and compare the models’ time series, ACFs and power spectra (Fig. 1b–d). The one-stage model responds too quickly to the step-function precipitation increase (Fig. 1a), reaching (1–1= e) of its equilibrium at t = τ = 6: 7 years, whereas for the flowline model this happens at t 15 years. However, the flowline model has the classic ‘tulip’, or sigmoidal, shape seen in many studies (e.g. JRW89; Oerlemans, 1997; Leysinger Vieli and Gudmundsson, 2004) and catches up quickly: by 20 years the flowline and one-stage models have converged quite closely and the final equilibrium lengths agree to within 5%. Roe (2011) showed similar agreement between the models for a wide range of step-function amplitudes, both melt-season temperature and precipitation, and also good agreement for the response to climate trends of the magnitude typical in observations.

Table 1. Parameters and geometry of control-case glacier. The first group of parameters are imposed, the second group are calculated from the flowline model (mean thickness is used for H) and used for the one-stage and three-stage model formulae. The simplified, pseudo-one-dimensional geometry means that not every aspect of the typical Mount Baker glacier can be matched at the same time. In particular, the standard glacier has a nominal length of 8 km, and the accumulation-area ratio is one-half, rather than two-thirds. Compare with values given by Roe and O’Neal (2009) for Mount Baker glaciers, and Oerlemans (2001) for flowline parameters

Fig. 1. Comparison of the one-stage model to the numerical flowline model. (a) Response to step-functions of ±0: 5 m a 1 in precipitation. (b) 500 year sample from a 10 000 year integration with stochastic, white-noise climate forcing of amplitude σ _T = 0: 8°C and _P = 1: 0 m a⁻¹. (c) Autocorrelation function, calculated from (b). (d) Power spectrum calculated from (b) using a modified periodogram (Hamming window with 16 segments overlapping by 50%, used for all spectra shown hereafter).

Fig. 2. (a, b) Power spectra and ACF for flowline-model glaciers of three different sizes, forced by stochastic white-noise climate variability. (c, d) Power spectra and ACF normalized by the τ and σ _L from the one-stage model for each glacier. The thick green curve is the ACF for the three-stage model, also normalized by τ.

In response to stochastic climate variability, the standard deviation in glacier length, σ _L, is 323 m in the flowline model and 361 m in the linear model – an overestimate of 12%. Figure 1b shows a 500 year segment of the 10 000 year integration. It is clear that the fluctuations of the flowline model are captured by the one-stage model, but that the one-stage model also has more high-frequency variability. These differences are manifest in the ACF (Fig. 1c), where the one-stage model is much less autocorrelated at short lags than the flowline model, and also in the power spectrum (Fig. 1d), where the one-stage model has more power at high frequencies. However the two models agree extremely well at low frequencies. As described by Roe (2011), in the low-frequency limit a glacier acts as a passive reservoir of ice in quasi-equilibrium with the climate forcing. The glacier's response is then dictated by its geometry, which is exactly what the one-stage model represents.

The step-function response, the ACF and the power spectrum of all one-stage models are constrained to the functional shapes seen in Figure 1a–d, and given in Eqns (5), (4) and (10). It is therefore clear that no one-stage model can characterize glacier variability on all timescales.

2.2. Glacier dynamics are governed by a single fundamental timescale

We have shown that the one-stage model gives good predictions for the equilibrium, the variance and the low-frequency response of glacier length, and is thus capturing some important aspects of the flowline behavior. Since these metrics are a function of τ in the one-stage model, we are motivated to explore whether the glacier behavior scales with τ more generally. We first vary glacier size in the flowline model by varying the basal slope, and take tan ϕ = 0.4, 0.2 and 0.1. This gives equilibrium lengths of = 8.0, 16.6 and 35.0 km and mean thicknesses of H = 44, 104 and 220 m. From Eqn (3c) for the one-stage model, this yields timescales of τ = 6: 7, 15.4 and 26.2 years. Figure 2a and b show the power spectra and ACF for these three glaciers forced by the same stochastic climate variability as in the previous section. The standard deviations of glacier length are σ _L = 323, 419 and 552 m, respectively. For the larger glaciers, the greater variance can be seen in Figure 2a, and the longer timescales are evident in Figure 2a and b.

We next normalize the power spectra and ACF using the τ and P ₀ predicted by the one-stage model (Eqns (3c) and (10c)), shown in Figure 2c and d. Strikingly, the spectra and ACFs now all plot on top of each other. The exception is at high frequencies (f > 1/τ), where the spectral power is very low and influenced by the model resolution.

The panels in Figure 3 show the same experiment as Figure 2, but for wide variations of the dynamics parameters in the flowline model (Eqns (1)). Again, when normalized by τ and P ₀ the curves plot almost on top of each other (the least agreement is found for a tenfold increase in f _s, likely due to the unphysical ice thickness of only 20 m that this produces). The same convergence of curves after normalizing by τ and P ₀ was obtained when varying the magnitude of the climate forcing by a factor of 40 (not shown).

Fig. 3. (a, b) Power spectra and ACF for four flowline-model glaciers with different dynamical parameters (see legend). (c, d) Spectra and ACF normalized by the for each glacier. The thick green curve is the ACF for the three-stage model, also normalized by τ.

These results are important. The fact that when the curves are normalized by the one-stage model (Eqn (3c)), they all collapse to approximately the same function demonstrates: (1) that the fundamental behavior of the flowline model remains the same for a very wide range of parameters; (2) that there is one underlying timescale for each glacier; and (3) that it is the τ of the one-stage model.

That there is only one timescale, and that it is defined by the glacier’s geometry, can be viewed as a consequence of scale invariance: the ice dynamics depends only on glacier shape (Eqn (1)), and the glacier shape varies self-similarly on a constant slope (e.g. Nye, 1951; JRW89). Therefore, with no differentiation in the glacier physics as a function of size, the functional shape of the temporal evolution of the glacier should also be independent of size. Consequently we are motivated to find a simple physical model that captures this evolution.

3. A Three-Stage Model

An indication of why the one-stage model fails at high frequencies can be seen from the response of the flowline model to a step-function increase in precipitation, ΔP = 0: 5 m a⁻¹, applied at t = 0. The dark curves in Figure 4 show the ice-thickness profile as a function of time, in increments of τ. After 1τ the portion of the flowline model upslope of the initial terminus is well on its way to equilibrium: the thickness in the middle of the glacier has reached 87% of its final value. In contrast, in that same time the terminus has only advanced 20% towards equilibrium (Fig. 1). After 2τ, the interior thickness is at 99.7%, and the terminus advance is 68% of equilibrium. After 3τ, the values are 100% and 92%, respectively. Thus, the interior approaches equilibrium relatively rapidly compared with the more slowly adjusting terminus. As noted above, one-stage models neglect this interior filling and presume that any mass imbalance acts immediately to change the terminus position. The lighter curves in Figure 4 show the profile for a decrease in precipitation of the same amount, and show an essentially symmetric behavior: rapid interior adjustment, slower length response.

Fig. 4. Flowline model glacier thickness, shown for time increments of τ, in response to a step-function increase (darker curves) and decrease (lighter curves) in precipitation, ΔP = ±0: 5 m a⁻¹.

The approach to equilibrium for a glacier advance can also be characterized in terms of the evolution of the mass fluxes beyond the initial terminus. There is a source flux due to the increased precipitation, and a sink flux due to extra melting as the glacier grows beyond its original terminus. The net flux, which is positive for this perturbation, allows the glacier to grow towards a new equilibrium, at which point the fluxes must ultimately come back into balance. Figure 5 compares the evolution of these fluxes for the one-stage and flowline models. In the one-stage model, the source flux instantly rises to its new equilibrium value at t = 0 (= ΔP · L ₀), whereas the sink flux is zero at t = 0 (since L’ = 0). Thus the net flux is initially large and positive, leading to rapid growth (Fig. 5). As time goes on, the sink flux gradually increases as the glacier terminus extends into higher temperatures, and so the fluxes return to a balance as the glacier exponentially asymptotes to equilibrium. By contrast, for the flowline model, at t = 0 both the source flux and the sink flux are zero. The source flux ramps up only slowly. The growing source flux does drive the glacier forward, and this growth generates a sink flux that lags slightly, leaving a small residual net flux that peaks around t = 1: 5 before declining towards zero at around 3τ. At 1, 2 and 3τ the source flux is at 53%, 96% and 100%, respectively, of its equilibrium value. In other words, the source flux lags the interior thickness but leads the terminus position. There is thus a sequence of three overlapping stages: interior thickening drives a flux past the original terminus, which in turn drives the glacier terminus forward.

Fig. 5. Anomalous fluxes past the initial glacier terminus for one-stage (thin curves) and flowline (thick curves) models, in response to the same step increase, ΔP = +0: 5 m a⁻¹.

3.1. Development of a three-stage model

A heuristic representation of the three stages of adjustment is suggested by considering the change in glacier geometry, illustrated schematically in Figure 6. This heuristic representation can be developed into a dynamical model which, we will show, effectively accounts for the response of the numerical flowline model. We loosely define a ‘terminus zone’ of length, with a (strongly negative) mass balance given by ḃ_term . In the interior, let the glacier have a characteristic thickness H. Support for the argument that the interior and terminus zones can be separated in this way comes from the linearized solutions to equations like Eqn (1) presented by Nye (1960) and JRW89: diffusion and advection are much more efficient in the interior than near the terminus, consistent with the profile development shown in Figure 4.

Fig. 6. Schematic illustration of the transition from an initial profile of length to a new larger profile of length + L ′. The three sectors (labelled 1, 2 and 3) that contribute volume changes, V ’_1,2,3, and the fluxes between them are shown. The length of the terminus zone is Λ. The volumes of ice within the terminus zones of the initial and new profiles are assumed to be the same.

Consider the transition from an initial profile to a new, larger profile in response to a small increase in mass balance, , assumed uniform over the glacier surface. Three volumes can be identified, V ′_{1, 2, 3} (m²), whose sum is the change in volume due to . From Figure 6, and with further assumptions that Λ ≪ and h ′₂, h ′₃′ H, we have V ′₁ ≃ h ′₁ , V ′₂ ′ h ′₂ ˄ and V ′₃ = (L ′ – ˄)H. We assume the profile of the terminus does not change substantially (e.g. JRW89), and V ′_1, 2, 3 get incorporated into the interior of the new glacier. The thickening of the glacier drives anomalous fluxes, F′ ₁ and F ′₂ (m² a⁻¹). In the spirit of Nye (1960) and JRW89, we seek first-order solutions assuming that the perturbations in mass balance, interior thickness and length are all much smaller in magnitude than the equilibrium values for the initial profile. Mass conservation in each section gives

(12a)

(12b)

(12c)

Adding the three equations confirms that . One-stage models typically neglect V ’₁ and V ’₂ and set V ’_tot = HL ’, which then gives the one-stage equation (i.e. Eqn (2)). However, the interior changes in thickness cannot be neglected. Suppose that the anomalous fluxes are proportional to the anomalous volume in each segment:

(13a)

(13b)

where τ ₁ and τ ₂ are coefficients of proportionality, and represent the timescales on which an excess of volume accumulating in either segment drives a downslope flux of ice. The assumption of proportionality is not rigorous, but could be loosely rationalized as a linearization of the ice flux with respect to thickness.

The above relationships can be written as a system of three coupled equations for the sequence of changes in interior thickness, terminus flux and length, which we term the three-stage model:

(14a)

(14b)

(14c)

Equations (14) qualitatively characterize the three overlapping stages of glacier adjustment. We established in the previous section that there is only one underlying timescale for a glacier on a constant bed, and therefore τ ₁ = τ ₂ and τ ₃ cannot be independent. We set τ ₁ = τ ₂ = τ ₃ = ϵτ. In part, this is for simplicity as it produces clean expressions for metrics such as the variance and the power spectrum of the glacier response; in part it is supported by JRW89 who demonstrate τ ₁ τ ₃, if the profile changes self-similarly; finally, it presumes the volume–flux relationship is the same in each segment (Eqns (13) with τ ₁ τ ₂). This choice also requires an ad hoc multiplication of the right-hand side of Eqn (14c) by 1/ϵ to retain the accurate equilibrium solution of the one-stage model, so Eqn (14c) becomes

In Section 4 we show that ϵ = 1/√3 is a natural choice. In the limit t → ∞ (and writing h’₁ → h’; F’₂ → F’),h’, F’ and L’ all asymptote to well-defined values: . These equilibrium changes are all fundamentally constrained by the glacier geometry: we saw above that the one-stage model gives a good prediction for L ′_eq; the self-similar shape of the glacier profile then sets h ′_eq; and F ′_eq must accommodate the imposed mass-balance anomaly over the initial length.

For a step-function forcing, we compare the evolution of h ′, F ′ and L ′ in the flowline model with that predicted by integrating the discretized versions of the three-stage model (i.e. , etc.). Results are shown in Figure 7. No simple geometric approximation can capture every detail of the evolution of a nonlinear diffusion equation; for instance, the flowline-model thickness asymptotes to equilibrium faster, and the flowline-terminus flux is slower to begin adjusting, and then equilibrates more rapidly, than is predicted by the three-stage model. However the overall sequencing of h ′, F ′ and L ′ is well captured, as are the details of the sigmoidal evolution of L ′. Figure 7 confirms the basic order of events suggested in Figure 5, and the ability of the three-stage model to represent it.

Fig. 7. The sequential response of H, F and L to a step increase in precipitation for the flowline (solid) and three-stage (dashed) models, as a function of t/τ. For the flowline model, H is the mean thickness and F is the flux past the original terminus. Variables are plotted normalized to their final equilibrium values.

Although the three-stage model has been developed considering a glacier advance, the evolution of the numerical flowline is very nearly symmetric with respect to retreat (e.g. Fig. 1a): reduced precipitation first thins the glacier, which then reduces the flux into the terminus zone, and the terminus subsequently pulls back. Temperature anomalies of either sign produce analogous sequences of events.

Eliminating h ′ and F ′ from Eqns (14) yields a single, third-order linear equation for L ′:

(15)

or, more concisely,

(16)

The discretized versions of Eqns (14) can also be combined into one equation for the dependence of L ′ _t on previous states and on climate forcing:

(17)

where k (1 – Δt/ϵτ). Equation (17) has the form of an autoregressive moving-average process (ARMA; e.g. Box and others, 2008), in which the state variable, here L ′, depends not only on its own values at previous times but also on the forcing at previous times. An ARMA(p, q) model is defined by

(18)

where a_i and c_j are coefficients, and C ’ _t is the stochastic forcing. Equation (17) is thus an ARMA(3,3) model with a_i = {3k, – 3k², k³g and c_j = {0, 0, 0, 1}. The statistical fitting of ARMA(p, q) models to time series of data is widely used in many fields to reveal the underlying dynamical processes. Here it has been derived from a physical model. Standard relationships have long been established for the dynamical properties of ARMA(p, q) models (e.g. Jenkins and Watts, 1968; Box and others, 2008), which we adapt for the three-stage model:

1. Power spectrum

For white-noise climate forcing, the spectrum of the discrete three-stage model can be written as

(19a)

(19b)

As before, the phase of the response can be derived from H(f), and P₀ is the spectral power in the limit f → 0. We set P₀, which is an overall scale, by demanding that the three-stage model give the same low-frequency behavior as the flowline model and the one-stage model: P₀ = 4τ (σ ² _L)_1s, where the subscript has now been added to make it clear it refers to the variance of the one-stage model (which is given by Eqn (9)). Eqns (19) simplify to

(20)

Fig. 8. Comparison of the flowline, three-stage and one-stage glacier models for (a) step-function forcings of ± 0.5 m a 1, (b) ACF, (c) power spectrum and (d) phase, for a 10⁴ year integration driven by stochastic climate variability.

Fig. 9. Time series of the flowline, three-stage and one-stage glacier models’ response to stochastic climate forcing with standard parameters, but varying the basal slope: (a) tan ϕ = 0. 4, (b) tan ϕ = 0. 2 and (c) tan ϕ = 0. 1. The equilibrium glacier profiles are shown in the inset panels.

2. Variance

The variance of a process can be obtained by integrating P(f) over all frequencies (0 f ≤ 1/2Δt). For Eqn (20) there is an analytic solution:

(21)

For the standard glacier and the ratio of the variances in the three-stage and one-stage models (i.e. Eqns (21) and (9)) is 0.76. Equivalently, σ_L|_3s is ~ 14% smaller than σ_L|_31s. This correction is only weakly dependent on τ and asymptotes to 19% in the limit of large τ.

3. Autocorrelation function

The ACF for a discrete ARMA(p, q) model can be calculated from a set of equations known as the modified Yule–Walker equations (e.g. Box and others, 2008), but no closed-form solutions exist. An explicit expression for the ACF can, however, be calculated for the continuous form of the three-stage model (Eqn (15)):

(22)

This function closely matches the ACF calculated from the discrete three-stage model: the areas underneath the ACF curves differ by <10%.

4. Degrees of freedom in a glacier record

It is important to know the number of degrees of freedom (i.e. independent pieces of information) in a time series when performing statistical tests to establish, for example, the significance of a trend. It is closely related to the ACF: more persistence equals fewer degrees of freedom. For a time series of n years of annual observations of glacier length, the degrees of freedom, n ’, is given by

(23)

based on matching the variance with a time series of independent observations (e.g. von Storch and Zwiers, 1999). Using Eqns (4) and (22), we find

(24)

(25)

For τ≫1 and there are ~ 35% fewer degrees of freedom in the three-stage model. For example, for our standard glacier with τ = 6.7 years, and a 100 year record, there are 6.9 degrees of freedom in the one-stage model, and only 4.6 in the three-stage model. Both numbers are small for the purposes of statistical detection, and this reinforces the point that century-scale glacier-length records have a very limited ability to resolve climate trends, compared with either glacier mass-balance records or instrumental climate data (e.g. Roe, 2011).

3.2. Performance of the three-stage model

If climate time series of T ^’ _t and P ^’ _t are imposed on the three-stage model, L ^’ _t may be found by directly integrating Eqn (17), or by applying H(f) to T ^’ _t and P ^’ _t using the filter function in MATLAB^®, for example, and using Eqn (21) to set the variance. For the same step-functions and 10⁴ year stochastic climate time series used to generate Figure 1, and for standard parameters (Table 1), we evaluate the ability of the three-stage model to reproduce the behavior of the numerical flowline model in Figures 8 and 9a. The results from the one-stage model are also shown for reference.

Fig. 10. For the three glacier parameter sets used in Figure 9: (a) average return time of an advance, as a function of the size of the advance beyond equilibrium; and (b) the chance of exceeding a given total excursion (i.e. maximum minus minimum) in a 1000 year period, as a function of the excursion size. The curves show the predictions from Eqns (29) and (30), and the symbols show results calculated from long (10⁵ year) simulations of the numerical flowline model forced by stochastic climate variability.

The three-stage model captures the sigmoidal evolution of the flowline model response to a step function very well (Fig. 8a), and the equilibrium amplitudes differ by <5% (by construction, equilibrium changes for the three-stage model are identical to the one-stage model). There is an indication of a slight nonlinearity of the flowline model, evident in the small difference between the step-increase and decrease in precipitation. The three-stage model also does an excellent job of emulating the flowline model response to stochastic climate variability. Figure 9a shows the glacier-length time series. The match is almost exact, and there is none of the excess high-frequency variability of the one-stage model, also shown for reference. The σ _L for the three-stage model is 314 m, which compares very well with the 323 m for the flowline model. The excellent match is also reflected in the ACF, the power spectrum and the phase (Fig. 8b–d). The convexity of the flowline-model ACF at short lags is also seen in the three-stage model, as is the rollover (i.e. decrease) in the power spectrum and phase near f 1/ 50 years. The inclusion of a lag in the climate forcing in Eqn (17) (i.e. using an ARMA(3,3) model) is important for capturing the phase at high frequencies. For an ARMA(3,0) model the phase would turn back towards zero, as is also seen for the one-stage model. The amplitude and phase of the three-stage model diverges from the flowline model at high frequency (f > 1= 10 a⁻¹), but the spectral power is so low as to not matter for the time series.

Figure 9b and c show that the three-stage model performs equally well for larger glaciers. For tan ϕ = 0.2, σ _L = 396 m for the three-stage model and 419 m for the flowline model, and for tan ϕ = 0. 1 these values are 481 m and 552 m, respectively.

4. Glacier Excursion Statistics

We have seen that the glacier length response to stochastic climate fluctuations is characterized by a probability distribution, with width σ _L. However as a practical matter, historical observations are never long enough to adequately sample the full distribution, and only the maximum extents of glaciers are recorded in moraine sequences. Thus the questions of practical relevance are: What is the expected return time of a given advance? How likely is it that in n years we would see a total excursion of x km?

Standard results exist to answer these questions (Rice, 1948; Vanmarcke, 1983), and were employed for the one-stage model by Roe (2011). The expected rate of upcrossings over a given length threshold, L ₀, is given by

(26)

where σ _L and are the standard deviations of the glacier length and its time-rate-of-change (i.e. dL= dt), respectively. The average return time of an L ₀ advance (R (L ₀)) is the reciprocal of λ. Roe (2011) derived an expression for /σ_L for the one-stage model, but required a timescale much longer than τ to produce agreement with flowline models. Here we improve on Roe (2011), using the statistical relation

(27)

and Eqn (22) for the analytical functional form for ACF(t), to give

(28)

Setting yields simply = σ_L/τ, and thus the average return time for an advance of size L ₀ is given by

(29)

Equation (29) supersedes Eqn (5) of Roe (2011). The average time between readvances past equilibrium is R (L ₀₎ = 0) = 2, or ~42 years for standard parameters. Figure 10a shows that Eqn (29) accurately captures the return times of upcrossings for the three glaciers in Figure 9. Note the logarithmic y –axis. For standard parameters the average return time of a 500 m advance is ~ 130 years; for a 1 km advance it balloons to 6000 years.

Roe (2011) also derived a formula for the probability of a glacier exceeding a maximum total excursion (i.e. maximum minus minimum positions) exceeding L in any given interval of time, T, for extreme events following a Poisson distribution:

(30)

Thus, based only on the glacier geometry, the likelihood of excursions can be characterized for a wide range of glacier 4sizes. Figure 10b shows the good agreement between Eqn (30) and calculations from the numerical model. For the standard parameters, in any 1000 year period there is a 95% chance of finding a total excursion exceeding 1400 m, but only a 5% chance of an excursion exceeding 2100 m. Larger excursions are more likely for the bigger glaciers. Roe (2011) demonstrated that the excursion probabilities are quite sensitive to parameter uncertainty, which must be accounted for when interpreting the climate represented by past moraines (e.g. Anderson and others, 2014).

5. A Case Study: Nigardsbreen

Finally, we present a case study for Nigardsbreen, Norway, arguably the most intensely observed and analyzed glacier anywhere. The results are not intended as an optimized simulation; rather, we want to determine how well the three-stage model performs relative to a flowline model when realistic basal geometry and variable valley widths are included. The geometry is shown in Figure 11a and b (Oerlemans, 1986). The numerical model is adapted for variable flowline width, following Oerlemans (1986). For the climate variability we take σ _T = 0.9°C, σ_p = 0.7 m a⁻¹. This gives summertime, wintertime and annual-mean standard deviations in mass balance of 0.59, 0.71 and 0.92 m a^-1, respectively, which compares well with 0.60, 0.61 and 1.02 m a⁻¹, respectively, from 49 years of observations, as reported by the World Glacier Monitoring Service database (WGMS, 2012). Although the mass-balance parameterization does not account for changing length of the melt season with altitude on the glacier, the modeled ELA is 1470 ± 210(1σ)m, compared with 1500 ± 164(1σ) m from observations (WGMS, 2012).

Fig. 11. Comparison of the flowline, three-stage and one-stage models for the geometry of Nigardsbreen, driven by 10 ka of stochastic climate variability with magnitude based on observations. (a) The basal topography and (b) the flowline width (Oerlemans, 1986); (a) also shows equilibrium flowline glacier profiles for = 11 and 14 km (dashed curve). (c–e) ACF, power spectrum and a 5 ka interval of the time series for the three models. (f) For completeness, the observed length variations of Nigardsbreen (Leclercq and others, 2013), scaled to the same axes as (e); tickmarks are 50 years. Note that (e) is not intended to be a direct simulation of (f).

The linearization leading to the one-stage and three-stage models requires that the formulae for α, β, γ, and τ are amended to account for variable valley width. These equations are presented in the Appendix. We linearize around a mean glacier position of 11 km, and from the flowline model we diagnose H = 180 m from the maximum thickness in the lower reach of the equilibrium glacier (JRW89). This then yields 227 m a^-1C⁻¹, β = 350 and τ = 44 years, with all other parameters as before.

We integrate the flowline, three-stage and one-stage models for 10 000 years, forced by stochastic white noise in P ^’ and T ^’, whose standard deviations are given above. From these integrations, we find the value ofσ _L for the flowline model is 1063 m, vs 1222 m for the three-stage model and 1501 m for the one-stage model. Figure 11d and e shows that, as expected, the one-stage model has excessive high-frequency variability.

The three-stage model does a good job of capturing the flowline model behavior, in particular the shapes of the ACF (Fig. 11c) and power spectrum (Fig. 11d), which suggests that τ is the correct timescale. The differences are also instructive. The 15% overestimate in σ _L is mainly because the three-stage model generates larger advances than the flowline model (Fig. 11e). The reason for this is a distinct shallowing of the basal slope beyond 11 km. Advances on these shallower slopes require a significant thickening of ice (Fig. 11a). This thickening is supplied by an ice flux that, in the three-stage model, would have fueled a greater advance. Thus Nigardsbreen is inherently asymmetric to advance and retreat, which no linear model can fully capture.

Although we again emphasize that we are not trying to optimize a simulation of the natural variability of Nigards-breen, the historical observations and reconstruction of terminus position since the 17th century (Leclercq and others, 2013) are shown in Figure 11f. The variable basal geometry means that the value of σ _L is quite sensitive to the assumed equilibrium position. For example, Oerlemans (2000) used ≃ 14 km, near the historical 18th-century maximum, and found σ _L 610 m. For the same we find σ_L = 650 m (not shown).

The choices made for and H in the three-stage model are reasonable, but are in part a calibration to the flowline model, and so were not made a priori, nor are they determined easily from observations. For a different basal geometry the mean, rather than the maximum, thickness might be more appropriate. The main point for this study is that the basic dynamical structure of the flowline model of Nigardsbreen is expressed in the ACF and power spectrum, and that these dynamics are well characterized by the three-stage model in a way that is not possible using a one-stage model. Some caution is warranted for highly variable basal slopes because self-similarity of the glacier profile can no longer be assumed for large fluctuations. Under such circumstances the response of any glacier model is sensitive to the assumed equilibrium position.

While not a focus of the present study, an area of future work in detailed simulations of individual glaciers is a comparison of different treatments of mass balance. For a mass balance based on positive degree-days (e.g. Braithwaite, 1984; Reeh, 1991) and daily observations from nearby stations (not shown), we found σ_L = 830 m compared with 1050 m for the melt-factor treatment used here. This is likely due to a difference in the pattern of how stochastic forcings in temperature and precipitation are felt on the glacier.

6. Summary and Discussion

At low frequencies and long timescales, a glacier responds in quasi-equilibrium with the climate forcing. The glacier’s length is dictated by the geometry it must attain to achieve a quasi-balanced mass budget, and the ice dynamics primarily affects glacier thickness. At high frequencies and short timescales, the ice dynamics does play an important role, the salient aspect of which is that mass redistribution is much more efficient in the interior than near the terminus. This gives rise to three overlapping stages in the adjustment: changes in interior thickness drive changes in terminus ice flux, which in turn drive changes in glacier length. On a constant slope, the temporal evolution of the length response is controlled by a single timescale that is a function of only the glacier geometry, the melt factor and the lapse rate. Fundamentally, this is the same timescale as that derived by JRW89. That there should only be one timescale can be viewed as a consequence of the self-similarity of glacier profiles on a bed of constant slope.

This timescale must not be interpreted as an e-folding timescale, however. The characteristic sigmoidal evolution of glacier length in response to a step change in climate, or equivalently, the shape of the power spectrum and ACF, cannot be matched by any one-stage model of the form of Eqn (2). The glacier has more short-term persistence, and is more damped at high frequencies, than Eqn (2) would imply. For example, if τ is estimated from the time it takes for a numerical glacier model to reach (1–1/e) of its new equilibrium length in response to a step change, the inferred timescale will be about twice that of the actual one (Fig. 7). This goes some way to explaining discrepancies in the literature between timescales calculated using different methods (e.g. Oerlemans, 2001; Leclercq and Oerlemans, 2012). Correct identification of the timescale is important for estimating glacier sensitivity to climate changes and trends (via Eqns (6) and (7)).

We have developed a three-stage model that explicitly represents the three phases of the glacier response, as a chain of three first-order equations for thickness, flux and length. Across a wide range of parameter space, this three-stage model closely emulates the temporal evolution of a numerical model for step-function climate changes, climate trends (not shown) and stochastic climate variability that occurs even in a constant climate. The three-stage model has the same number of parameters (three) as the one-stage model. By construction, the equilibrium-length changes and the long-term response to a climate trend are the same as the one-stage model, but the three-stage model is better able to emulate the high-frequency/short-timescale aspects of the numerical model. It might seem surprising that dynamic parameters do not enter into the model, but, as discussed by JRW89, the dynamical parameters are in fact implicit in the glacier’s geometry (e.g. greater sliding leads to a thinner glacier).

The three-stage model has the form of a third-order autoregressive moving-average process (ARMA(3,3)), for which standard formulae exist. Some important metrics of glacier fluctuations (e.g. sensitivity, variance, ACF, power spectrum and excursion probabilities) are analytic functions of the glacier geometry, which can be estimated from field observations. Uncertainties in parameters can be efficiently propagated into uncertainty in these metrics (e.g. Anderson and others, 2014). Moreover, implementing the standard numerical algorithms takes only a single line of MATLAB® code (or other scientific programming language) to filter a known climate time series to give the glacier length, or to invert the glacier-length record to recover the undamped portion of the climate record.

Harrison and others (2003) modified a one-stage model by relaxing the coupling between glacier volume and glacier length variations, thus deriving, in our terminology, a two-stage model; they then calibrated the coefficients to South Cascade Glacier, Washington, USA. Lüthi (2009) derived a similar formulation. The high-frequency limit of a two-stage model is a 1808 phase lag of the response relative to the forcing (e.g. Box and others, 2008), and thus such a model cannot emulate the phase behavior of the numerical model, which reaches 270° lag (Fig. 8d). A two-stage model also cannot match the rollover in the power spectrum (Fig. 8c). Similarly, a four-stage model produces too large a phase lag, and too steep a rollover. Therefore, in addition to the clear physical interpretation (i.e. thickness changes → flux changes → length changes), the flowline model output supports a three-stage model as the correct description. These three stages are fundamental to the adjustment of a glacier, and the literature suggests that using higher-order ice dynamics in a flowline model does not change this picture (e.g. Leysinger Vieli and Gudmundsson, 2004). In this study each stage was assigned the same timescale, but it is possible to imagine that for sufficiently variable bed topography or variable bed conditions, that assumption might need to be relaxed. It may also be possible to emulate the behavior of quasi-cyclical surging glaciers with a different set of coefficients (Budd, 1975; Raymond, 1987; personal communication from A. Robel, 2013).

The three-stage model is, of course, not an exact solution to the nonlinear diffusion equation (Eqns (1)). It was not rigorously derived from the dynamical equations but rather developed heuristically, in order to represent the essential features of the glacier adjustment. There are discrepancies: not every aspect of the flowline evolution is captured (e.g. Fig. 7); and the collapse to a single ACF and spectrum when normalized by τ is only an approximation (Figs 2 and 3), albeit a good one. The three-stage model offers an efficient approach to calculating some useful properties of glacier response. Its discrepancies from the flowline model (and indeed the purpose of including more complex ice dynamics) should be gauged relative to other uncertainties in the problem, among the most consequential of which are: the relationship between glacier mass balance and atmospheric variables (here distilled into a melt factor, µ, a lapse rate and an assumption that accumulation equals annual precipitation); basal hydrology and substrate, and their impact on sliding; bed topography; and, importantly, where the real equilibrium length lies for a given mean climate.

We did not find it necessary to incorporate a height/mass-balance feedback (e.g. Harrison and others, 2001) into the three-stage model, even though the mechanism operates in the flowline model. It is not important for the scales of glaciers we have considered (down to glaciers with a 5% slope, not shown). One measure of when height/mass-balance interactions become important is when the surface slope of the glacier becomes significantly different from the bed slope, averaged over the ablation zone, at which point the geometry of the three-stage model is violated. Other modeling aspects omitted are longitudinal stresses and the detailed pattern of glacier mass balance. Several studies have concluded such details are not important on constant-slope beds (JRW89; Boudreaux and Raymond, 1997; Leysinger Vieli and Gudmundsson, 2004).

We calibrated the three-stage model to a flowline model of Nigardsbreen after adapting it for variable flowline width and bed topography. For realistic white-noise interannual climate variability, even without any climatic persistence, the numerical flowline model produces substantial glacier fluctuations (σ_L = 1 km), which the three-stage model reproduces to within 15%. There is a significant break in the bed slope around the position of the modern terminus, which creates an asymmetry between advance and retreat that the linear, three-stage model cannot capture. Where such circumstances apply, the response of all models depends sensitively on the assumed equilibrium position, and past glacier variations may not be a good guide to the future, in the face of a changing climate.

Simplified models of a glacier’s response to climate have long been a focus in glaciology. The three-stage model presented here represents a substantial improvement over prior models in emulating the ice dynamics that governs the flow of alpine glaciers, and which controls the glacier response on a range of timescales from annual to multi-decadal. It also offers analytic expressions that have a clear physical basis for how some important metrics of glacier response depend on the parameters and geometry of a given setting. These metrics include: a glacier’s sensitivity to climate change and climate trends; a glacier’s variance, autocorrelation, power spectrum and excursion probabilities in response to stochastic climate variability; and the statistical degrees of freedom in a glacier record. It is therefore an efficient tool for interpreting the centennial and millennial climate variability expressed in moraine records; for inverting continuous glacier-length records to recover climate history; and for providing insight into the varying degrees of glacier response in the era of the instrumental climate record and onwards into the future.