Interferogram generation

Readers interested in background information on radar interferometry may consult Reference Zebker, and Goldstein,Zebker and Goldstein (1986), Reference Goldstein,, Zebker, and Werner.Goldstein and others (1988), Reference Gabriel,, Goldstein, and Zebker,Gabriel and others (1989) and Reference Zebker,, Rosen,, Goldstein,, Gabriel and Werner,Zebker and others (1994). The basic principles of radar interferometry will not be repeated here.

Two interferograms were formed using image 2 (orbit 3248) as the reference image. The complex amplitude radar images were first co-registered with sub-pixel accuracy, including additional pixel offsets over the fast-moving part of the ice. The registered images were then cross-correlated. The normalized correlation or phase coherence of the cross-products, denoted ρ. taking values between 0 (no coherence) and 1 (perfect temporal coherence), is high (ρ > 0.8) over most of the scene, yielding high-quality interferometric fringers ((Fig. 2)). Phase unwrapping was performed using Reference Goldstein,, Zebker, and Werner.Goldstein and others’ (1988) unwrapping algorithm after smoothing of the data using a two-dimensional spectral filter.

Upstream of the grounding line, the interferograms exhibit a complex pattern of closely spaced fringes (360° variations in phase) with pairs of concentric circles where phase unwrapping is difficult to perform. A similar fringe pattern is seen in radar interferograms of the southwestern flank of the Greenland ice sheet. The pairs of concentric circles are attributed to variations in the vertical component of the ice-motion vector as ice flows past bumps and hollows in surface topography, several meters in height and several kilometers in diameter, created by faster ice-sheet flow over the bedrock topograph) near the ice margin (Reference Joughin,, Winebrener, and Fahnestock,Joughin and others, 1995b; Reference Rignot,, Jezek, and Sohn,Rignot and others, 1995). In those regions of more rapidly varying surface slope, phase coherence is red med compared to that of the surrounding ice. because the illumination angle of the ice blocks changes slightly as they move past the bumps and hollows in surface topography.

Fig. 2. Flattened interferogram combining orbits 3205 and 3248. with a color intensity modulated by phase coherence. Dark areas indicate low phase coherence.

The baseline or distance separation between the successive positions of the ERS radar antenna was estimated by least-square fitting using 1400 tie-points selected from a digital elevation model (DEM) of the glacier. 0.005° in latitude spacing and 0.025° in longitude spacing, provided by S. Ekholm and R. Forsberg of KMS and referenced to sea level. The KMS DEM data were projected into the radar-imaging geometry, interpolated using a bilinear interpolation, and registered to the radar scene within 1-2 image pixels using one tie-point.

Interferometric products

I denote V_x, V_y, V_z the components of the steady-motion vector of the ice along the x, y and z axes. The term “steady” here refers to the ice motion over a time-scale much larger than the tidal cycle. The x axis is in the cross-track direction, pointing north. The y axis is in the along-track direction, pointing west. The z axis is the vertical axis. Under tidal influence, the ice tongue undergoes upward and downward motion along the z axis of amplitude Z. 1 use the sign convention that the phase, ϕ measured by the radar is equal to −4πR/ λ, where R is the range distance between a point at the surface of the glacier and the center of the synthetic aperture, and λ is the radar Wavelength (5.66 cm for ERS-1 radar). The phase difference, ϕ ₁₂ = ϕ ₂ − ϕ ₁, measured between antennae 1 (orbit 3205) and 2 (3248), may then be expressed as

where B ₁₂ is the baseline or distance separating antennae 1 and 2, θ _z is the illumination angle with the horizontal for a point at elevation z, δθ _z is equal to θ _z − θ ₀, θ ₀ is the illumination angle with the horizontal at the center of the scene for a point at a reference elevation z = 0, α₁₂ is the baseline angle with the horizontal, B ₁₂⊥ = B ₁₂Sin [θ ₀+ α₁₂) is the component of the baseline perpendicular to the direction of the radar illumination, B _12|| = B ₁₂ cos(θ ₀ + α ₁₂) is the component of the baseline parallel to the direction of the radar illumination, (t ₂ − t ₁) is the time lag between the two images, i is the local incidence angle of the radar illumination with the vertical, and θ ₁₂° is an absolute phase offset.

