Introduction

Berendon Glacier (lat. 56° 15’ N., long. 130° 05’ W.) is situated in the Coast Mountains of British Columbia, Canada. Reference Untersteiner and NyeUntersteiner and Nye (1968) (referred to in this paper as U-N) calculated the frequency response of this glacier and made predictions about its possible future movement. They were particularly concerned to see whether, in the next 40 years or so, the glacier was likely to advance sufficiently to be a danger to the mining installations of the Granduc Operating Company, which are situated near the snout of the glacier. The data they used, however, "necessarily involved a good deal of guesswork". Since the publication of U-N, detailed field measurements have become available (Reference Rogerson and StanleyRogerson, in press; Reference StanleyStanley, in press). It was decided, therefore, to repeat the earlier calculations using these improved data, to see how they affected the predictions and also to see how sensitive the theory is to changes in the input data.

Procedure for Setting up the Datum State

The theory used, which is concerned only with changes in the glacier that are small compared with the overall dimensions, is developed in five papers (Nye, Reference Nye1960, Reference Nye1963[a], Reference Nye1963[b],1965[c]). The nomenclature used here is the same as that used in those papers. The datum glacier is characterized by three functions B_o(x), c _o(x), D _o(x) where x is a coordinate denoting distance down the glacier, the functions being, respectively, the width of the glacier in a datum state, the velocity of kinematic waves multiplied by B_o , and the diffusivity of kinematic waves multiplied by B_o . The suffix zero signifies the datum state.

To estimate these functions for the Berendon Glacier we used the same procedure as that described in Nye (Reference Nye1963[b] and Reference Nye1965[b]). Since the essential difference between this paper and U N is the data used in determining B_o(x), c _o(x) and D _o (x), we now consider in detail the data used in the two cases.

The data used by U-N were as follows:

1. Three 1:12 000 topographical maps for 1949, 1956 and 1963. None of these maps, however, showed the upper reaches of the accumulation zone, and here U-N had to resort to a 1:250 000 map.

2. The cross-sectional area, S_o(x), of the glacier. The bedrock topography was not known and had to be inferred from cross-sections of the valleys in the lower part of the glacier and estimated in the upper part.

3. The variation of the net mass balance, a, with x, estimated by comparing Berendon with eight other glaciers along the coast between Washington and Alaska.

The data used in this work were:

1. A new 1: 10 000 map made in 1968, which has closer contour intervals (10 m) than any previous map and which gives complete coverage of the upper accumulation areas. From this, B(x) is determined and B(x) is assumed to be B_ox).

2. Net mass-balance data for 1967, 1968, and 1969 (Rogerson, in press; Stanley, in press). These are shown as a function of x in Figure 1. For 1967/1968 the data are extrapolated for x < 1.8km. From these measurements, a datum state mass balance is derived by taking the mean of the two years’ results and displacing it downwards by 0.72 m to satisfy the condition (Reference NyeNye, 1963[b]), and this is also shown in Figure 1.

3. Measured velocities for the bottom half of the glacier for 1967/68 and 1968/69.

Fig. 1. Mass-balance data for Berendon Glacier. The data for 1967/68 and 1968/69 are shown, together with the mean of these and also the datum state chosen.

A slightly different method was used to obtain the datum velocities u _o(x). U-N inferred their u_o(x)’s from the equation q _o = S_ou_o , using their estimate of S_o and a value of q _o given by . The present authors used measured values of u(x) for the lower part of the glacier and then multiplied each value by u_o(L)/u(L) where x = L is the snout of the glacier, in a manner suggested by Reference NyeNye (1963[b]). In this procedure, u_o(L) is determined from the equation

(1)

which holds for a wedge-shaped terminus of wedge-angle θ_o. If α_o is the surface slope at the terminus and β_o is the bedrock slope, then θ_o = α_o — β_o . In our case a_o(L) is fixed, α_o = 11° was taken from the 1968 map and β_o = 2.5° was estimated by first determining the depth of the ice, h, at x = 8 400 m from the equation for a parabolic profile, S = ⅔hB_o , using U-N’s value of S, the cross-sectional area of the glacier. Hence u_o(L) was determined from Equation 1. For our data, U_o(L)/U(L) = 1.0, by coincidence, which means that the snout is currently in equilibrium. To determine the datum velocities for the upper half of the glacier, the values used by U-N were taken and moved up the velocity axis so that they joined up smoothly with the values found above for the lower half of the glacier. The resultant values of u_o(x) are given in Table I.

