1. Introduction

A theoretical treatment has been given in other papers (Reference NyeNye, 1960, Reference Nye1963[a], Reference Nye1963[b], Reference Nye1965) referred to as Nye [I], Nye [II], Nye [III] and Nye [IV] respectively, of the relation between the fluctuations of accumulation and ablation on a glacier and the associated fluctuations of the glacier itself. As a matter of history, while the variations in accumulation and ablation on glaciers have only been systematically measured since about 1945 the records of the changes in length go back in some cases for 300 yr. The geomorphological record of course goes back, with less detail, into the Pleistocene. The possibility therefore presents itself of taking the known records of changes in length of the glaciers and trying to infer the changes of accumulation and ablation that caused them, thereby considerably extending the record of accumulation—ablation history into the past. The present paper takes the differential equations of the theory as its starting point, and considers how to apply them numerically to deduce the input (accumulation—ablation record) that gave a given output (advance—retreat record). The calculation scheme thus developed is applied to the terminus records of South Cascade Glacier, Washington, U.S.A., and Storglaciären, Kebnekaise, Sweden, since the functions that characterize these two glaciers are known from previous work.

The method described in Nye [III] for solving the same problem, which uses a series of coefficients λ _n , entails first filtering out the high frequencies in the terminus record and then computing time derivatives. By contrast, the method to be described here works directly from the annually spaced observations of the terminus without filtering; it is comparatively simple to apply and is, in general, much better than the earlier method. However, the method using the λ _n . series is still useful in theoretical studies where the terminus record is taken as a polynomial in the time t.

2. The Differential Equations

The notation follows that used in the previous papers. x is the distance down the glacier, t is the time. A datum glacier (suffixes 0) is considered in which the accumulation rate (thickness of ice added per unit time) is a ₀(x). a ₀(x) is negative in the ablation area. A linear perturbation treatment [Nye III] leads to the two simultaneous partial differential equations

(1)

(2)

where q ₁(x, t) is the perturbation in discharge from the steady datum rate, h ₁(x, t) is the perturbation in thickness, and a ₁(x, t) is the perturbation in accumulation and ablation rate. B ₀(x), c ₀(x) and D ₀(x) are fixed functions that characterize the particular glacier. B ₀(x) is the datum-state width of the channel, c ₀(x)/B ₀(x) is the velocity of kinematic waves of q ₁ down the glacier, and D ₀(x)/B ₀(x) is their diffusion coefficient. (1) arises from the equation of continuity, while (2) is a consequence of the assumed nature of the flow process. If q ₁ is eliminated between (1) and (2) we obtain a type of diffusion equation

(3)

where the primes denote differentiation with respect to x.

We shall mostly be concerned with the case where a ₁(x, t) is independent of x, so we simply write it as a ₁(t) Thus a ₁(t) is the function representing the variation in accumulation rate that we are going to try to deduce. It represents the increase in the accumulation rate or the decrease in the ablation rate from the steady state value a ₀(x) at any point on the glacier. Other assumptions about the x dependence of a ₁(x, t) could be made, but some assumption must be made if the basic problem of deducing budget history from advance—retreat history is to have a solution. Otherwise we should be trying to deduce a function of two variables from a function of only one variable.

The given data, besides the coefficient functions B ₀(x), c ₀(x), D ₀(x) could be the time variations of thickness, that is, of h ₁ at some fixed point on the glacier. What is more commonly observed, however, is the variation in the position of the end of the glacier. We denote this by l ₁(t), measuring l ₁ from the end of the datum glacier. l ₁ is geometrically related to h ₁ at the wedge-shaped end of the datum glacier (fixed at x = L) by the equation

(4)

where θ is the angle of the wedge. The problem is then, given l ₁(t), to find h ₁(L, t) from (4), and then find what source function a ₁(t) in (3) will produce the required h ₁(L, t).

3. The use of Influence Coefficients

The record of l ₁(t) and therefore of h ₁(L, t) is not, in practice, a continuous one. Observations are commonly made at the end of the ablation season at intervals of one year. One way of filling in the gaps in the record would be to make some assumption about the winter-summer variation between the observed points. If we did this our final result for a ₁(t) would include the winter-summer variations. We are not usually interested in these, and in any case they would merely reflect what was assumed for the winter-summer variation of l ₁(t). To pursue the same line of thought, the true function l ₁(t) includes also very high-frequency detail due to individual snowfalls, hot days and so on, and all this is faithfully reflected in the true variation of a ₁(t) by the simple relation

