1. Introduction
1.1. Nivlisen
The ice shelves around Antarctica are important elements in the climate coupling of the Antarctic ice masses. The grounding line is the output gate for mass loss of the grounded ice sheet. The ice shelves themselves have intensive interactions with the atmosphere (through ice accumulation and ablation, radiation and albedo) and with the ocean (through bottom melting and freezing and iceberg calving). Hence, ice shelves are particularly sensitive to climatic influences and changes.
Here we investigate Nivlisen, one of a number of small ice shelves bordering Antarctica’s Atlantic sector. Figure 1 shows its geographic situation. The ice shelf extends over about 80 km in the south–north and 130 km in the west–east direction. Ice thicknesses vary from about 700 m in the southeast to 150m at the front. At some locations the ice rests on the bedrock, forming ice rises and ice rumples. The ice shelf is nourished by glaciers draining the inland ice of Wegenerisen through the mountains of Wohlthatmassivet. Typical flow velocities are in the order of 100 ma–1.
The surface regimes, and, accordingly, firn and ice structures in the upper layers, vary significantly over Nivlisen. A bare ice region extends over about 100 km to the southeast. This is predominantly a large ablation zone preserved by dry katabatic winds and summer melting due to the reduced albedo. To the north, there is a transition to an extended accumulation zone. On the border between the inland ice and the ice shelf, there is an ice-free area of bedrock (a so-called oasis) maintained due to the ablation regime of the surrounding ice. This area is called Schirmacheroasen.
The ice shelf can be considered in hydrostatic balance at positions a few kilometres away from grounded ice. This property is supported both by mechanical modelling (Reference Rabus and LangRabus and Lang, 2002) and by observational evidence. Synthetic aperture radar (SAR) interferometry from the European Remote-sensing Satellite 1 and 2 (ERS-1/-2) tandem mission shows distinct zones of vertical deformation along the grounding lines, indicating that the ice shelf adapts to tidal and other sea-level displacements along these zones and keeps its bulk in hydrostatic balance (Reference Metzig, Dietrich, Korth, Perlt, Hartmann and WinzerMetzig and others, 2000; see also section 2.3). Reference Korth, Perlt, Dietrich and DietrichKorth and others (2000) used several kinematic global positioning system (GPS) tracks across the grounding zone to show in detail how the reaction to tidal displacements evolves from zero on the grounded side to the full tidal displacement over a distance of 2–4 km.
Not least due to the logistic basis provided by the Soviet/Russian, Indian and former East German (GDR) research stations in Schirmacheroasen, the region has been the subject of extensive geoscientific investigations for decades. Reference Bormann and FritzscheBormann and Fritzsche (1995) give a comprehensive review of research up to the beginning of the 1990s. Since that time, geodetic–glaciological field observations have been widely extended. New types of satellite-aided observations have become available, such as GPS measurements on the ground, and satellite altimetry and SAR interferometry from space.
This paper analyzes some of these newer observations which allow a more extensive, more detailed and more accurate description of Nivlisen’s glaciological regime than presented thus far. After a description of our basic observations and their processing (section 2), we present in section 3 the additional inferences from combining the individual types of data under two aspects: hydrostatic equilibrium and mass balance. (These fundamental links between major ice-shelf parameters are briefly reviewed in section 1.2 and 1.3.) In section 4 we discuss our findings, addressing the spatially varying firn conditions visible from the hydrostatic balance relation, basal processes and possible temporal changes, the interplay between snow accumulation and compaction and dynamical ice-flow conditions.
1.2. Hydrostatic balance of ice shelves
Free floating of an ice shelf implies the hydrostatic balance relation
where Z is the total ice thickness, h the surface height above sea level (freeboard height), ρ w the mean density of the displaced sea-water column, and the mean density of the ice-shelf column.
This relation has been widely exploited in ice-shelf studies. See, for example, Reference RennerRenner (1969) for an early review, and Reference Coslett, Guyatt and ThomasCoslett and others (1975) and Reference Shabtaie and BentleyShabtaie and Bentley (1982) for the methodology adopted here. Through hydrostatic balance, surface heights were derived from ice thicknesses (Reference Bamber and BentleyBamber and Bentley, 1994), and ice thicknesses from surface heights (Reference VaughanVaughan and others, 1995), or geopotential models were locally evaluated by the hydrostatic balance requirement (Reference Shabtaie and BentleyShabtaie and Bentley, 1987; Reference Jenkins and DoakeJenkins and Doake, 1991). Apparent violations of hydrostatic balance gave indications for grounding (Reference Shabtaie and BentleyShabtaie and Bentley, 1987), for varying firn and ice structures (Reference Budd, Corry and JackaBudd and others, 1982; Reference Bamber and BentleyBamber and Bentley, 1994), for sediments transported by the ice (Reference Thyssen and GrosfeldThyssen and Grosfeld, 1988) or for marine ice frozen to the shelf ice bottom (Reference ThyssenThyssen, 1988; Reference Jenkins and DoakeJenkins and Doake, 1991).
Equation (1) may be reformulated in the form
where ρ i is the density of pure ice and δ is defined by
That is, δ is the relative deviation of the ice-shelf density ρ from the pure ice density ρ i, integrated along the ice column. Resolved for its three main constituents, Equation (2) reads
In the common case that an ice shelf consists of a mixture of pure fresh-water ice and air, the meaning of δ is plausible: a (fictive) rearrangement of the ice shelf into a body of solid ice and a layer of air would result in the solid ice thickness
and the air layer thickness δ. Since most of the air is contained in the upper firn layers, δ depends on the accumulation and densification regime but has little or no dependence on total ice thickness. The mean density , in contrast, depends on both. Values of air layer thicknesses δ reported in the literature cited in this section vary between 10 and 20 m
Unless the ice shelf consists only of fresh-water ice and air, δ is not the air layer thickness. Therefore, δ is called the ‘apparent air layer thickness’ in the following. If, for example, a marine ice layer of thickness Z mar and density ρ mar > ρ i is frozen to the bottom of a meteoric ice part of thickness Z met, Equation (3) reads
δ is now the sum of the air layer thickness δ air and a negative correction δ mar. However, since ρ mar ≈ 925kgm–3 exceeds ρ i only slightly (see section 2.4), δ mar ≈ —0.02 Zmar tends to be small. Sediments enclosed in the lower ice-shelf layers are another possible phenomenon that adds a negative contribution to δ and may even cause δ to be negative in total. The apparent air layer thickness δ when computed from surface height and ice thickness using Equation (4) is a useful indication of density structures and of the validity of the free-float assumption. At the same time, it is sensitive to errors in the input data, in particular in Z and h. Gross errors in Z may arise if a marine ice layer is not recorded in the measurement, as may occur with radar techniques. Biases in the freeboard height h may arise from uncertainties in the geoid and sea surface topography (SST) data needed to reduce ellipsoidal heights to freeboard heights. Therefore, δ is frequently treated as a combination of the air layer thickness and an unknown height bias.
Once δ is known for a certain region of free floating, the relation between ice thickness and freeboard height is established and can be utilized.
1.3. Mass flux and mass balance on ice shelves
The mass conservation equation for a vertical column of an ice shelf (Reference PatersonPaterson, 1994) can be written as
Here, b s and b b (in units of mass per surface area per time, positive for mass gain) are the surface mass balance and the basal mass balance, respectively, Z i is the solid ice thickness (Equation (7)), ρ i δZ i/δt is the local ice mass change with time, and u is the depth-independent horizontal flow velocity vector with components u and v in the x and y coordinate directions, respectively. According to Reference Bamber and BentleyBamber and Payne (2003), we call b s + b b the specific mass balance and ρ i δZ i/δt the local mass balance. Their integral effect over a certain region A can be computed through Gauss’ theorem by balancing the ice in-flux and out-flux across the region boundaries:
where K is the boundary curve oriented anticlockwise, and β is the (anticlockwise) angle between K and the velocity u . Conveniently, some sections of K follow flowlines where the integral vanishes, and other sections form the flux gates.
The basal balance b b is commonly the most uncertain mass-balance term. The above equations may constrain it based on the geometric magnitudes of ice thickness, freeboard height and ice-flow velocity.
2. Data
2.1. Ice thickness
Figure 2 gives an overview of ice-thickness observations on Nivlisen. In this study, we use data from the three sources shown by tracks in Figure 2. A radar system operating from a helicopter was used by the German Federal Institute for Geosciences and Natural Resources (BGR) during the GeoMaud expedition in 1995/96 (Reference Damm and EisenburgerDamm and Eisenburger, in press). A similar instrumentation, but mounted on an airplane, was used by the Alfred Wegener Institute (AWI), Bremerhaven, in the same summer season (Reference Steinhage, Nixdorf, Meyer and MillerSteinhage and others, 1999; Reference Meyer, Steinhage, Nixdorf and MillerMeyer and others, 2005). Finally, radar measurements along a ground traverse crossing the east of the ice shelf were carried out by the AWI in 1992 (Reference FritzscheFritzsche, 2005).
The relative precision of the ice-thickness measurements was estimated to be 1–2%. In a crossover analysis for all radar tracks (including grounded ice), crossover differences at 61 positions, divided by the absolute thicknesses, showed an rms value of 4.5%. This value is affected by errors in the measurement positions and in the interpolation to crossover points which both, in turn, depend on along-track variations of the measurements. The ice-thickness tracks on Nivlisen seem to have two fluctuation regimes, with a distinct border along approximately –70.35° latitude (see Fig. 8 for an example). In the south, fluctuations are high, with typically 10–20 m differences between track points 300 m apart. In the north, the thickness curves are smooth. Accordingly, the rms of nine relative crossover differences in the northern part is 1.9%, indicating a measurement precision of 1.9%/ i.e. typically 1.4% × 350 m = 4.9 m. In what follows, we will rely on thickness data with a reduced resolution of 3.7 km. For these data a relative precision of 1.4% seems realistic over the entire ice shelf.
Despite this precision, biases in the thickness data must be considered. The along-track fluctuations in the southern part may indicate a complex (e.g. crevassed) structure of the ice-shelf bottom, with associated ambiguities in the radar echo. The estimated inaccuracy resulting from these effects is up to 4.4% of ice thickness. Variance combination with the above 1.4% precision leads to a total uncertainty of in the southern part, i.e. typically 4.6% x 450m = 21m. Another bias in the initial ice-thickness data results from using the signal propagation velocity for ice (168 × 106ms–1) in the radar processing. In section 3.2 we apply the correction for the higher velocities in firn.
Ice-thickness data from other measurements (see Fig. 2) were analyzed, but not included in the further computations. During the 1995/96 GeoMaud expedition, ice thickness was measured at ten ground positions on Nivlisen by radar sounding and partly also by seismic sounding (Reference ReitmayrReitmayr, 1996). For part of the data the interpretation of the signal return from only a single point measurement is ambiguous. In 1975/76, Soviet Antarctic Expedition groups performed two ice-core drillings through the ice shelf, and a third through the nearby grounded ice (Reference Korotkevich, Savatyugin and MorevKorotkevich and others, 1978). In the same year, radio-echo soundings were performed from airplanes (Reference Kozlovskii and FedorovKozlovskii and Fedorov, 1983) and along ground traces (Reference Boyarskii, Nikitin and StepanovBoyarskii and others, 1983; Reference Eskin and BoyarskiiEskin and Boyarskii, 1985). Differences between these historical data and the models derived hereafter are further addressed in section 4.3.
2.2. Ice surface height
2.2.1. Ellipsoidal height
Ellipsoidal heights were observed by satellite altimetry and by GPS measurements. Figure 3 gives an overview.
We performed a dedicated analysis of ERS-1 data of the geodetic mission phase E/F (April 1994–March 1995) processed by the Centre Nationale d’Etudes Spatiales, France (Reference Rémy, Shaeffer and Legrésy.Rémy and others, 1999). The data contain geophysically corrected and retracked ellipsoidal heights at satellite subtrack positions, together with the waveform parameters of the leading edge, the trailing edge and the backscatter obtained in the retracking procedure.
At the steep transition between floating and grounded ice, the altimeter tracker fails to follow the surface. This leads to an interruption of tracks along the grounding zone. In contrast to continental-scale processings, no attempt was made to interpolate heights across this gap, and only measurements on the ice shelf were used. Moreover, to avoid slope errors due to the poorly known slopes of the adjacent grounded ice, positions near the grounding line were excluded. The minimum distance was set to 8 km, allowing for a 0.02 rad slope beyond the grounding line.
Strongly reflecting targets such as meltwater surfaces can divert the on-board tracker from the nearest surface points. Such instances are indicated by exceptionally high back-scattering coefficients σº and steep trailing-edge slopes Fl of the altimeter waveform. Hence, measurements exceeding σº = 30dB or Fl = –0.024dB per gate (1 gate = 3 ns) were rejected and their neighbouring positions were checked for obvious corruptions.
Static GPS measurements were carried out during the GeoMaud expedition in 1995/96 (Reference ReitmayrReitmayr, 1996). Their vertical precision is better than 10 cm. Kinematic GPS measurements along a ground traverse through the eastern part of Nivlisen were carried out in 1998 by Technische Universität Dresden (Reference Korth, Perlt, Dietrich and DietrichKorth and others, 2000). The accuracy of the processed heights is better than 5 cm.
2.2.2. Reductions to freeboard height: geoid, sea surface topography and tides
The transformation from ellipsoidal heights H to freeboard heights h required in the hydrostatic balance relation reads
where the geoid height N is the height of a (hypothetical) ocean at rest, ζSST is the deviation from the geoid due to ocean dynamics and density variations, and ζTides is the tidal displacement of the sea surface and, correspondingly, of the freely floating ice.
The geoid height N is the largest correction term. Thanks to the new satellite gravity missions CHAMP (CHAllenging Mini-Satellite Payload, launched in 2000), GRACE (Gravity Recovery and Climate Experiment, launched in 2002) and GOCE (Gravity Field and Steady-State Ocean Circulation Explorer, to be launched in 2006), the accuracy of global geoid models is rapidly improving, in particular over Antarctica. Further refinements may be achieved with airborne and ground observations. For Nivlisen, such a refinement is available by means of a regional geoid model (Reference Korth, Dietrich, Reitmayr, Damm, Forsberg, Feissel and DietrichKorth and others, 1998) derived from surface gravimetry together with airborne and ground-based ice-thickness measurements. Its estimated accuracy is 50 cm for the absolute geoid level and 7 cm for the geoid variations within the model domain. Figure 4 shows a regional comparison between the geoids according to the Earth Geopotential Model 1996 (EGM96, a state-of-the-art geopotential model before the new missions), preliminary CHAMP and GRACE models and the regional model. All the models show an average level of 15–18m over the region, with EGM96 deviating most from the regional model. The larger-scale variations of the regional model are comparably well reproduced by the GRACE model and by EGM96, although the preliminary GRACE model uses only 39 days of satellite observations while EGM96 includes virtually all information available by 1996. However, even with its full performance, GRACE will not resolve small-scale geoid variations, due to their attenuation at satellite height. In the following, we use the regional geoid model as the best available source. Adding the permanent tide effect of –0.19m transforms it from the tide-free to the mean-tide system to be consistent with the ellipsoidal height data.
For ζ SST, we use the large-scale model estimated simultaneously with the EGM96 geopotential model from a combination of satellite altimetry, satellite tracking data and gravimetric data (Reference Lemoine, Segawa, Fujimoto and OkuboLemoine and others, 1997), which gives a value of –1.7 m for Nivlisen. The geoid to which the EGM96 SST refers has a different reference gravity potential (W 0 = 6.263685857 × 107m2s–2) than our regional geoid model (W 0 = 6.263686085 × 107m2s–2, according to the Geodetic Reference System 1980). Hence, the SST height was corrected by the respective height difference (23 cm) to a value of –1.47 m (see, e.g., Reference SmithSmith, 1998, for the subtleties of geoid heights). Smaller-scale features of the SST are not corrected but do enter the error budget. The uncertainties assigned to the SST model are 0.3 m for its absolute level and 0.14 m for its spatial variations.
Temporal variations of SST are not corrected. Their largest part is the ocean’s inverse barometric effect. This was verified from tide gauge observations and had an rms in the order of 13 cm (Reference Dietrich, Liebsch, Dittfeld and NoackDietrich and others, 1995, Reference Dietrich, Dach, Korth, Polzin and Scheinert1998). The estimated total rms of temporal variations is 0.2 m.
The tidal displacements ζ Tides of the freely floating ice shelf are predicted from a precise tidal model. It was derived from 15 months of tide gauge observations in an epishelf lake in Schirmacheroasen which is linked to the ocean (Reference Dietrich, Liebsch, Dittfeld and NoackDietrich and others, 1995, Reference Dietrich, Dach, Korth, Polzin and Scheinert1998). The applicability of this model to the whole of Nivlisen was validated by kinematic GPS observations at several positions on the ice shelf (Reference Korth, Perlt, Dietrich and DietrichKorth and others, 2000). After reduction of the inverse barometric effect, the residual rms error of predicted tides was 4 cm.
2.2.3. Compiled elevation models
From the altimetry data that passed the selection described in section 2.2.1 (see Fig. 3) a model of the temporal mean ellipsoidal heights on Nivlisen was generated with a continuous curvature spline algorithm (Reference Wessel and SmithWessel and Smith, 1991). Heights on floating ice were corrected in advance for their instantaneous tidal displacements. The slope error was corrected by the simple slope correction method (Reference Brenner, Bindschadler, Thomas and ZwallyBrenner and others, 1983). Since computation of surface slopes already requires a topography model, the model generation was iterated. The resulting topography model has a resolution of 6’ × 2’, i.e., about 3.7 × 3.7 km2. It thus reflects heights averaged in a 3.7km scale. Values at arbitrary positions within the model domain are obtained by bicubic interpolation.
A rough error estimation of the model is summarized in Table 1. We distinguish several kinds of errors. First, errors in the altimetric data arise from orbit errors, errors in the tropospheric delay correction, instrument noise and, most importantly, the complex effects of surface roughness and penetration on the retracked heights (Reference Rémy, Shaeffer and Legrésy.Rémy and others, 1999). The rms of these errors is estimated to be 0.4 m. Second, uncorrected temporal SST variations may not fully average out. They were assigned with an rms of 0.2 m in section 2.2.2. Third, the limited spatial resolution and coverage of the measurements affects the model’s accuracy. The resolution results not only from the sampling but also from the 3 km scale pulse-limited footprint of the ERS altimeter. The applied slope correction does not consider surface curvature effects, but an additional curvature correction (Reference Brisset and Rémy.Brisset and Rémy, 1996) would exceed 10 cm only in the area north of the ice rumple about –70.4º S, 10.7° E, where the topography is probably too complex for this correction to apply. Near the boundary of the model domain, interpolation of heights turns into extrapolation, providing an additional error source. Based on the information available on smaller-scale slopes and curvatures (Reference Wiehl and DietrichWiehl and Dietrich, 2000), the estimated rms error due to the limited spatial sampling is 0.8 m for point heights and 0.5 m for the 3.7 km scale averages represented by the model, except for positions near the domain boundary and the above-mentioned ice rumple where our estimate is 1.5 m. Combining the three estimated errors in the sense of variance combination, the estimated rms error of the model is 0.67 m (0.92 m for point heights) in general and 1.6 m near the domain boundary and the ice rumple.
For internal validation, a crossover analysis of the selected data revealed an rms crossover difference of 0.51 m, supporting the error estimate for the altimetric data plus uncorrected temporal variations. For external validation, the altimetric model was compared to the GPS heights. The rms difference to the nine static GPS measurements in the altimetric model domain was 0.82 m, supporting the estimated 0.92 m error for point heights. A comparison with the kinematic GPS measurements along the traverse is shown in Figure 5. It shows agreement within 1 m except for the two margins. In the northwestern part near the ice front, the altimetric model smooths the topographic details of the ice rumples. In the southern part, extrapolation near the altimetric domain boundary causes larger errors.
The GPS heights (after their tidal correction at floating ice positions) were then added to the altimetric data to yield the ellipsoidal height model with extended data support and coverage.
Finally, the modelled ellipsoidal heights at floating positions were reduced by the geoid height and the constant part of the SST to obtain freeboard heights. The errors of this freeboard height model (Table 1) arise from variance combination of the ellipsoidal height model errors and the errors in the freeboard height corrections discussed in section 2.2.2.
2.3. Ice-flow velocity and grounding zone
2.3.1. Observations
Interferometric processing of satellite SAR data was applied to infer a two-dimensional high-resolution dataset of horizontal flow velocities as well as the limits of the grounding zones. Two ascending, partly overlapping ERS-1/-2 stripes were acquired over Nivlisen during the tandem mission in May 1996 (Table 2). The temporal baseline of exactly 1 day between subsequent acquisitions as well as a relatively short interferometric baseline yields good coherence in both interferograms. In addition to the SAR data, ice-flow velocities were observed by repeated GPS observations in 1992 and 1993 in the northeast of the ice shelf (arrows in Fig. 7). The velocities and flow directions have an accuracy of 2ma–1 and 2.5°, respectively, (Reference Korth and DietrichKorth and Dietrich, 1996) which is suitable for calibration purposes.
2.3.2. SAR processing
The SAR sensors of the ERS-1/-2 satellites peer perpendicular to the flight direction with an approximate incidence angle of 23°. Due to the sensor’s operating principle, only changes in its line of sight can be observed, so additional information is needed to infer the three-dimensional flow vector or even its two horizontal components. Using the different views of ascending and descending passes is a typical way to overcome this problem. However, even if interferograms from descending passes were available, this method would be difficult to apply over floating ice because different vertical displacements between acquisitions cause significant phase changes. Although a precise tide model exists for Nivlisen (see section 2.2.2), air-pressure differences and, hence, the inverse barometer reaction are unknown. However, if flow directions are known they can be deployed in conjunction with the assumption of surface-parallel flow to add the required information. Of course, this technique is also sensitive to ocean-induced vertical changes, but, assuming that these are uniform over the ice shelf, their effect cancels out when tying the velocity field to an absolute reference. Phase gradients of the interferogram can be used to verify this assumption and, moreover, can be exploited to discover regions of stronger non-uniform deformation.
For most parts of Nivlisen, flowlines were identified in the low-pass filtered geocoded SAR amplitude image (see Fig. 6a). They were digitized and their directional information was interpolated to a grid of flow azimuths. By combination with surface slopes calculated from the digital elevation model (DEM) derived in this study (sections 2.2.3 and 3.2), three-dimensional flow directions were obtained.
The interferograms were calculated and their phase gradients were determined (Fig. 6b). They support the assumption of a nearly uniform vertical displacement of the floating ice. Moreover, grounding zones can be clearly detected by the larger gradients induced by inhomogeneous vertical deformations. Upper and lower limits of flexure were mapped and used to mask out areas of non-uniform vertical deformation.
The effects of Earth curvature and topography were corrected by subtracting a synthetic interferogram calculated using precise ERS orbits (Reference Scharroo and VisserScharroo and Visser, 1998) and our DEM. The remaining interferogram (Fig. 6c) should include deformation effects only. After masking out grounded areas, grounding zones and regions with disrupted phase values, the modulo 2π phase was unwrapped by the branch cut method (Reference Goldstein, Zebker and WernerGoldstein and others, 1988) to obtain absolute phase differences which represent the projection of flow velocity differences onto the SAR line of sight. Then, using the three-dimensional flow direction, all three components of the displacement vectors were calculated. Each interferogram was processed separately in this way. In the region of overlap, the final product was generated by an average with linearly changing weights. Then the relative velocity field obtained so far was referenced to absolute values by the GPS ground-truth observations.
2.3.3. Results
Figure 7 presents the horizontal velocity field obtained for Nivlisen.
To assess the accuracy of the applied technique, we consider that the horizontal velocity |u| is a function of three observations: the range change in the satellite’s line of sight (divided by the temporal baseline), r; the angle between the flow direction and the satellite’s line of sight azimuth, α, which is inferred from the digitized flow directions; and the surface slope, β, inferred from the DEM employed. Neglecting the comparatively small effect of β, the functional relation is
where the radar incidence angle θ is precisely known. Those regions where 74° < α < 106°, i.e. where the factor cos¯1 α leads to >350% error magnification, were masked out a priori. An error propagation from r and α to |u| was then computed by
The error σα of the digitized flowlines was estimated to be 5°. An evaluation of σr is more difficult since any errors in the interferometric analysis (e.g. co-registration and interpolation errors, DEM and baseline errors or atmospheric or ionospheric influences) propagate into r. With reference to the precise orbits of Reference Scharroo and VisserScharroo and Visser (1998) and to the stochastic model for radar interferometry by Reference HanssenHanssen (2001), an estimated range change error of one-quarter of a fringe, leading to σr = 2.6 m a–1, seems reasonable. Figure 6d depicts the resulting velocity error σ|μ|. Its variations are predominantly induced by variations in α.
These error estimates were supported by a comparison between the two SAR tracks in their region of overlap, which indicated even better precision. However, the only absolute reference for the entire velocity field is the three GPS points in the northeast, matching with ±3.1 ma–1. Since the assumption of uniform vertical displacement of the ice shelf is only an approximation, residual deformations wrongly interpreted as ice flow may cause additional errors.
The detected grounding-zone limits (see section 2.3.2; Fig. 6) are used as criteria for free float later in this study. Seaward from the flexure zones the ice is considered to be in hydrostatic balance. Low tide limits are used where they differ from high tide limits. In our maps, however, we depict the most landward flexure limits to indicate the grounding line.
2.4. Surface mass-balance and density data
As mentioned in section 1.1, Nivlisen is divided into an ablation zone in the southeast and an accumulation zone in its remaining parts. Accumulation is due to both precipitation and drift snow. Part of the accumulated snow melts in summer and percolates into deeper horizons, forming ice layers and lenses (Reference Bormann and FritzscheBormann and Fritzsche, 1995).
Ablation rates were observed over up to 5 years (1988–93) at stakes along two traverses through the grounded part of the ablation area (Reference Korth and DietrichKorth and Dietrich, 1996). The mean ablation at 47 stakes was 230 kgm–2 a–1 with a 100 kg m–2 a–1 standard deviation of the spatial variations. For the close vicinity of Schirmacheroasen, Reference SimonovSimonov (1971) reported 340 kg m–2 a–1 ablation. Accumulation rates along the traverse crossing the ice shelf (cf. Fig. 2) were measured at seven stakes in 1992/93. The mean value is 316 kgm–2 a–1. Reference SimonovSimonov (1971) reported a mean accumulation rate of 230–240 kg m–2 a–1 in the years 1959–65. Summing up these references, the surface mass balance is typically –250± 100 kgm–2 a–1 in the ablation zone and +300± 100 kg m–2 a–1 in the accumulation zone.
Density information is available down to 610m depth from a deep core drilling in the grounded accumulation area south of Schirmacheroasen (Reference Kislov, Manevskii and SavatyuginKislov and others, 1983). These data together with densities from a nearby shallow ice core (Reference FritzscheFritzsche, 1996) may be fitted by the exponential relation (Reference PatersonPaterson, 1994)
with the depth z, the surface snow density ρ s = 420 kg m–3 and the densification parameter C = 0.029 m–1. The resulting air layer thickness from this fit is δ≈ (ρ i–ρ s)/(ρ i C) = 18.7 m. On the ice shelf itself, densities down to 15 m depth were observed from three ice cores (see Fig. 2 for locations) by the 1982/83 Soviet Antarctic Expedition (Reference KhokhlovKhokhlov, 1984). To the south, in the ablation area, vertically homogeneous ice with a large amount of air bubbles and a mean density of 879 kg m–3 was found. In the transition zone from the ablation to the accumulation regime, the ice below 12 m resembled the ice in the ablation zone, while the upper 12 m consisted of both ice and consolidated firn. The third core near the ice front contained frozen firn of 640 kg m–3 density at 10–15 m depth. In the upper 10 m, ice layers dominated in the interchange with firn layers. These cores may illustrate the complexity and spatial variability of meteorological influences on the firn and ice structure.
For the densities of ocean water, ρ w, and pure fresh-water ice, ρ i, required for hydrostatic balance considerations, we use the values ρ w = 1029 kg m–3 and ρ i = 917 kg m–3. Ocean water densities applied, for example, by Reference Shabtaie and BentleyShabtaie and Bentley (1982), Reference Jenkins and DoakeJenkins and Doake (1991) and Reference VaughanVaughan and others (1995) at different ice shelves vary slightly between 1027.5 and 1030 kgm–3, indicating an uncertainty of 1.5 kg m–3 in our value. The pure ice density, ρ i, applied in the works cited above, varies between 915 and 920 kg m–3. We assign an uncertainty of 3 kg m–3 to our pure ice density value. In hypothetical computations assuming a marine ice layer (sections 1.2 and 4.3) we use the marine ice density ρ mar = 925 kgm–3 according to Reference Jenkins and DoakeJenkins and Doake (1991).
3. Derived Glaciological Information
For the combined analysis of the various kinds of data compiled in section 2 we assume that all the observations included represent one and the same epoch (say, 1 January 1996). Although they are obtained in a certain time-span, mainly 1992–98, temporal variations of the observables are considered to be small compared to errors in the observations themselves. The issue of temporal changes is further discussed in section 4.2.
3.1. Spatial variations of the hydrostatic balance parameter δ
To exploit the hydrostatic balance relation between ice thickness and freeboard height, we selected those thickness measurements lying in the domain of our freeboard height model and on the freely floating ice shelf according to the limits explained in section 2.3.3. For these positions we applied Equation (4) to compute the apparent air layer thickness δ.
Figure 8 illustrates the result for the data along one radar track. From the irregularly distributed values of δ a gridded model was derived by a continuous curvature spline algorithm (Reference Wessel and SmithWessel and Smith, 1991). This model, shown in Figure 9, is intended to represent δ in the sense of 3.7 km scale averages.
Clearly, δ shows a large spatial variation with a distinct pattern. In the southeast of the ice shelf, values generally vary between 1 and 4 m. To the northwest of this region, δ increases from 4 to 14 m within a distance of about 22 km. Northwest of this zone of increase, values are generally 14–20 m.
Table 3 summarizes an error assessment for δ obtained by variance propagation from the errors of the individual inputs of Equation (4). We consider the total error as a sum of the error of the spatial variations and the error of the absolute level. In particular, errors in ρ w and ρ i induce errors in δ which are proportional to Z–h, i.e. positive everywhere. Therefore, the error propagated from ρ w and ρ i is the sum of a constant part associated with the mean value of Z–h and a spatially varying part associated with the spatial variations of Z–h. Errors of interpolating and slightly extrapolating δ are roughly estimated from the local variability of δ. The resulting estimated errors in the spatial variations of δ are 2.8 m or less. The discovered variation of δ is, hence, much larger than its uncertainty and must reflect a geophysical signal.
3.2. Improved ice-thickness and height information
In section 4.1 we argue that the obtained model for δ can in fact be considered to represent the air layer thickness. In anticipation of this outcome, we use the information on firn density inherent in δ to determine the firn correction ΔZ for the ice-thickness measurements (see section 2.1). It can be given as
(Reference Jenkins and DoakeJenkins and Doake, 1991), with the pure ice refractive index n i = 1.78. Strictly speaking, the values for δ computed with the uncorrected ice thicknesses must also be corrected. With δ m denoting the values computed with the uncorrected ice thicknesses, we have
where the correction Δδ is given through Equation (4) by
From the last three equations we find
That is, in the northwestern part of Nivlisen, ΔZ is about 7 m. Corrected values for Z and δ are used in the further analyses as well as in Figures 8 and 9.
Now, by the parameter δ, we know the spatially varying relationship between ice thickness and ice surface height and can use this knowledge to transform heights into ice thicknesses and vice versa, thus improving our information on both.
Thus, an ice-thickness model (Fig. 10) was generated merging direct observations and height data at floating positions converted to ice thicknesses (Equation (6)) using the established model of δ. The model domain was extended to the grounded ice east of Nivlisen, where a satisfyingly dense coverage with original thickness measurements is available. The resolution is again about 3.7 km.
For an error estimation of the obtained ice thicknesses, a simple error propagation from all inputs of Equation (6) would mean an overestimation because, by the genesis of the δ model, systematic errors of δ are correlated with errors in ρ w, ρ i and the absolute level of the freeboard correction (N + ζ SST + ζ Tides) such that they largely cancel out in Equation (6). Therefore, only the errors in the spatial variations of h and δ (see Tables 1 and 3) were propagated. The resulting estimated error in ice thicknesses derived from heights is 23 m in general and 27 m near the model margin and the ice rumple at about –70.4º S, 10.7° E. The latter value is also reasonable for the eastern ice-shelf margin where thicknesses were interpolated over larger distances. For a model of the solid ice thickness Z i = Z–δ the accuracy is practically the same as for Z.
Analogously, an ellipsoidal height model (Fig. 11) was generated from both direct height measurements and ice-thickness data converted to heights using Equation (5). In this way the model described in section 2.2.3 was extended, in particular at the eastern ice-shelf margin. At gridcells without original thickness observations, we used the ice-thickness model.
Concerning the errors of heights converted from thicknesses, an analogous argumentation applies as for thicknesses obtained from heights. Propagating only the errors of Zand the errors in the spatial variations of δ, we estimate the error of the obtained heights to be 3.8 m
3.3. Mass flux and mass balance
Figure 12 shows the horizontal ice mass flux on Nivlisen computed from the models of solid ice thickness and of horizontal flow velocity. Visual inspection suggests that the ice body loses mass while flowing northwards.
For a quantitative assessment of this loss, the net effect of specific mass balance and local mass balance, b s + b b – ρ i ∂Z i/∂t, was evaluated using Equation (9). Figure 13a shows the result. Its spatially varying accuracy is indicated by dashed brown contours.
For this accuracy estimation, we expressed the velocity components u and v of Equation (9) through the two original observables r (velocity in the radar line of sight) and α (flow direction; see section 2.3.3). Then we performed an error propagation from the errors in r, α and Z i (23 m, 2.6 ma–1 and 5°, respectively, according to sections 2.3.3 and 3.2) and from the errors in their spatial derivatives. Those errors were roughly estimated to be 2 m km–1 for ∂Z i/∂x, ∂Z i/∂y, 0.2 m a–1 km–1 for ∂r/∂x, ∂r/∂y, and 0.7° km–1 for ∂α/∂x, ∂α/∂y. The largest contributions to the error in b s + b b – ρ i ∂Z/∂t originate from the spatial derivatives of flow direction and of velocity.
The mass-balance magnitudes in Figure 13a are negative in most places, but show a variable pattern and have large uncertainties. To obtain more significant information, we performed a mass-flux evaluation (Equation (10)) over areas A-C shown in Figure 13b. Table 4 lists the results. Accounting for the known surface mass balance (section 2.4), we find that, on average over the three areas, the combined effect of basal mass balance and temporal change is –654 ± 170 kg m–2 a–1.
A rigorous error propagation for the integration in Equation (10) would require spatial correlation information on the error of the integrand |ρ i Z i u|sin β. Our error assessment for each area assumes uncorrelated errors of the fluxes through its four boundary sections: upstream, downstream and the two lateral flowline boundaries where zero cross-flux is expected but subject to flow direction errors (‘lateral in- and out-flux’ in Table 4). Along each section, we compute the cross-flux error as the integral of the one-sigma error thus allowing for full error correlation within the section. equals for flowline sections and for the perpendicular flux gates.) Area C covers part of the flow across the grounding line. Surface velocities of the grounded ice are known with 5ma–1 accuracy from the interferometric analysis of SAR observations from the ERS-1/-2 tandem mission (Reference Baessler, Dietrich, Shum, Jekeli and ShumBaessler and others, 2003). Here, the depth dependence introduces an additional uncertainty in the depth-averaged velocity. Since likely values are 80–100% of the surface velocity (Reference PatersonPaterson, 1994, p.252) we applied 90% ± 10%.
4. Interpretation and Discussion
4.1. Firn structure
In section 1.2 we explained that the magnitude δ, derived from Equation (4) at free-floating positions, represents the amount of air contained in the ice shelf (the ‘air layer thickness’) if, first, the ice shelf consists of a mixture of pure ice and air only, and second, there are no gross errors in the observations used. There is no indication of a violation of these assumptions. On the contrary, the values of δ and their spatial variations as shown in Figure 9 are in remarkable agreement with external information.
For comparison, we estimate δ from density information via Equation (3). In the bare ice region, the upper 15 m contribute only 0.6 m to S, based on the 879 kg m–3 mean density cited in section 2.4. For the density below 15 m we apply two different models: a density profile of the Ross Ice Shelf, Antarctica, (Reference GowGow, 1963) starting from 879 kg m–3, and an exponential densification analogous to Equation (14): with Z 0 = 15 m, ρ 0 = 879 kg m–3 and C = 0.029 as obtained in section 1.2. The resulting total air layer thicknesses are 3.6 and 2.0 m, respectively, and compare well with our model.
The increase of δ to the north corresponds to the reported northward-growing firn layer (section 1.1). Moreover, the isolines of δ appear parallel to the blue-ice area boundary as seen in satellite imagery (e.g. Figs 1 and 2). Our δ values north of the increasing zone are in the 10–20 m range reported by other authors (cf. section 1.2) and they compare with the 18.7 m value derived from grounded ice-core data from the region (section 2.4).
We conclude that δ can be interpreted as air layer thickness to reflect the variations of firn conditions over Nivlisen.
4.2. Mass balance and basal processes
Although the spatial distribution of the mass balance calculated in section 3.3 (Fig. 13a) is uncertain so that its interpretation would be speculative, we have found that on average over a 2400 km2 subregion of Nivlisen (34% of its total area) the net effect of basal mass balance and temporal ice mass change is –654 ± 170 kg m–2 a–1 or, equivalently, millimetres water equivalent per year. Under the steady-state assumption ∂Z i/∂t = 0, this would mean average basal melting of 654 kg m–2 a–1. Basal melt rates at other Antarctic ice shelves reach several metres per year (Reference RignotRignot, 2002), so our value is within the expected range. An additional temporal change would change the basal melt rate. For example, a thinning by 100 kgm–2 a–1 would enhance the melt rate by another 100kgm–2 a–1, while a thickening would accordingly reduce it.
A dominance of bottom melting appears likely from our results. Melting is also reported from the two Soviet ice cores on Nivlisen (see Fig. 2). Texture analysis of the basal ice at drilling site 2 showed that no bottom freezing took place and that the ice was of meteoric origin (Reference Korotkevich, Savatyugin and MorevKorotkevich and others, 1978). For site 3 also, texture and isotope analyses indicated lasting thawing processes (Reference Gordienko and SavatyuginGordienko and Savatyugin, 1980; Borman and Fritzsche, 1995, ch.6). On the other hand, basal processes may vary considerably over small scales. It has been estimated (Borman and Fritzsche, 1995, ch. 4) that 13.4 × 106m3a–1 of summer meltwater runs off into the lakes in Schirmacheroasen which are linked to the ocean water beneath the ice shelf. A similar meltwater entry may occur at the tidal fractures along the grounding zone. This fresh water, being lighter than sea water, is likely to refreeze quickly on the ice-shelf bottom. Bottom melting in other zones may then overcompensate this local freezing and lead to the negative net effect.
Independent observations also indicate a possible ice thinning with time. The bare ice area on its grounded part next to the ice shelf showed a negative mass imbalance of –100 to –200 kg m–2 a–1 in the years 1991–98 (Reference Korth, Perlt, Dietrich and DietrichKorth and others, 2000). In Schirmacheroasen, a decrease in snow-patches on a decadal scale has been observed (Reference Bormann and FritzscheBormann and Fritzsche, 1995, ch.6). Extrapolation of these signals to the ice shelf, with its peculiar regimes of ice dynamics, meteorology and ice–ocean interaction, is of course not justified. On the other hand, ice shelves may respond even more rapidly to mass imbalances than grounded ice (Reference VaughanVan der Veen, 1986). If temporal changes in the ice-shelf thickness occurred, their implications for the overall evolution of the ice shelf would be an interesting question. For example, Reference Shepherd, Wingham, Payne and SkvarcaShepherd and others (2003) observed a progressive thinning of Larsen Ice Shelf, Antarctica, prior to its partial collapse. A comparison of the present ice-thickness data to historical observations is tempting but very questionable (see section 4.3).
In conclusion, no quantification can be given from the present data on how bottom melting and temporal ice-mass changes contribute to their negative net effect. A dominance of bottom freezing, however, seems unlikely, at least for central and eastern Nivlisen.
4.3. Discrepancies with historical ice-thickness observations
In 1975/76, Soviet Antarctic Expedition groups performed three ice-core drillings (see Fig. 2) to the bottom of the ice (Reference Korotkevich, Savatyugin and MorevKorotkevich and others, 1978). The ice thickness reported at core 1 (374 m) agrees well with our nearby radar measurements, but the reported thicknesses at cores 2 (357 m) and 3 (447 m) are 65 and 61 m larger, respectively, than in our ice-thickness model. This discrepancy is striking. On the one hand, drilling ought to be a direct and accurate thickness observation, but on the other, the differences far exceed the uncertainties estimated for our data.
Geophysical variations alone are an unlikely explanation. Temporal variations of >60m in 20 years would be very large and are not supported by comparison of our data to radar observations during the ice drilling (Reference Kozlovskii and FedorovKozlovskii and Fedorov, 1983). Spatial variations at smaller scales than our model resolution (as assessed from variations along our radar tracks) rarely exceed 20 m in the vicinity of core 3 and are much smaller near core 2.
While we have tried to thoroughly assess the uncertainties in our data, it is difficult to assess the reliability of the historical ice-drilling observations. Tilted boreholes may have led to core lengths exceeding the ice thickness (e.g. a 65 m excess length at core 2 could result from a 35° tilt). Some discrepancies within the reports on the 1975/76 ice drilling and radar activities remain unexplained. Drilling-site positions reported by Reference Korotkevich, Savatyugin and MorevKorotkevich and others (1978) and by Reference Eskin and BoyarskiiEskin and Boyarskii (1985) differ by about 6 and 4 km for cores 1 and 3, respectively. These two publications agree on the position of core 2, but Reference Kozlovskii and FedorovKozlovskii and Fedorov (1983), in their ice-thickness map derived from radar observations, mark it about 14 km away. Otherwise they would have revealed a difference in the order of –100m between their mapped thickness and the drilled thickness.
Accepting, nonetheless, the reported ice thicknesses at drilling sites 2 and 3 leads us to two questions which we cannot plausibly answer: how are those large ice thicknesses compatible with our observed freeboard heights in view of the hydrostatic balance? and what is the horizon observed by various radar observations? A hypothesis on a marine ice bottom layer does not answer either question. First, the meteoric origin of the bottom ice was reported from ice-core analyses (see section 4.2). Second, with a marine ice layer of, say, 65 m (corresponding to the difference between ice-core and radar observations), the resulting higher total ice thickness would be incompatible with the observed freeboard height because the marine ice reduces δ only slightly (cf. Equation (8)) and, accordingly, only slightly changes the hydrostatic thickness-to-height relation (Equation (5)). Denser material transported by the ice (e.g. sediments) could resolve the problem of the hydrostatic balance but would not explain the internal horizon seen by radar observations.
We conclude that if no gross error is responsible for the discrepancy described, it is most likely a peculiar coincidence of several effects such as borehole tilts, position errors, errors in our ice-thickness model, small-scale variations not included in this model, and temporal changes. In any case, we are sceptical about including the reported borehole thicknesses in a geophysical interpretation.
4.4. Snow accumulation and compaction
The evolution of the air layer thickness δ along the ice flow can be used to constrain snow accumulation and firn compaction parameters. Air contained in the ice shelf is added in accumulating snow. For example, with the accumulation rate b s = 300 kgm–2 a–1 (section 2.4) and a typical fresh snow density ρ s = 300 kg m–3, the air layer thickness added by snow accumulation can be derived from Equation (3) to be b s (ρ —1 – ρ i –1) = 0.67 ma–1. On the other hand, contained air escapes from the ice shelf when snow and firn are compacted, either by the continuous transformation of snow to ice or by melting with subsequent refreezing. Much of the compaction occurs shortly after accumulation. For example, Reference Korth and DietrichKorth and Dietrich (1996) use a 600 kgm–3 mean density of the previous year’s accumulation layer. With b s and ρ s as above, this implies that at the end of a year the accumulated snow has already been compacted by 0.5 m and retains no more than 0.17 m of air.
In analogy to the ice mass-flux and -balance methodology (section 1.3), the net effect of contained air gain and loss is described by the divergence of the horizontal flux of contained air, div(δ u). We computed the integral effect of this divergence for three areas by balancing the in- and outflow of contained air through their boundaries. We chose areas I–III marked in Figure 13b since they represent different glacio-logical regimes (cf. section 4.1). The obtained average divergence of contained air was 0.01, 0.04 and 0.02 ma–1, respectively. That is, the effect of accumulation slightly predominates over the effect of compaction, and this dominance is largest in the central area II. The compaction effect in area II is in the order of 0.04–0.67 = –0.63 ma–1, or, excluding the first-year compaction, 0.04–0.17 = –0.13ma–1.
We speculate that the obtained qualitative pattern, irrespective of its uncertainties, reveals a geophysical signal induced by spatial variations of accumulation and compaction. A higher accumulation rate in area II could result from drift snow transported by katabatic winds from further south across the ablation area. A more rapid compaction in areas I and III could result from enhanced melting and refreezing in the vicinity of the ablation area (area I) and the coast (area III). In addition, area III has a thicker firn body with a larger air content, so there is more firn to be compacted.
4.5. Ice-dynamical features
The maps of flow velocity, ice thickness and ice mass flux (Figs 7, 10 and 12) show clearly how the ice flow across Nivlisen is channelled between areas of local grounding or, conversely, how these grounding areas form obstacles to the ice flow, leading to lower ice thicknesses in their lee. This underlines the complexity of the ice-dynamical regime of Nivlisen which may be sensitive to changing grounding conditions induced by changes in either ice thickness or sea level.
5. Conclusions
Among the ice shelves along Antarctica’s Atlantic coast, Nivlisen is remarkable for having a peculiar glaciological regime, a long history of exploration and a wealth of observations. The compilation and combined analysis of these data yield a consistent geometrical and glaciological description of Nivlisen that is more comprehensive, accurate and detailed than previous descriptions. Among other things, we constrained the spatially varying firn conditions by utilizing the hydrostatic balance relation of floating ice, and we estimated from mass-flux considerations that, under a steady-state assumption, basal melting with an average rate of 654± 170mma–1 w.e. occurs in a large part of the ice shelf.
Our results may be valuable for ice-dynamics modelling efforts like those of Reference Paschke and LangePaschke and Lange (2003) either for validating models or for providing new and improved constraints.
Our study may also serve as a starting point for investigating temporal variations. Whether and how the negative local mass balance observed upstream on the grounded ice continues on the ice shelf remains an open question relevant for understanding the dynamics of the whole Nivlisen drainage system. Additional observations and analyses could address this question. Apart from the possibility of repeated ice-thickness measurements subject to the logistic constraints of Antarctic fieldwork, a further exploitation of satellite altimetry should certainly be considered. Repeated observations will soon be available from a whole sequence of altimeter missions – ERS-1, ERS-2, Envisat, ICESat and, from 2005, CryoSat – spanning about 15 years. Advanced methods of radar altimeter analysis may use the full along-track resolution and include several altimetric waveform parameters, thus facilitating the separation of real surface height changes from changes in firn conditions (Reference Legrésy, Rémy and VincentLegrésy and others, 2000). Repeated GPS measurements could serve as ground truth. In such analyses of originally ellipsoidal height data, exact knowledge of vertical displacements due to ocean tides and the inverse barometric effect is of great value. New information on bottom processes could be gained from more sophisticated radar applications (Reference Corr, Jenkins, Nicholls and DoakeCorr and others, 2002). An extension of mass-flux and mass-balance observations at western Nivlisen, including a mass-flux estimation across the southwestern grounding line, is another open issue.
In view of the extensive knowledge already available, and the favourable logistic conditions, it seems worth pursuing the open questions on Nivlisen and its drainage basin as a case study of the complex interactions of land ice, bedrock, ice shelf, ocean and atmosphere.
Acknowledgements
The main parts of this research were supported by the German Research Foundation (DFG). We acknowledge B. Legrésy and G. Vinay for supplying processed altimeter data. We also thank B. Legrésy for valuable discussions on radar altimeter data analysis. C. Walter is acknowledged for the elaboration of initial studies about the hydrostatic equilibrium on Nivlisen performed for her diploma thesis. Comments by B. Breuer, H. Rott and A. Wendt helped to improve the manuscript. CHAMP and GRACE geopotential models were provided by the GeoForschungsZentrum Potsdam.