I. Introduction

This is a review of the ray and wave theory required to infer the form of glacier beds B, and to measure glacier motion relative to B, by analysing observations of radio echoes from B. The radio source-receiver P is situated at a point R ₀ in a horizontal plane above B (Fig, 1), and emits quasimonochromatic pulses with predominant wavelength λ (typically 5 m in ice) and spatial length σ (typically 50 m in ice), which travel through ice with speed c. The echo received back at P usually consists of several broadened pulses, arising from different parts of B, and also from features in the ice between B and P. Here we ignore the reflections and scattering in the ice, and also refraction at the ice-air interface; these effects can be taken into account as explained by Robin and others (1969) and Reference HarrisonHarrison (1970, Reference Harrison1971)· We do, however, consider absorption in the ice, and the effects of the finite angular beamwidth of the source and receiver.

Fig. 1. Notation for source-receiver, pulse and glacier bed (showing two specular points).

The analysis falls naturally into three parts, corresponding to increasing fineness of distance resolution. First, in section 2, we discuss briefly the use of geometrical optics in deducing the geography of B; this is the topography averaged over distances of the order of the average depth h of B below the plane in which P lies (typically h is several kilometres). Then in section 3 we discuss the use of diffraction theory in deducing the roughness of B; this is the term used for the statistical properties of the topography, considered on scales down to about λ/10. Finally in section 4 we discuss the use of diffraction theory to understand the detection and measurement of displacements of P relative to B on a scale which may be as small as λ/500.

The emphasis throughout will be on describing and explaining the various methods and formulae, rather than on deriving these from first principles (derivations can usually be found in the papers cited).

2. Inferring Bed Geography From Specular Echoes

Specular echoes are those from the “specular points” on B, that is, from the points where the normal to B passes through the source-receiver position P (Fig. 1). For P fixed at R ₀, let the specular echoes arrive with delay times measured to be t_i(R₀ ), where i = 1, 2,..., N(R₀ ), N, being the total number of specular echoes received at P. Construct spheres centred on P with radii ct_i (R ₀/2) (Fig. 1). Then we know that B must touch each of these spheres at least once, but we do not know where each point of contact is. To resolve the ambiguity, let P wander over the plane R ₀ (so that the specular points, and their ranges, change) and repeat the sphere construction at each point. The glacier bed B is then the (lower) envelope of all these spheres, because B must touch every sphere.

We emphasize that this simple method yields a unique reconstruction of the “accessible” parts of B -that is, those parts of B for which the outward normal does not intersect B again — only if P makes a full two-dimensional traverse of R ₀. Usually this is not done: P traverses just a line in R ₀ (e.g. the track of an aeroplane carrying the radio-echo sounder). In this case the envelope of spheres is a “worm” with rotational symmetry about the locus of P and with varying radius (Fig. 2), and all we can infer with certainty is that B touches the worm in some unknown line. In practice, this ambiguity is often resolved by assuming that B varies effectively one-dimensionally, i.e. that the bed is a series of ridges normal to the traverse direction; then the cross section of B normal to the ridges is simply the lower line in which the vertical plane through P intersects the worm (in Fig. 2 this is the line ℒ). Another technique occasionally used for roughly locating the line where B touches the worm is to oscillate the aerial at P and use the fact that the beam is somewhat dircctional (private communication from J. F. Nye and M. E. R. Walford).

Fig. 2. “ Worm” envelope of spheres corresponding to ranges of specular echoes received at each point as P traverses a line.