We do not have this detailed record of h ₁(L, t) nor do we wish to know the corresponding detail in a ₁(t) In this situation we propose to adopt a different way of filling in the gaps in the record.

We consider an artificial function a ₁(t) that is constant throughout any one (budget) year. It thus consists of a number of steps (Fig. 1a). The corresponding function l ₁(t) consists of a number of nearly straight (more accurately, exponential) segments, with discontinuities in slope (Fig. 1b). The observations are taken to be the points p₁, p₂, p₃, .… The method of deducing a ₁(t) that we describe is valid when the time Δt = 1 yr. is sufficiently short to be taken as the time interval in a certain finite difference solution of the equations. Now it could happen that a time interval of 1 yr. was loo long for this to be true, but in fact the time constants of glaciers are long enough for the approximation to be a good one. (We verify this later for the individual glaciers treated.) This is a fortunate circumstance for it means that annual observations may be used to infer the step curve of a ₁(t). Thus we may say that if the real a ₁(t) actually followed the step curve (although we know it does not), then the resulting l ₁(t) curve would pass through the observed positions. That is the precise meaning of the a ₁(t) step curve that we deduce.

Fig. 1. (a) Assumed step function a ₁(t). (b) Resulting response funetion l ₁(t)

To a good approximation this step curve may be taken as the mean over the year of the actual a ₁(t) function. This mean is closely connected with what is usually, but wrongly, Footnote * called the net budget. The net budget of a glacier at the point x is, in our notation,

where Δt = 1 yr. The second term in the right-hand bracket is what we deduce. When added to a ₀(x) and multiplied by Δt it gives the net budget of the glacier at the point in question. In short, the step function a ₁(t) that we deduce gives the time variation in the net budget of the glacier.

The numerical method to be described comes from a finite difference approximation to the original differential equations. The time derivative ∂/∂t is replaced by a centred difference using an interval Δt = 1 yr. This means that the equations contain values of the dependent variables at one-year intervals; in particular, they contain the h ₁ values at the instants of observation. As we have just noted, Δt = 1 yr. is short enough for this approach to be feasible.

Let the values of h ₁ at given x, at intervals Δt, be denoted h(x,  1), h(x, 2), …, back into the past, beginning at a certain instant shown as r in Figure 2. h(x,  n) is thus a continuous function of x but a discrete function of time. ‘l’he values of q ₁ are similarly denoted q(x, 1), q(x, 2), …. a ₁(t), which is now taken as constant over each interval, has values a(1), a(2), …. Note in Figure 2 that these values are placed at the centres of the intervals. We now assume plausibly that

(5)

(6)

(7)

where e(x, 1), e(x, 2), …, are influence coefficients to be determined. Each a(n) in the past, but not in the future, will have an influence on h(x, 1) and since the equations are linear their effects will be additive. This gives equation (5). Equations (6), (7), etc., make the same statement for h(x, 2), h(x, 3), etc., and, since the theory is concerned with small perturbations, they use the same set of coefficients e(x, n).

Fig. 2

We shall show in a moment how the e(x, n) may be computed numerically. For a particular x, which could be x = L we write the equations simply as

(8)

If these equations are now solved for a(1), a(2), …, in terms of h(1), h(2), …, we find

(9)

where the g(n) are a new set of coefficients which, by substitution in (8), satisfy

(10)

Notice that the matrix of g’s in (9) has the same form as that of the es in (8), each row being repeated but displaced by one position.

Equations (10) allow the g’s for a particular x to be found in succession once the es are known, by a simple algebraic process. Equations (9) are the result we need for our problem, for they give the sequence a(1), a(2), …, directly in terms of the observed sequence h(1), h(2), …

The g’s are thus our objective, but to calculate them we must first find the e’s. For this purpose first assume a set of equations for the q’s similar to (5), (6), (7), etc.:

(11)(12)(13)

