1. Introduction

Modeling climatic change on time-scales longer than 100 years must take into consideration the dynamics of ice sheets, which in turn alter the boundary conditions of the other components of the climatic system. During the last 20 years, data on the long-term variations in global ice volume and in many other climate variables have been derived from deep-sea sediment records, which initiated the development of a hierarchy of ice-sheet models, including coupling with models of the atmosphere, the ocean, and the continents. The most sophisticated ice-sheet model has been developed by Reference JenssenJenssen (1977) and was applied to the Greenland ice sheet.

For climate studies, the relevant ice-sheet variable, feeding back into the climate system, is ice thickness as a function of geographical position and time (Reference Oerlemans and van der VeenOerlemans and van der Veen 1984) and there is now growing interest among climate modelers in using three-dimensional ice-sheet models of the Jenssen type to describe the evolution of the (ice-age) ice sheets during the Pleistocene. Meaningful shorter interval integrations may be carried out for the Northern Hemisphere ice sheets, with reference to their initial build-up, starting about 120 000 years ago, or by modeling ice-sheet decay, where the initial conditions are set by data from 18 000 B.P. (CLIMAP Project Members 1976). In both cases, the initial conditions and the mass balance can probably be specified.

In this paper an attempt is therefore made to take up Jenssen’s idea once more and to construct a high-resolution ice-sheet model for testing on the Antarctic ice sheet. The present version still does not include basal sliding and excludes the ice-shelf regions. It is planned to include ice shelves at a later stage by coupling this ice-sheet model to the ice-shelf model ofReference MacAyeal and Thomas MacAyeal and Thomas (1982). In section 2 the equations, approximations, and boundary conditions are described. Preliminary results showing the (stationary-state) velocity and temperature field for the Antarctic ice sheet are presented in section 3 and discussed in section 4.

2. Model Description

In the usual notation (see Reference PatersonPaterson 1981, which is also used to assign reasonable values to the material constants of ice), the general set of dynamic equations is given by the force balance

and the flow law

In Equation (1), g is the vector of the Earth’s acceleration (

), σ is the stress tensor and ρ is the density of ice taken as constant. A value of ρ = 0.9 × 10³ kg m⁻³ (which is the density reached below a depth of about 100 m in ice sheets) has been chosen. Equation (2) relates the rate of change in the components of the deformation tensor έ_ik to the components of the stress deviator tensor defined by:

and

for the non-diagonal elements, where i, k denote the axes of a Cartesian coordinate system, with x, y as the horizontal axes and z directed vertically upwards. In Equation (2), a non-linearity enters through the second invariant σ raised to the power n - 1, where n - 3 and the invariant σ is defined by:

. The temperature-dependent coefficient Α(Τ^’), with T^’ measured above pressure-melting point T _Μ, is based on measurements supplied by Paterson (1981).

Since the rate of change in deformation can be expressed in terms of velocity gradients

, the set of nine equations (Equations (equ.1) and (equ.2)) can be solved in principle for the nine unknowns, the three velocity components u _x, u _y, u _z, and the six components of the symmetrical stress tensor σ.

In the case of an ice sheet (following Reference MahaffyMahaffy 1976), some approximations are introduced. In the force balance (Equation (equ.1)), the shear stresses σ_ik are small compared to the normal stresses σ_ii; the latter are assumed to be isotropic (σ_xx = σ_yy = σ_zz): that is, longitudinal deviatoric stresses are ignored. Recent theoretical and modeling results (Reference Veen van der, Veen van der and Oerlemansvan der Veen 1987) show that, in the presence of basal sliding, the longitudinal deviatoric stresses are not negligible over a distance from a few kilometers to a hundred kilometers up-stream of the grounding line. However, the horizontal resolution of this model is just 100 km and no attempt is made to model details below the grid size. If the grounding line is allowed to move (it is fixed in this version of the model), the effect of longitudinal deviatoric stresses will probably have to be included in some way.

The small aspect ratio of an ice sheet (where the thickness is small compared to horizontal extension) implies that the horizontal derivative of any ice-sheet variable is small compared to its vertical derivative. This is essentially the “shallow ice approximation” described by Reference HutterHutter (1983). It can be derived formally by expanding the set (Equations (equ.1) and (equ.2)) to zero order with respect to the aspect ratio as the expansion parameter. If we discount atmospheric pressure, the set (Equations (equ.1) and (equ.2)) is then easily integrated, yielding the two horizontal velocity components as a function of the vertical coordinate z:

with

In Equation (3), u _x(B) and u _y(B) are the horizontal components of the basal sliding velocity (which is zero in the present version of the model) and in Equation (4), h _s and h _B are the heights of the surface and the bottom of the ice sheet respectively. The vertical component u _z follows from the incompressibility condition ∇․u=0;

