1. Introduction
The problem of predicting the response of a glacier to a given climate scenario has received a great deal of attention over the years. The conventional approach is to use a numerical model capable of estimating the full longitudinal profiles of surface elevation and velocity, the terminus position and their evolution with time. An analytical approximation, which we call the “traditional” approach, is also useful because it shows the general nature of the results to be expected (e.g. Reference NyeNye, 1960, Reference Nye1963). These models contain some of the physics of flow dynamics, but are semi-empirical because they also contain major assumptions, approximations, or adjustable parameters to account for the incompletely known geometry of the glacier and processes of ice deformation, basal motion, longitudinal-force transmission and terminus dynamics. Even though these models are semi-empirical, they are complex enough that simple controls on the amplitude and time-dependence of the response, what their orders of magnitude are, and how accurately they can be determined, tend to be buried in detail.These issues can be addressed with another, much simpler, semi-empirical approach which started with papers by Reference Jóhannesson, Raymond, Waddington and OerlemansJóhannesson and others (1989a, Reference Jóhannesson, Raymond and Waddingtonb).
The idea behind the alternative approach is that a simple approximate-response theory can be obtained by ignoring the spatial distributions of map area, thickness and mass-balance rate. One considers only the total map area, volume and glacier-wide balance rate. It would be more conventional to use length instead of area, but because the termini of some glaciers are complicated, area is more fundamental, especially in the case of small glaciers (Reference Ye, Yunjian, Fengjing and CaohaiYe and others, 2003). Area and length are closely related because most area changes occur near the terminus. The resulting theory is contained in two equations for the two unknowns, area and volume. One equation expresses the dynamics of flow in terms of a relationship between area and volume and the other expresses mass continuity in terms of the same variables.
This type of “macroscopic” approach has a long history in physics. It is sometimes advantageous to ignore the fine details of physical processes in favor of simple equations of state or constitutive laws which may have some physical motivation, but whose calibration and validation depend primarily upon observation. The motivation in this case comes first from the fact that our existing version of macroscopic theory (Reference Harrison, Elsberg, Echelmeyer and KrimmelHarrison and others, 2001) would be exact for a perfectly plastic glacier, which actual ice masses approach to some extent. Second, it has been shown with a conventional numerical model that the shape of a glacier is mainly controlled by the glacier-wide balance rate, while the details of balance-rate distribution have a relatively minor effect (Reference Boudreaux and RaymondBoudreau and Raymond, 1997), as is assumed in macroscopic theory. Here we generalize our existing version of macroscopic theory to account for departure from perfect plasticity. It is noteworthy that our approach is a first-order one limited to small changes.
2. A Relationship Between Area and Volume
2.1. Formulation
Our first task is to formulate an appropriate relationship between glacier area A and volume V. We shall motivate a simple relationship, and parameterize it in a way that allows its range of usefulness to be tested with observations. We adopt the point of view that in this relationship area is the dependent variable; it “responds” to volume. Volume is taken as the independent variable because mass continuity (considered in the next section) puts a constraint on the rate of volume change.
With area the dependent variable, departure from perfect plasticity is interpreted as introducing a delay between volume change and area response to it. The limiting case, in which volume changes so slowly that area effectively stays in adjustment with it, is of special interest. The states corresponding to this situation lie on a curve in the A−V plane which we sketch in Figure 1 and label (A a, V a). We call these the “area-adjusted states”, realizing that although area and volume are in adjustment with each other, they may not be in adjustment with climate. When they are, the term “steady state” will be used.
Relationships in the plane defined by area and volume. The area-adjusted states lie along the curve Aa(V), whose tangent at V′ is shown by the broken line. (A′, V′) are the initial area and volume. ΔA0 is the difference between the initial area and that required for adjustment with the initial volume.
For the real situation in which the volume changes at a finite rate, we take the rate of adjustment of area dA/dt to be proportional to the amount by which area is out of adjustment with volume.We express this as the difference between the instantaneous area A and the area in a particular areaadjusted state, the one A a(V) which would correspond to the instantaneous volume V . Thus we write
in which 1/τ A is a proportionality factor and the negative sign makes the response stable for positive τ A. We call τ A the “area time-scale”. It may depend upon climate, but if the glacier is not far from steady state (a requirement of first-order theory), one can show that climate dependence enters via second-order terms. τ A will therefore be taken to be constant.
Equation (1) is our tentative relationship between area and volume, but it is not in the most useful form. The most obvious problem is that our approach is a first-order one, and therefore its scope is limited to the calculation of small changes in area and volume from some reference state, rather than area and volumes themselves. (This is mitigated by the fact that the change in volume is often a measured quantity or of more interest than the volume itself, especially in hydrologic or sea-level applications.) The other problem with Equation (1) is that the curve representing the area-adjusted states A a(V) is usually unknown. Both problems can be addressed by changing variables from area and volume to their changes, ΔA and ΔV , as measured with respect to the initial state at time t = 0. Primes denote quantities associated with the initial state. (A′, V′) is usually not an area-adjusted state. Formally,
ΔV is measured by geodetic methods or by accumulating annual balance measurements in ice-equivalent units, at least if the near-surface density structure of the glacier is approximately constant in time (see, e.g., Reference BaderBader, 1954). ΔV would then be the cumulative balance in ice-equivalent units, which we use exclusively. ΔA can be measured by a variety of mapping methods. Elimination of A in Equation (1) and minor rearrangement give
in which we have used the fact that dA/dt = dΔA/dt, valid because A′ is a constant.
Equation (4) still contains the unknown function A a(V) with its dependence upon the usually unknown volume V. A power-series expansion of A a(V) about the point V′ expresses it in terms of ΔV and parameters that can be determined from observations:
where
evaluated at V′. H is defined to be the “thickness scale”. We neglect the higher-order terms in accordance with our first-order approach. This means that we approximate the curve of area-adjusted states by its tangent at point V′ (Fig. 1), which is valid for small changes from the initial state. Combination of Equations (4) and (5) gives
where
ΔA 0 specifies the condition of the initial state, or more precisely, the amount by which area exceeds that needed for adjustment with volume initially (Fig. 1). It is determined by the history of the glacier prior to t = 0. In the case of no delay in area response (τ A = 0), Equation (7) gives
which is the tangent line approximating the curve of area-adjusted states in the (ΔA, ΔV ) plane.
Equation (7), or equivalently its general solution for ΔA = 0 at t = 0,
is our tentative relationship between area and volume, in which ζ is an integration variable. Equation (10) contains three parameters: the area time-scale τA and the thickness scale H, which characterize the dynamics of the glacier, and the area ΔA 0, which characterizes the condition of the initial state. It is a generalization of the simplest example of a relationship between area and volume (Reference Jóhannesson, Raymond, Waddington and OerlemansJóhannesson and others, 1989a, Reference Jóhannesson, Raymond and Waddingtonb; Reference Harrison, Elsberg, Echelmeyer and KrimmelHarrison and others, 2001) in which the glacier is plastic. Then τ A = 0 and ΔA 0 = 0 because area responds instantly to volume, all states are area-adjusted and there is only the single parameter H.
Assuming for the moment that Equation (7) is valid, the determination of its three parameters τ A, H and ΔA 0 requires at least three sets of measurements of ΔA, dA/dt and ΔV . The time required to collect the data will depend upon the measurement accuracy and how rapidly the glacier is changing, but a time resolution significantly < τ A is needed. Once the parameters are known, future measurements of ΔA and dA/dt, which can be made relatively simply, can be used in Equation (7) to determine the cumulative balance ΔV . This idea, that balance may be estimated from observations near the terminus, is an old one (Reference NyeNye,1965). Our approach is simpler and focuses on the cumulative balance, which is more directly determined than the annual one.
Finally, it is worth pointing out that there is an alternative to our Equation (7) or (10) for describing the relationship between area and volume. It uses scaling relationships (Reference Raper, Briffa and WigleyRaper and others,1996; Reference Bahr, Meier and PeckhamBahr and others, 1997), but its accuracy is insufficient for our purposes.
2.2. Application to South Cascade Glacier
We now test our area–volume relationship with data from South Cascade Glacier, a 2 km2 valley glacier in the North Cascade Mountains, Washington, U.S.A. It has a long record of annual balance measurements made by conventional glaciological methods, a continuous and reasonably accurate record of area, and an extensive set of measurements of volume change made by geodetic methods, especially since 1970 (Reference KrimmelKrimmel, 1999). By then the glacier had retreated out of a lake into which it had calved earlier, a process which could have complicated the relationship between area and volume. For these reasons we use a cumulative balance series that starts in 1970.The accuracy of the uncorrected version depends upon how the estimated error in the annual balances (0.2 m a−1 for the area average) is partitioned between random and systematic components, and thus how it accumulates when the annual balances are summed. This is unknown, but the unusually abundant geodetic data at South Cascade Glacier permit corrections to be made to the cumulative balance series (Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others, 2001). ΔA is shown as a function of the resulting ΔV in Figure 2.
Measured ice-equivalent cumulative balance ΔV as a function of change in area ΔA for South Cascade Glacier (filled circles). The best fit with Equation (10) is shown by the solid curve.The tangent to the curve of area-adjusted states (Equation (9)) is shown by the broken line. The scale on the right side is for the fractional area change ΔA/A′, where A′ is the area in autumn 1970.The scale on the top is for the average thickness change.
The best values for the parameters τ A, H and ΔA 0 can be determined easily by fitting Equation (7) or (10) to the data in Figure 2 by the method of least squares. We used Equation (10) with the integral approximated by the trapezoidal rule, and commercial software employing the Levenberg–Marquardt algorithm. The data were weighted equally. The fit is shown in Figure 2, together with the tangent to the curve of area-adjusted states (Equation (9)). A complete evaluation of the success of the fit requires knowledge of the errors in ΔA and the cumulative balance ΔV ; the latter are probably the more important. Considering the small ambiguities which we encountered in making the geodetic corrections to the ΔV series, we can only judge that the deviations between the fit and the data in Figure 2 are comparable to the uncertainties. Thus our three parameter area–volume relationship seems to fit the data reasonably well, including its complex fine structure. This provides support for its usefulness.
The parameter values and their errors as estimated from the quality of the fit are summarized in Table 1. These errors are significant despite the reasonable quality of the fit and do not necessarily include all the effects of systematic errors in the data. This illustrates the importance of measurement accuracy. The area time-scale τ A is roughly 8 years. The thickness scale H (∼123 m) is less than the value obtained under the assumption τ A = 0, as discussed by Reference Harrison, Elsberg, Echelmeyer and KrimmelHarrison and others (2001). It is seen from the value of ΔA 0 in Table 1 that the initial area, in autumn 1970, was roughly 4% too large to be in adjustment with the initial volume.The fit is significantly poorer if the initial state is assumed to be area-adjusted.
Parameters and 1-σ errors in the area–volume relationship for South Cascade Glacier
3. The Response to Climate
3.1. Mass continuity
The macroscopic theory is completed by a statement of mass continuity:
is the “effective” specific balance rate at the terminus (a negative number), 
is the “effective” gradient of the balance rate with elevation, and 
is the glacier-wide “reference-surface” balance rate, the balance rate as it would be on the surface of the glacier if it remained in its initial t = 0 “reference” state. Details are given by Reference Harrison, Elsberg, Echelmeyer and KrimmelHarrison and others (2001; equation (2)) and Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others (2001; equations (9) and (10)). We have replaced their 
with 
, thereby consistently using lower-case symbols for specific or “local” balance quantities (m a−1) and their gradients (a−1), and upper-case symbols for glacier-wide balance quantities (m3 a−1). All balance quantities are expressed in ice-equivalent units. 
is a continuous function of time but is usually adequately approximated by the reference-surface annual balance expressed as a rate 
and 
can be determined from suitable field observations as done by Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others (2001) for South Cascade Glacier. The terms 
and 
account for the effects on the glacier-wide balance rate of changing map area and surface elevation, respectively. Equation (11), like Equation (7), is valid for small changes from the initial state.
Equation (11) is useful because the two variables, ΔV and ΔA, are the same as those in our area–volume relationship (Equation (7)), and because the reference-surface balance rate 
, unlike the conventional one, is unaffected by changes in the surface of the glacier as it responds. It is thus the appropriate purely climatic forcing term for macroscopic theory; the response to particular climatic variables can be found by expressing 
in terms of them. Notice that although 
is not affected by changes in the surface, it does depend upon the initial configuration of the surface.
3.2. Separate equations for area and volume
Equations (7) and (11) are two simultaneous first-order differential equations for the two variables (ΔA, ΔV ), with initial conditions (0, 0). They constitute the complete macroscopic theory for the response of area and volume to climate as specified by the glacier-wide reference-surface balance rate 
. Since numerical solution is easy, the discussion could be ended here but our purpose is to exploit the simplicity of the macroscopic theory with an analytical approach which, like the traditional one, indicates some general properties of the response.
Although it is not essential, we have found it instructive to decouple ΔA and ΔV in Equations (7) and (11) by producing two new differential equations in which they occur separately. Solving Equation (7) for ΔV and substituting in Equation (11) gives
while solving Equation (11) for ΔA and substituting in Equation (7) gives
where τ V is defined by
Equation (12) requires H and τ
A to be constant, while Equation (13) requires the same of 
and 
; we ignore their interannual variations. τ
V is the “volume time-scale” of Reference Jóhannesson, Raymond, Waddington and OerlemansJóhannesson and others (1989a) as modified by Reference Harrison, Elsberg, Echelmeyer and KrimmelHarrison and others (2001) to account explicitly for the effect of changing surface elevation on balance rate. A small mathematical detail is that since we now have second-order differential equations, we need two more initial conditions. These are obtained with no new assumptions by setting (ΔA, ΔV ) = (0, 0) at t = 0 in Equations (7) and (11), giving 
at t = 0.
3.3. Solutions
Equations (12) and (13) can be solved for ΔA and ΔV once the time dependence of 
is specified. This is done in the Appendix for three cases: 
a constant; 
an impulse; and 
any function of time. Several properties of the solutions merit consideration.
3.3.1. Decomposition of the response
The solutions for all 
cases have the common property that both ΔA and ΔV can be separated into two terms. The first is the response as if the initial state were an area-adjusted one (ΔA
0 = 0), and the second as if ΔA
0 ≠ 0 but the climate were neutral 
for t > 0. We shall sometimes refer to these as the “direct” and the “transient” terms. The direct term is determined by the time dependence of the climate after t = 0, and the transient term by the climate history prior to t = 0.Thus the transient term is independent of 
, which we specify only for t > 0. Significantly, the transient term approaches a constant different from zero at large time. We shall use the notation ΔA
T and ΔV
T for the transient part of the response. If 
is constant for t > 0, the direct term is the response of a glacier initially in a steady state to a step-change in climate.
3.3.2. Damped spring and mass analogy
The lefthand sides of Equations (12) and (13) contain the terms which can be thought of as characterizing the glacier: the righthand sides, the climate forcing and the condition of the initial state. With this in mind, there is a useful, well-studied analogy. Each of these equations is similar to that describing the response of a damped spring and mass system to an external force (Fig. 3). ΔA in Equation (12) and ΔV in Equation (13) are analogous to the displacement of the mass from the initial position. Examples of analogous parameters on the lefthand sides of these equations are τ A → mass and 1 = 1/τ V → spring constant. A list of analogs is given inTable 2. As an example of the use of these analogs, Reference Harrison, Elsberg, Echelmeyer and KrimmelHarrison and others (2001) pointed out that at least in the special case τ A = 0, ΔA and ΔV are unstable when τ V is negative.This is also true in the more general case considered here because 1 = 1/τ V is analogous to the spring constant, which if negative would mean that the spring would push when it would normally pull, making the mechanical system unstable.
Corresponding quantities in damped spring and mass analogy. ΔA in Equation (12) and ΔV in Equation (13) are analogous to the displacement of the mass from its initial position
Damped spring and mass system with governing equation. g is the gravitational acceleration.
3.3.3. Damping
Both the damped spring and mass system and the glacier have quite different responses depending on the strength of the “damping”. The damping of the glacier is characterized by the non-dimensional parameter p:
The origin of p is motivated by a simplification of Equations (12) and (13) given in the Appendix. For p < 1 the glacier is “under-damped”, which means that if 
becomes constant, ΔA and ΔV overshoot their final values and undergo decaying oscillations, or “ringing”, about them. If p > 1, the glacier is “over-damped”, and the approach of ΔA and ΔV to their final values is monotonic. This behavior may be more complicated when the transient terms are significant. For p = 1 the glacier is “critically damped”. In the example considered below, p ≈ 1 so the damping is close to critical. It is shown in the Appendix that the time dependence of the response is most naturally expressed in terms of the time 
, the geometric mean of τ
A and τ
V.
3.3.4. Steady climate
If the climate is steady (
constant for t > 0), area and volume will eventually approach the steady-state values (ΔA
∞, ΔV
∞), at least if τ
V > 0.These can be found by setting the derivatives in Equations (12) and (13) equal to zero while holding 
constant:
These equations are important because they are a statement about the ultimate evolution of a glacier under steady-climate conditions. The first equation, when divided by the initial area A′, is a generalization of the approximation of Reference Harrison, Elsberg, Echelmeyer and KrimmelHarrison and others (2001; equation (10)), which included the effect of surface elevation on balance rate through the definition of τ V, but not the effect of the condition of the initial state as characterized by ΔA 0. That approximation, in turn, was a generalization of the original and simplest one given by Reference NyeNye (1960; equation (40)).
It is of interest to consider what determines the relative importance of the direct 
and transient (ΔA
0) terms. For the area response (Equation (16)), it depends upon the relative magnitudes of 
and 
. For the volume response (Equation (17)) it depends upon the relative magnitudes of 
and 
. In the damped spring and mass analogy, the issue would be the relative magnitudes of the external force and the weight (Table 2). The transient term must be more important in the volume than in the area response, because 
for τ
V > 0 by Equation (14), which is the condition for stable response.
In Equations (16) and (17) ΔA ∞ and ΔV ∞ are proportional to the volume time-scale τ V. In other words, τ V scales the ultimate response, or “sensitivity”, of glacier area and volume to the combined effects of long-term climate and the condition of the initial state. This is an example of how τV specifies an amplitude as well as a time-scale.
4. An Example of Climate Response: South Cascade Glacier
4.1. Interpretation with the macroscopic approach
We now consider the response characteristics of South Cascade Glacier. The input values are as follows. For τ
A, H and ΔA
0 we use the values already determined (Table 1): 8 years, 123 m and 0.094 × 106 m2 respectively. We use the measured spatial pattern of mass-balance rate to estimate 
and 
. The former value is slightly different from that used by Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others (2001) and Reference Harrison, Elsberg, Echelmeyer and KrimmelHarrison and others (2001) and was found by a slightly more accurate method.
These five input parameters produce the following results (Table 3). A value of (48 ± 17) years for the volume time-scale τ
V follows from Equation (14), which is 6 × τ
A. The geometric mean time 
, which most naturally characterizes the time-dependence of the response, is 20 years. By Equation (15) the damping parameter p = 1.0 ± 0.3, which indicates that our best value for it is the critical one. The errors in 
and H are the main contributors to the error in τ
V, while the error in τ
A is the main contributor to the error in p. The errors are not well-known, but they are substantial.
Parameters and 1-σ errors characterizing the response of South Cascade Glacier to climate
We next consider how the area and volume of this glacier should respond to some climate scenarios. The simplest example is steady climate, 
constant for t > 0, where t = 0 is autumn 1970, when the reference surface is defined and when A′ is the area (Table 1). We take a value of 
such that 
, which is its mean value for 1970–97 (Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others, 2001). This simple scenario will approximately reproduce the behavior of the glacier starting in 1970.
The response of ΔA follows from Equations (A5), (A6) and (A7) and is shown in Figure 4. As discussed above, it is the sum of the direct term (due to non-zero 
for t > 0, labeled ΔA
S), and the transient term (due to the initial misadjustment of area and volume, labeled ΔA
T). ΔA
S starts off with zero slope and finally settles down to a value of ΔA
S/A′ of about −39%. ΔA
T initially decreases and continues to do so for several years. In the damped spring and mass analogy, this would be because of the initial non-zero momentum, which is analogous to −ΔA
0 (Table 2). ΔA
T/A′ finally changes sign and settles down to a value of about 5%. In this scenario the glacier stabilizes after losing about one-third of its area relative to 1970, a process that is about two-thirds complete after the time τ
V, 48 years. ΔA for the special case τ
A = 0 (Equation (A11)) is shown by the broken curve in Figure 4. In this case ΔA
T = 0 because ΔA
0 = 0 for τ
A = 0, as discussed above.
Response of area ΔA to a constant glacier-wide average reference-surface balance rate 
of −1 m a−1 (ice). ΔA is the sum of the direct and transient components ΔAS and ΔAT.The response for the special case τA = 0 is shown by the broken curve.
The response of ΔV follows from Equations (A8), (A9) and (A10) and is shown in Figure 5. ΔV for the special case τ
A = 0 (Equation (A12)) is shown by the broken curve.The most notable difference from the ΔA case is the magnitude of the transient term, which decreases the magnitude of the final response by 22%, a value which would have been even larger had 
not been so large. This effect arises from an initial state in which area is out of adjustment with volume by only 4%, illustrating the importance of the condition of the initial state, particularly in the response of volume. The volume stabilizes after a loss of 38 m in average thickness, about a third of the thickness scale H.
Response of volume ΔV to a constant glacier-wide average reference-surface balance rate 
of −1 m a−1 (ice). ΔV is the sum of the direct and transient components ΔVS and ΔVT.The response for the special case τA = 0 is shown by the broken curve.
Finally, we consider the responses of area and volume to an impulse in climate, which means that 
is proportional to a δ function in time.The direct responses, ΔA
I for area and ΔV
I for volume, are calculated from Equations (A14) and (A17) and shown in Figures 6 and 7. The complete responses will also include the same transient terms, ΔA
T and ΔV
T, as the steady-climate case, but only ΔA
I and ΔV
I are shown.
Response of area ΔAI of an initially area-adjusted state to an impulse in glacier-wide average reference-surface balance rate 
sufficient to add an average thickness of 1 m (ice) to the glacier.The normalization of the “traditional” curve is arbitrary.
Response of volume ΔVI of an initially area-adjusted state to an impulse in glacier-wide average reference-surface balance rate 
sufficient to add an average thickness of 1 m (ice) to the glacier.
4.2. Comparison with the traditional approach
The response characteristics of South Cascade Glacier were also analyzed by Reference NyeNye (1963, Reference Nye1965), who used the “traditional” approach. In the latter paper he represented them with a function that describes the evolution of the terminus in response to a uniform unit thickness instantaneously added to an initially steady state, one which approximated the actual surface in 1959 (his fig. 4a and table 1). This can be compared with our response function ΔA I (Equation (A14) with unit B′). Agreement would mean that the theories give similar results, because the response to any climate scenario can be expressed in terms of this response function (Equation (A18)). One problem is that his function was expressed in terms of the thickness change at the initial position of the terminus, rather than the area change ΔA. The relation between the two depends upon the details of the geometry of the glacier near the terminus, and to convert one to the other quantitatively would require a detailed geometric analysis, which we have not attempted. Instead, we merely multiply Nye’s thickness-change function with a factor which best normalizes it to our ΔA I function, as determined by the method of least squares.
The result is shown in Figure 6. Because of our normalization procedure, there is nothing to be learned by comparing the amplitudes of the functions, but it is seen that the time-dependencies are almost identical, the macroscopic result being slightly but not significantly faster. This means, by Equations (A1) and (A14), that the geometric mean of τ
A and τ
V, 
, is the same as its undefined equivalent in Nye’s theory, close to 20 years. The similarity of the two functions indicates that at least the time-dependence of Nye’s approach can be duplicated very simply by treating the glacier as a critically damped spring and mass even though his approach is much more complicated, because it involves the distribution of several quantities along the length of the glacier. The agreement may lie in part with his assumption that the change in reference-surface specific-balance rate (our terminology) is the same over the full length of the glacier. This input does not create any complex spatial detail that cannot be tracked by the macroscopic approach, which considers only the integrated, or glacierwide, balance rate.
The agreement in time-dependence is remarkable when one considers the different methods for parameter determination. Nye characterized the dynamics in terms of kinematic wave speed and diffusion coefficient along the glacier, giving the sensitivity of ice flux to changes in thickness and slope. These were derived from the observed glacier geometry and distribution of speed, using conventional flow theory, some of the limitations of which were discussed in our introduction. Our approach, which did not use measured speed at all, seems more direct and therefore is probably more accurate, but it does have the disadvantage that it requires measurement of area and volume changes over a period of years.
Nevertheless, the agreement in the time-dependence must be partly fortuitous. First, the large uncertainty in the numerical values deduced by both the macroscopic and traditional theories would by itself make the near-agreement in the time dependence surprising. Second, the glacier was larger and terminated in a lake during the period of the traditional study. Third, the traditional analysis does not take into account the effect of surface elevation on balance rate, which is potentially the most important difference between the approaches. This should affect both the amplitude and time-dependence of the response. For example, if we were to simulate the effect by setting our 
, τ
V would be decreased by roughly a factor of two, and the ultimate amplitude of the response to steady climate by the same factor.
5. Discussion
It is worth reviewing the origins of the five parameters which occur in this version of macroscopic theory. ΔA
0 is the amount by which the initial area is out of adjustment with the initial volume, τ
A and H characterize the dynamics of flow, and 
and 
are mass-balance rate quantities. τ
V and 
are defined in terms of these five and thus are not independent parameters. If the damping is known to be critical, there is a constraint among τ
A, τ
V and 
(Equation (15) with p = 1), which would leave a total of four parameters.
The parameters ΔA 0 and τ A account for the breakdown of perfect plasticity, at least approximately. This is the essential difference between this theory and its predecessor (Reference Harrison, Elsberg, Echelmeyer and KrimmelHarrison and others, 2001). The breakdown of perfect plasticity affects the response in two ways. First, τ A changes the time dependence of the “direct” part of the response, which would be represented by the ΔA S and ΔV S terms in the example of steady climate. It has been pointed out to us (personal communication from T. Jóhannesson, 2002) that in this case the non-zero τ A tends to decrease the time to complete the response, basically because the resulting lag of area with respect to volume increases the negative feedback represented by the area term in Equation (11).
The second effect of the breakdown of perfect plasticity is the appearance of the “transient” terms ΔA
T and ΔV
T in the response. These are proportional to ΔA
0, and approach constants different from zero at large time. It is significant that this second effect, although perhaps less obvious than the first, may cause the larger departure from the predictions of the perfect-plasticity theory. This will occur when the climate remains close to neutral 
for t > 0, but area is out of adjustment with volume at t = 0 so that ΔA
0 ≠ 0. A partial discussion of the importance of the transient terms for the case of steady climate was given in connection with Equations (16) and (17). In this case ΔA
0 affects the ultimate response but τ
A does not.
The ultimate response to steady climate is scaled by τ V. Thus τ V is more than a time-scale; it controls the amplitude of the response as well. Our definition of τ V is essentially the same as that of Reference Harrison, Elsberg, Echelmeyer and KrimmelHarrison and others (2001). It accounts explicitly for the effect of changing surface elevation on balance rate, which increases its value by roughly a factor of two for South Cascade Glacier, and probably much more for low-slope glaciers. This effect has been noted in numerical models by several authors (e.g. Reference OerlemansOerlemans, 1997). Reference Harrison, Elsberg, Echelmeyer and KrimmelHarrison and others (2001) point out that there are basic problems in finding an accurate value for τ V, even when the effects of surface elevation are included, and that this implies a fundamental difficulty in making an accurate prediction of response with any model, even a detailed numerical one. Our best value for τ V, which resulted from one of the best datasets in the world, has an uncertainty of at least 35%.
The balance between the achievements and limitations of our approach remains to be summarized. On the positive side the approach is essentially simple, and thus leads to new insights. To the points just discussed, we may add the ability to reproduce the time-dependence of the response of South Cascade Glacier with a simple critically damped spring and mass analogy, and the potential to determine the cumulative mass-balance history from the more easily measured history of area.
On the negative side there are several limitations. First, the approach suffers from a limited range of validation by observation, which ideally would be carried out for larger glaciers and for a period of more complicated trends in climate. If the recent, relatively steady climate at South Cascade Glacier persists, a second limitation may appear because the predicted ultimate changes in both area and volume are large, roughly a third of their initial (1970) values. These are about twice as large as the changes during the period in which we determined parameter values. Our theory may fail under such extreme conditions, because in common with the traditional one it is a first-order theory and its range of validity is unknown.
A third limitation is that our approach is essentially empirical. Our reasoning has been that testing our theory with data gives a more direct indication of its usefulness than would comparison with conventional theory, which has its own limitations. The main disadvantage with an empirical approach is that it does not provide the physics which determines the value of τ A. Thus we cannot estimate τ A when sufficient data for a direct determination are unavailable, nor address questions such as why the damping of South Cascade Glacier (which depends upon τ A) is approximately critical. A comparison with the conventional theory at a basic level may provide the missing physics, and is an obvious next step.
Acknowledgements
We are grateful for the significant improvements to the paper which resulted from the comments of D. H. Elsberg, J. F. Nye, A. Rasmussen, M. Truffer and the referees, G. H. Guðmundsson and C. S. Hvidberg. The input from the Scientific Editor, T. Jóhannesson, was especially important. The primary support was from the U.S. National Science Foundation, grant No. OPP-9707515.
Appendix
Solutions
A.1. Non-dimensionalization
Equations (12) and (13) can be simplified by writing them in terms of the non-dimensional time t* defined by
This choice is the most useful because it reduces to one the number of parameters on the left hand sides of these equations. The result is
where p is a non-dimensional parameter defined by
p characterizes the strength of the damping and determines the mathematical form of the solution.We focus on the special case of critical damping, p = 1, but also give an example of the situation p → ∞. The general case is slightly more complicated.
A.2. Solution for steady climate
The simplest case is 
a constant for t > 0. The critically damped solution for ΔA is
Where
and
The corresponding solution for ΔV is
where
and
Because these results apply for critical damping only, τ V and τ A are not independent but related by Equation (A4) with p = 1.
The special case τ A → 0 (p → ∞) can be shown to give
(e.g. Reference Harrison, Elsberg, Echelmeyer and KrimmelHarrison and others, 2001).
A.3. Solution for an impulse
Another simple case is an impulse in climate, which means that 
is very large at t = 0 for a short interval. The strength of the impulse, which we designate by B′, is 
integrated over this interval. The solution for ΔA can be found by superposing the solutions due to two step-changes of opposite sign, slightly offset in time. Using the same decomposition of the solution into two terms as before,
one finds
ΔA T is given by Equation (A7).To obtain the solution for ΔV one first needs the appropriate initial conditions, which are
and
The first arises because the impulse gives an immediate discontinuity in ΔV of B′. The second follows from Equation (11). The resulting ΔV I is
ΔA T and ΔV T are given by Equations (A7) and (A10). As before, τ V and τ A are related by Equation (A4) with p = 1.
A.4. General solution
The general solution for a climate scenario specified by any 
can be written in terms of the impulse solution:
where (Δa
I, Δv
I ) ≡ (ΔA
I/B′, ΔV
I/B′), the unit-impulse responses obtained from Equations (A14) and (A15). As before, τ
V and τ
A are related. The response for harmonic variation of 
, for example, follows.