where the f (x, n) are influence coefficients for q. Now write the original differential equation (1) for the point p in Figure 2 and replace ∂h ₁/∂t by {h(x, 1)−h(x, 2)}/Δt. Approximate ∂q ₁/∂x at p by

where the primes denote differentiation with respect to x. If we now substitute into equation (1), using (5) and (6), and equate coefficients of a(1), a(2), …, we obtainFootnote *

(14)

where K(x) = 2B ₀(x)/Δt. (14) is thus a finite difference approximation to the original differential equation (1).

Now write the original differential equation (2) for the point r (not p) in Figure 2, substitute from (5) and (11), and equate coefficients of a(1), a(2), …, to obtain

(15)

Remembering that all quantities are functions of x and dropping the x’s for brevity, the systems of equations (14) and (15) may be written

(16)

where the Z(n), which are functions of x, are given by the recurrence formulae

(17)

More compactly we may write (16) as

(18)

Eliminating f(n) gives the alternative form

(19)

Starting with n = 1 we therefore have a second-order ordinary differential equation to solve for e(x, 1). When this has been done Z(2) may be found from (17), and then (19) with n = 2 gives an equation to solve for e(x, 2), and so on. Thus the influence coefficients e(x, 1), e(x, 2), …, may be computed in succession.

The nature of the above process is made clearer if we notice the physical interpretation of the e(x, n) coefficients. If all the a’s are zero except for a single one, say a(r), which is equal to 1, equation (5) shows that h(x, 1) = e(x, r). Thus e(x, r) is the h ₁, response at time 1 to a unit pulse a(r). If a unit pulse of a ₁, occurs as shown in Figure 3, the response h ₁(t) at given x will be a certain curve. The height of the curve at time rΔt from the beginning of the pulse is e(x, r). In short, the sequence e(x, n) represents the h ₁ response at intervals Δt to a unit pulse of a ₁. Similarly the sequence f(x, n) represents the impulse response of q ₁.

Fig. 3. Interpretation of the influence coefficients e(x, r) as the response to a unit pulse

The same interpretation follows in the continuous case (Δt→0). Here, by analogy with (5) we should write the integral equation

(20)

If we apply a unit pulse (Dirac δ-function) of a ₁ at time t−t ₀, the only contribution to the integral comes from τ = t ₀, and we obtain

Thus e(x, t ₀) is the response h ₁ at time t to a unit pulse that occurred at time t−t ₀. Similarly for f.

Accordingly, if we wished to find the continuous function e(x, τ) we should try to solve equations (1) and (2) for an impulse of a, occurring at time t = 0, say. Hence we should wish to solve the equations

for initial condition e = f = 0, where we have simply put e for h ₁ and f for q ₁ in (1) and (2).

Now our equations (18) are merely a finite difference form of this pair of equations. If a rectangular pulse of unit height is placed in the first interval (Fig. 3) and if ∂f/∂x is approximated by

So far only t has been put in discrete steps, not x. If, however, we write the simplest type of central difference for the x derivatives in (18) we find we have precisely the Crank—Nicholson scheme for integration of the diffusion equation (Reference Crank and NicholsonCrank and Nicholson, 1947; Richtm)er, 1957, p. 93). This is in fact the scheme we have used.

4. Boundary Conditions

It was shown [Nye II, III] that the functions B ₀(x), c ₀(x), D ₀(x) that characterize the glacier may be taken to have the following behaviour at x = 0 and x = L, the end points of the datum glacier: near x = 0, B ₀ non-zero, c ₀ ~ x, D ₀ ~ x ²; near x = L,  B ₀ non-zero, c ₀ non-zero, D ₀ ~ L−x. The simultaneous equations (18) for a given n then have the same form as those studied in appendix I of Nye [IV], with e and f substituted for H and Q. From this previous work we know that a unique solution is determined if we require as boundary conditions that e be not infinite at either x = 0 or x = L. This condition at x = 0 is equivalent to the condition f = 0. The boundary condition at x = 0 gives a one-parameter set of solutions which all have the same leading term, namely

(21)

Z(n),

(22)

c ₀,

To summarize, at x = 0 we have a simple condition on e, namely (21); while at x = L we have a condition involving a linear combination of e and e′, namely (22).

5. Solving the finite Difference Equation

Points on the x axis are numbered backwards from j = 1  (x = L) to j = J (x = o) at intervals Δx. In finite difference form equation (19) at any given n is then