For the sphere construction to work, it is essential to be able to identify the specular reflections t_i . To study this question it is necessary first to state what is meant by roughness and smoothness; we define B as rough (or smooth) on a given horizontal scale x if B possesses (or does not possess) irregularities whose average horizontal dimension is x. Identifying t_i is easy if B is smooth on all scales smaller than h, because then geometrical optics is applicable (i.e. the reflection from B is wholly specular) and, moreover, the number N (R ₀) of specular points is small (indeed N is unity if all centres of curvature of concave parts of B lie above the plane R ₀); thus the arrival times t_i of different reflections are usually separated by much more than the pulse duration σ/c, and the reflections are easily identified. The situation is more complicated if B is smooth on the scale of λ but rough on scales between λ and h, which is always the case if h is large enough. Then many specular reflections arrive at P; these are separated by times short in comparison with σ/c and so are difficult to identify from the echo, which is their superposition. A possible deconvolution procedure for identifying the t_i , in these circumstances has been suggested by Reference BerryBerry (1972), but is likely to be difficult to apply in practice. The simplest method that can be applied in this case, and also in the case where B is rough on the scale of λ—so that most of the scattering is non-specular—is to identify delay times when the echo is large (especially the first returns from B), and use the fact that t_i are continuous functions of R ₀. The times t_i thus obtained correspond to the few specular echoes that would be received if the actual rough surface B were smoothed on all scales smaller than h. This method works quite well in practice (see the paper by Harrison (1970), which also shows how to include the effects of refraction from air into ice).

3. Inferring Bed Roughness from the Echo Envelope

The roughness is the statistics of the small-scale departures of B from its “geography”, this being defined as the surface B_h smooth on scales smaller than h which is determined by the methods of the last section. If points on B are specified by their height B(R) above the point R (= x, y) on horizontal reference plane h below the R ₀ plane containing P, then the departure of B from its geography B_h (R) may be specified by a height function f (R) = B(R)—B_h (R) (Fig. 3), and regarded as a perturbation of the plane R.

Fig. 3. Definition of height function (R) as departure of bed B its geography B_N, and polar diagram of source-receiver P.

Statistics of f (R) are defined by averages over R, these being denoted by brackets < >. By definition the mean height <f(R)> is zero. The r.m.s. height S is defined by

the autocorrelation function C(R_I ) between two points separated by R_I on the surface is defined as

and does not depend on R. Simple forms for C(R_I ) are

both of which correspond to isotropic roughness with “correlation length” L; L and S are the average horizontal and vertical dimensions of the irregularities of B (Fig. 3). (We are assuming here that there is a single horizontal scale L for the roughness. If there is not— that is, if the surface is irregular on all scales (see Nye, 1970, p. 385 ff )-then the separation into roughness and geography seems artificial. It is not clear how to deal with this problem, but the short duration (relative to h/c) of most observed echoes suggests that actual beds might not be rough on all scales.)

Qualitatively different surfaces may have the same roughness statistics S and L; it is difficult to give an exhaustive description of all types of surface, so Reference BerryBerry (1973) introduces two extreme cases. The first (“surface α”), is .smoothly undulating, with “Gaussian” statistics, and C(R_I ) may take the first form in Equation (3) but not the second. The other case (“surface β”), consists of steps, that is vertical walls separating flat facets and C(R_I ) may take the second form in Equation (3) but not the first. Surface a might correspond to an undulating rock bed or to the bottom of rugged floating sea ice, while surface ß might correspond to the ice-water interface beneath flat ice-floes.

We shall also consider a third type of bed (“surface γ”) where the roughness is not due to departures from geography but arises from variations in the (complex) reflectivity Z(R); over a smooth surface. Thus f (R) = o and the roughness is the statistics of Z(R); the auto-correlation function of Z(R) may be modelled using either of the forms in Equation (3). Surface γ might correspond to a smooth bed made up of varying rock types.

We wish to use echo data to infer S, L (and possibly also C(R_I )) and distinguish between roughness types α, ß and γ. We represent the echo signal received with P at R ₀ by the voltage ψ (τR ₀) at the aerial terminals; τ is time delay, measured from the instant that a specified part (say the centre) of a specularly reflected pulse would arrive at R ₀ from the smooth “geographical” surface. It is usually convenient to smooth away the echo oscillations with the frequency of the carrier wave ω₀ (= 2πc/λ), and record the more slowly-varying echo envelope p, defined as

