Combining SAR interferometry and the equation of continuity to estimate the three-dimensional glacier surface-velocity vector

Abstract

Until now, an assumption of surface-parallel glacier flow has been used to express the vertical velocity component in terms of the horizontal velocity vector, permitting all three velocity components to be determined from synthetic aperture radar interferometry. We discuss this assumption, which neglects the influence of the local mass balance and a possible contribution to the vertical velocity arising if the glacier is not in steady state. We find that the mass-balance contribution to the vertical surface velocity is not always negligible as compared to the surface-slope contribution. Moreover, the vertical velocity contribution arising if the ice sheet is not in steady state can be significant. We apply the principle of mass conservation to derive an equation relating the vertical surface velocity to the horizontal velocity vector. This equation, valid for both steady-state and non-steady-state conditions, depends on the ice-thickness distribution. Replacing the surface-parallel-flow assumption with a correct relationship between the surface velocity components requires knowledge of additional quantities such as surface mass balance or ice thickness.

Introduction

Until the recent implementation of satellite interferometric synthetic aperture radar (InSAR) techniques, the density of surface velocity observations on most glaciers was low because of high manpower and logistics costs related to the traditional methods of measuring glacier velocities. Although repeated differential global positioning system (GPS) positioning has revolutionised in situ velocity observations on remote, large ice masses such as the Greenland and Antarctic ice sheets, high cost is still an obstacle to applying this technique for a general regional-scale mapping of glacier surface velocities.

The use of InSAR techniques to estimate ice motion and surface elevation has developed rapidly in recent years. Using digital elevation models, slant-range velocity has been separated from topographic effects (Joughin and others, 1995; Rignot and others, 1995), and using multiple interferograms, surface elevation and one slant-range velocity have been determined simultaneously (Joughin and others, 1996a, b; Kwok and Fahnestock, 1996: Rignot, 1996). One velocity component is sufficient only if the direction of flow can be estimated by other means, such as in the interior of an ice sheet where ice flow is in the direction of maximum surface slope. Usually, however, knowledge of the full three-dimensional (3-D) flow pattern is essential for validation of ice-flow models. For mass-balance assessments by flux-divergence calculations (Reeh and Gundestrup, 1985), at least the two-dimensional surface velocity pattern must be known.

It has been demonstrated that the 3-D surface flow pattern and surface elevation can be obtained by combining InSAR measurements from descending and ascending satellite tracks having different look directions (Joughin and others, 1998; Mohr and others, 1998). Three look directions, though, are required in order to directly measure the full 3-D velocity vector. The satellite systems presently available have used only two look directions. The InSAR measurements must therefore be combined with additional ice-flow information. Joughin and others (1998) and Mohr and others (1998), for example, used the assumption of surface-parallel flow to express the vertical velocity component in terms of the horizontal velocity vector. This assumption neglects the contribution to the vertical velocity of the local mass balance as well as a contribution arising if the ice sheet is not in steady state.

One purpose of the present paper is to discuss this assumption. We find that the mass-balance contribution to the vertical surface velocity is not always negligible as compared to the surface-slope contribution. It may even be the dominant term over large areas of the ice-sheet surface. Moreover, the vertical velocity contribution arising if the ice sheet is not in steady state can be significant. If tie points of known location and velocity are used to generate InSAR velocities, and the vertical velocity component of the tie-point information is ignored, then an inversion which assumes a zero vertical velocity would, of course, provide the correct horizontal velocities at the tie points. However, in areas away from the tie points where the vertical velocities are different, increasing horizontal errors would be observed. In Greenland there are many tie points on the coast, particularly on exposed rock, where the velocity can be assumed to be zero. When using such tie points, the compensation procedure described above would not account for the vertical motion of the adjacent ice, and erroneous ice velocity vectors would be the result In the interior of Greenland, tie points include the GPS measurements performed along the 2000 m elevation contour as part of the Program for Arctic Regional Climate Assessment (PARCA) project (Thomas and others, 1997) and measurements in the Summit area (Hvidberg and others, 1997). Here, the 3-D position and 3-D velocity of the tie point are used in an "absolute sense" to determine the InSAR measurement parameters (e.g. satellite baselines and the absolute phase of the interferogram). It is important to stress that the measured vertical velocity is used in determining the InSAR parameters. If the InSAR data were inverted without considering the correct vertical velocity, the horizontal velocity of the InSAR data would not be correct, even at the tie points.