(23)

where

• A_j = −a_j −γ_j ,

• B_j = −2a_j −N_j ,

• C_j = −a_j +γ_j ,

• a_j = (D ₀) _j /(Δx)²,

• γ_j = V_j /2Δx,

• 
  ,

• 
  ,

Thus the coefficients A _j , B _j and C _j in (23) may be computed once and for all; they do not depend on n.

The boundary conditions give, from (21),

(24)

and from (22)

(25)

To solve the second-order difference equation (23) directly we should need to know two starting values at one end. Instead we have one condition at each end. This boundary value problem is easily dealt with by the method described by Richtmycr (1957, p. 103–04).

We discuss the one-parameter family of solutions S satisfying the boundary condition at x = L (j = 1). They will satisfy a certain first-order difference equation, say

(26)

and our problem is to find E _j and F _j so that (26) is true for any member of S. Since both (25) and (26) are true for all members of S, it follows that

(27)

By substituting E _j−1 e _j +F _j−1 for e _j−1 in (23) and comparing with (26) we find

Since E ₁ and F ₁ are now known from (27) these equations allow E ₂, E ₃, … and F ₂, F ₃, … to be computed in succession.

Equation (26) may now be used to compute the e _j in succession, starting with j = J−1 and using the value of e _J given by (24). The method is very efficient.

Accordingly, the numerical procedure is to solve equation (9) by the above method for n = 1. Z(2) is then found from (17) and equation (19) is solved for n = 2, and so on. This determines the e(x, n).

Some methods of this sort for solving the diffusion equation by forward integration in time suffer from instability when the time step is too long in comparison with Δx (Reference RichtmyerRichtmyer, 1957, p. 9). However, the Crank-Nicholson method, which we are using, is unconditionally stable. This important result was secured essentially by taking central differences at p in Figure 2 rather than, say, one-sided time differences at r.

We have now described how, from the functions B ₀(x),  c ₀(x),  D ₀(x), the influence coefficients (ex, n) in equations (5), (6), (7), etc., may be obtained numerically—essentially by forward integration of the diffusion equation in time steps of 1 yr. The next step is to invert the e array by equations (10) and thus obtain the g’s. This need only be done for the point x that is of interest, normally x = L, the end of the datum glacier. Having found g(1), g(2), …, they may be applied directly to a record of the changes in terminus position to deduce the time variation of the net budget.

A check on the values obtained for the e(n) and g(n) may be made by considering the case of a perturbation steady in time. Thus if a ₁(t) = A say, h ₁ , will come to the steady state H(x) say. ln the notation of Nye [III] we have

In the notation of this paper, we put

in equations (8) and (9) to obtain

Thus

(28)

If μ ₀(x) and λ ₀(x) are known by some other method (for example by integration of equations (1) and (2) with ∂/∂t = 0, or from the frequency response at zero frequency) these formulae give a check on the values of e(n, x) and g(n, x).

Interpretations may also be given of the other coefficients μ ₁, μ ₂, … and λ ₁, λ ₂, … introduced in Nye [III]. For example, take the series

and let h ₁ = pt, where p is a constant. Then, at t = 0,

But by (9) we find for a ₁ the expression

Hence

In general, by putting h ₁ = pt^m we find

A similar equation relates μ _m with the e(n).

A remark about the convergence of the two series (8) and (9) is in place before we proceed to the numerical work. Experience with the response of the special model [Nye II] leads us to expect that the impulse response of a glacier will persist for several hundred years. Thus several hundred e(n,  x) will be needed if (8) is to be used. On the other hand, we found [Nye II] that only a much shorter past history of h ₁ was needed in calculating a ₁. So we expect that the g series in (9) will converge much more rapidly than the e series in (8). To calculate, say, the first 50 g’s we do not have to calculate more than the first 50 e’s. Even though at n = 50 the e’s are still very appreciable, we expect that the g’s will be quite small.

Let us now sec how these procedures work out in practice.

6. Application to South Cascade Glacier

(i) Computing the influence coefficients