The first line of Equation (1) depends only on the glacier topography and is scaled by both the baseline and the radar wavelength. The second line is the term of ice motion along the radar line of sight caused by the steady motion of the ice. The component V_y is absent from Equation (1) because the y axis is parallel to the flight direction, and surface displacements are not measured in that direction. The third line corresponds to changes in surface elevation along the radar line of sight caused by the downward motion of the ice under tidal forcing. with a second interferogram combining images 2 and 3. I obtain

with a different baseline separation, B ₃₂, angle, α ₃₂, and relative tidal displacement, Z ₂−Z ₃; but with the same ice motion vector, V_x and V_z .

If radar scenes 1, 2 and 3 are acquired in sequence and exactly 3d apart, adding Equations (1) and (2) eliminates the term of steady ice motion

Using tie-points from the KMS DEM on both rock and ice, at the exclusion of the floating section, I estimate the baseline parameters of Equation (3) and remove the topography term to obtain

which depends only on the tidal displacements. The subscript “flat” designates a phase value for which the effect of the baseline and of surface topography has been removed. The map of the relative tidal displacement, Z ₂ − 0.5(Z ₃ − Z ₁). between scene 2 and the average of scenes 1 and 3, is shown in Figure. 3.

Model predictions from the clastic-beam theory indicate that tidal displacements at a given point along an elastic beam vary linearly with the tidal amplitude. Several experimental studies have shown that the elastic-beam model matches observations official displacements well (Reference Holdsworth,Holdsworth. 1969). If we assume that tidal forcing is the same everywhere along the beam as in Reference Holdsworth,Holdsworth’s (1969) study and that the clastic damping factor of the ice does not change with tidal amplitude (Reference Holdsworth,Holdsworth. 1977), a different realization of tidal forcing should exhibit the same pattern of tidal displacement as that given in Equation (4), scaled by a different relative tidal amplitude. Under these circumstances, I rewrite Equation (1) after removal of the topography term as in

Which yields the ice velocities, V_x and V_z , once γ₁₂ is known. To determine γ₁₂ either the tidal amplitude or the ice velocity must be known at one location, otherwise there is an infinite number of solutions for γ₁₂.

Here, I estimated 11 control velocities by tracking a set of crevasses below the grounding line in two ERS-1 radar images separated by 1 year. The rms error in the velocity estimates is 50 ma⁻¹. The least-square estimate of γ₁₂ is 1.8 ± 0.2.

The x velocity. V_x , is deduced from Equation (5) assuming V_z = 0. in effect, the vertical motion associated with glacier thinning is negligible compared to the horizontal motion, and ice flows nearly in the horizontal plane since the glacier slope is less than 1%. The x velocity was subsequently transformed into a two-dimensional velocity by assuming a flow direction parallel to the dynamic center line of the glacier (dotted line in Figure. 1), The dynamic center line was drawn based on intense surface crevassing at the center of the floating part of the glacier and using the line of maximum x velocity on grounded ice. The result is shown in Figure. 4.

Once the tidal amplitude is known, it is also possible to estimate the glacier topography, at an enhanced spatial resolution compared to the KMS DEM. A color composite image of the glacier topography is shown in Figure. 5. Holes correspond to areas where phase unwrapping failed because of low phase coherence.

Fig. 3. Tidal displacement, (color-coded between −50 and 180 mm and modulated by the radar brightness for display purposes), and hinge line (continuous white line) of Pelermann Gletscher, Dark patches indicate areas with no interferometric data.

Measurement uncertainties

The rms error of the phase values calculated during the baseline estimation process was 3 rad for pair 3248-3205 (Equation (1)), 0.6 for pair 3248-3291 (Equation (2)) and 5 for the two pairs combined (Equation (3)). These phase errors translate into uncertainties in surface topography of. respectively, 75, 400 and 100 m. The errors are large because the perpendicular baselines are short (respectively 58, 2 and 60 m). The velocity errors are conversely small because they do not depend on the baseline separation and are, respectively, equal to 4 and 1 mm d⁻¹ for the individual pairs, which means about 1 ma⁻¹ uncertainty in ice velocity; and 10 mm for the two pairs combined, which means 5 mm uncertainty in relative tidal displacement. These errors are, respectively, three orders of magnitude less than the velocity of the glacier (1000 ma⁻¹) and two orders of magnitude less than the largest relative tidal displacement which can be recorded over a complete cycle (800 mm; see below).