In Mahaffy’s (1976) work, the coefficient A(T^’ ) was assumed to be constant. Thus Equation (4) can be integrated analytically. A further integration with respect to z yields the vertically integrated ice flux actually calculated in her model to describe the evolution of Barnes Ice Cap. A vertically integrated form was applied also by Reference OerlemansOerlemans (1983) to build up the Antarctic ice sheet; it included a simple (quasi-stationary) description of ice shelves. Jenssen (1977) used a temperature-dependent empirical relation between the horizontal ice-velocity component and the calculated shear stresses (instead of the flow law in Equation (2)). Furthermore, a “σ-coordinate” scheme with a variable grid size in the vertical was chosen. Although this allows a finer vertical resolution at the ice-sheet edge, it was not used here, because at the ice-sheet edge the “shallow ice approximation” breaks down and additional numerical errors are expected, due to a coarse horizontal resolution.

The rate of change of temperature Τ within the ice sheet is given by the energy-balance equation

where k = 1.15 × 10⁻⁶ m² s⁻¹ is the thermal diffusivity and d is the contribution of deformational heat given by

with c = 2009 J kg⁻¹ K⁻¹, the specific heat capacity of ice. Equation (6) yields the temperature evolution, starting from an initial temperature distribution within the ice sheet (−10°C throughout) and prescribed boundary conditions (see below).

At the ice-bedrock interface, the geothermal heat flux G enters the ice. For temperatures below freezing point, this fixes the vertical temperature gradient at the ice bottom:

where λ = 2.1 W m⁻¹ K⁻¹ is the thermal conductivity of ice and G = 5 × 10⁻² W m⁻². If the bottom temperature reaches the melting temperature T_M , Τ = T_M replaces the boundary condition (Equation (equ.8)) and the geothermal heat flux is used for melting. The melting temperature is corrected for pressure: T_M = —apgh, where h is the thickness of the ice sheet and a = 7.42 × 10⁻⁵ Κ (kPa)⁻¹.

At the ice-atmosphere interface, the surface temperature T_s is prescribed. For the preliminary calculations described in section 3, a temperature of TSL = − 7°C was chosen at sea-level hss, and an atmospheric lapse rate of Y = 13 Κ km⁻¹ was assumed. From these numbers, a surface temperature is defined for a given height hs of the ice surface: T_s = T_SL - Y(h_S - h_ss), which yields an approximate fit to the observations.

Fig. 1.

Definition of the staggered grid and volume element ΔV = Δx Δy Δz used for numerical calculations (see also section 2).

Equations (3), (4), and (5) for the ice flow, and Equation (6) for the temperature evolution are discretized using centered differences on a staggered grid (Fig. 1). The temperature is defined on grid points (i,j,k), which mark the centers of volume elements ΔV = Δx Δy Δz, for which the energy balance (Equation (equ.6)) has been formulated. To calculate the advective input and output heat fluxes, it is convenient to define the normal components of the ice flow at the centers (i ± 1/2, j ± 1/2, k ± 1/2) of the six planes, forming the volume element ΔV. The resolution of the grid is Δx = Δy = 101 km in the horizontal, and Δz = 285 m in the vertical, and the time step is Δt = 10 years. To suppress the evolution of numerical instabilities during the (explicit) time integration of the temperature equation (Equation (equ.6)), numerical diffusion is added in the form of an up-stream scheme. This scheme introduces artificial diffusion, of the order of the horizontal advective heat transport. In a later version of this model the up-stream scheme will probably be replaced by implicit time-stepping. Calculations are run on a Cyber 205 (CDC) computer. Using the above resolution, computing time is 15 min for 10 000 model years of the Antarctic ice sheet.

3. Results

All the calculations shown describe the model ice sheet in a stationary state. The assumed stationary shape of the Antarctic ice sheet has been derived by digitizing the maps in Reference DrewryDrewry (1983). Ice-sheet heights and thicknesses have been smoothed by allowing the ice sheet to deform for 1000 years according to the mass-balance equation, assuming zero accumulation. Calculations showed that small errors in the surface gradient produce large errors in the vertical velocity component u _z, yielding unrealistic spatial fluctuations in the sign of u _z. The resulting smoothed topography is still close to the unsmoothed topography. During the integration of the temperature equation (Equation (equ.6)) the (smoothed) shape was retained.