The above method for finding the e(n,  x) and g(n,  x) was programmed for an I.B.M. 709 digital computer. The functions B ₀(x), c ₀(x), D ₀(x) could be read in in tabular form and the programme made linear interpolation as necessary. The range of x was divided into 500 intervals and Δt was taken as 1 yr. The B ₀(x), c ₀(x), D ₀(x) functions for South Cascade Glacier, Washington, U.S.A. have been derived [Nye III] from the field observations of Reference Meier and TangbornMeier and Tangborn (1965). They gave the values of the e(n) and g(n) at x = L shown in Table I and plotted in Figures 4a and b. The impulse response e(n) behaves as expected (Fig. 4a). We recall that e(n) is the response of h ₁ to a pulse of a ₁ of unit amplitude, say 1 m./yr., lasting for 1 yr. from n = 0 to n = 1. One year after the start of the pulse h ₁ = 1.084 m. The pulse now ends. After a further year h ₁ = 1.275 m. and so on down the table. The fact that h ₁ grows at first, rather than decays, is the instability phenomenon treated in Nye [I]. The annual increment in h ₁ also grows at first, becomes a maximum between n = 6 and 7, and then diminishes. h ₁ itself reaches a flat maximum of 4.967 m. at n = 23—the maximum of the flood wave—and then slowly diminishes. Its fastest rate of fall is at n = 43 to 44. After too years, h ₁ is still quite appreciable, 0.622 m. This, remember, is the result of a mere 1 m. of ice added to the glacier 100 yr. previously.

Table I Values of e(n) and g(n) for Soute Cascade Glacier and Storolaciaren

Fig. 4a. Influence coefficients e(n) for the termini of South Cascade Glacier and Storglaciaren. The curves also represent the response to a unit pulse

Fig. 4b. Coefficients g(n) for South Cascade Glacier. The coefficients are shown as a step curve to emphasize their interpretation as a series of a(n). Scale at left refers to filled-in curve. The unfilled-in curve (scale at right) shows the same data with scale enlarged 20 times

Fig. 4c. Same as Figure 4b, but for Storglacidren. Scale at left refers to filled-in curve. The unfilled-in curve (scale at right) .shows the same data with scale enlarged 4o times

The behaviour of g(n) is more complicated (Fig. 4b). g(1) = 0.922 and g(2) = −1.085. After this the variation settles down to an oscillation of decreasing amplitude. The g’s diminish much faster than the e’s, as forecast; after n = 46, |g(n)| < 0.00001. The values of g(1) and g(2) may be understood by considering the artificial situation where the glacier was a perfect integrator. Then all the e(n) would be simply Δt( = t). a, would then be ∂ _h1/∂t or a(1) = {h(1)−h(2)}/Δt. Comparing with equations (9) we see that

Thus the departures of g(1) and g(2) from the values t and —1, and the subsequent oscillations of g(n), measure the extent to which the glacier is not a perfect integrator.

It is also illuminating to notice that after a pulse of unit height and length Δt the thickness of new material added to the glacier is Δt. Hence as Δt→0, e(1) and e(2) both approach Δt. It then follows from the first two equations of (10) that g(1)→(Δt)⁻¹ and g(2)→(Δt)⁻¹ (In the continuum limit the g’s in fact involve the first derivative of a s-function.)

The physical meaning of the g(n) sequence is found by regarding the g(n) as a series of a1 values moving forward in time. The g(n) sequence represents that time sequence of a ₁ which would be necessary to produce a time sequence of h ₁in which the first value was 1 and all others were zero.

As a check on the numbers obtained for the two sequences we find

compared with

[Nye IV, table I, col. (2)]. The first too terms thus give excellent agreement with equation (28) for the g’s; but, as is clear from their slow rate of decay, more than ion e’s would be needed before their sum attained the limiting value μ ₀ to 3 figures.

As a check to see whether Δt = 1 yr. was small enough the time interval was halved, but the pulse length was kept the same, now running from n = 0 to n = 2. This showed that the e(n) computed with Δt = 1 yr. do indeed represent the true response to a unit rectangular pulse of duration 1 yr. to an accuracy of 0.3 per cent. We conclude that the basic premise, that we can work solely with annual values, is justified.

(ii) Application of the g(n) to the terminus record