Fig. 4. Ice velocity of Petermann Gletscher, color-coded between 0 and 1200 ma⁻¹. and modulated by the radar brightness. Dark patches indicate areas with no interferometric data.

Fig. 5. Surface topography of Petermann Gletscher color-coded between 0 and 900 m. Dark patches indicate areas with no interferometric data.

Errors in tidal displacement are visible in several areas outside of the main ice stream of Petermann Gletscher. For instance, in the center top of the scene, running almost east west, a 10 km wide segment shows a relative tidal displacement of −10 to −50 mm (colored white in Figure. 3) outside of ihe glacier area. This anomaly coincides with the ice-covered areas of Washington Land and Kane Plateau (see Reference Higgins,Higgins, 1991, (Fig. 2)). This error is probably not due to baseline uncertainties, since its magnitude does not increase with surface elevation. It is probably due to a combination of atmospheric and ionospheric propagation delays (Reference Goldstein,Goldstein, 1995), and surface effects including for instance a change in snow thickness of the ice caps. To average out these errors, additional interferograms are necessary.

Tidal displacements

The pattern of relative tidal displacements ((Fig. 3)) delineates the part of the glacier that is afloat. Nearly the entire section of the glacier below the grounding line undergoes tidal motion, with a sharp discontinuity between the rock margin and the ice tongue. This observation suggests there is little mechanical coupling between the ice tongue and the rock margin. The tidal displacements increase rapidly from zero to a maximum value about 6 km downstream, and subsequently decrease slowly toward an asymptotic value.

On the eastern side of the glacier, where Porsild Gletscher merges with the main stream of Petermann Gletscher, the pattern of tidal displacement is more complex. Phase unwrapping failed at the junction between the two glaciers, but Porsild Gletscher is likely to undergo tidal motion as well. The pattern of tidal motion probably reflects the interplay of the grounding zones from both glaciers.

To explain the pattern of tidal motion derived from the interferometric data and to characterize the flexural rigidity of the ice, I compared a tidal profile extracted along the western half of the ice tongue ((Fig. 1)) with model predictions from an elastic beam of infinite length, with one end rigidly clamped on bedrock (Reference Holdsworth,Holdsworth. 1977). The predicted tidal amplitude at time t and abscissa x along the beam is

where Z_t is the asymptotic value of the tidal displacement at time t referenced t mean sea level, and β is the elastic damping factor of the ice given by

where E is Young’s elastic modulus of ice, ρw = 1030 kg m⁻³ is the density of sea water, g = 9.81 ms⁻² is the acceleration of gravity, v = 0.3 is the Poisson coefficient for ice, and h is the glacier thickness. The best fit is obtained for β = (4.7 ± 0.1) × 10⁻⁴m⁻¹. with a rms fit error of 0.8 mm ((Fig. 6)).

Judging from the low rms error and the high number of points used in the comparison, the model predictions fit the measurements very well and explain the pattern of tidal displacements measured by radar interferometry. The inferred value of β is in reasonable agreement with the curve of Reference Vaughan,Vaughan (1995) relating α to the ice thickness at the grounding line. To obtain a measurement point lying exactly on his curve, the glacier thickness would have to be 863 m, or 288 m thicker than the ice thickness measured by an ice-sounding radar (see below). In an equivalent fashion, the value of E calculated from Equation (7) is 3 ± 0.2 GPa. or three times larger than the value derived by Reference Vaughan,Vaughan (1995) using data from several glaciers. Hence, the method of using Equation (7) is probably not a sensitive predictor of the elastic modulus of the ice.

Fig. 6. Tidal displacements (dots) along a 160 in wide profile shown in Figure. 1. compared to model predictions (solid line) from an elastic-beam theory. The arrow points to the area of maximum bending stress, or hinge line.

Residual errors in model fitting are present downstream from the point of maximum tidal deflection ((Fig. 6)), where we earlier noticed the presence of residual phase errors running across the scene. These errors remain small compared to the tidal displacements recorded on the ice tongue.

Hinge line and grounding line