In practice p is measured using the second expression, that is by rectifying and smoothing, but the first expression is more useful for theoretical purposes. The approximation in Equation (4) introduces a fractional error of at most λ/2σ, which is small for the quasimonochromatic pulses commonly used. The averages <ρ(τ)> or <ρ²(τ)> over R ₀ give the mean shape of the echo, from which the bed roughness is to be determined.

To relate the average echo shape <P²(τ)> to the bed roughness, we neglect, polarization effects and use the scalar-wave Kirchhoff diffraction theory, together with the Fresnel approximation, which consists of replacing cos θ by 1 —½ θ² for all incident or scattered waves making angles θ with the vertical. The Fresnel approximation greatly simplifies the mathematics, which would otherwise be intractable. The neglect of large-angle scattering (beyond θ ≈ 45°) is justified, because

• (1) the angular beam width of the source and receiver while often large, is finite;

• (2) absorption in the ice reduces the intensity of waves scattered through large angles, because of the greater path lengths involved;

• (3) back-scattering of large-angle radiation is insignificant unless the surface slopes are precipitous;

• (4) Kirchhoff theory neglects shadowing and multiple scattering, which inevitably occur at large angles.

The emitted radio pulse is specified by the voltage a(t) cos w₀t at the terminals of the transmitting aerial, where the duration of the smoothly-varying pulse envelope a(t) is σ/c, which is assumed large in comparison with 2π/w₀ (“quasimonochromaticity”). Absorption in the ice is specified by the attenuation distance d at which the intensity of a plane wave is reduced by a factor e; it is the large value of d (typically ≈ 1 km) that makes radio-glaciology possible. Finally, the finite angular beam width of the source-receiver (Fig. 3) is modelled by the “polar diagram” intensity function exp (—tan² θ/tan² θ₀). This form, circularly symmetric about the vertical, is chosen for mathematical convenience; beams used in practice are usually rather broader along the direction of motion of P than across it, but this effect, and the fact that the beam may have side lobes, rarely affects the interpretation of roughness scattering because θ₀ is often large—say 25° to 50° (M. E. R. Walford, in a private communication, states that the shapes of beams used in practice are rarely, if ever, precisely known).

An elaborate analysis (Berry, 1973), after a slight generalization to include beam-width and absorption, then shows that the mean echo shape from surfaces α, β or γ can be written in the form

where g_t is the transmitting aerial power gain; S_r is the receiving aerial absorption cross section; δ is a coherence parameter, defined as

Q,(T) is an echo-tail function, defined as

where for surface γ we have written the formula corresponding to the first form forC(R_I) in Equation (3); finally,. ℐ is a decay time that includes the effects of absorption and beam-width, and is defined as

The first and second terms in the ( } in Equation (5) describe respectively the coherent and incoherent contributions to the mean echo envelope, the relative strengths being determined by δ. When δ ≈ 1 the echo is mostly coherent, and has the same shape a²(τ) as the original pulse intensity, with duration σ/c. If δ ≪ I the echo is mostly incoherent, and is characterized by a tail for long delays τ; this tail arises from scattering by roughness not directly below P—the roughness on the ring |R—R ₀| = (chr’) ½ gives the contribution at τ’ to the integral in Equation (5). Very slight irregularities produce incoherence: Equation (6) shows that δ is reduced to 0.2 if S ≈ λ/10.

A convenient measure of the length of the echo is the half-power time T, defined as the time (measured from τ = o) when half the echo energy has arrived, that is the value of τ bisecting the area under the graph of <ρ²(τ)>. Berry (1973) shows that Tis given by

(Remembering the definition of the instant τ = o, we see that the time over which half the energy may be considered to arrive is actually T+ σ/2c.) Equation (9) ignores the effects of beam-width and absorption; this is justified with the wide beams commonly used, because echo tails are rarely observed to extend with significant strength out to, say, τ ≈ h/c. However, if the surface were so rough that absorption or beam-width dominated the decay of <p² > at large τ, then Equation (9) would have to be replaced by