The principles behind the use of InSAR to derive glacier surface elevation and motion have been described in several papers (e.g. Joughin and others, 1996a). Currently, key limitations of the InSAR technique include atmospheric perturbations, differential snow layers (a snow layer in only one image) and radar-instrument limitations including clock drift, etc. All of these errors can propagate to erroneous baseline parameters if tie points are used to determine these parameters. However, all of these effects are random in character and will tend to average out when larger areas and many datasets are analyzed. In particular, it should be noted that the basic noise limitations of InSAR (e.g. 10 m elevation and 3 m a ^-1 to 5 m a ^-1 velocity accuracy) are reduced by spatial averaging. However, the issue discussed in this paper is systematic in nature and will not average out.

We apply the principle of mass conservation to derive an equation relating the vertical surface velocity to the horizontal velocity vector. This equation, valid for both steady-stale and non-steady-state conditions, depends on the ice thickness distribution. Thus, replacing the surface-parallel-flow assumption with a correct relationship between the surface velocity components requires knowledge of additional quantities such as surface mass balance or ice thickness.

Surface Boundary Condition

The boundary condition at the surface of an ice sheet relating the rate of change of surface elevation ∂S/∂t to the vertical ice-particle (pole) velocity w_s, the horizontal velocity vector the surface slope and the specific mass balance b_s is (see Fig. 1)

(1)

Equation (1) is the local kinematic boundary condition which is valid at any instant, b_s is the instantaneous specific surface balance (snow accumulation or snow/ice melt) measured in m snow/ice per unit time.

Fig. 1. Kinematic boundary condition at the surface.

In order to avoid problems with decorrelation, synthetic aperture radar (SAR) images for interferometry must be acquired from periods of low, preferably zero, accumulation/ablation, i.e. b_s = 0.

Consequently, under preferred InSAR conditions, the correct surface boundary condition is

(2)

The assumption that flow is parallel to the surface is expressed by the equation

.

Equation (2) shows that the assumption of surface-parallel flow requires that ∂S/∂t must be small compared to . For an ice sheet in steady stale, it is straightforward to estimate the relative magnitudes of ∂S/∂t and .

Steady-State Ice Sheet

It follows from Equation (1) that the surface boundary condition for a steady-state glacier is

(3)

f(t) is a seasonal variable subject to the condition lim

, where ∆t = t ₂ – t ₁. b _S is the mean annual specific mass balance, i.e. the submergence/ emergence velocity measured in m snow/ice of surface density per year. Observations with our airborne SAR indicate that the penetration depth of radar waves into Greenland firn may be several meters even at C-band (5.66 cm wavelength). To account for the fact that the radar apparently does not record the motion at the top surface, but rather the motion at some depth below the surface, the density at this depth should be used in the evaluation of the term -f(t)b_s. Another problem is that firn compaction may not progress at a uniform rate. Firn-settling observations over a 3 month period at the Jarl-Joset station on the Greenland ice sheet, performed during the Expédition Glaciologique Internationale au Groenland (EGIG) expedition 1957-60, indicated a uniform firn-settling rale, except for one event, when a ”firnstoss" gave rise to an instantaneous settling (Haefeli and Brandenberger, 1968). The magnitude of this abrupt surface-elevation change is unknown, since only relative measurements were performed. Further research is needed to investigate the importance of radar-wave penetration into firn and ice and of discontinuities in the firn-settling rate. For the moment, we accept Equation (3) as the surface boundary condition valid for steady-state conditions.

For planar glacier flow, Equation (3) is reduced to

where u_s is the surface velocity in the x direction. The relative importance of the contributions to the vertical particle velocity w_s from the u_s ∂S/∂x and b_s terms is evaluated by using detailed observations from the Greenland ice sheet of surface elevations (Mälzer, 1964), surface velocity (Bauer and others, 1968; Hofmann, 1975), mass balance (Benson, 1962; Ambach, 1963) and surface-layer density (Benson, 1962; De Quervain and others, 1969) along the EGIG profile in central West Greenland (see Figs 2 and 3). Figure 3 shows that, although short-scale fluctuations of u_s ∂S/∂x, over a section of the profile from approximately km 350 to km 480,numerically reach values that often exceed both the large-scale average value of u_s ∂S/∂x and the snow/ice-equivalent specific balance term b_s , the latter term is in general not negligible, but is the dominant term over the main section of the profile.

