Introduction
The cryosphere plays a distinctive role in the global climate system. The high albedo of snow and ice surfaces relative to land and ocean strongly influences the global radiation budget of the Earth. Thus, changes in the spatial extent of the cryosphere will have a direct bearing on radiation balances and thus on climate. The ice sheets of Greenland and Antarctica comprise the major storage of fresh water on our planet. Any volume change of these bodies will have, and has had in the geologic past, an immediate effect on global sea level. The proximity of ambient snow and ice temperatures to the melting point of ice leads to a high sensitivity of snow and ice surfaces to even moderate climate warming. In addition, positive feedback loops between temperature changes, ice and snow extent, decreasing albedo and enhanced absorption of solar radiation constitute the particular sensitivity of the cryosphere to climate change. Similarly, changes in snow and ice masses have a direct bearing on climate change.
However, there are additional links between the cryosphere and the global system that require consideration. The interaction between major ice shelves and the underlying ocean ultimately leads to the formation of Antarctic Bottom Water (AABW), a water mass that influences oceanic processes on global scales far beyond the Antarctic realm. The Weddell Sea and the adjacent ice shelves in the Atlantic sector of Antarctica play a distinct role in this regard.
These issues are central to our work at the University of Münster, Germany. We strive to describe and better understand the present and possible future behavior of the cryosphere in this part of Antarctica, mainly by means of numerical model studies. The choice of numerical models as our primary method is motivated by (i) the relative inaccessibility to in situ measurements of the regions studied, (ii) the possibility of joining modeling and remote-sensing observations, both in combined investigations and as a means for mutual verification, (iii) the potential to identify areas that deserve particular attention in field studies through numerical modeling, and (iv) the possibility of deriving projections of the future development of ice masses by means of numerical models.
In this paper, we briefly describe some of the methodologies that we have developed and provide insight into the dynamics of the Antarctic cryosphere. More specifically, we address the dynamics of the temperate ice cap on King George Island, which requires modifications to numerical methodologies geared towards cold ice bodies; the breakup of the Larsen ice shelves; and the dynamics of the eastern Weddell Sea ice shelves (EWIS), Ekströmisen and the inland-ice/ice-shelf system of Nivlisen (Fig. 1). While we cannot deal with all of these issues comprehensively, we highlight some of the major results of our studies that underline the importance of these ice bodies in the global (climate) system.
Fig. 1. Geographic overview of the study areas in the Atlantic sector of Antarctica.
Methodology
General
The dynamics of ice bodies can be described by conservation equations of mass, energy and momentum and a constitutive equation (for more details and explanations on modeling strategy, see Reference GlenGlen, 1955; Reference Paschke and LangePaterson, 1994). In addition, accumulation and ablation rates and other meteorological and oceanographic boundary conditions must be specified. Being highly non-linear, the partial differential conservation equations and the constitutive equation must be solved numerically. Since conditions at an ice–bedrock (inland ice) and at an ice–ocean interface (ice shelf) differ drastically, the numerical treatment of these two regimes requires separate considerations.
Inland ice
Stress and velocity fields as principal parameters to be derived from modeling are computed through a three-dimensional higher-order flow model developed by Reference PatersonSandhäger (2000). This takes into account the role of longitudinal and transverse stress gradients, yielding for stress equilibrium:
where Tij (with i, j = x, y, z) represent the components of the deviatoric stress tensor, the mean density of the overlying ice column, g the gravitational acceleration, h the surface elevation above a reference surface (usually sea level) and Rzz the vertical resistive longitudinal stress:
σxz , σyz represent the vertical shear stresses in x and y direction, respectively. The value of Rzz is often negligibly small, but in regions with significant subglacial bedrock topography, near the grounding line and close to ice divides and mountain tops, the stress field is mainly affected by so-called bridging effects, especially in the lower part of the ice body (Reference ThomasVan der Veen and Whillans, 1989).
Regarding the energy balance, ice temperatures are derived from a slightly simplified version of the general heat-transfer equation (e.g. Reference Paschke and LangePaterson, 1994):
with u, v and w being the x, y and z components of the ice velocity vector, T temperature, t time, k i thermal conductivity, c p the specific heat capacity of ice, and T the effective stress. A and n are parameters that need to be specified in Glen’s flow law as described below.
The numerical integration scheme proposed by Reference BlatterBlatter (1995) is used to solve the resulting equations after a transformation into sigma coordinates.
We focus on the description of glacier movement in our studies for which Reference Paschke and LangePaterson (1994) distinguishes three components. The plastic deformation of polycrystalline ice (first component) is usually described by the empirically derived Glen-type flow law, which relates strain rate żij to the stress deviators Tij and the effective stress T:
where m and n are adjustable parameters, m often called the enhancement factor. The Arrhenius factor A depends on different conditions in the ice body, such as the potential temperature T *. If ice temperatures approach the pressure-melting point, the water content W cannot be neglected. The latter is especially important for the simulation of ice dynamics of the temperate ice cap of King George Island.
The sliding of ice over its bed (second component) and the deformation of the bed itself (third component) both describe basal motions which largely occur when the temperature at the base of a grounded ice body reaches pressure-melting-point conditions. Since both components are hard to distinguish, an empirical relation between the horizontal flow velocity , the basal vertical shear stress and the effective pressure is implemented in the numerical flow model:
Z represents the height above buoyancy that replaces the effective pressure following the suggestion of, for example, Reference HuybrechtsHuybrechts (1992). The basal flow parameter C b depends on the inverse bedrock roughness and on the type of underlying sediments. Reference Oerter, Graf, Schlosser and OerterPaschke and Lange (2003) discuss the importance, and the difficulties, of setting this value focusing on the drainage area of Nivlisen (see below).
Ice shelves
The diagnostic simulation of ice-shelf dynamics is based on numerical calculations of an approximate solution of two coupled differential equations which describe the distribution of the horizontal ice-flow velocity. The differential equations are obtained by combining the continuum-mechanical momentum and mass-balance equations with Glen’s flow law, the ice-shelf approximation (depth independency of the horizontal velocity field) and the incompressibility condition. With regard to a regular Cartesian x, y, z coordinate system, these governing equations of ice-shelf flow may be written as (e.g. Reference Lythe and VaughanMacAyeal and others, 1986):
and
with
where F is a function of the effective viscosity, gGx,y describe the pressure gradient due to gravity, is the horizontal velocity vector, H the ice thickness, h ice surface elevation above sea level and ρ c the density of completely consolidated ice (915 kgm–3). The parameters A(T), n and m of Glen’s flow law are predetermined as follows: the flow factor A(T) as recommended by Reference Paschke and LangePaterson (1994); the exponent n = 3; the enhancement factor m = 1, or m 6≠ 1 if model adjustment is needed. The temperature distribution in the ice body is calculated using Equation (4).
The stress conditions in the ice body, i.e. the distributions of the horizontal deviatoric stress components T xx, T yy and T xy, are determined by means of Glen’s flow law from the strain rates (gradients of the velocity field). In combination with boundary values (forat the ice front; for T at the surface and bottom of the ice body), this set of model equations describes the distributions of the relevant ice-dynamic quantities u, v, T, T xx, T yy and T xy, subject to the respective ice body geometry.
The model grid used for the numerical simulations is composed of block-centered cells. Each cell represents a volume element of the ice body with a constant base Δx×Δy and a height varying with depth level between 0.01H and 0.2H. The specific numerical method incorporated into the flow model is a relaxation procedure as described by Reference Herterich, Van der Veen and OerlemansHerterich (1987). This procedure facilitates an iterative calculation of an approximate solution of the governing equations of ice-shelf flow (see above) to their finite-difference forms.
The results of a diagnostic flow model yield the steady-state mass exchange between the ice body and the atmosphere and ocean by means of the vertically integrated mass-balance equation in its time-independent formulation
where a s is the surface accumulation rate and a b the melting/ freezing rate at the base of the ice body. a s and a b are measured in mice a–1.
Ocean circulation
Ocean circulation in sub-ice-shelf cavities is modeled using a three-dimensional general ocean circulation model derived from the Bryan (1969) and Reference CoxCox (1984) primitive equation formulation. The model numerically solves a set of differential equations, including the non-linear partial differential Navier–Stokes equation, describing the motion of a fluid in Eulerian notation,
the continuity equation which ensures mass conservation
and Laplace equations describing the diffusion of the conservative properties salinity and potential temperature,
In Equations (10–12), the following notations are used: u, v are horizontal velocities, λ and ϕ are longitude and latitude, respectively, r E is the Earth’s radius, f is the Coriolis parameter or planetary vorticity (f = 2Ωsin ϕ), with Ω being the angular velocity of the Earth), ρ 0 is the mean density of the ocean, ρ the varying density and p the pressure.
In spherical coordinates, d/dt is given as
and the frictional terms Fλ and Fϕ are in their simplest form describable by the Laplace equation. The molecular diffusivities for temperature and salinity have to be replaced by eddy diffusivities A T and A S as a consequence of the turbulent nature of oceanic currents.
The utilized equations are simplified by applying a number of approximations, namely incompressibility, hydrostatic conditions, the ‘shallow-water’, Boussinesque, traditional and rigid-lid approximations (e.g. Reference Haidvogel and BeckmannHaidvogel and Beckmann, 1999). The model is formulated in vertically scaled sigma coordinates, i.e. terrain-following coordinates (Reference GerdesGerdes, 1993; Reference Haidvogel and BeckmannHaidvogel and Beckmann, 1999), that allow for a better representation of kinematic boundary conditions at the ice-shelf base and of the ocean bottom in the cavity. Density is calculated through an empirical equation of state (Reference Markus, Kottmeier, Fahrbach and JeffriesMellor, 1991),
where hydrostatic instabilities in the water column are removed by a convective adjustment scheme that ensures neutral stability at the end of each time-step. The present numerical sigma-coordinate code ROMBAX (Reference SimmonsThoma and others, 2005) is a rigorously restructured version of previous numerical ocean models used in regional studies (e.g. Reference Grosfeld, Gerdes and DetermannGrosfeld and others, 1997).
The thermodynamic interaction between ocean and ice-shelf base is mathematically formulated by an empirical function for the freezing-point temperature (Reference Foldvik and KvingeFoldvik and Kvinge, 1974) and conservation equations for heat and salt at the ice–ocean interface, described in detail by Reference Holland and JenkinsHolland and Jenkins (1999). The amount of melting and freezing at the ice-shelf base influences the density stratification of the ocean in terms of fresh-water (salt) fluxes. It also changes the geometry of the ice shelf, which becomes important when modeling ice-shelf dynamics in a coupled ice-shelf–ocean system (Reference Grosfeld and SandhaägerGrosfeld and Sandhäger, 2004).
Case studies
The methodology described above has been applied to a number of case studies in specific regions of the Weddell Sea sector of Antarctica (Fig. 1). In the following, we present a selection of major results of these investigations.
The temperate ice cap of King George Island, South Shetland Islands
Cold ice sheets with temperatures below the pressure-melting point and negligible amounts of meltwater (not affecting the ice flow) are distinguished from usually smaller and temperate ice caps (temperatures at pressure-melting point) where meltwater is present within the ice body (Reference Paschke and LangePaterson, 1994). An example of the latter is the ice caps of the sub-Antarctic islands which have been subject to increasing surface temperatures over the last few decades (cf. Reference Van der Veen and WhillansVaughan and Doake, 1996). King George Island (KGI), the largest of the South Shetland Islands, is located in the marginal sea-ice zone north of the Antarctic Peninsula at 62˚ S, 58˚W and is influenced by subpolar maritime climate conditions. The temperate ice cap on KGI contains a significant amount of water, considerably affecting its dynamics. This has to be taken into account through modifications of the three-dimensional numerical flow model describing its dynamics. The modifications mainly affect Glen’s flow law, namely the Arrhenius factor A and the enhancement factor m (Equation (5)). In order to account for the ice cap’s water content W within the model equations, an approach similar to that proposed by Reference Greve, Weis and HutterGreve and others (1998) has been applied. We hereby explicitly account for the effect on ice rheology of water present in the ice body and determine the time-varying water content of the ice cap.
Comparison of modeled (grey arrows and contour plot) and measured (black arrows) horizontal surface ice velocities for parts of the KGI ice cap (for additional details, see text).
Boundary conditions describing the ice geometry have been derived from field measurements during austral summer 1997/98. In the course of this campaign, an area of 250km2 in the northwestern part of the ice cap on KGI was investigated using differential global positioning system (GPS) and radio-echo sounding methods. This enabled the determination of surface and bedrock topographies and the ice-thickness distribution.
Figure 2 represents first results of our modeling study. Shown here are modeled vs measured ice velocities for parts of the KGI ice cap. The modeled ’surface’ ice velocities at stations 3, 7, 8 and at the Basis are particularly noteworthy, as they are based on measured ice thicknesses, while the thicknesses east of the ice divide were obtained from extrapolations. Ice thicknesses in the region west of the ice divide have been surveyed extensively, but only one velocity measurement was obtained (station 8). The comparison nevertheless yields good agreement between modeled and measured ice velocities and encourages us to continue our modeling efforts with regard to the KGI ice cap.
Larsen C ice shelf
Map of the LIS C region showing the distribution of water-column thickness beneath the ice shelf and in the adjacent Weddell Sea area. A shaded relief map of the inland is added. (b) Mass balance at the ice-shelf base resulting from the ocean modeling (B b– indicates melting, B b+ freezing). B b is the net basal volume balance, āb the spatial mean of the rate. (c) Steady-state mass balance calculated using the ice-shelf modeling results and an estimate for the surface accumulation rate. Symbols as in (b).
The major section of the Larsen C ice shelf (LIS C; Fig. 3a) covers an area of about 57 000km2 at the eastern Antarctic Peninsula between 66˚ S and 69.5˚ S. Ice surface elevation and ice-shelf thickness data for this region are available from several sources (e.g. Reference Bamber and HuybrechtsBamber and Huybrechts, 1996; BAS and others, 1998; Reference Lucchitta, Mullins, Allison and FerrignoLythe and others, 2001) and can be combined into digital geometric models suitable for detailed ice-dynamic modeling (Reference SandhägerSandhäger, 2003). In contrast, the seabed relief is largely unknown. The first information on bedrock topography underneath the Larsen Ice Shelf was published by Reference DomackDomack and others (2000) after the decay of the Larsen A ice shelf. Based on this information, the bedrock is highly structured into troughs and ridges. We constructed a new bedrock elevation map for the entire Larsen Ice Shelf region, compiling bathymetric data in front of the ice shelf (IOC and others, 1997) and information on bedrock depths at the grounding line (Reference Lucchitta, Mullins, Allison and FerrignoLythe and others, 2001). Underneath the ice shelf, the unknown relief is extrapolated based on considerations of glacio-geological evidence derived from data at the ice-shelf boundary. The resulting complete set of digital geometric models, including a gridded water-column thickness dataset (Fig. 3a), facilitates dual estimates of essential mass-balance parameters of LIS C by means of ocean and ice modeling.
The described ocean model is applied at 0.1˚×0.1˚ horizontal resolution and 14 vertical layers, ranging from 2% to 23% of the water-column thickness. The model is forced with a climatological wind field after Reference Korth and DietrichKottmeier and Sellmann (1996) and is initialized with monthly-mean temperature and salinity values from the large-scaleWeddell Sea model study by Reference Sandhäger and BlindowSchodlok and others (2002). The general flow regime (not shown) indicates a strong northward coastal current along the continental-shelf break. Due to topographic steering in a deep trough of the continental shelf, warm waters from the deep sea enter the region south of Kenyon Peninsula, leading to high melt rates along the grounding line (Fig. 3b). This region is separated from the central part through a ridge system connecting Kenyon Peninsula and Gipps Ice Rise, preventing the penetration of warm waters into the deep cavity. In the central part of LIS C, a separate flow regime establishes, leading to low melt rates. Increased melt rates occur only along the grounding line, where various glaciers enter the ice shelf. Since the parameterization of the melt rate is proportional to the difference between the pressure-dependent freezing point and the ocean temperature at the ice-shelf base, it is directly linked to the ice thickness, which is largest where drainage glaciers enter the ice shelf. Therefore, the melt rate strongly depends on both ice-shelf and seabed topography and is possibly overestimated in our model results. In the northern cavity, basal freezing occurs adjacent to the melting zone, as a result of fast upwelling of supercooled waters along the rapidly descending ice-shelf base. Since the freezing rates are very low, we do not expect any accretion of marine ice beneath LIS C, and identify this feature as a model artifact rather than a realistic depiction of the true melt/freeze pattern. The total net basal melt rate amounts to about 75km3 of ice per year, indicating a potentially substantial contribution to the Weddell Sea fresh-water budget from this comparatively small ice-shelf region.
The diagnostic simulation of ice-shelf dynamics was performed with an extended version of the ice-flow model to account for the possible effects of weakness zones (shear fractures, crevasse bands) on the flow regime of the Larsen Ice Shelf (Reference SandhägerSandhäger, 2003). Even though a relatively simple approach was used to consider fracture-induced flow effects and prescribed glacier inflow velocities (based on rough estimates and a modification of the flow law (Equation (5)), plausible model results were obtained for the horizontal ice-shelf velocity. This dataset can be combined with the digital ice-thickness model to quantify the vertical mass balance a s + a b by means of the continuity equation. To estimate the basal mass balance a b of LIS C (Fig. 3c), we considered a distribution of surface accumulation rates similar to those mapped by Reference Vaughan and DoakeVaughan and others (1999). The accumulation rates vary between 0.40 and 0.67m ice a–1, with a mean value of 0.52mice a–1. However, uncertainties of ice velocity, ice thickness and surface accumulation rates cause correspondingly large errors of basal melt/freeze rates given in Figure 3c.
Both model studies for LIS C indicate a high net basal melt rate. Although results differ by about a factor of two in their estimates, the general pattern indicates increased melting along the grounding line, while only minor melting occurs in the central part. The discrepancies are mainly due to the almost unknown ice thickness and seabed topography distributions in the inflow regions of the glacier systems, which ultimately govern the dynamic flow regime of the ice shelf and the oceanic circulation pattern underneath.
The eastern Weddell Sea ice-shelf region
The eastern Weddell ice-shelf (EWIS) region represents an important water source area for the entire Weddell Sea. Warm Deep Water originating from the northern branch of the Weddell Gyre and the westward-flowing Antarctic Coastal Current enter the southern Weddell Sea in this region, which is known to be the most efficient region for Antarctic Deep and Bottom Water formation. The EWIS region is of particular interest as it is separated from the Antarctic Coastal Current by a narrow, steep continental-shelf margin. Hence, direct interaction of the strong coastal current with the ice-shelf cavity provides for high freshwater input rates from glacial melting.
The EWIS is located at 71–76˚ S, 10–28˚W and consists of Riiser-Larsenisen in the north and the Brunt Ice Shelf in the south, separated by Lyddan Island (Fig. 4). Two additional ice rumples near the ice front of Riiser-Larsenisen act as anchor points for the ice-shelf flow. The Brunt Ice Shelf is mainly fed by discharge from Stancomb-Wills Ice Stream, while the major ice stream draining toward Riiser-Larsenisen is Veststraumen. The largest ice thicknesses in the EWIS are up to 600m where the ice streams from the hinterland enter the ice shelves. The total ice-shelf area and volume of the EWIS are 75 000 km2 and 16 000 km3, respectively.
Stream function (a) and basal mass balance (b) in the EWIS region simulated with the ocean model. The northern Riiser-Larsenisen yields the highest melt rates due to the direct inflow of warm deep waters. Increased melt rates are modeled where ice-shelf thickness is largest, namely along Veststraumen (at about 74˚ S) and for Stancomb-Wills Ice Stream (entering the Brunt Ice Shelf at 75.3˚ S). (c) Ice velocities as obtained by the ice-shelf model.
In order to investigate the principal flow regime of the EWIS region and the fresh-water production due to glacial melting, we apply ROMBAX in a medium-resolution version of a 0.3˚ longitude 0.1˚ latitude grid. 14 sigma layers in the vertical (i.e. layers in the terrain-following coordinate system used in the simulation) provide for high resolution near the surface (20 m), whereas coarser resolution is chosen in the abyss (23% of water-column thickness). Model initializations for temperature and salinity are taken from the hydrographical dataset of Reference GouretskiGouretski and others (1999). In the inflow region, ‘Newtonian damping’ is applied to World Ocean Circulation Experiment (WOCE) A12 transect data along the prime meridian (Reference Fahrbach, Hoppema, Rohardt, Schröder and WisotzkiFahrbach and others, 2004). Newtonian damping is a numerical procedure designed to avoid heterogeneities/discontinuities between fixed boundary conditions and computed results in the vicinity of that boundary (e.g. Reference Barnier, Marchesiello, de Miranda, Molines and CoulibalyBarnier and others, 1998). The computational cells thus modified are sometimes referred to as the ‘sponge layers’. The ocean surface is forced with monthly climatological wind fields (Reference Korth and DietrichKottmeier and Sellmann, 1996) and is Newtonially damped to the initial or winter-adjusted T – S values, depending on the season.
Figure 4a presents the mean state for the vertically integrated mass-transport stream function. It depicts a coastal current of about 4 Sv flowing along the continental-shelf break, where the bathymetry reaches the abyssal plain, while the transport in the central Weddell Sea is much larger. The coastal current in the EWIS region directly interacts with the ice-shelf region. The flow beneath the EWIS is characterized by a cyclonic through-flow from northeast to southwest. Only in two areas near Lyddan Island do small, weak anticyclonic gyres establish. Simulated basal mass balances in the EWIS region (Fig. 4b) show that northern Riiser-Larsenisen yields the highest melt rates due to the direct inflow of warm deep waters. Increased melt rates are modeled where ice-shelf thicknesses reach their highest values, namely along the two major ice streams entering the ice sheet at about 74˚ S (Veststraumen) and 75.3˚ S (Stancomb-Wills Ice Stream).
Our ice-shelf model was also applied to the EWIS region in order to explore possible ice-shelf–ocean interactions. In so doing, we restrict our model domain to the ice-shelf area, i.e. we deal with stress-free bottom conditions. This implies that we must specify depth-invariant ice inflow velocities along the grounding line, thus introducing considerable potential uncertainties in our calculations. However, we are able to constrain the inflow velocities through results from satellite remote-sensing studies for the Brunt Ice Shelf (Reference Schodlok, Hellmer and BeckmannSimmons, 1986; Reference Kottmeier and SellmannLucchitta and others, 1993). For Riiser-Larsenisen such data do not exist, except for individual ice velocities for the Veststraumen and Plogbreen ice streams. However, our results indicate that a choice of different but reasonable inflow velocities has negligible effects on the final results in diagnostic (i.e. ‘snapshot’) models, whereas the specific choice of inflow velocities will have a major bearing on results in prognostic (i.e. time-dependent) models.
The flow regime of the Brunt Ice Shelf is dominated by Stancomb-Wills Ice Stream, with surface velocities of up to 1200ma–1 near the grounding line (Fig. 4c). Riiser-Larsenisen derives most of its input from Veststraumen and Plogbreen in the south of the ice shelf, with surface velocities of about 100ma–1. Lyddan Island, separating the two ice shelves, is most pertinent to ice flow. The flow velocity of Riiser-Larsenisen increases to about 700ma–1 near the ice front. The values derived here are in the same order of magnitude as established by Reference Thoma, Grosfeld, Mohrholz, Lange and SmedsrudThomas (1973) for the 1960s.
Ekströmisen
Ekströmisen and its catchment area extend over about 8700 and 20 700km2, respectively, and represent a comparatively small drainage system in the Atlantic sector of East Antarctica. The ice shelf is the northeastward extension of the EWIS. The Southern Ocean in this region is characterized by an extremely narrow continental shelf about 35km wide. A reliable description of ice surface topography (Fig. 5a), ice-thickness distribution and subglacial bedrock relief was obtained by determining detailed digital geometric models from airborne radio-echo sounding and altimetry data (Reference PatersonSandhäger and Blindow, 2000). Additionally, the sea-floor topography is taken from the General Bathymetric Chart of the Ocean (IOC and others, 1997).
Relief map of Ekströmisen and adjacent areas north of 72˚ S. The flat ice shelf is bounded by several ice domes and the Ritscherflya ice-sheet region to the south. Mass discharge from inland into the ice shelf is mainly concentrated on four active grounding zones (stars). (b) Modeled distribution of depth-averaged ice velocity and associated flowlines indicating the direction of horizontal ice flow. (c) Annual mean mass balance from the ocean model (negative values represent melting).
Ice dynamics of the drainage system were quantified by means of a flow model in time-independent/diagnostic mode (for more details end explanations on modeling strategy, see Reference PatersonSandhäger, 2000). To compute velocity, stress, strain rate and temperature fields of the ice body, the ice-shelf–ice-sheet system was assumed to be in a steady state, and the geometric datasets as well as estimates for the geothermal heat flux (G = 0.0546Wm–2) and the distribution of mean annual surface temperature T s were used as basic model input. T s depends on freeboard height h (in m) as follows: T s = (–18˚C– 0.004˚ Cm–1)h. The horizontal resolution of the model grid is 1 km.
Figure 5b shows the modeled distribution of depth-averaged ice velocities, which are mostly <25ma–1 inland, but exceed 250ma–1 at the front of Ekströmisen. The mass discharge from inland into the ice-shelf area is concentrated into four active grounding zones where mass fluxes of >100 000 t a–1m–1 occur. Passive grounding-line sections, however, are characterized by mass fluxes of typically <10 000 t a–1m–1. A thorough validation of the model results reveals overall good agreement with direct observations (Reference PatersonSandhäger, 2000).
The ice-shelf–ocean interactions were analyzed by applying a three-dimensional thermohaline circulation model based on a version of an Ocean General Circulation Model (OGCM) (Reference BryanBryan, 1969; Reference CoxCox, 1984). One of the main advantages of this model approach is that it permits sufficiently high resolution (3.4–4.2 km zonally and 5.6km meridionally) to resolve ice-shelf cavity and shelf processes. A detailed model description and geometric setting as well as initialization and restoring are given by Reference GerdesGerdes (1993), Reference Grosfeld, Gerdes and DetermannGrosfeld and others (1997) and Reference MellorNicolaus and Grosfeld (2002).
Based on the steady-state assumption, the fundamental mass-balance quantities of Ekströmisen can easily be calculated from the digital ice-thickness model, the simulated ice-velocity distribution and an estimate of a s = 0.34mice a–1 for the mean surface accumulation rate (Reference Nicolaus and GrosfeldOerter and others, 1997). The mass flux over the grounding line (~4.3 Gt a–1) and the mass gain due to surface accumulation (~2.7 Gt a–1) are compensated to ~46% (~3.2 Gt a–1) by calving and melting at the ice-shelf front and to ~54% (~3.8 Gt a–1) by melting at the ice-shelf base. Thus, the mean basal mass balance of Ekströmisen is a b–0.48mice a–1 (Reference PatersonSandhäger, 2000). Corresponding results from the ocean model range from –0.95mice a–1 (winter) to –1.02mice a–1 (summer), which is about twice the value obtained from the ice model. In addition, Reference MellorNicolaus and Grosfeld (2002) showed that most of the basal melting occurs where the ice streams enter the ice shelf, and that little or no accumulation of marine ice takes place (Fig. 5c). All of these basal mass-balance studies are in good agreement with a two-dimensional ice pump model applied by Kipfstuhl (1991) and field measurements by Reference MacAyeal, Shabtaie, Bentley and KingMarkus and others (1998) yielding melt rates of 0.5–1.0mice a–1.
Nivlisen: coupled ice-sheet/ice-shelf regime
Following the application of the numerical flow model to Ekströmisen by Reference PatersonSandhäger (2000), Reference Oerter, Graf, Schlosser and OerterPaschke and Lange (2003) applied the model to Nivlisen, East Antarctica, and its drainage area (70.5–73˚ S, 9–14˚ E). The ice-shelf and catchment areas comprise about 7200 and 50 000km2, respectively. The catchment is characterized by a complex ice geometry and consequently a complicated ice-flow regime that had to be taken into consideration. This includes extended parts of the inland ice sheet; 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. As a consequence of the interactions between all of these different features and the non-linear character of the higher-order model equations, Reference Oerter, Graf, Schlosser and OerterPaschke and Lange (2003) introduce additional constraints for the diagnostic ice-dynamics simulation. They demonstrate that setting one central flowline of each glacier system to pressure-melting-point conditions results in a steady-state solution that agrees significantly better with the glaciological situation than the results obtained without this modification. Comparisons between modeled and measured ice velocities suggest the introduction of spatially varying boundary conditions in contrast to the usual uniform characteristics (e.g. with regard to the basal flow parameter Cb (Equation (6)). The value of this parameter depends on the inverse of the bedrock roughness and on the type of underlying sediments. However, due to the scarcity of sufficient field data, a constant value for Cb is often used.
Implementing these modifications in the numerical flow model results in satisfactory agreement between computed and measured ice velocities, i.e. deviations between these velocities do not imply significant differences in velocity magnitude or in the flow direction. Figure 6 shows the horizontal velocity at the ice surface. The grey shading indicates modeled ice velocity, and the white arrows flow direction. For comparison, black arrows denote the direction of surface velocities obtained through in situ measurements by Reference KipfstuhlKorth and Dietrich (1996). Most of the points show good agreement concerning flow direction. More details are given in Reference Oerter, Graf, Schlosser and OerterPaschke and Lange (2003). They conclude that complex flow regimes may be better treated by employing a nested grid solution with higher spatial resolution for particularly complicated flow situations. In addition, they suggest that decoupling of the velocity fields of fast-flowing glacier systems and near-stagnant portions of the ice sheet may yield further improvements of the numerical results.
Surface ice velocities for Nivlisen and its drainage system (white arrows and shading) as derived from the numerical flow model of Reference Oerter, Graf, Schlosser and OerterPaschke and Lange (2003). Unscaled black arrows give the direction of measured surface velocities of Reference KipfstuhlKorth and Dietrich (1996).
Discussion and Conclusions
Given the practical and logistical difficulties of field operations in Antarctica, numerical models of ice-dynamics processes may be considered the most feasible and cost-effective way of gaining a better understanding of these processes. Moreover, numerical models are the only alternative when it comes to predictive studies that aim to address the impacts of climate change on the cryosphere. While the basic methodology needed to construct such models is well known, the complexity of ice flow in a natural setting and the scarcity of necessary boundary conditions continues to present significant challenges to the ice-dynamics modeler.
The case studies described above illustrate some of these challenges. The modeling chiefly rests on well-established mathematical descriptions of conservation relations for mass, energy and momentum in combination with constitutive equations that describe the rheological behavior of an ice body. These formulations are usually applied in a straightforward manner, whereas we strive to apply them to situations usually not considered for Antarctic ice dynamics.
The case study of the ice cap on KGI is the first case in point. Here we modify Glen’s flow law in order to account for the temperate nature of the ice cap and the presence of significant amounts of meltwater. The resulting ice velocities largely agree with measured values or those indirectly derived from measurements. However, some discrepancies remain. They most likely reflect both insufficient observational data and deficiencies in model formulations. More work on this subject is currently underway including the acquisition of additional field data.
The work carried out for Nivlisen and its drainage area represents possibly the most challenging investigation carried out so far. The transition between inland ice and an ice shelf is still open to several questions. However, the transition considered here does not take place through a ‘regular’ grounding zone, but through relatively narrow outlet glaciers between steep mountain tops. Consequently, we had to pursue this problem by applying new and somewhat arbitrary parameterizations of ice flow across the bedrock. While we consider the results satisfactory, we cannot exclude that they are fortuitous. Clearly, more work is needed to come up with a more satisfactory conclusion.
In comparison, the study on Ekströmisen represents a more ‘conventional’ type of ice-dynamics investigation. The uniqueness of the work presented here lies in the relatively well-documented characteristics of this ice shelf and parts of its drainage area. This is due to the fact that Ekströmisen has been and continues to be the ‘home base’ of German Antarctic research and has consequently been well studied. This also implies that we have a solid foundation for verifying our model results and for fine-tuning our methodologies. Thus it is not surprising to learn that the methodologies developed in the context of the Ekströmisen model have been the starting point for most of our ice-dynamics modeling. Moreover, the Ekströmisen work was also the starting point for our work on ice–ocean interactions, and led to very promising first results.
Consideration of the interactions between an ice shelf and the underlying ocean has been a major advance in our work. It implies that we have to couple two models that operate on very different timescales. While the flow of ice shelves, even though fast in comparison to the inland ice, takes place on a scale of hundreds of years, the ocean circulation underneath an ice shelf requires months or years at most. The two case studies presented, for the LIC and the EWIS region, illustrate the significant increase in information gained by following such an approach. The models not only result in more realistic ice-dynamics representations, but offer insight into the consequences of melting and freezing for the water masses adjacent to the ice shelves. Both case studies, though representing very different settings, demonstrate that even smaller ice shelves must be taken into account when aiming to address the contribution of ice– ocean interactions to the formation of major Antarctic water masses.
Thus, while offering an overview of some of our work, our studies demonstrate that significant progress has been made in ice-dynamics modeling in Antarctica during the last decade. While several of the investigations reported here are still in progress, the results obtained so far are truly encouraging.
Acknowledgements
We gratefully acknowledge financial support for some of our projects by the Deutsche Forschungsgemeinschaft and the Federal Research Ministry (bmb+f contract 03F0377C). We also appreciate the many helpful and constructive comments and suggestions provided by numerous colleagues in the course of our projects. Finally, we express our gratitude to the two reviewers of the paper for their positive and constructive comments and suggestions and to the scientific editor, G. Hamilton, for providing valuable help.