Reference OswaldOswald (1975) uses the time delay at which the echo maximum arrives, rather than the half-power time, as a measure of echo duration.

To infer the surface statistics it is sufficient to determine δ and T, and use Equations (6) and (9). For surfaces α or β this enables S and L to be found, while for surface γ this fixes |<Z>/²|<|Z|>² and L. Any uncertainty as to which type of surface (α, β or γ) lies beneath the ice can in principle be eliminated by varying the carrier frequency ω₀, keeping σ constant: if T does not change, the surface is type α; if δ does not change, the surface is type γ; if both T and δ change, the surface is type β. Direct measurement of T is not difficult unless the echo is almost coherent (δ ≈ I) (in which case all we can say is that B is smooth between scales h and λ). To measure δ it is necessary to split <p² > up into its two terms (Equation (5)), and use the fact that the ratio of the areas under their graphs is δ^-1-1. Carrying out the splitting would not be easy (especially in the “incoherent” case δ ≈): probably the best procedure is to vary σ this would leave the incoherent term largely unaffected, but would change the coherent term, which could thereby be identified. To facilitate the splitting, the instant τ = o could be precisely located by measuring the spatial fading rate (see section 4). (Another method for determining δ is described by Oswald (1975); he measures the variance of p² as it fluctuates about its mean value <p² >.)

On existing records there is no possibility of varying ω₀ and σ, and the best that can be done is to measure T, assume an undulating surface α, and infer the r.m.s. slope S/L of the roughness from Equation (9). The surface type α. is perhaps most likely to occur in practice; types β and γ, moreover, would probably not occur as pure cases. Therefore any further sophistication of the admittedly crude deconvolution procedures suggested here is probably unwarranted at the present time (an example of such sophistication would be an attempt to deduce the surface autocorrelation function—which need not be given by Equation (3)— from measurements on the echo tail, as suggested by Berry (1973); a knowledge of the form of C(R_I) would be useful in testing theories of glacier sliding, see e.g. Nye (1970)).

4. Inferring Ice Motion From Echo Fading

As the position R ₀ of P is changed, the echo wave function and envelope for any given τ fluctuate about their mean values. For small displacements the fluctuation is quasi periodic, and is called spatial Jading. This arises by scattering from the roughness of B in a manner first explained by Robin and others (1969): the echo at τ comes from irregularities lying in a certain annulus on B, whose centre is R = R ₀; the new annulus about a slightly displaced position R ₀+ΔR ₀ contains essentially the same irregularities, but the phases of their contributions at P will now be different. A fading “period” corresponds to a distance Δ R ₀; for which the average phase of the contributions has changed by 2π. The fading rate is most easily measured as the average number N_P(τ) of times that the envelope p passes through a given datum level pD as P moves along unit length of a straight line. If pD is very large or very small, N_p is small. For some intermediate value (actually pD = p _rms/√2),N _p is a maximum, and a formal analysis (Berry, 1973) shows that in this case

where τ( τ) is a modified delay time defined for the case of incoherent echoes by

this function is plotted in Figure 4 for the Gaussian pulse envelope

a(t) = Λεχρ (—cWja¹).

We see that τ ≈ τ whenever τ > σ/c. Thus the fading is fastest in the tail, and slowest near the first return. The formula for N_P(τ) has been confirmed experimentally by Walford (1972).

Fig. 4. Modified delay time τ(τ) required to calculate the spatial fading rate, N_p .

Equation (12) is valid whenever the echo is incoherent, that is, always in the tail (τ > σ/2c), and for all τ if δ is small. We emphasize thai the formula for N_p(τ) does not involve the roughness parameters S and L ; the only requirement for this type of spatial fading is that a large number of independently scattered waves reach P from B, and the condition for this is

which is not very restrictive (if h = 1 km and σ = 50 m it gives h < 400 m). We conclude that no information about bed roughness can be obtained from the spatial fading of incoherent echoes. However, measurements of N_p could be employed in conjunction with Equation (11) for precise measurements of time delay; this might be useful in view of the fact that it is not easy simply by inspection of echo curves to identify the instant τ = 0.