The record of the position of the terminus of South Cascade Glacier (Reference Meier and TangbornMeier and Tangborn, 1965) (Fig. 5a) consists of observations in 1928,Footnote * 1953, and annually from 1955–63. The total length is thus 35 yr. The record of a ₁(t) (shown shaded in Figure 5b)Footnote † extends from 1952 to1963 with a gap for 1953/54. It is not unreasonable to apply the foregoing theory and see to what extent the theoretical a ₁(t) agrees with that observed. (Some comparison has already been made in Nye [III].)

Fig. 5a. Record of terminus position of South Cascade Glacier.● observed positions (Reference Meier and TangbornMeier and Tangborn, 1965), ○ interpolated positions

Fig. 5b. South Cascade Glacier. Time variation of net budget. Shaded step curve: observed, inferred from data of Reference Meier and TangbornMeier and Tangborn (1965). Unshaded step curve: computed from unsmoothed terminus record

Fig. 5c. South Cascade Glacier. Time variation of net budget. Dots show 10 yr. running means of values computed from terminus record; open circle, observed 9.yr. mean

We must first fill in the gaps in the terminus record l ₁(t); this was done arbitrarily by linear interpolation for the years 1929–52 and for 1954. The record was also extrapolated back linearly before 1928. The angle θ in equation (4) was taken as the current value of 6.7 degrees. Then the g(n) gave the record of a ₁ by equations (9); it is shown in Figure 5b (unshaded). There is some correspondence with the observed a ₁ record but it is not at all good; peaks occur in the right places, but they are all accentuated. The best agreement would be expected for the most recent points; detailed agreement ought to fall off into the past as the amount of information in the terminus record diminishes. On this basis the 1962/63 value is discordant. (The theoretical values for 1953/54 and 1954155 are replaced by their mean because separately they depend heavily on the missing l ₁ value for 1954.) The computed a ₁(t) before 1953 is linear and reflects the assumed linear retreat.

We have given the direct comparison between theory and observation. Now we must look a little more closely at possible sources of difficulty. To begin the discussion let us take note of the following well-known relation which must hold at the terminus of a glacier ending on land:

(29)

where

is the rate of advance measured along the bed, —a ₁ is the rate of ablation measured perpendicular to the ice surface, u _t is the ice velocity parallel to the bed, and θ is the angle between the upper surface of the ice and the bed. (If ablation occurs at the bottom surface of the ice an extra term will appear.) This relation is a geometrical requirement and is quite independent of any theory of the sort we have been dealing with. The reason this simple formula cannot be used to relate the changing ablation-rate to the retreat of the glacier is that changes higher up the glacier will continually be propagated down and will change the values of θ and u _t . It could be said that our purpose in formulating the theory with kinematic waves and diffusion is precisely to deal with this difficulty. Nevertheless, over short enough periods θ and u _t will stay relatively constant, and then

will simply be linearly related by (29) to a _t . Detailed consideration shows that over a short period like 6 yr. on South Cascade Glacier it is reasonable to take u _t and θ as constant. Let us then do so and plot

against a _t for the various years (Fig. 6). The points show a good deal of scatter. The straight line shown corresponds to the measured values θ = 6.7 degrees and u _t = 31.4 ft./yr. (9.6 m./yr.). The errors shown for

are obtained from the errors in the individual

values given by Meier and Tangborn. The average error in the ablation at the snout is perhaps ±1 ft./yr. (0.3 m./yr.), so one sees at once that it is the error in measuring

that is the source of the trouble. The ablation is measured with ample precision in most years, but the retreat rate is not.

Fig. 6. South Cascade Glacier. Showing the poor correlation between

, the measured annual retreat of the terminus, and −a _t , the measured annual ablation at the terminus

Now South Cascade Glacier ends in a small lake, and it may be that this is the main cause of the scatter. For, in the first place, the occasional calving of the glacier into the lake causes abrupt changes in the terminus position; secondly, the geometry of the terminus does not strictly fulfil the conditions assumed in the derivation of equation (29). In these circumstances the proper comparison to make would be with the changing height of the ice surface at a transverse section near the terminus, rather than with the changing terminus position, if these data were available for a sufficient number of years. Unfortunately they are not. Thus, the relevance of the theory to the existing data on year-to-year changes of this particular terminus is questionable. The calculation of a ₁ with the g coefficients depends heavily on the current retreat rate; it is therefore not surprising that the correspondence in Figure 5b is rather poor.