TABLE I. Datum state values calculated for berendon glacier

With these data, a datum-state glacier was set up and the result is shown in Table 1. This may be compared with table I of U-N (p. 207). The major difference is in the velocity column u_o where our velocities are in general, lower than U-N’s. This causes a corresponding fall in c _o, because c _o = 3B_ou_o (except near the terminus) for a parabolic channel (Reference NyeNye, 1965[a], table IIIc, p. 678).

Comparison of the Predictions

TABLE II. Summary of the results of this paper compared to U-N’S

The theory given by Reference NyeNye (1963[b]) connects the time variation in the mean mass balance a_i(t) , with the time variation in the thickness of the glacier, h_i(t), at the particular value of x through the equation

(2)

where the dots denote differentiation with respect to time, and the λ’s are numerical coefficients that may be computed for the particular x in question. The values of the λ’s at x = 8 550 m (corresponding to the location of the portal to the mine) were computed by the method described by Reference NyeNye (1963[b]) using the functions B_o(x), c_o(x), D_o(x) as input data. The values obtained for λ₀, λ₁, and λ₂ are shown in Table II compared with those found by U -N. The problem is to assess the likelihood that h ₁, at the portal will reach a critical level within a specified time. U-N took 1963 as their reference time, t = 0, and specified times of 20 and 40 years. We took 1968 as our zero time and specified times of 15 and 35 years to compensate for the five-year difference. Moreover, whereas U-N considered that in 1963 the ice had to rise 19 m to overwhelm the portal, we considered a rise of 35 m was necessary from its 1968 level. Our procedure was the same as U-N, namely, to assume a polynomial in t for h₁(t) that is (a) consistent with what is known about h_l in 1968 and (b) will make h₁(t) have the critical value 35 m at a time 35 (or 15) years later. We then use Equation (2) to compute the corresponding function a₁(t), and assess the likelihood of such a variation in the light of the mass-balance data.

Fig. 2. Curves L₁ and L₂: postulated change in ice level near tunnel portal. Ice reaches parlai in 35 and 15 years respectively. Curves B₁, and B₂.: calculated mean mass balances corresponding to curves L₁ and L₂.

From the topographic maps we deduced that the ice thickness in the region under consideration was diminishing by about 4.5 m/year. (U-N used 2.4 m/year.) Thus, taking t = 0 in 1968, we have h₁ (0) = 0 (by definition), h₁ (0) = —4.5 m/year, h₁ (35) = 35 m. The simplest polynomial function h₁(t) consistent with this information is (curve L, in Figure 2)

where o ₁, and p₂ are constants that can be determined. Then

and all higher derivatives are zero. So, substituting this in Equation (2) and using the values of the λ’s and p’s we find

in units of metres and years. This is an almost linear variation in a, from the value +0.24 m/year at t = 0 to +0.99 m/year at t = 35 years (curve B₁ in Figure 2). When the calculation is repeated with the postulate that the ice reaches the critical level in 15 years instead of 35 years, the resulting curves are shown in L₂ and B₂ of Figure 2. Table II summarizes the results of this paper and compares them with U-N’s. It can be seen that the required values that a₁ (t) should attain by 1983 and 2003 to cause the portal to be buried are somewhat greater in the present calculation although we did require the glacier to respond slightly more quickly, and by a somewhat greater amount, than did U-N. However, U-N thought such a deterioration in climate unlikely whereas our data give a mean mass balance for 1967/68 and 1968/69 of +0.72 m/year, and there has been a steady thickening in the accumulation basin since observations were begun in 1963 (W. H. Mathews, private communication in 1969). Stanley (in press) gives a mass balance of +0.78 m/year but this was calculated in a slightly different way and so the two values are considered to be in agreement, If our value of +0.72 m/year persists, we calculate that the ice level will reach the portal in about 1994.Footnote *

Sensitivity of the Predicted Response to the Datum State Chosen

