Introduction
For the present study, we focused on Nivlisen, an East Antarctic ice shelf, and its drainage area. This region comprises four sub-regions, each characterized by a different ice geometry and consequently a unique flow regime. Extended inland ice-sheet parts; laterally small, but deep, outlet glaciers located between high, nearly ice-free mountain tops; meandering ice streams; blue-ice areas with non-negligible mean annual ablation rates; and an ice-shelf area have to be taken into consideration. The sum and the interaction of all these different features represent a challenging task for the adaptation of the three-dimensional numerical flow model developed by Reference K., Sandhäger and LangeSandhäger (2000) and previously applied to a more homogeneous flow regime.
The numerical model needs as initial and boundary condition reliable gridded ice-geometry data. Up to now, only less detailed information has been available for the study region, based on larger-scale datasets for the entire Antarctic ice sheet, a few ground-based in situ measurements and limited coverage by remote sensing. Therefore, we combined several available datasets to obtain the respective information on a 625 m grid. This horizontal resolution enables a numerical simulation of the outlet glaciers in the mountainous region and other small-scale features like ice rises and ice rumples within the ice-shelf part. Whereas the exact derivation of the geometrical models is described in Reference PaschkePaschke (2002), we will focus in this paper on the numerical flow model, its necessary modifications and the simulated results.
The Study Region
The region under consideration, Nivlisen and its drainage area, is located in East Antarctica (69–75° S, 9–14° E; Fig. 1) and represents a relatively small grounded-ice/ice-shelf system. Following Reference Bormann and FritzscheBormann and Fritzsche (1995), this system can be subdivided as follows:
Fig. 1. Shaded relief map of the digital elevation model derived in the present study (illumination from northwest). Contours for 50 m intervals (up to 200 m), 100 m intervals (up to 1000 m) and 200 m intervals (above 1000 m) are added. The position of the grounding line is marked. In addition, the four sub-regions (A–D) are indicated as described in the text. The positions of the two profile lines (e.g. Fig. 2) are given, as well as the horizontal spacing of the gridpoints for the numerical simulation.
Fig. 2. Vertical cross-sections of the ice body from south north (top) and west to east (bottom). The positions of profiles are indicated in Figure 1.
The ice shelf Nivlisen (A, Fig. 1) represents a relatively flat region with surface elevations mostly below 50 m. Some research stations (e.g. the Forster station of the former German Democratic Republic at 70°46’ S, 11°51’E) which are located at the southern edge of the ice shelf serve as bases for numerous airborne and ground-based glaciological investigations. In addition, satellite data are available for the entire study region, which provide insight into surface topography and ice flow.
The region of increasing surface elevations between the ice shelf and the Wohlthat Massif (B) is mainly characterized by a relatively smooth surface topography of the ice body, with a snow cover textured by sastrugi. A notable feature of this region is the presence of nunataks with elevations of up to 250 m above the ice surface. Their presence has distinct effects on the ice flow, especially on the northward-meandering main glacier.
The Wohlthat Massif (C) is part of the alpine mountain range in Dronning Maud Land and reaches altitudes of 1000 to > 2000 m a.s.l. The outlet glaciers are located in between these high, nearly ice-free mountain tops. This situation leads to a complex flow regime which is discussed in detail further below.
The southernmost portion of the drainage area (D) is part of the Wegener Inland Ice. The maximum ice thickness observed in the study region (43300m) is found here.
To quantify the ice flow and mass balance of such a complex morphology by means of numerical modelling, a detailed description of ice and bedrock geometry is needed. Since no high-resolution digital elevation and ice-thickness models at intervals of 5 2 km spacing are available for this region, we combined several datasets of different horizontal resolution by utilizing a suitable gridding algorithm (kriging). Through this a new database has been constructed which finally represents the aforementioned four zones on a 625 m grid. The exact derivation of the datasets has been described in Reference PaschkePaschke (2002).
The surface elevation grid (Fig.1) is based on the contour lines of the Antarctic Digital Database (ADD) published by the British Antarctic Survey (1998), the 5 km grid provided by the BEDMAP Consortium and published by Reference Liu, Jezek and LiLiu and others (1999) and some unpublished airborne radar-altimetry measurements obtained during the GeoMaud-Expedition 1995/96 by the Federal Institute of Geosciences and Natural Resources (BGR) in Hannover, Germany (personal communication from V. Damm, 2001). The surface elevation over the ocean was set to be zero regardless of the presence or absence of sea ice, as this is not treated in the numerical flow model.
The results obtained agree well with the flow structures identifiable on the RADARSAT mosaic (Reference JezekJezek, 1999). Furthermore, a comparison between our digital elevation model (DEM) and the dataset published by Reference Bamber and HuybrechtsBamber and Huybrechts (1996) reveals a high level of consistency. Eighty per cent of the points within the ice-shelf area show a mismatch of 510 m. Regarding the locally varying thickness of the firn layer, discrepancies of up to 10 m can be explained by, for example, the use of horizontal smoothing and/or neglected reduction of positioning uncertainties. Within the grounded-ice region, 73% of the points differ 5100 m. Those points with >100 m mismatch are located in the alpine region of the Wohlthat Massif which is partly covered by radar-altimetry measurements and is therefore possibly better represented than in the 10 km grid by Reference Bamber and HuybrechtsBamber and Huybrechts (1996).
Concerning the ice thickness, there is only sparse information available for the study region. We used the data providedby the BEDMAP Consortiumand published by Reference Lythe and VaughanandLythe and others (2001), unpublished radio-echo sounding (RES) data from BGR Hannover obtained concurrently with the surface elevation data mentioned above (personal communication from V. Damm, 2001) and RES data from the Alfred Wegener Institute (AWI) in Bremerhaven, Germany, obtained as part of the European Project for Ice Coring in Antarctica (EPICA) 1995/96 and published by Reference Sclater, Jaupart and GalsonSteinhage and others (1999). These two airborne datasets are characterized by a high horizontal resolution along the flight-lines, with a focus onto the region between the Wohlthat Massif
and the southern edge of Nivlisen. However, south of the alpine mountain belt, data of less than 10 flights are available for the study area. This leads to an ice-thickness model with high accuracy north of the Wohlthat Massif and significant uncertainties south of it.
The subglacial bedrock topography for the grounded-ice also region is calculated based on the surface elevations and ice used thicknesses. To specify the positions of the grounding line, ice rumples and ice rises, we employed the ADD in combination with a personal communication from J. Perlt (2001) and the RADARSAT mosaic published by Reference JezekJezek (1999). The identifiable flowlines and the derived surface topography model serve to define the boundaries of the drainage area.
To illustrate the aforementioned complicated ice geometry of the study region, two profile lines in the north–south and east–west directions, respectively, are shown in Figure 2. One easily recognizes the complex situation: ice masses from the thick Wegener Inland Ice have to pass the Wohlthat Massif through laterally small, but thick, outlet glaciers in the middle of high, nearly ice-free mountain tops before meandering towards Nivlisen. A steep gradient in ice thickness is visible in the transition zone from the grounded-ice part to Nivlisen, which is also represented in the dataset published by Reference Bamber and HuybrechtsBamber and Huybrechts (1996) as well as in some RES flight-lines.
Based on the ice-thickness model, we calculated several basic parameters for the ice body in the study region (Table 1). Although Nivlisen is a relatively small ice shelf whose drainage area comprises about 0.43% of the East Antarctic ice sheet (Reference Drewry, Jordan and JankowskiDrewry and others, 1982), the complexity of its flow regime is comparable with more extended regions of the Antarctic ice sheet.
Table 1. Ice volume (Vice), ice-covered area (Aice), mean ice thickness (Hice) and ice mass (mice) for the specified subregions of the study area (labels as in Fig. 1) as derived from the ice-thickness model
The Numerical Simulation
The three-dimensional, numerical flow model
The major tool used to quantitatively describe the ice dynamics and mass balance in our study region is a three-dimensional numerical flow model developed by Reference K., Sandhäger and LangeSandhäger (2000). The structure of the model is schematically depicted as a flow chart in Figure 3. It is mainly based on two approaches for the grounded and the floating ice portions, respectively.
Fig. 3. Flow chart of the numerical flow model (Sandhaäger, 2000) used for the ice-dynamics simulation.
Within the floating ice shelf, the main assumptions are that the strain rates are independent of depth and that the shear strain rates in a vertical plane are zero. This leads to the commonly used ice-shelf equations (e.g. MacAyeal and others, 1986). Reference K., Sandhäger and LangeSandhäger (2000) solved the unscaled version of these equations using the relaxation approach proposed by Reference Herterich, Van der Veen and OerlemansHerterich (1987). The ice-flow velocities are thereby calculated, depending on the simulated temperature distribution. The latter is derived by solving the steady-state heat equation under the assumption that the heat production results mainly from internal ice deformation (Reference PatersonPaterson, 1994).
This approach is also utilized for the temperature distribution within the grounded-ice part of our model domain.
where m and n are adjustable parameters, R is the gas constant and Q the activation energy for creep. The following values are chosen for our study: m=1, n=3, A0= 1.14 × 1 0 –5 Pa– 3 a– 1 and Q = 60 kJ mol –1 if the absolute temperature corrected for the dependence on the pressure-melting point satisfies T* < 263.15 K. If T * > 263.15 K then A0 = 5.47 6 10 10 Pa– 3 a– 1 and Q = 139 kJ mol –1 are used (Reference PatersonPaterson, 1994).
The second and the third component both describe basal motion which can occur if the temperature at the base of a grounded-ice body reaches pressure-melting-point conditions. These two processes are presently not completely understood and depend on several parameters related to the bedrock and the substrate. Thus, they are hard to distinguish. An empirical expression describing the relation between the horizontal flow velocity 
the basal vertical shear stress 
and the effective pressure pe
ff has been proposed by Reference BindschadlerBindschadler (1983) who modified a Weertman-type sliding law for the effect of subglacial water pressure. Reference K., Sandhäger and LangeSandhäger (2000) integrates this relation together with the approach of, for example, Reference HuybrechtsHuybrechts (1992) which implies that the effective pressure can be replaced by the height above buoyancy Z (with Z ≥1m for the grounded-ice region):
The basal flow parameter Cb depends on the inverse bedrock roughness and on the type of underlying sediments. The scarcity of sufficient field data may imply the use of constant value(s) for Cb, despite the obvious simplification of such an approach given the true nature of the glacier bed. Alternatively, if surface velocities are available, they may be used to infer spatially dependent values of Cb. Initially, the values of the three parameters are chosen in the same range as found by Reference K., Sandhäger and LangeSandhäger (2000) for Ekstromisen: Cb = 2.625×10–7 m7/4 Pa–2 a–1, p = 2 and q = 0.75.
Since Equations (1) and (2) are derived empirically, further investigations have to be undertaken to establish suitable parameters. In the present study, we focus on the basal motion, namely on the basal flow parameter Cb and its aforementioned possible spatial variation.
Application to the study region
The ice-geometry models serve as initial and boundary condition for the numerical model. However, due to limited computer power, we were forced to minimize the computational requirements of the model. This constraint together with the desire to maintain the necessary accuracy led to a Cartesian grid with a horizontal resolution of 1.25 km (Fig. 1) and 13 normalized depth layers such that the glaciers transporting the ice from the Wegener Inland Ice through the Wohlthat Massif to the northern part of the study region are represented with at least approximately 10 gridpoints in the east–west direction. Furthermore, we exclude the southernmost part of the study region in our model runs to save computational resources. This can be justified because this part appears to be characterized by minor and/or homogeneous lateral ice flow and has thus almost negligible influence on the numerical results. However, the diversion into separate glacier systems directly south of the mountain belt (RADARSAT mosaic (Reference JezekJezek, 1999)) is maintained and calculated using the present grid. The output of the diagnostic simulation described in the present paper includes ice velocities, internal stress and temperature distributions.
A major challenge addressed in our study lies in the difficulty of correctly describing the ice dynamics of a complicated flow regime through a numerical model initially tested for a less complex situation. Up to now, the flow model was used by Reference K., Sandhäger and LangeSandhäger (2000) to analyze the flow regime of Ekstromisen and its drainage system, an ice-shelf/ice-sheet region in the western part of Dronning Maud Land.
Applying the model to Nivlisen, we found that, using the same boundary conditions and parameterization, the calculated magnitudes of the vertically averaged ice velocity, especially within the alpine mountain belt, were too small by almost one order of magnitude (Fig. 4, dotted line). We first tried to improve the results by varying some initial parameters, such as the geothermal heat flux (interval [0.050;0.065] W m–2 (Reference SandhägerSclater and others, 1980)), but without any positive influence on the simulated velocity field. This changed when we introduced two principal modifications describing the basal motion of glaciers:
Fig. 4. Horizontal flow velocity (vertically averaged) across the west–east transect presented in Figure 2. The dotted line represents the model run without any modifications; the dash-dotted line depicts the model run where pressure-melting conditions at the base of the glacier systems have been implemented. The solid line shows the effect of enhanced and reduced basal flow parameter for central and marginal parts of the glacier systems, respectively.
Due to the non-linear character of the model equations, different mathematically correct solutions, highly sensitive to initial and boundary conditions, are possible. Since the simulations using the aforementioned adjustable parameters did not succeed, we introduced an additional constraint to obtain a glaciologically reasonable solution. We used the RADARSAT mosaic (Reference JezekJezek, 1999) to identify the positions of the main glacier systems (i.e. glaciers and ice streams) and set the basal temperature for one central flowline per glacier to the pressure-melting point, which is essential for basal sliding. The steady-state solution thus derived provides the initial condition for time-dependent and self-consistent simulations in the future. This modification improved the resulting velocity field significantly (Fig. 4, dash-dotted line), but still not sufficiently when compared to expected values.
Therefore, we focused on the basal sliding law (Equation (2)), for which no general theories are available (Reference Bamber and HuybrechtsBamber and Huybrechts, 1996). Reference HarborHarbor (1992) examined the effects of varying parameters of the basal sliding law and of Glen’s flow law (Equation (1)). He found that for a glacier cross-section, best results were obtained with scaled values for the bedrock roughness using lowest values for the centre line and increasing ones towards the margins. This assumption is reasonable since relatively high concentrations of basal and englacial debris from subglacial and subaerial sources are observed there. In the present study, we utilized the suggestion of Reference HarborHarbor (1992) and enlarged the basal flow parameter Cb, which is inversely proportional to the bedrock roughness, for the aforementioned central flowline of each glacier by a factor of two and halved its value for the remaining grounded-ice area affected by basal sliding. The additional improvement of the resulting velocity distribution, especially within the Wohlthat Massif, as a result of this procedure is also shown in Figure 4 (solid line).
Using the two modifications, the magnitudes of the horizontal velocity are five times larger than those obtained with the constant basal flow parameter and without predetermined basal temperatures for the central flowlines of the glacier systems. Unfortunately, no comparison with measurements is possible within this region due to the lack of field data. Nevertheless, the modelled horizontal ice-velocity values within the glacier systems are 530 m a–1.This is certainly still too small since these glacier systems are identifiable as outlet glaciers based on the RADARSAT mosaic (Jezek, 1999).
The explanation for this mismatch probably lies in a combination of at least two different aspects. First, the vertically averaged dynamic velocities depend on the ice thickness to the fourth power, i.e. a change of 10% in ice thickness can influence the ice velocity up to 46% (Reference Bamber and HuybrechtsBamber and Huybrechts, 1996). Since the alpine mountain region is characterized by a complicated ice geometry and due to the rather sparse measurements south of the Wohlthat Massif, there is still an uncertainty left in the accuracy of the derived ice-thickness model.
Another possible difficulty is due to the present situation where fast-moving ice of glacier systems lies next to the more stagnant parts of the inland ice. This results in a discontinuity in slip resistance between these two units. Reference RaymondRaymond (1996) notes that differentially moving portions of an ice sheet also cause an equally discontinuous stress distribution at their interface. This situation can be compared to the flow regime in the vicinity of ice rumples and ice rises within an ice shelf. In order to adequately represent this situation in a numerical ice-shelf model, Reference RaymondSaheicha and others (in press) propose a mechanism that decouples the horizontal velocity fields of that part of the ice body belonging to the free-floating ice shelf from the one moving across the ice rise. Such a treatment would cause effects comparable to those obtained by spatial variation of the enhancement factor m (Equation (1) A similar approach in our study may allow the coexistence frozen-on parts of the ice sheet next to parts belonging to the fast-moving glacier systems without having the latter slowed down by excessive horizontal stresses at their interface. However, it remains to be verified if such a decoupling approach will indeed minimize the current mismatch between the calculated ice velocities within the Wohlthat Massif and the more typical flow velocities for glacier systems.
Results of the numerical flow model
Adapting the three-dimensional numerical flow model to our study region and using the aforementioned modifications, we obtained results for the ice velocities, internal stress and temperature distributions. However, in this paper we will only analyze the velocity field.
Figure 5 shows the calculated horizontal ice velocities for the entire model domain. The spatial velocity distribution appears to be glaciologically reasonable: The highest values occur along the glacier systems. Within ice-sheet areas, the ice is characterized by low velocities. Furthermore, one easily recognizes the increase of the inland ice velocities towards the grounding line. The flow regime within the ice shelf is dominated by the effect of the ice rises and ice rumples and shows the expected pattern, with velocities increasing strongly towards the ice-shelf edge. The directions of the horizontal ice flow are indicated by unscaled arrows and agree closely with observed flowlines (e.g. on the RADARSAT mosaic from Reference JezekJezek (1999)).
Fig. 5. Ice-flow velocities (vertically averaged) for the entire model domain as derived from the numerical flow model. The grey shading indicates velocity magnitudes, and the arrows denote the flow direction.
To validate our model results with observations or field measurements, the only available dataset within the study region comprises velocity measurements located between the Wohlthat Massif and the southern edge of Nivlisen (Reference Korth and DietrichKorth and Dietrich, 1996) (Fig. 6b). Comparing these values with the modelled surface velocities in the same positions (Fig. 6a), we note a broad consistency regarding the flow directions, indicated by the arrows. Concerning the velocity magnitude, i.e. the scaled length of the arrows, the a modelled results are on average about 35% too small. This is). certainly a kind of after-effect of the aforementioned low of velocities within the mountainous region. To tackle this problem, we propose to utilize laterally varying boundary conditions (e.g. with respect to the basal flow parameter). Moreover, in order to adequately describe the complex geometry in the area of the Wohlthat Massif, a nested grid seems advisable. This would enable higher spatial resolution at specific points of the study area while maintaining reasonable computer times.
Fig. 6. Flow directions and the horizontal speed at the surface derived from the numerical model (a) compared with in situ measurements published by Reference Korth and DietrichKorth and Dietrich (1996) (b).
Regarding the depth dependence of the speed within the grounded-ice region, one has to distinguish between glacier systems and ice-sheet regions: due to the cold base, the horizontal velocity of the latter increases from low values at the bottom where the ice is frozen to the bedrock to higher values at the top of the ice sheet. But these magnitudes are still significantly smaller than those of glacier systems, which, in contrast, are dominated by low vertical velocity gradients. This can be seen in Figure 7 which shows the lateral distribution of the local mass flux across the grounding line in a vertical transect. One easily recognizes the difference between so-called active grounding lines (glacier-systems–ice-shelf) characterized by a high local mass flux and passive grounding lines (ice-sheet–ice-shelf) with a low transition of ice mass.
Fig. 7. Top: Local mass flux across the grounding line. Bottom: Vertical distribution of the horizontal ice velocity across the grounding line of Nivlisen.
Mass balance
Comparing mass-balance values calculated based on the numerical model with those obtained using a map of net accumulation rate for the whole of Antarctica (Reference Giovinetto and BentleyGiovinetto and Bentley, 1985) represents a rather crude approach. The following brief discussion explains several occurring uncertainties in detail and should therefore be considered more as a classification of the mass-balance quantities than as a comparison or validation. Furthermore, the arguments given below justify large error tolerances with respective signs.
As mentioned above, the modelled magnitudes of the ice velocity, especially in the mountainous region, are presumably still too small by a factor of two, which also has impacts on the ice-flow regime in the northern part of the study area. Since the simulated horizontal velocities are taken to calculate the mass flow across the grounding line, the explained after-effect leads to a positive error assumption ((1800+ 1000) Mt a–1).
Contrariwise, the map of the net accumulation rate published by Reference Giovinetto and BentleyGiovinetto and Bentley (1985) also bears a large inaccuracy considering a small-scale region. First, the isolines are given in intervals of 50 kg m–2 a–1 such that an uncertainty of 50% occurs in most cases. The second problem due to the large-scale map is that smaller ablation regions like blue-ice areas might be neglected, which seems reasonable considering the entire Antarctic ice sheet. But for a drainage area the size of our study region, this would lead to higher values of accumulation than those observed locally. This is certainly the case for the considered region since several blue-ice areas are observed between the Wohlthat Massif and the southern part of Nivlisen (Reference Bormann and FritzscheBormann and Fritzsche, 1995) such that the error of the net surface accumulation rate within the drainage system is assumed to be negative ((4500–2000) Mt a–1). Keeping this in mind, the mass flow across the grounding line and the accumulation rate within the grounded-ice area agree within the error tolerances.
Within the relatively small ice-shelf area, the accumulation rate from Reference Giovinetto and BentleyGiovinetto and Bentley (1985) served as input data contributing to the vertical velocity component at the surface. Thus, the expected agreement between modelled values and those calculated from the map (net surface accumulation rate of (1450±600) Mt a–1) is obtained. Integrating the continuity equation for the steady-state mass flux vertically, the basal melting rate is derived which comprises 20 cm a–1 for Nivlisen. This leads to a mass loss of –1970 Mt a–1 at the bottom of the ice shelf. The total mass balance of Nivlisen can be calculated (210Mt a–1) based on the mass contribution across the grounding line and due to the net surface accumulation, on the one hand. On the other hand, the mass loss at the ice front derived from the simulated horizontal velocity (–1070Mt a–1) and due to basal melting is considered. Given the large aforementioned uncertainties, a contribution of 210 Mt a–1 appears almost negligible. Thus, it can be concluded that Nivlisen is currently close to steady-state conditions.
Conclusion
We present basic ice dynamics for Nivlisen and its drainage area, East Antarctica, resulting from a three-dimensional numerical flow model. This region was chosen as it is characterized by a complex ice geometry and therefore a complicated ice-flow regime. Four sub-regions can be distinguished which all have their own flow characteristics. We adapted an existing numerical flow model (Reference K., Sandhäger and LangeSandhäger, 2000) to the study area, introducing some modifications in order to describe these four sub-regions collectively using only one numerical model.
Due to the non-linear character of the higher-order model equations, it is reasonable to use an additional constraint for the diagnostic simulation. It is shown that setting one central flowline of each glacier system to pressure-melt-ing-pointconditions, a steady-state solution is reached which agrees better with the glaciological situation than the results obtained without this modification. Moreover, comparisons between the numerical model results and measured ice velocities thus derived suggest the introduction of spatially varying boundary conditions in contrast to the usually assumed uniform characteristics (e.g. with regard to the basal flow parameter). Implementing these modifications in the model that has been derived by Reference K., Sandhäger and LangeSandhäger (2000) results in satisfactory agreement between computed and measured ice velocities. The still existing mismatch is mainly due to the complicated flow regime, especially within the alpine mountain belt. This also has obvious impacts on the mass balance. Nevertheless, the present study can be considered a significant step on the way to simulating regional effects such as deep outlet glaciers in the middle of nearly ice-free mountain tops while treating an entire drainage system. In this context, the usefulness of spatially varying boundary conditions such as the basal flow parameter has been shown.
In future studies, the employment of a nested grid for specific parts of the study region and the decoupling of the velocity fields of fast-flowing glacier systems and near-stagnant portions of the ice sheet may yield further improvement of the numerical results.
Acknowledgements
Special thanks are due to H. Sandhäger. He not only provided the numerical model that was used in this study, but also gave invaluable advice and encouragement throughout the course of this project. We are grateful to K. Grosfeld and K. Saheicha for many constructive and helpful comments and suggestions. We thank V. Damm of BGR Hannover who provided unpublished datasets of ice thickness and surface elevation for the ice-geometry models used. We are also grateful to J. Perlt for contributions to the geometric models. Last but not least, we wish to thank R. C.Warner, C. B. Ritz and an anonymous reviewer for valuable comments on an earlier version of this paper.