In the elastic-beam theory (Reference Holdsworth,Holdsworth. 1969), the bending stress of the ice. , where z is the vertical distance to the neutral axis of the beam, and is the bending strain rate, is maximum at the ice surface for z=±h/2 and at the hinge line for x = 0. Locating the maximum of the bending stress, however, involves second-order derivatives of the phases, which increase the noise level of the data. Instead, I propose to define the interferometric grounding line as the location of the minimum relative tidal displacement measured along a lidal profile extracted in the glacier-flow direction ((Fig. 6)), In the elastic-beam theory, the point of minimum deflection also coincides with the hinge line.

The hinge line may not necessarily coincide with the grounding line or with the line of hydrostatic equilibrium of the glacier (Reference Smith.Smith. 1991). The grounding line is typically downstream from the hinge line and upstream from the line of hydrostatic equilibrium. On Petermann Gletscher these three zones are separated by 1-2 km. as discussed later.

The hinge line of Petermann Gletscher is shown in Figure. 3. overlaid on die tidal displacements. Based on the phase noise of the tidal signal ((Fig. 6)), I estimate the hinge line can lie detected within I pixel or 20 m. Near the center of the glacier, the precision is less, because the hinge line shifts in the cross-track direction by several pixels over an across-flow distance of about 500 m, and the tidal profile no longer exhibits a sharp minimum. At the glacier margin, the hinge line cannot be detected accurately, since the fringe rate is too high, phase coherence is much lower and the phase values cannot be unwrapped.

The achieved mapping precision still remains more than one order of magnitude superior to that quoted by Reference Goldstein,, Engelhardt, Kamb, and Frolich,Goldstein and others (1993) who utilized a single radar interferogram. The reason for the lower precision of the single interferogram technique is that it includes the longitudinal gradients in ice velocity which tend to smooth out the local minimum in tidal displacement and bias its location. Here, the bias in absolute location is of the order of several hundred meters. In order to map the hinge line of a floating glacier precisely, it is therefore essential to utilize multiple interferograms and eliminate the longitudinal gradients in ice velocity.

Ice velocities

The ice velocity varies from 400 ma⁻¹ at 900 m elevation to 1100 ma⁻¹ at the grounding line, decreasing thereafter to about 900 ma⁻¹ toward the edge of the scene ((Fig. 7) ). Removal of the tidal signal clearly reduces the variations in ice velocity across the grounding line, yielding a more reasonable velocity profile.

Fig. 7. Ice velocities in the flow direction along the glacier center line before (dotted, grey line) and after (solid line) tidal corrections. The width of the profile is 80 m.

In several parts of the floating section of the glacier, large discontinuities in ice velocity are detected ((Fig. 4)). Along the center line, the eastern slab of the floating ice tongue moves about 30-50 ma⁻¹ faster than the western slab. About 20 km downstream from the grounding line, where the velocity difference between the two slabs reaches 50 ma⁻¹, the velocity of the eastern slab abruptly decreases by 40 ma⁻¹. The discontinuity in velocity and apparent surface rupture reveals an overriding of the northern slab by the faster-moving southern slab. A further 20 km downstream, a similar discontinuity in velocity occurs on the western side of the glacier. The two slabs subsequently move at comparable speeds.