Although it is evident that the details of the retreat curve cannot be used with much confidence, the average retreat rate over say 10 yr. is surely still significant. With this thought we have taken 10 yr. running means of the computed a ₁ values. They are plotted on Figure 5c. For comparison we show the measured 9 yr. mean for 1954/55–1962/63. The agreement is very good. With the short terminus record of South Cascade Glacier we cannot go much further than this, except to deduce that on average the net budget was only very slightly negative in the 1930’s and 1940’s. It seems that for year-to-year fluctuations the geometry of the terminus and calving prevent effective use of the theory. But for longer-period trends the comparison between observation and calculation is favourable for the theory.

7. Application to Storglaciaren

(i) Computation of the influence coefficients

In Nye [IV] the observations of Dr. V. Schytt and Dr. E. Woxnerud on Storglaciären, Kebnekaise, Sweden, were used to deduce B ₀(x), c ₀(x) and D ₀(x) for this glacier. These tabulated functions were then used to find the e(n) and g(n) at the snout of Storglaciären shown in Table I and in Figures 4a and c. The general behaviour of e(n) is the same as we have already seen but the time scale is about twice as long. The maximum of e(n), the peak of the flood wave, is reached at n = 48 (n(48) = 7.100) instead of n = 23 for South Cascade Glacier. The check whether Δt = 1 yr. was small enough for e(n) to be a true representation of the impulse response showed that the maximum error was 0.08 per cent; the figure for South Cascade Glacier was 0.3 per cent; the smaller value for Storglaciären again reflects its longer time-scale.

The g’s (Fig. 4c) do not fall off as fast as one would wish:

The error incurred by not using enough g’s depends, of course, on the record of h ₁; but if the record were steady we should obtain 0.0014 h ₁ instead of 0.0017 h ₁ as the value of a ₁, an error of 18 per cent.

(ii) Application to the terminus record

