The above introduction simply summarizes the relevant parts of Nye [III] without adding anything new.

valid for high frequencies. A scheme for computing the *λ*’s and *ν*’s from known functions *B*
_{0}(*x*), **c**
_{0}(*x*), **D**
_{0}(*x*) was developed and applied to data derived from Dr. Mark F. Meier’s work on South Cascade Glacier, Washington, U.S.A. In this way the frequency response of South Cascade Glacier was computed. The work was done on a desk calculator.

When automatic digital computers became available to the author two possibilities presented themselves. First, to write a programme for the calculation of *μ* and *λ* coefficients previously done by hand, which could then be used for data on any glacier. The results from this programme will be described later. Second, instead of obtaining the frequency response in two stages, by first computing the *λ*’s and *v*’s and then using the two series (7) and (8), it should be possible to compute solutions of (3) and (4) directly for a number of different *w*g’s, thereby finding the response curves. This would be a more direct way of finding the frequency response, and at the same time would avoid any difficulties in the intermediate range of frequencies where both (7) and (8) might break down.

### (i) Results for analytical model

As a start I investigated solutions of (3) and (4) with *B*
_{0}, **c**
_{0}, **D**
_{0} having the special forms (5). *σ* is taken as the unit of time and *l* as the unit of length. *A* is also put equal to 1 without loss of generality since both *Q* and *H* are directly proportional to *A*. The equations are then

with

and with the datum glacier running from *x* = 0 to *L*, where*L* = 1−*δ*. Thus 0 ≤ *x* ≤ *L* < 1.

The boundary conditions at *x* = 0 and *x* = *L* need special consideration. **D**
_{0} = 0 at both end points, which are regular singularities of the equations. The problem is fully discussed in Appendix I (see also Nye [II]) and here we need only quote certain results. A unique solution for any given w is determined if we simply require that H be not infinite at *x* = 0 or *x* = *L* The requirement at *x* = 0 is equivalent to putting *Q* = 0 at this point. Near *x* = 0 there is a one-parameter family of solutions satisfying the boundary condition, in the sense that a particular solution is fixed by the value of one complex constant; but the leading term is the same for all members of the family, namely

where the dots stand for terms of higher order in *x*.

Near *x* = *L* there is also a one-complex-parameter family of solutions satisfying the boundary condition. The leading term in this family of solutions, when they are expanded about the point *x* = *L*, is an arbitrary complex constant. For all this family of solutions, at *x* = *L*,

and

where the primes denote differentiation with respect to *x*. (14) is derived by putting **D**
_{0}
*H*′ = 0 in (10) (more strictly in (14)); the second may be obtained by differentiating (10), putting **D**
_{0}
*H*″ = 0, and substituting for *Q*′ from (9). The jusgification for putting **D**
_{0}
*H*′ = 0 and **D**
_{0}
*H*″ = 0 at *x* = *L* is that **D**
_{0}(*L*) = 0 and, as shown in Appendix 1, *H*′(*L*) and *H*″(*L*) are not infinite in the family of solutions satisfying the boundary condition.

The method of solution is to integrate the simultaneous equations backwards from *x* = *L*, using (9) and (10) at all points except *x* = *L*; at *x* = *L*
(10) does not suffice to fix *H*′ because **D**
_{0} is zero; we therefore use (15) instead. The starting values of H and 0 are not known, and must be chosen by trial so that the solution satisfies the boundary condition at *x* = 0, that is, so that *Q* = 0 or H is not infinite. To avoid the trial solutions becoming infinite the testing is done not at *x* = 0 but at *x* = *ϵ*, *ϵ* being small. The value we have to achieve is

from (13). The trials are made as follows. Write *H* = *u*+*iv* and *Q* = *r*+*is* where *u*, *v*, *r*, *s* are real. Choose arbitrary starting values for *u*, *v*, say (0, 0), at *x* = *L*. Find *r* , *s* from (14). Then integrate back to *x* = *ϵ*. Let the resulting values of *r*, *s* at *x* = *ϵ* be *r*
_{1}, *s*
_{1}. Now choose a new pair of starting values for *u*, *v*, say (1, 0), and repeat the integration, obtaining at *x* = *ϵ* the values *r*
_{2}, *s*
_{2}. Since *Q*(*ϵ*) depends linearly on the starting value *H*(*L*) we have, in general,