In Figure 2, the velocity and temperature field is plotted for a vertical plane (indicated by the line A-Β in Figs 3 and 4) after integrating the heat equation (Equation(equ.6)) for 150 000 years, using present-day boundary conditions of the Antarctic ice sheet (see section 2). The ice velocity is increasing monotonically from the bottom of the ice to the surface. Most of the increase with height occurs within the lower half, whereas the ice velocity is almost constant within the upper half. The temperature is increasing monotonically from the surface down to the bottom, reaching the pressure-melting point over large areas (see also Fig. 4). The temperature gradients at the surface are generally smaller in relation to bottom gradients. Where the velocity is small, the temperature gradient is almost constant, as is to be expected for a stationary state if heat conduction is the only mechanism which can distribute heat. Since the bottom of the ice is at or near the melting point, the isolines of temperature within the ice follow the bottom topography to some extent.

Fig. 2.

Vertical structure of velocity and temperature (°C) for the profile A—Β as indicated in Figures 3 and 4. Heights a.s.l. and the horizontal extension are given in kilometers. (Flow vectors give correct directional information. For plotting reasons, however, the magnitude of the flow is related to the length of the vector non-linearly.)

An evaluation of the ice-flow model is obtained by calculating the divergence of the vertically integrated flow, which should be equal to the observed accumulation rate, provided that the Antarctic ice sheet is in a stationary state. In Figure 3, the sign of the (stationary-state) divergence in the model is shown (it is negative in the hatched areas). In the central areas of the ice sheet, the magnitude of the divergence in the model (5 cm a⁻¹) is near the observed accumulation. However, some areas have a negative mass balance. These areas can generally be identified as the convergence areas of the vertically integrated mass flux, and are probably due to ignoring basal sliding. The negative values of the mass balance at the very edge of the Antarctic ice sheet may be attributed in part to underestimating the surface gradient there (due to coarse horizontal resolution).

Fig. 3.

Hatched areas show negative model divergence (convergence) of the vertically integrated flow.

The resulting bottom temperature is of special interest, because it may control the occurrence of basal sliding, which can contribute considerably to the ice flux. Just small changes in the bottom temperature determine whether there is melting at the bottom or whether the ice is frozen to the bedrock (the hatched area in Fig. 4). The bottom temperature of the model is quite consistent with the location of sub-ice lakes (the dots in Fig. 4) as observed by Reference Oswald and Robin.deOswald and Robin (1973). A comparison of the patterns in Figures 3 and 4 shows that areas of negative mass balance generally lie in regions where the pressure-melting point has been reached. This opens up the possibility of correcting the ice flow by adding some parameterized form of basal sliding, which also depends on the basal water flow. The production of melt water can be calculated directly from this ice-sheet model.

Fig. 4.

Hatched areas show where the model ice is frozen to the bedrock, and dots give the locations of observed sub-ice lakes.

4. Conclusion

In this paper, a preliminary version of a three-dimensional model of the Antarctic ice sheet is presented; it includes the fully coupled fields of velocity and temperature. The present version of the model does not include basal sliding and excludes the ice-shelf regions. The results shown are for an almost stationary state of the model, which is reached after 150 000 years. The calculations for the stationary state cannot yet be used to prove the usefulness of the model for describing the future evolution of the Antarctic ice sheet. However, a lower limit for the necessary degree of complexity of such a model can be given:

(i) For time-dependent calculations, starting with some initial conditions, the surface topography has to be known fairly accurately. Otherwise the vertical velocity component will include a large degree of error, at least within the first 1000 years of integration time. For longer time integrations, these initial disturbances will even out.

(ii) The model has to include the temperature dependence of the flow. This is already clear from the flow law (Equation (equ.2)), where the temperature-dependent factor A(T^’) varies by two orders of magnitude within the temperature range −30° to 0°C. Only those calculations which take into account the vertical structure of the temperature field will yield the observed accumulation rate of about 5 cm a⁻¹ in the interior of the Antarctic ice sheet that is necessary for a stationary state.

(iii) The three-dimensional temperature and velocity field is crucial also for estimating the temperature at the ice bottom, which decides whether basal sliding will occur. Since the bottom temperature also depends on the geothermal heat flux (which is only poorly known), it may also be necessary to calculate the temperature distribution within the bedrock, especially in transient situations (Reference MacAyeal and ThomasMacAyeal and Thomas 1980).

(iv) Some parameterized form of temperature-dependent basal sliding has to be added. The omission of sliding probably accounts for the areas of negative mass balance in Figure 3. In these regions, the bottom temperature is also at the pressure-melting point. Thus the next modeling step will be the use of a parameterized form of basal sliding (Reference Budd, Keage and BlundyBudd and others 1979), including a continuity equation for the basal water flow.

Acknowledgements

I should like to thank R. Calov for programming assistance, and E. Maier-Reimer for helpful discussions on the numerics.