Introduction

In numerical simulations of the glacial cycles of the Pleistocene with climate models, the ice sheets form the central prognostie component, while the atmosphere and the ocean follow the boundary conditions (in part set by the ice sheet) in a quasi-stationary way. If one is interested only in time scales longer than 1000 years, even parts of the ice sheet may be simulated by a quasi-stationary description. In particular, the ice shelf and the transition zone between ice sheet and ice shelf have short time scales compared to the ice sheet, and may thus be modelled as quasi-stationary. Reference OerlemansOerlemans (1982) developed a three-dimensional ice-sheet/ ice-shelf model but without including an explicit formulation of the transition zone. The influence of such a transition zone on the evolution of the coupled system was studied by Van der Veen (Reference Van der Veen1985) for the two-dimensional case, demonstrating its dynamic (stabilizing) effect. A detailed three-dimensional treatment of the transition zone in models of global ice volume to simulate glacial cycles is still outside presently available computing facilities. Thus, either the dynamics of the transition zone is taken into account at only a few grid points close to the grounding line, or, as will be described in this paper, the effect of the transition zone is parameterized. Neither approach, however, provides an adequate description of ice streams.

A three-dimensional ice-sheet model is used here with a fully coupled flow-temperature calculation but without basal sliding (Reference HutterHerterich, 1988). The model was further developed to include a quasi-stationary description of the ice shelf. The two regimes (ice sheet and ice shelf) are coupled by a parameterized form of the transition zone.

Model Description

The model consists of three components: the ice sheet, the ice shelf and the transition zone in between. The ice-sheet model uses the shallow-ice approximation (Reference HutterHutter, 1983), which alows us to derive an analytical expression for the horizontal flow field

for a given ice-thickness distribution and temperature field (Reference HutterHerterich 1988):

where, x, y, z form a Cartesian coordinate system and h_sf,h_B are the heights of the ice surface and the ice base, respectively. The density of ice is ρ = 910 kg m⁻³, g = 9.81 m s⁻² is the gravitational acceleration, and n = 3. The coefficient A(T′) depends on the temperature Τ′ above pressure melting point, derived from measurements compiled by Reference PatersonPaterson (1981). Given the flow, the further evolution of both the ice thickness and the temperature field can be calculated by integrating the well-known conservation laws for mass and energy.

The ice-shelf model, which calculates the stationary-state thicknesses within the ice shelf is a physical empirical model. The ice-shelf thickness h at some point P(x,y) is calculated from

where h_g,i is the thickness at the grounding line, d_i the distance from some point P(x,y) in the ice shelf to the grounding line, and m the number of points by which the grounding line is represented on the coarse grid of the model. The parameter d₀ is fitted to observations. The form in Equation (2) is the solution for a stationary (two-dimensional) ice-shelf profile (m = 1) near the grounding line, with zero accumulation rate, and is supposed to be an approximate generalization to three dimensions, by averaging along the grounding line (i = 1, m). With a parameter value d ₀ = 250 km, which controls the thinning length scale of the ice shelf along a flow line, the ice-shelf model matches the observed positions of the ice-shelf front to some extent, where the simulated front is defined by the model’s 200 m thickness contour.

Fig. 1. Vertical profile of the transition zone in the vicinity of the grounding line and definition of variables used in the text.

The width of the transition zone may be small compared to the coarse grid resolution (≈100 km) of our model (Reference Herterich, Veen and OerlemansHerterich, 1987), and the transition-zone physics may therefore be treated as a sub-grid scale process. The properties of the transition zone of main interest here are the location of the grounding line and the thickness at the grounding line. These quantities have to be expressed in terms of variables defined on the coarse grid of the model. To proceed further, we formulate the integral force balance along a closed loop in a vertical plane perpendicular to the grounding line. The loop is defined by the numbers 1 to 4 indicated in Figure 1, connecting the two neighbouring grid points of the (coarse) grid, just upstream and downstream the grounding line. The integral force balance then reads:

where h_i , h_s are the thicknesses at x = x_s and x = x_s , respectively, σ′_xx is the horizontal stress deviator at x = x_s, and τ the mean vertical shear stress along the path from x_g to x_i . At x=x_i , the shallow-ice approximation holds, the ice-shelf approximation (no vertical shear) is valid at x = x_s. In addition, the slope of the ice base is assumed to be small. In the following, we will neglect the stress deviator σ′_xx on the ice-shelf side at x = x_s, mainly because we cannot calculate σ′_xx from the ice-shelf model (Equation (2)). Thus the model is restricted to situations where σ′_xx is small compared to both the pressure difference

across the transition zone and the bottom friction τ(x _g - x _i ). With the definition h =(h _i+h _s)/2, and assuming h_s ≈ h_g , we obtain:

The non-dimensional ratio μ = tan α (see Fig. 1) is formally a friction coefficient (ratio between friction force and normal load) which depends on the physical properties at the ice bottom. It will be treated as a free parameter in the following. Equation (4) relates the parameter μ to the slope of the ice surface at the margin of the ice sheet and thus determines the outflow from the ice sheet (Equation (1)). The position of the grounding line x_g and the thickness at the grounding line h _g follow from Equation (4) and from the floating condition:

where ρ, ρ _w are the densities of ice and water, respectively, h_s1 is the sea-level height and h_T is the height of the bottom topography.

The numerical calculations for the coupled ice-sheet/ice-shelf system proceed as follows. First, the thickness distribution and the temperature field within the ice sheet are given at some initial time t₀ (initial state). Thicknesses were taken from the maps of Reference DrewryDrewry (1983) and the temperature field resulted from a 150 000 years integration of the temperature equation (energy balance) for a fixed shape and constant boundary conditions (presently observed surface temperature and a geothermal heat flux of 5 × 10⁻² W m⁻²). For a fixed value of μ in Equation (4), the position of the grounding line and the thicknesses at the grounding line can then be determined (using Equations (4) and (5)), together with the ice-shelf thicknesses (using Equation (2)). With the flow resulting from Equation (1), the rate of change of ice-sheet thicknesses follows from the mass-balance equation, thus defining the new position of the grounding line, and hence the new distribution of ice-shelf thicknesses.

Fig. 2. Thickness distribution of ice sheet and ice shelf (mainly white areas) after 10 000 model years integration, starting from (modelled) present conditions, for three different parameter values (a) μ = 0.05, (b) μ = 0.01, and (c) μ = 0.25.

Results

We attempted to simulate the presently-observed areal distribution of ice sheet and ice shelf by performing sensitivity studies with respect to the parameter μ, assumed to be independent of position. This was achieved with μ = 5 × 10⁻² (see Fig. 2a). During a model integration over 10 000 years, with currently-observed accumulation rates (Reference Budd, Jenssen and RadokBudd and others, 1971), the position of the grounding line did not change noticeably compared to the initial state. There was, however, a slight increase in ice-sheet thickness (especially for the East Antarctic ice sheet), indicating that the initial state was not exactly stationary.

The following two 10 000 year simulations show that the model is sensitive to changes in the parameter μ. With μ = 0.01 (Fig. 2b), the present ice-shelf areas were replaced by an ice sheet, while with μ = 0.25 (Fig. 2c) those areas of the West Antarctic ice sheet where the ice bottom is presently below sea level disappeared. Reducing the parameter μ by a factor 1/5 compared to the stationary state (μ = 0.05) effectively cuts off the ice-sheet discharge at locations where the model identifies a grounding line Since the flow at the ice margin is proportional to μ³ (Equation (1)), the outflow from the ice sheet is reduced by a factor (1/5)³. Thus, with the same accumulation rate as used in the experiment with μ = 0.05, the ice-sheet volume grows and the ice-sheet margin is shifted to the boundary of the continental shelf. In the experiment with the parameter μ increased by a factor 5, the outflow increased by a factor 5³, which resulted in a decay of the West Antarctic ice sheet.

Discussion and Conclusions

The ice-sheet model used here (together with the extensions to include an ice shelf and a transition zone), was designed to allow paleoclimatic simulations on time scales of the order of 100 000 years (the length of a typical ice-age cycle) and longer. Due to limitations in computer resources, the model needs to be simple and should have a rather coarse resolution (≈100 km in the horizontal in our model). The possible degree of simplicity will be constrained by the performance of the model in simulating the paleoclimate. We should be able to make a realistic prediction of the thickness distribution of the ice sheet in time and space, which is the first order variable feeding back into the climate system. Although this rather simple model has still not been used in paleoclimatic simulations, we believe that such simulations are feasible. Our optimism is based on two results:

(1) The present distribution of ice shelves for a fixed grounding line was adequately simulated by a simple and stationary model (Equation (2)), with a constant parameter value d₀ = 250 km.

(2) In a dynamically coupled experiment, the grounding line maintained its presently observed position, using a spatially constant parameter value μ = 0.05.

It is significant that we were able to find a mean value for the parameter μ, since the model is rather sensitive to changes in μ, depending on the physical conditions at the ice bottom, which may vary with position. When reducing complicated physics to only one parameter (μ), the fitted parameter value will generally deviate from corresponding observed quantities. Indeed, the value μ = 0.05, which is equal to the surface slope at the ice-sheet margin in the model, is rather high. This high value may partly be due to the neglect of basal sliding. In the model, the flow due to sliding can only be modelled by pure deformational flow which then affords a higher surface slope. The successful simulation of ice shelves with a constant parameter d₀ is a less surprising result, since Equation (2) partly reflects the proper dynamics along a flow line.