Introduction
It is a challenging problem to predict how a glacier responds to climate, because the result is sensitive to the input data, and because the physical processes acting in the terminal region and at the bed are poorly understood. Even with gross approximations for these processes, the resulting numerical or analytical theories tend to be so complicated that it is easy to overlook or fail to exploit the simplifying effect of a fundamental constraint, mass continuity. This is why a simple approach by Reference Jóhannesson, Raymond, Waddington and OerlemansJóhannesson and others (1989a, b) is potentially so useful and why it has attracted so much attention. The basic idea is that mass continuity, together with a minimum of input from theoretical glacier dynamics, should determine a “volume” time-scale (the time to “fill” or “empty” the glacier volume to a new value after a climate change) that characterizes the gross response of a glacier. Applications for a simple approach like this include several related problems: the effect of a glacier’s shape and size on its response characteristics, the regional extrapolation of mass-balance measurements, and the effect of glaciers on sea level.
Here we develop a modified version of this simple theory which, unlike the original, accounts explicitly for the effect of surface elevation on mass balance. We start with mass continuity, identify the mathematical requirement for glacier response to be characterized by a single timescale, and demonstrate a method which uses field data to test whether this is a reasonable approximation in a given situation. We could also make this test within the framework of the conventional theory of glacier dynamics, but because of its inherent limitations we prefer the comparison with data, from which the time-scale, or limits on it, can be determined directly in suitable situations. Data from South Cascade Glacier, Washington, U.S.A., are used for illustration.
Theory
Formulation
The basis of our approach is the concept of a reference-surface balance rate (Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others, 2001). This is simply the balance rate that would occur on some fixed reference surface, which in this paper is taken to be the surface of the glacier as it was at the beginning of a balance time series, defined to be at time t = 0. The usefulness of this balance is that, unlike the conventional one, it is a function of climate alone and is not influenced by changes in the shape and the size of the glacier. For this reason it could be called the “climatic” balance. The concept is applicable both to specific balance (defined at a point) and to the glacier-wide integrated balance. It is important to distinguish between the instantaneous balance rate and the balance accumulated over some finite interval, which we emphasize is usually >1 year in this paper. Both balance rate and balance are considered to be continuous functions of time in our approach. Specific and glacier-wide balances are designated by lower and upper-case symbols, respectively, and reference-surface quantities are primed. The concept of the reference-surface balance rate and the notation are summarized by Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others (2001, fig. 1, table 1).
The starting point is a relationship between the conventional and reference-surface glacier-wide balance rates, and , respectively (Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others, 2001, equation 8):
For our applications ΔV is the difference in the volume between the instantaneous and the initial states, ΔA is the corresponding quantity for map area, is the “effective” specific balance rate in the vicinity of the terminus, and is the “effective” gradient of the specific balance rate with elevation. If the actual gradient is independent of elevation, then takes on that value; otherwise it is a weighted average of the gradient. The weighting factor is the elevation difference between the instantaneous and the initial states. Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others (2001) give methods for evaluating and from field observations. The term in Equation (1) arises when one chooses any point in the map plane, and considers the difference between specific balance rates on the initial and instantaneous surfaces. This difference is the product of the gradient with the elevation difference between the surfaces; integration over the map area to obtain the glacier-wide balance rates gives the ΔV factor. The two balance rates, and , differ because the initial and instantaneous surfaces are different in both area and elevation. The effect of the difference in their areas is accounted for by the term , and the effect of the difference in their elevation by the term .
Equation (1) was used by Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others (2001) to transform between the reference-surface balance rate and the conventional one, but we use it here for a theory of glacier response to climate. We take the reference surface to be the initial surface, the surface as it was at time t = 0, as noted above. It is simplest to use ice-equivalent units for all balance quantities, and to assume that the near-surface density structure of a glacier does not change significantly with time (Sorge, 1935; Reference BaderBader, 1954). Then ΔV is equal to the ice-equivalent conventional balance at time t, accumulated since the initial surface was defined (possibly an interval of many years). This is often called the cumulative balance. Also, , in which V(t) is the instantaneous volume. This is the same as d(ΔV)/dt because ΔV is the instantaneous volume minus the initial volume, the latter of which is a constant. Equation (1) can therefore be rewritten as
Equations (1) and (2) are based on mass continuity and contain no assumptions about glacier dynamics.
Equation (2) is a simple differential equation for the two dependent variables ΔV and ΔA in terms of the independent variable t. The climate-forcing function is , which varies with time but is unaffected by changes in the glacier. Although is really a continuous function of time, in a given year it can be approximated by its time-averaged value. Recalling that , the term in Equation (2) can be interpreted as a negative feedback term; if the climate becomes less favorable for growth, the glacier reduces its rate of loss by reducing its ablation area (ΔA< 0). The term is a positive feedback term; as thinning lowers the surface, the rate of loss tends to become larger because the balance rate becomes more negative. The relative importance of these terms determines whether the glacier has a large or even unstable response to climate.
A simple theory of glacier response could be based on Equation (2) if one could assume a functional relationship between area A and volume V. This is a simple condition from a mathematical point of view, but except for a perfectly plastic glacier, it is empirical and a potentially drastic approximation to the dynamics of a real glacier. It can be tested with data from a balance program that also measures area on a routine basis. We expect a relationship between A and V to be simplest and best defined when the glacier is changing very slowly. Then time lags are negligible and area tends to remain “in phase” with volume. The simplest plausible relationship is that the rate of change of volume with respect to area is constant, or
where H is the “thickness scale” (see Reference Jóhannesson, Raymond, Waddington and OerlemansJóhannesson and others, 1989a, p. 350). A useful integrated version of Equation (3) is
This implies that ΔV = 0 when ΔA = 0, which is consistent with our definition of the reference surface. Equation (4) can be used to eliminate ΔA (or ΔV) in Equation (2), which becomes
where τ v is a time-scale given by
Equation (5) describes the change in the volume of a glacier in response to the climate forcing , the reference-surface balance rate. Equation (3) is the essential condition for the validity of Equation (5); the single parameter H is the constraint on glacier dynamics that leads to a simple response characterized by the single parameter τ v defined by Equation (6). Of course and also occur in Equation (6), but they are climatic parameters. Note that τ v is constant only to the extent that H, and are constant. We use long-term averages of and in the application discussed below.
Our time-scale τ v has the same meaning as the “volume time-scale” or the “memory length” of Jóhannesson and others (Reference Jóhannesson, Raymond, Waddington and Oerlemans1989a, Reference Jóhannesson, Raymond and Waddingtonb), but there is a significant difference because ours accounts explicitly for the effect on balance rate of the changing surface elevation of the glacier, via the term in the denominator of Equation (6). Their τ v was defined to be (our notation), in which characterizes the balance rate at the elevation of the ice surface in the terminus region, while our characterizes it at the elevation of the bed there (Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others, 2001). One may therefore expect to be less than . This accounts partially for the effect of changing surface elevation on balance rate in their approach, but their time-scale estimates still will be shorter than ours. The term is important, as discussed below.
Steady state
Equation (5) can be applied to the important question of the sensitivity of a glacier to climate change. First we compare our approach to that of traditional theory (Reference NyeNye, 1960). Suppose that the climate, its effect characterized by the reference-surface balance rate since t = 0, has been constant for a time t ≈ τ v. Then if the response is stable, a steady state is approached which can be found by setting d(ΔV)/dt = 0 in Equation (5), giving
and
by Equation (4). ΔV ∞ and ΔA ∞ are the final changes in volume and area measured relative to their values at t = 0. To facilitate comparison with the traditional approach, it is useful to characterize the climate by a slightly different version of the reference-surface balance rate, its glacier-wide average :
where A′ is the map area of the glacier in its initial state. Then if Equations (6), (8) and (9) are combined, one finds that
Equation (10) is different from the corresponding result from traditional theory (Reference NyeNye, 1960, equation 40), which in our notation is
where ΔL ∞ is the change in length of the glacier from an initial steady-state condition, and L is the total length. One obvious difference between Equations (10) and (11) is that the variable on the lefthand side is area in the first case and length in the second. This is because Equation (11) is restricted to glaciers of uniform width. Another difference is that the balance rate at the terminus occurs in Equation (11), rather than the effective value , which is roughly 25% smaller in magnitude in the example discussed below (Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others, 2001). However, the most important difference is the absence of the term in the denominator of Equation (11), which means that it does not account for the effect of the surface elevation of the glacier on balance rate, as was recognized by Reference NyeNye (1960). The dependence of Equation (10) upon H means that it depends upon glacier dynamics, unlike Equation (11). This is a complication but it is realistic because the final adjustment of the glacier must depend upon its mechanical properties. Finally, Equation (11) requires the initial condition of the glacier to be a steady-state one, while Equation (10) does not. This advantage is more apparent than significant, because we expect Equation (3), and therefore Equation (10), to be most accurate when the initial area has had time to adjust to the initial volume; in other words, when the initial condition is near steady state.
Solution for different reference-surface balance-rate scenarios
It is instructive to consider solution of Equation (5) for two different scenarios of the reference-surface balance rate . If the climate is steady for t > 0, and is therefore constant, the solution is
Because of Equation (4), the solution for ΔA is the same except for the factor H. In this simple case the volume and area follow an exponential approach to a new steady-state geometry, where the time constant is τ v and the final state is that given by Equations (7) and (8). If the climate were to undergo a steady rate of change from time t = 0, with
where C is a constant, the solution would be
Because below we wish to consider the responses of glaciers of different sizes, it is useful to rewrite these solutions in terms of the glacier-wide averages of and ΔV. We use the reference-surface glacier-wide average balance rate , defined by Equation (9), and the glacier-wide average balance , defined by
is the spatially averaged thickness change of the glacier since the reference surface was defined at t = 0, as might be determined by geodetic or traditional glaciological methods. It is useful only if changes in area are small enough that A ≈ A; otherwise one must keep the variable ΔV. With this approximation, division of Equations (12) and (14) by A ′ gives
and
where c = C/A′. If we denote as the value of for t ≈ τ v, then Equations (16) and (17) both give
In the former case is constant, and in the latter it is proportional to time. Equation (18) is equivalent to Equation (7) for small area changes.
The general solution to Equation (5) for any form of the reference-surface balance rate is
which indicates a fading memory to previous climate. The usefulness of this relation is affected by several factors, most notably the limited validity of Equation (3).
Meaning of the time-scale
The time-scale τ v defined by Equation (6) is positive if , which is the case if the negative feedback term in Equation (2) is larger than the positive one. This is the usual case; the response is stable and a steady state is eventually attained if climate remains constant. The critical value of τ v occurs when . In this special case the positive and negative feedback terms cancel, |τ v| is very large, and the reference-surface and conventional balance rates remain equal even though the initial and instantaneous surfaces are different. τ v is negative if , which is the case if the positive feedback term is the larger of the two. In this case the response is unstable. This idea of instability or near-instability was studied by Reference BöðvarssonBöðvarsson (1955), Reference NyeNye (1960) and Reference WeertmanWeertman (1961) with characterizations of glacier dynamics different from our Equation (3).
The time-scale has a dual role. Equation (18) is an example of how it determines the amplitude of the final response for t ≫ τ v (Reference Jóhannesson, Raymond, Waddington and OerlemansJóhannesson and others, 1989a). In this sense it can be said to determine the sensitivity of a glacier to climate. It also determines the rate of response, as its name implies, but not at small time. Because ΔV ≈ 0 for small time, ΔV/τ v ≈ 0 also, so by Equation (5) , which is independent of τ v. This is easily confirmed for the two examples above by expanding the exponentials in Equations (16) and (17). This result is not surprising, because the initial and instantaneous surfaces are the same at small time. Thus there are characteristically different responses at large and small time, as is reasonably well known (e.g. Reference Warrick, le Provost, Meier, Oerlemans, Woodworth, Houghton, Filho, Callander, Harris, Kattenberg and MaskellWarrick and others, 1996).
One can make a more general statement about the response at small and large time. If two glaciers are subjected to the same reference-surface average balance rate , their average thickness change will be the same for small time. In contrast, at large time the response is proportional to the time-scale. The result is that given two ice masses with the same , the one with the larger time-scale will have the larger response (larger ) for all t > 0, at least within our approximate version of glacier dynamics.
This result brings up a point that is sometimes confusing in the literature. Reference Oerlemans and FortuinOerlemans and Fortuin (1992), for example, stated that for times on the order of a century or a decade the large ice sheets should be relatively unimportant contributors to sea-level change because of their relatively long time-scales compared with those of glaciers and small ice caps. There are two issues concerning this assertion. The first, not discussed here, is whether the implied relation between size and time-scale is correct (see, e.g., Reference Bahr, Pfeffer, Sassolas and MeierBahr and others, 1998). The second is whether the implied relation between time-scale and magnitude of the response is correct. In the light of our discussion, the average thickness changes of Oerlemans and Fortuin’s two classes of ice masses would be the same at small time for a given reference-surface average balance rate. Then because Greenland, for example, has such a large area it would have a much greater effect on sea level than the mountain glaciers and ice caps, ignoring the possible effects of calving. As time proceeds, Greenland would become even more important if it had a larger timescale because glaciers with the largest time-scales have the largest response for t > 0, as noted above. Of course these comparisons require the different ice masses to be subjected to the same reference-surface average balance rate, which is not likely. Nevertheless, it is not correct to ascribe differences in short-term response to differences in the volume time-scale.
Application
The ideas described above can be illustrated with data from South Cascade Glacier, a small glacier in the North Cascade Mountains, northwest U.S.A. It has a record of annual balances measured by traditional glaciological methods beginning in 1958 (Reference KrimmelKrimmel, 1999b), an extensive set of volume measurements made by airborne photogrammetry (Reference KrimmelKrimmel, 1999a), and a continuous record of area. We use the geodetically corrected cumulative balance series constructed from these data by Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others (2001), which is about 50% more negative than the uncorrected one. This series, in ice-equivalent units, is our AV, measured relative to September 1970. We do not use the earlier data since we are uncertain about their accuracy.
ΔV is shown as a function of ΔA in Figure 1, both measured relative to autumn 1970. There is a reasonably linear relationship between the two, which suggests that Equation (4) and therefore Equation (3) are satisfied approximately. The least-squares value for the slope dV/dA is 171 m, which approximates the thickness scale H. The corresponding value for the time-scale τ v can be estimated from Equation (6) using values of the effective balance rate near the terminus and the effective balance-rate gradient determined by Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others (2001, fig. 5). These are −6.2 m a−1 (1970–97 average) and 0. 024 a−1 (1985–97 average), respectively. The resulting value for τ v is 82 years, although this value is probably too large as discussed below. Reference Schwitter and RaymondSchwitter and Raymond (1993) also discussed the time-scale for this glacier.
The role of the time constant in determining the sensitivity to climate is easily illustrated with this value of the time-scale. The climate between 1970 and 1997 can be roughly characterized by a reference-surface balance rate of −1 m a−1 (Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others, 2001). If this climate were to persist, and if the 82 year value for the time-scale were correct, one would expect the glacier to stabilize after losing an average thickness of 82 m relative to 1970 (Equation (18)). This example illustrates the simplicity of the approach, but it is doubtful that the simple theory would be valid for such large changes. It is worth pointing out that as a steady state is approached in this scenario, the conventional average balance rate approaches zero, but the reference-surface average balance rate remains at −1 m a−1.
Discussion
Meaning of the thickness scale
As noted above, the formal definition of the thickness scale H is the derivative of volume with respect to area when the glacier changes sufficiently slowly. The name “thickness scale” could be misinterpreted. Our value for H is comparable to depths measured in the thickest part of South Cascade Glacier by Reference HodgeHodge (1979) from 1973 to 1977, but it is easy to imagine glacier shapes in which the derivative is quite different from the maximum thickness. An example would be a glacier whose upper part is flat and thick, and whose lower part is steep and thin. Our derivative would be representative of the thin lower part of such a glacier.
The variation of H with the size of a glacier would affect the range of climate change for which our results are valid. Our approach is insensitive to this variation for South Cascade Glacier over the duration of the study period, during which the area change was about 15%. Theoretical analyses of the effect have been carried out by Jóhannesson and others (Reference Jóhannesson, Raymond, Waddington and Oerlemans1989a, Reference Jóhannesson, Raymond and Waddingtonb) and by E. D. Waddington (personal communication, 2000).
Uncertainties
A discussion of some of the uncertainties illustrates the difficulty in making accurate estimates of H and τ v. The most obvious uncertainty is in the effective balance rate near the terminus , which is not well known because its determination requires extrapolation of balance rates measured at a few points higher on the glacier to the terminus or beyond, while also taking into account the three-dimensional nature of the terminus. An additional complication is that the balance rate may vary rapidly near the terminus. also has a weak trend with time and area (Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others, 2001). By Equation (6) this means that our simple approach, in which the time-scale τ v is assumed to be constant, is probably valid only for a small range of area change.
A less obvious and probably more serious source of uncertainty is in the value of 171 m for the thickness scale H, which came from the slope of the line in Figure 1. The problem is rooted in our fundamental assumption about the dynamics of the glacier, Equation (3), which implies that area is able to react instantaneously to volume, and therefore to remain “in phase” with it. If part of the slope of the line in Figure 1 arises from the failure of area to do so, then 171 m would be an upper limit on H. Considerations beyond the scope of this paper suggest that a more appropriate value, one associated with truly slow changes in the glacier, would be about 30% smaller. This would decrease our estimate of the time-scale, and therefore of the sensitivity to climate, by about a factor of two.
Importance of the surface elevation term
The dependence of the balance rate on surface elevation has a major effect on the response. Taking a best value for H about 30% lower than our 171 m as just discussed, we find that the ratio of the two terms in the denominator of Equation (6), , is about 0.5. In other words, the effect on the balance rate of surface elevation (the origin of the second term), although smaller than that of area (the origin of the first), is significant. The result is that at South Cascade Glacier the surface elevation term causes a doubling of the time-scale and of the sensitivity to climate.
Qualitative statements can be made about the effect of surface elevation on the responses of other ice masses. We write the ratio of the two terms in the denominator of Equation (6) in a different way, using the fact that the quantity is a rough measure of the elevation of the equilibrium line relative to the terminus on South Cascade, and probably on most other glaciers. Then the ratio, , is roughly H/Z t→eq. For low-slope glaciers, which are usually thick, one would expect this ratio to be relatively large, but as we just saw, it is on the order of 0.5 even for South Cascade Glacier, which with a characteristic surface slope of roughly 10° is relatively steep. For lower-slope ice masses H/Z t→eq may approach unity, which by Equation (6) would make τ v, and therefore the sensitivity, very large. Thus the effect of surface elevation seems to be critical in the response of low-slope ice masses, and significant in that of most others.
If the magnitude of , or equivalently of H/Z t→eq, is close to unity, there is a fundamental problem in calculating the response of a glacier, because of the near-cancellation of the surface elevation and area effects. This is brought out clearly by our model, because the evaluation of the timescale, and therefore the sensitivity to climate, requires the subtraction in the denominator of Equation (6) of two numbers with the same order of magnitude, neither of which is ever known accurately. Although this problem may be less obvious in other models, it must exist in all, whether analytical or numerical, because it is due to the fundamental interplay between positive and negative feedback. For example, numerical models which specify balance rate as a function of elevation automatically account for the effect of elevation on balance rate, but the possibility of near-cancellation of the surface elevation and area effects remains, and it amplifies the effects of errors in the model or the input data. This means that the error in the calculated response could be very large, and that it needs to be investigated in individual cases.
Our discussion points out the limitations of traditional theory, the best-known versions of which do not consider the effect of surface elevation on balance rate at all. This is worth keeping in mind because traditional theory still tends to be the basis of our perception of what is important in the response of ice masses to climate. Nevertheless, the limitations of traditional theory may sometimes be less significant than the uncertainty in the balance rate near the terminus, which has long been known to be a key to the prediction of response (Reference NyeNye, 1960), but is rarely measured.
Acknowledgements
We are grateful for the important suggestions of T. Jóhannesson (Scientific Editor), C. F. Raymond, E. D. Waddington and M. Kuhn. This work was supported by U.S. National Science Foundation grant OPP9707515.