or, taking real and imaginary parts,

and

where *j* = *j*
_{1}+*ij*
_{2} and *k* = *k*
_{1}+*ik*
_{2} are complex constants. *j*
_{1}, *j*
_{2}, *k*
_{1}, *k*
_{2} are fixed by the definitions of *r*
_{1}, *s*
_{1} and *r*
_{2}, *s*
_{2} so that we have

Thus the constants in the linear relation are determined by the known values *r*
_{1}, *r*
_{2}, *s*
_{1}, *s*
_{2}. If *r*(*ϵ*) and *s*(*ϵ*) are now assigned the values required by (16), (17) and (18) may be solved to give the correct starting values *u*(*L*), *v*(*L*). A final integration is now performed using these starting values, and it should strike the target values *r*(*ϵ*), *s*(*ϵ*) required by (16). This is the required solution.

This method* was programmed for an I.B.M. 709 digital computer, using a Runge-Kutta integration method, with *ω*, *E*, *δ* and *ϵ* as variable parameters. It worked successfully for values of *ω* from 0 to about 30 (with *δ* = 0.01), and was checked against the known analytical solution for *E* = 1. The precision was checked by changing *ϵ* and the interval. At the highest frequencies *ω* ≥ 18 the fit with the boundary condition at *x* = *ϵ* deteriorated. This is because, at these frequencies, a very small change in the starting values at *x* = *L* has a large effect at *x* = *ϵ* (or, put the other way round, a large input *Q* at *x* = 0 would be needed to produce a sensible effect at *x* = *L*, because of attenuation by diffusion); thus for a reasonable fit at *x* = *ϵ* very many significant figures would be needed in the starting values at *x* = *L*. But this effect need not trouble us, because the main purpose was to find the response at *x* = *L*, and this is still determined quite precisely even though the fit at *x* = *ϵ* is bad. At still higher frequencies (say *ω* > 50), where overflow may occur in the trial integrations, the result is best obtained by using the high-frequency series (8), of which the leading term is simply *H*(*x*)/*A* = (*iω*)^{−1}.

The results at *x* = *L* for *E* = 0.1 and 1 are shown in Figures 1a, b, c, d. (Results for *E* = 2 did not differ much from those for *E* = 1.) If we look first at the curve for *E* = 1 in Figure 1a (already known from Nye [III), the amplitude |*H*| of the response to unit amplitude *A* approaches the value *δ* = 100 at low frequencies, and approaches zero as *ω*
^{−1} at high frequencies. On the plot of log |*H*|: log *ω* in Figure 1b the curve for *E* = 1 is below the straight line |*H*|−*ω*
^{−1} (broken line) at low frequencies, above it at medium frequencies, and approaches it asymptotically at high frequencies. This may be explained by saying that at high frequencies the glacier acts as a perfect integrator of the changes in accumulation rate, that is ∂_{
h1}/∂*t* = *a*
_{1}, giving the response |*H*| ~ *ω*
^{−1}. At medium frequencies the response at the terminus is greater than this because it is reinforced by kinematic waves from higher up the glacier. But at low frequencies the glacier no longer integrates, because the leakage rate is too fast, and we have *h*
_{1} simply proportional to *a*
_{1}

Fig. 1a. Frequency response at the terminus of an artificial glacier. The response amplitude |H| to a unit amplitude variation of accumulation and ablation rate is plotted against angular frequency ω. ω is measured in natural units (σ = 1). For illustration a scale of periods, taking σ = 6 yr., is shown at the top

Fig. 1b. Same as Figure 1a, but showing |H| on a logarithmic scale

Fig. 1c. Same as Figure 1a, but showing the phase lag ϕ of H on A

Fig. 1d. Frequency response at the terminus of an artificial glacier, with E = 0.1. A polar plot of the end of the vector H as ω changes. Numbers against the curve are values of ω in natural units (σ = 1)

The curves of |*H*| for *E* = 0.1 in Figures 1a, b are similar in general form to those for *E* = 1, but they cross at *ω* = 1.2 (owing presumably to more destructive interference between the waves).

The phase lag *ϕ* of the response at *x* = *L*, defined by

is plotted against log *ω* in Figure 1c. The curve for *E* = 1 goes through a maximum (at

) as found in Nye [II]. The curve for

*E* = 0.1 is similar but more sharply peaked.

Figure 1d is a polar plot for

*E* = 0 of the end of the vector

*H* as it changes with

*ω*.

It seemed desirable for completeness to calculate the response for *E* = 0, which corresponds to no diffusion and therefore pure kinematic wave propagation down the glacier. The differential equations become of first order and hyperbolic, instead of second order and parabolic. The numerical method used for non-zero *E* fails for this case. It is shown in Appendix II that, if *B*
_{0}(*x*) ≡ 1 and **D**
_{0}(*x*) ≡ 0, the response at *x* = *L* may be expressed quite generally as

where

*T*(*ξ*) is the travel time for a wave from *x* = *ξ* to *x* = *L*; in the case under consideration where **c**
_{o} is given by (11)

The integral in (19) was evaluated numerically. At small *ξ*, *ωT*(*ξ*) in the integrand is large and rapidly changing, so the integral was approximated in the range *ξ* = 0 to *ϵ* by the first two terms in the power series in *ϵ*. Over the rest of the range *ξ* = *ϵ* to *L*, Simpson’s Rule was used. With *ϵ* between 10^{−5} and 10^{−2}, and with the remainder of the range divided into 400 intervals, this procedure was satisfactory for values of *ω* up to about 2. The accuracy was checked by changing *ϵ* and the interval.

For *ω* > 2 a very small interval would have been needed near *ξ* = *ϵ* to cope with the rapid change of *ωT*. It seemed preferable to deal with this range of *ω* by direct integration of the original differential equations. At the same time this gave a check on the previous results. Equation (4) simply reduces to *Q* = **c**
_{0}(*x*)*H*, with **c**
_{0}(*x*) = *x*(1−*x*), and equation (3) therefore becomes

The boundary condition is *Q* = 0 at *x* = 0. At *x* = 0, which is a singularity, (21) does not suffice to determine *Q*′ directly. In the neighbourhood of *x* = 0 the equation is sufficiently well represented by

which has the general solution

where *P* is an arbitrary constant. *P* = 0 by the boundary condition, and hence

(21) and (22) were integrated from *x* = 0 to *L* (*L* = 1−*δ*, *δ* = 0.01) by a Runge-Kutta method. The previous results from (19) were checked and the range of *ω* extended up to *ω* = 5.6 (where the ranges 0 ≤ *x* ≤ 0.891 and 0.891 ≤ *x* ≤ 0.99 were divided into 540 and 600 intervals respectively).

The results for *E* = 0 are seen in Figures 1a, b, c, e, f. Figure 1a shows |*H*|:log *ω*. The drop from maximum response takes place in a rather narrower frequency band when there is no diffusion. In Figure 1b, which shows log IHI: log w, the response amplitude is seen to oscillate before settling down to fall off as *ω*
^{−1}. The phase lag *ϕ*, is no longer confined within a single 360° range. This is seen in the polar plot of Figures 1e, f; as *w* rises from 0 to 2.95 the vector *H* rotates clockwise through

revolutions; it then rotates a further 22° before coming in with minor oscillations to

*ϕ* = 90°. On the plot in

Figure 1c
*ϕ*. is taken to be between 0° and 360°, and, for clarity, only the two end branches of the curve for

*E* = 0 are shown. Appendix II shows that a good approximation for

*ϕ*. except at high frequencies is

for *δ* = 0.01. This gives values of *ϕ* correct to about 1 per cent up to *ω* = 2.

A polar plot of *H*, such as Figures 1d, e, f, must start at *ϕ* = 0 and finish at *ϕ* = 90°. For *E* = 0.1 it does not encircle the origin (Figure 1d), but for *E* = 0, as we have seen, it encircles the origin twice. This poses a continuity problem. As *E* changes continuously between 0 and 0.1 the curve must change the number of times it encircles the origin first from 2 to 1 and then from 1 to 0. How it can do this in a continuous fashion is an interesting question which has not been studied in detail. Probably, for certain small values of *E*, the curve actually traverses the origin for finite *ω*. This would mean that *H* = 0 for certain finite *w*—with an infinite downward oscillation on a plot of log |*H*|: *ω*. The terminus would be a node for these frequencies.