These discontinuities in ice velocity occur in the turning section of the glacier. Surface rupturing could be due to the differential velocity between the two sides of the glacier associated with flow turning. Overriding of the ice indicates that the resistance to ice (low is larger downstream from the grounding line, during the glacier turn, which is consistent with additional friction at the side margins during flow turning.

Ice thickness

I estimated the glacier thickness using the KMS DEM data by assuming that the floating glacier ice tongue is in hydrostatic equilibrium. To perform this calculation. I used an ice density of 917 kg m⁻³ and a density of sea water of 1030 kg m⁻³. The corresponding values of the ice thickness averaged across the glacier width, h , are shown in Figure. 8, with error bars corresponding to one standard deviation in ice thickness across the glacier width.

Fig. 8. Ice thickness measured by the University of Kansas’s airborne ice-sounding radar (ISR, grey line) along the glacier center line compared with the ice thickness deduced from the KMS data (DEM. solid line) assuming hydrostatic equilibrium of the ice and averaged across tat glacier width.

Ice thickness was coincidentally measured by an airborne ice-sounding radar, designed and operated by the University of Kansas (Reference Chuah,, Gogineni,, Allen and WohletzChuah and others, 1996) and flying on board a NASA P-3 aircraft, along the center line of Petermann Gletscher, on 26 May 1995. The ice-sounding radar operates a coherent radar system at a center frequency of 150 MHz.

At the hinge line, ice thickness is slightly over-estimated by the DEM data ((Fig. 8)). The two thickness estimates become equal at about 2.7 ± 0.5 km from the hinge line, which indicates the approximate location of the line of hydrostatic equilibrium of the glacier tongue. Between km 25 and 45. the ice-sounding radar indicates a larger ice thickness than calculated from hydrostatic equilibrium of the ice. Conversely, between km 5 and 25. ice thickness is slightly overestimated by the DFM data. If the KMS DEM is accurate, this result suggests that either a large part of the ice tongue is not in hydrostatic equilibrium, or the ice-sounding radar data include spatial irregularities in ice thickness that do not correspond to steady-state conditions. Multiple ice-sounding radar profiles are needed to interpret these differences more completely.

About 1 ± 0.5 km from the hinge line, the ice-sounding radar data transition to a regime of larger variations in ice thickness (Fig- 8) and the basal reflections produce hyperbolic radio-echo records (Jezeck and others. 1995). This feature is likely caused by bottom crevassing of the glacier and indicates the approximate position of the grounding line (Reference Jezek,, Gogineni and Rignot.Jezek and others. 1995). The grounding line, the hinge line and the line of hydrostatic equilibrium of Peterniann Gletscher are therefore not coincident and are probably separated by about 1-2 km.

Ice discharge

The ice flux, Q. from Petermann Gletscher at and below the hinge line is calculated using

where V _x is the ice velocity along the x direction averaged across the glacier width, and W_y is the glacier width measured in the y direction (Reference Paterson,Paterson, 1994). The glacier thickness used in the calculation is that derived from the KMS data. Where the glacier is afloat. basal velocities should equal the surface velocities, so the surface velocities measured by radar interferometry arc equivalent to vertically integrated ice velocities.

Ice discharge is 12.3 ± 1 kin³ a⁻¹ at the hinge line and 12.1 ± I km³ a⁻¹ about 1 km downstream where the glacier is more likely to he in hydrostatic equilibrium Fig. (9). The measurement error is associated with uncertainties in ice velocity at the ice margin, and with uncertainties in ice thickness. About 30 km upstream from the hinge line, Joughlin and Others (1995b) estimated an ice flux of 12.7 km³ a⁻¹, consistent with our results. Downstream from the hinge line, the ice discharge decreases rapidly. At km52, the ice discharge is only 2.1 km³ a⁻¹. These results are consistent with those obtained by Reference Higgins,Higgins (1991) near the ice front ((Fig. 9)). At the glacier front, Reference Higgins,Higgins (1991) estimated a calf-ice production of only 0.59 km³ a⁻¹, or 20 times smaller than the ice discharge at the hinge line. Ablation processes therefore melt more than 95% of the ice that crosses the hinge line. Calf-ice production plays a minor role in the mass discharge to the ocean from Petermann Gletscher.

Fig. 9. Ice discharge (dots) of Petermann Gletscher at and below the hinge line, and average discharge (solid line) averaged over 1.5 km segments. Diamond symbols indicate the ice-flux estimates by Reference Higgins,Higgins (1991).

Melt rates

Assuming that the ice density is everywhere the same, the equation of mass conservation integrated vertically and across the glacier width is

where is the glacier net balance, positive if the glacier accumulates mass (Reference Paterson,Paterson. 1994), I now assume that the glacier is in steady-state conditions, meaning ∂ h /∂t = 0, and calculate the glacier net balance, , from the gradient in ice flux divided by die glacier width.

The largest source of error is the uncertainty in thickness gradient. To reduce that error. I calculate the glacier net balance at a discrete number of locations, over 5-10 km long segments. The results are shown in Figure. 10. The largest value, about 24 ± 5 ma⁻¹. is recorded close to the grounding line. Near the ice front, the glacier net balance is only 2-3 m a¹ (Reference Higgins,Higgins, 1991). The ice flux decreases from 12.1 to 1.5 km³a⁻¹ over an ablation area of 886 Km², corresponding to an average melt rate of 12 ± 2 m a¹.

Fig. 10. Melt rate of the floating ice tongue vs the downstream distance from the hinge line with 20% error bars.