In the case where the echo is partially coherent, Equation (11) is still valid, but Equation (12) for the modified delay time τ(τ) is no longer applicable. When the surface is known (or assumed) to be of type α, the value of τ( τ) for any degree of echo coherence can be obtained from Figure 4 if the labelling of the abscissa is changed by the replacement

It follows from this that for very gentle slopes and/or narrow-angle beams the echo fades much more slowly than in the incoherent case; this is the regime investigated by Oswald (1975)·

If P is fixed with respect to the upper surface of a glacier, then R ₀ will change if the ice flows. Spatial fading of the echo is a very sensitive indicator of this motion. It was suggested by Nye and others (1972[b]) that it might be possible to measure the motion of a glacier in a few days, and this was confirmed in a field experiment in Antarctica by Walford (1972), who detected a motion of 0.38±0.03 m d ^-1, which corresponds to an accuracy of better than λ/100. To achieve maximum sensitivity with this method two conflicting requirements must be balanced: the fastest fading rates occur in the tail of the echo, but the most sensitive measurements can be made on the main body of the echo where the signal is strongest. Let us suppose that the maximum of the echo envelope as a function of τ is p_m , and that the instrumentation is such that it is possible to detect a change Δp in the echo envelope. Then Berry (1973) shows that the smallest detectable displacement |ΔR ₀|_min corresponds to choosing a delay time τ ≈ 4T/3, where T is the half-power time defined in section 3 (of Equations (9) and (10)). If we choose a place R ₀ where p is fluctuating n times as fast as its r.m.s. rate of change, then |ΔR ₀|_min for the case of an incoherent echo is given by

As a numerical example, suppose cT is observed to be h/4, n = 2 and Δp/p_m = 0.01; then Equation (15) gives |ΔR ₀|_min = λ/500, so that an ice displacement of 1 cm would be delectable in this case. This method depends, of course, on choosing R₀ so that the echo comes from a fixed rock bed and not, for example, from a moving moraine.

Spatial fading of the echo envelope then reveals very small horizontal displacements Δ R ₀ of P, It does not, however, reveal small vertical displacements Δ h, because the minimum detectable Δh is proportional to λ h/σ, and this is generally larger than | Δ R ₀|_min (Equation (16)) by one or two orders of magnitude. In order to measure Δ h (and so begin to discover whether, for example, the Antarctic ice sheet is getting thinner or thicker) it is necessary to study the phase χ of the echo rather than its envelope. The phase fades horizontally as well as vertically, and Nye and others (1972[a]) suggest using the fading of p to recover horizontal position and then using χ to detect vertical movement. It is possible that in the future a new level of precision of phase-sensitive techniques may be achieved, based on the recently discovered “wave-front dislocations” (Reference Nye and BerryNye and Berry, 1974).

5. Discussion

We see that in principle radio-echo sounding can be used to ascertain the geography and roughness of subglacial terrain, and to monitor horizontal and vertical motion of the source-receiver relative to a frame of reference fixed with respect to the rock bed. In practice the methods of section 2 for reconstructing geography have been used to study essentially one-dimensional relief (Harrison, 1970), and the methods of section 4 for measuring horizontal motion have been successfully tested by Walford (1972) and C. S. M. Doake (private communication). The Bristol group propose shortly to begin a field test of the method suggested for measuring vertical ice motion. There is a need for more experiments, both in the field and in the laboratory (using ultrasonic simulations of radio-echo sounding pioneered by Nye and others (1972[b])), to test the various techniques whose principles are outlined here, especially the proposed methods for measuring roughness (section 3).

Concerning the theory, we suggest three lines of development. First, a careful study should be made of the separation of bed relief into geography and roughness, to discover how far the methods suggested here apply when the rock surface is rough on all scales between λ and h. Second, the statistical methods of sections 3 and 4 should be refined to include the case where the bed roughness is strongly anisotropic. Third, the information about the bed contained in the polarization of the echo should be investigated.