Fig. 2. The EGIG profile in central West Greenland. (a) Surface elevation profile (Mälzer, 1964) and distribution of horizontal velocity ( Bauer and others, 1968; Hofmann,1975). (b) Distribution of mass balance (Benson, 1962; Ambach, 1963) and surface-layer density (Benson, 1962; De Quervain and others, 1969).

Fig. 3. Comparison of the contributions to the vertical velocity from the specific mass balance ḇ_s and the slope term u_s ∂S/∂x along the EGIG line.

The EGIG-profile data represent the conditions of the central Greenland ice sheet, but probably have general validity for the entire Greenland ice sheet because the u_s ∂S/∂x and ḇ_s terms follow similar trends elsewhere in Greenland. Both terms increase when moving southwards from central Greenland, and decrease when moving northwards. Ofcourse, local deviations from these trends occur that may cause a change in the relative importance of the two terms.

The variation along the EGIG profile of the vertical surface velocity is found by adding the u_s ∂S/∂x and ḇ_s terms, possibly after modifying the b_s term by a seasonally varying factor. However, this gives the correct vertical velocity only for steady-state conditions. By comparing surface elevation profiles obtained by levelling in 1959 and 1968, Seckel (1977) documented that, in this period, the surface elevations of the central part of the EGIG profile increased by an average of 9 cm a^-1. In the period 1968-92, however, the surface elevations of the same section showed an opposite trend, decreasing by an average of 15 cm a^-1. (Kock, 1993; Möller, 1996). These changes are small as compared to u_s ∂S/∂x — ḇ_s (see Fig. 3).

Non-Steady-State Ice Sheet

Observations on Storstrømmen, a large outlet glacier from the ice sheet in northeast Greenland, show that, at present, this glacier is out of balance by up to several meters per year (Reeh and others, 1995). This shows that the vertical velocity associated with non-steady-state glacier flow may reach a magnitude that, if neglected, will impose significant errors on the InSAR-derived horizontal velocity field, in addition to the fact that a non-steady-state vertical velocity component is in itself an important quantity to measure.

In the case of non-steady state, it is not possible to derive an expression for ∂S/∂t and hence a relationship between the components of the surface-velocity vector, by using the surface boundary condition alone. In this case, an additional equation is needed that can be derived by means of the principle of mass conservation (see Appendix A ). The equation is

(4)

where h is ice thickness, F = u /u_S and subscripts "S" and "B" refer to the upper and lower glacier surfaces, respectively. Underlined quantities denote column mean values. Primes denote ice-equivalent quantities.

With typical values of F ≈ 0.9 and ρ _B/ ρ ≈ 1.01, we find (ρ _B/ ρ – 1)(F – 1) ≈ –0.001, and hence, to a very good approximation,

which is identical to the equation used by Rech and Gundestrup (1985) to estimate the local mass change at Dye 3 on the South Greenland ice sheet (see also Paterson, 1994, p. 257). The rigorous derivation of Equation (4) given in Appendix A exposes the approximations behind using the concept of ice-equivalent thickness in mass-flux calculations.

In the case of a non-steady-state glacier, ∂S/∂t must be derived by translating the local mass change as given by Equation (4) into surface elevation change. The outcome of this procedure depends on how the density distribution in the vertical column changes with time. In the ablation zone, where the column consists of ice from surface lo bottom, the density distribution will not change. In the accumulation zone, on the other hand, snow densification processes might result in time-changing depth-density profiles.

As stated by Paterson (1994, p. 14), most discussions of snow densification assume that, at a given place and depth in the glacier, density does not change with time, except near the surface, where newly fallen snow undergoes rapid transformation. This assumption is plausible as long as accumulation rate and snow temperature remain constant. It may also be plausible under non-steady-state conditions as long as the rates of change are moderate. The assumption might, however, be violated during rapid thickness changes such as occur during a surge.