The terminus record of Storglaciären (Fig. 7a) goes back to 1897 with annual observations since 1944 (personal communication from V. Schytt; [Reference AhlmannAhlmann and others], [1950]. The record of annual net budget obtained by Dr. Schytt and Dr. Woxnerud (personal communication from V. Schytt) extends from 1945 onwards.

Fig. 7a. Record of terminus position of Storglaciären. Dots show observed positions (personal communication from V. Schytt, [Reference AhlmannAhlmann and others, ] [1950]). Curve shows smoothed position used for computation

Fig. 7b. Slorglaciären. Time variation of net budget. Shaded step curve: observed, inferred from data of Stkytl and Woxnerud (personal communication from V. Schvtt). Unshaded step curve: computed from unsmoothed terminus record

Fig. 7c. Storglaciären. Time variation of net budget. Step curve is computed from the smoothed terminus record. Open circles show 10 yr. running means of observed values

The correlation between a _t and

implied by equation (29) is again poor (Fig. 8). As with South Cascade Glacier the scatter is presumably due to difficulties with

rather than with a _t ; but there is no lake to blame in this case. The difficulty seems to be one of sampling and averaging. Schytt points out (personal communication) that the retreat varies appreciably across the ice front and that the terminus shows relief in the transverse direction—that is, a variation both in position and in angle θ. During some summers parts of the front may stay snow-covered, and again during other summers moraine-covered. Nevertheless, in spite of these complications, if all quantities in equation (29) are measured at strictly the same point on the terminus, the equation should hold precisely—apart from bottom melting. It should also hold precisely if all quantities are averaged, in an appropriate way, across the terminus. It is presumably the difficulty of finding precise averages from a limited number of points that prevents a better year-to-year correlation between a _t and

, interpreted as averages across the glacier. A further reason for discrepancy is that the survey of the front is very seldom simultaneous with the end of the budget year; the error in l ₁ from this source could be about 2 m. (personal communication from V. Schytt).

Fig. 8. Storglacicären. Showing the poor correlation between

the measured annual retreat, and −a _t , the measured annual ablation at the terminus

In this situation it is plain that we must be very cautious in attaching significance to the observed annual fluctuations in

. Bearing this clearly in mind, let us nevertheless see what happens if, just for the moment, we take the observed l ₁ values literally and fill in gaps by linear interpolation and extrapolation (not shown in Figure 7a). The result for a ₁ (using θ = 18.2 degrees) is shown in Figure 7b by the unshaded step curve. The sharp changes at 1908 and 1922 merely reflect the arbitrary assumption made to fill in the l ₁record.

We now wish to compare this theoretical a ₁(t) curve with observation. As illustrated in Nye [IV], the observed a ₁(x, t) shows some dependence on x from year to year on Storglaciären, although when averaged over several years the approximation a ₁(x, t) = a ₁(t) is roughly obeyed. The question arises therefore of how best to interpret the a ₁(t) of the theory in this case. One way is to take it as an average down the glacier. However, it seems better to take it as the a ₁(t) measured at the terminus. By doing this we preserve the theoretical correlation between rate of retreat and ablation rate at the terminus for short-period changes—and for long-term changes it does not matter what point on the glacier is chosen for measuring a ₁(t). In this way we may hope to make use of the short-period information in the computed a ₁—which we should otherwise have to discard as not meaningful. Accordingly, the a, values shown in Figure 7b as observed (shaded step curve) are in fact a ₁ at the current glacier terminus. There is some relation between observed and calculated values. That the relation is not closer we believe to be due primarily to the difficulty of measuring

(averaged across the terminus), as evidenced by Figure 8.Footnote *

Again, as with South Cascade Glacier, it seems we must attach little weight to the annual fluctuations of

, but we should be entitled to trust the general trend. Therefore we have drawn a smooth curve (Fig. 7a) through the points for 1897, 1908, 1922 and 1929, and passing evenly among the 1944–62 points. The a ₁ step curve calculated from this smooth curve is shown in Figure 7c. We believe the trends to be significant, but of course the values must become less reliable as they go into the past, and as the information in the l ₁ record diminishes. For example, the linear trend before 1908 is purely the result of the linear interpolation between 1908 and 1897 and the quite arbitrary linear extrapolation before 1897. The sharp down-turn in the l ₁ record about 1920 is reflected in the sharp decrease of a ₁ at this time (note, no time lag). But, interestingly, even though a ₁ returned to zero about 1930 the retreat still continued, because of effects from the early 1920’s propagated down the glacier. Looked at in another way, if a ₁ had not risen in the 1930’s the retreat caused by the pulse in the 1920’S would have accelerated, by the instability effect—for the e(n) show that the greatest rate of retreat after a negative pulse occurs 10 yr. after the pulse. However, a reduction in a ₁ in the late 1930’s was needed to sustain the retreat. The theory finds that a ₁ has been fairly constant between −0.3 and −0.7 m. of ice/yr. since 1940. The 10 yr. running mean of the observed a ₁ (at the terminus) is also shown in Figure 7c. It agrees quite well. By taking other smooth curves through the l ₁ data the agreement could probably be made even better.

One complication should be mentioned (personal communication from V. Schytt). Our theory assumes that the slope of the bed is much the same at the successive positions of the terminus. Before about 1935–40 the slope of the bed was in fact opposite in sign to the slope of the ice surface; subsequently it was the same. The main effect, presumably, will be to make θ larger for the earlier years than the value we have taken. Thus, for these years, we have underestimated the absolute magnitude of a ₁.

Our conclusion is that, just as for South Cascade Glacier, the observed year-to-year fluctuations of the terminus cannot be used with any certainty to obtain the annual budget; but the general trend of advance and retreat over longer periods gives a mean a ₁ that agrees well with recent observation. We therefore think that the curve in Figure 7c, with its maximum in the 1930’s, has genuine significance, and extends the budget record back into the period before it was measured. In saying this we must repeat the proviso that the curve necessarily becomes less trustworthy as it goes into the past.

8. Acknowledgements

I must again express warm thanks to Dr. V. Schytt, Dr. E. Woxnerud, Dr. M. F. Meier and Mr. W. Tangborn for making available their observations before publication. Dr. Schytt and Dr. Meier have also reviewed the manuscript most helpfully—but, naturally, they must not be held responsible for its shortcomings. Most of the computations were done on the I.B.M. 709 digital computer at Yale University; a few were done on an I.B.M. 7094. digital computer at the University of California at Los Angeles.