The major difference between the datum states used by UN and the present authors was in the velocity values used, as mentioned earlier. On average, our velocities were 1.5 times smaller than U-N’s although u_o(L) was 2.0 m/year in this calculation against U-N’so.7 m/year. The response of the terminus is directly related to u_o(L) (via c _o(L)) and an error in c _o(L) leads directly to an equal percentage error in the calculated response. But we were not interested in the response of the terminus (x = 8 900) but at x = 8 550 m where the value chosen for c _o(L) is not so important. Figure 3 shows the computed response of the glacier to a unit pulse in a₁(t) and can be compared with U-N’s figure 3a, p. 211. The response at x/L = 0.96 is close to the mine portal position and the main difference between U-N’s value and Figure 3 is the time of the maximum peak which is 60 years in our work and about 35 years in UN. This is almost entirely due to the smaller velocity figure that we used, because the velocity of a kinematic wave is directly proportional to the ice velocity.

Fig. 3. Computed response of Berendon Glacier to a unit pulse, h₁ is the increase of thickness in metres. The mean balance was taken as + 1 m/year from t = 0 to t = 1 year and zero at other times. Numbers against the curves refer to the values of x/L. For points higher up glacier than X/L ≈ 0.15 the initial increment in thickness begins to decay at once. For points lower down the glacier a flood wave builds up and then decays.

In order to see how sensitive the predicted response is to changes in the input data, the programs were run with experimental datum states constructed with the 1968 datum as a basis. In calculating the λ’s the 1968 datum stale data were altered in a number of ways:

1. By putting in U-N’s D _o (x) function, or equivalently using their a_o(x) function.

2. By putting in U-N’s c _o (x) or equivalently using their u_o(x) function, but keeping c _o(L) = 2.0.

3. By putting in values of c _o (x) equal to the 1968 values near x = 0 and x = L but equal to a constant value for all other values of x. The constant value was the average of the 1968 c _o (x) function.

4. By putting in an a_o(x) different from the U-N and the 1968 function. The U-N and 1968 a_o(x) resulted in very nearly the same ice discharge q _o (x), while this trial mass balance function gave a q _o (x) approximately twice that of U-N and the 1968 datum state.

Of all these experimental datum states only (2) gave λ’s significantly different from those of the 1968 datum state and even then the required future values of a ₁ in 2003 were all between 0.5 and 1.1 m year^-1.

The same experimental data were used to calculate the response to a unit pulse. Only (2) gave a response markedly different from that represented in figure 3 for the 1968 datum state. The maxima at X/L — 0.96 all occurred close to t = 60 years except (2) where the maximum was 47 years after the end of the pulse. This was because, as mentioned earlier, U-N’s velocities were on average 1.5 times greater than ours.

From using these different experimental datum states it was concluded that the results were most sensitive to the u_o(x) function, but also it seemed that detailed knowledge of the velocities was not necessary for this type of calculation. Knowing u₀(x) near x = 0, x = L and having an idea of its average value in between seemed to be adequate. Similarly, details of a₀(x) did not appear necessary for these calculations. The results were rather insensitive to variations in a_o(x) and probably a straight-line approximation of a_o(x) would produce all die details reported above. However, a detailed knowledge of a(x) is still needed to calculate the net balance.

It follows, of course, that the main cause for uncertainty in making predictions was the u_o(x) function. The method used to find it from the measured values of a u(x) was indirect and although the values were reasonable, it would have been reassuring to have had some measured ice depths from which we could obtain valley profiles S(x), and obtain u_o(x) from u_o(x) = q _o (x)/S_o(x). making the assumption that S(x) and _o(x) did not greatly differ.

Conclusions

The predicted response of Berendon Glacier calculated using a considerable amount of field measurements is similar to that found by U-N using very little data. For the mine portal to be buried by 2003 A.D., U-N predicts that the net balance must rise to 0.70 m/year, while we predict it has to rise to 0.99 m/year. Our field measurements, however, show that the 1968 mean mass balance is +0.72 m/year, rather than close to zero (for 1963) as was thought by U-N. It is concluded, therefore, that there is a much greater chance of the portal being buried in the mid-1990’s, than was thought by U-N.

The theory is shown to be most sensitive to changes in the velocity figures used, u_o(x), although a detailed knowledge of u_o(x) is probably unnecessary. Similarly a detailed knowledge of mass balance is not necessary to set up the datum state, although the mean mass balance is needed to compare with the final prediction.

Acknowledgements

We would like to thank the many people who collected the necessary field data, without which this paper could not have been written. In particular we thank Mr R. J. Rogerson and Dr A. D. Stanley for allowing us to use their unpublished results. We are grateful to Professor J. F. Nye, University of Bristol, England, for his comments on an early draft of this manuscript, and for kindly supplying a copy of his computer programs.

MS. received 21 April 1970 and in revised form 17 July 1970