Here, we follow the standard assumption that density depends only on the depth below the surface. As shown in Appendix B, the general equation relating the three components of the surface-velocity vector can then be expressed as

(5)

Equation (5) is not dependent on assuming the glacier to be in steady state. But contrary to Equation (3), the relation involves quantities not at the glacier surface, since it also depends on ice thickness, the depth distributions of density and horizontal velocity, and the mass balance at the ice-sheet base. Except for the ice thickness, the quantities involved generally can be estimated with sufficient accuracy by using glaciological information obtained from field observations and/or theoretical studies.

For a grounded glacier, ∂B/∂t = 0, since isostatic motion of the glacier bed is negligible on the short time-scales involved in SAR interferometry. For a glacier frozen to the base, b_B = 0, but even when there is melting at the base, b_B is at most a few centimeters per year for a grounded glacier. Hence, for a grounded glacier, Equation (5) is approximated by

(6)

In the ablation zone, the glacier consists of ice of constant density from bed to surface. Therefore, ρ = ρS = ρB and Equation (6) is further reduced to

(7)

Besides surface velocity, surface gradient and ice thickness, this equation contains only the shape factor of the horizontal velocity depth profile F, which for glacier flow is typically in the range 0.9 to 1.

Conclusions

Data from the EGIG profile in central West Greenland show that, in the central section of the profile, the decadal rates of surface-elevation change are small compared to the combined mass-balance and slope terms . This indicates that the ice sheet along the EGIG profile is close enough to steady state to justify the use of Equation (3) as an approximation to the relationship between the surface velocity components.

Moreover, we have shown that, although short-scale fluctuations of the surface-slope term S over a section of the EGIG profile numerically reach values that often exceed both the large-scale average value of the surface-slope term and the snow/ice-equivalent specific balance term ḇ_s , it is often not justified to neglect the ḇ_s term compared to the S term and assume surface-parallel flow.

On the other hand, over the large interior accumulation region of the Greenland ice sheet, the error imposed on the horizontal velocity components due to neglecting the ḇ_s term(on the order of magnitude of l m a^-1) is limited. The 230 incidence angle of the ERS-1/2 SAR amplifies a 1 m a^-1 error in vertical velocity to an error of the horizontal velocity of approximately 3 m a^-1. This error is of the same order of magnitude as optimistic estimates of the precision of InSAR-derived horizontal ice velocities that can presently be achieved (Mohr and others, 1998). On the other hand, it is a systematic error which can be avoided. In the interior accumulation region, the error introduced by neglecting the slope term S is even smaller (see Fig. 3).

Nearer to the ice-sheet margin, where both the b_s term and the slope term S may reach values of several meters per year, neglecting any of these terms might impose an error of the horizontal velocity component of some tens of meters per year, i.e. an error that significantly exceeds the noise level. Neglect of a possible non-steady-state vertical velocity component might even considerably enhance the error. Nonetheless, the surface-parallel-flow assumption is, in itself, a significant improvement compared to completely neglecting the influence of the vertical velocity. The surface-parallel-flow assumption also reduces the influence on the SAR measurements of short-scale local surface undulations, thereby easing adjustment of interferograms.

The principle of mass conservation has been used to derive an equation relating the vertical surface velocity to the horizontal velocity vector. Equation (5), valid for both steady-state and non-steady-state conditions, depends on the ice-thickness distribution, and is also influenced, to some extent, by the horizontal velocity-depth and density-depth profiles. We show that replacing the surface-parallel-flow assumption with a correct relationship requires measurement of additional quantities.

The local rate of change of ice thickness or ice mass expressed in Equation (4) is a corollary of the derived surface-velocity condition. Application of this equation for calculating mass changes of ice-sheet regions has hitherto been hampered by the lack of detailed, accurate surface-velocity observations. With the development of the InSAR method for simultaneous determination of ice-sheet surface elevation and velocity, this obstacle to applying the method is now eliminated. By means of InSAR-derived detailed surface elevations and velocities, the local change of ice thickness or ice mass can be estimated over extended regions of the ice sheet provided that ice thickness and specific net balance are also known. Such estimates can be used to calculate regional ice-sheet mass changes which are important for assessing the contribution from the ice sheet to global sea-level change. Also, regional variations of surface elevation change can be studied, which will give important new information about ice-sheet dynamics.