The Dynamics of the Antarctic Ice Sheet

Abstract

A three-dimensional dynamic, thermodynamic ice-sheet model has been developed to simulate the past, present, and future behaviour of the Antarctic ice sheet. The present ice velocities depend on the deep ice temperatures which in turn depend on the past changes of the ice sheet, including surface temperature, accumulation rate, and ice thickness. The basal temperatures are also strongly dependent on the geothermal heat flux. The model has therefore been used to study the effect on the basal temperatures, of changes to the geothermal heat flux, as well as the past changes of surface temperature and accumulation rate based on results obtained from the Vostok deep ice core. The model is also used to compute the distribution of surface velocity required to balance the present accumulation rate and the dynamics velocity based on the stress, temperature, and flow properties of ice, for the internal deformation, plus a component due to ice sliding. These velocities are compared to observed surface velocities in East Antarctica to assess the state of balance and the performance of the dynamics formulation.

Introduction

Prognostic dynamic models of the whole Antarctic ice sheet are needed for comprehensive global atmosphere-ice-ocean models to study long-term changes of climate and sea-level. A large part of the coastal ice flux from the Antarctic is channelled into narrow deep channels forming relatively fast-flowing glaciers or ice streams. Earlier dynamic modelling of the Antarctic ice sheet by Budd and Smith (1982) used a resolution of 100 km which was not adequate to simulate the differentiation into the high-speed flow of the outlet glaciers. It was therefore decided to increase the resolution to 20 km based on the digitization of the surface and bedrock from the maps of the SPRI Map folio by Drewry (1983). Although these maps do not have accurate data over the whole of the continent, they provide a useful basis to start from, and the regions where the data are reliable are indicated by the coverage of flight lines and traverses on the maps. A first use of this digitized data set for computing the ice velocities which would be required to balance the present surface snow-accumulation regime (the “balance velocities”) was described by Budd and Smith (1985). Further uses of the high-resolution data for thermodynamic and dynamic modelling have been given by Budd and others (1984, 1987) and Radok and others (1986). This modelling computes the “dynamics velocities” from the stress and temperature distributions and the flow properties of ice. To allow more rapid computations of the full three-dimensional representation of the temperature and velocity distribution, the model has been reformulated by Jenssen and has been used to simulate the present ice sheet and its past and future changes.

A useful test for the performance of the model is the degree to which it can simulate the large number of surface velocities of the ice sheet now becoming available, e.g. from Hamley and others (1985) and Hamley (unpublished). The model computations of surface-dynamics velocity require deep ice temperatures which depend on the past history of the ice sheet. The model has therefore been run forward from 500 000 B.P. with a variety of past histories to study the impact on the basal temperatures. The results from the deep Vostok ice core have been used as a guide to the past changes in the interior of the ice sheet over the last 160 000 years for such variables as the surface temperature or accumulation rate, cf. Lorius and others (1985), Jouzel and others (1987), and Yiou and others (1985).

The purpose of the present paper is to examine the computed surface velocities in relation to the velocities required for balance (with the present accumulation) and also to the observed velocities, in order to provide an indication of the present state of balance and the degree to which the observed velocities can be simulated by the dynamics. The sensitivity of the dynamics to the major unknowns such as the geothermal flux, the past accumulation rates, and the ice-flow relation is also examined.

Ice sliding also becomes important, particularly near the coast, and the observed flow of the outlet glaciers provides a check on the model ice-sliding formulation as given by Mclnnes and Budd (1984) and Budd and others (1987) based on the empirical studies of Budd and others (1979).

The details of earlier versions of the model have been described by Radok and others (1982, 1986). Some of the essential components are described in the Appendices. Only a brief outline is given here to indicate the basic principles of the model and the major changes from previous versions.

Outline of the Model

The model covers a domain of 281 = 281 grid points with 20 km spacing over the Antarctic region with 31 points in the vertical through the ice and the possibility for an additional 20 points in the bedrock reaching several kilometres depth. The equations for heat conduction and ice velocity are solved together, as given in Appendix I, allowing for internal heating varying with depth in the ice sheet. The horizontal shear stress (t _xz) is taken as the only significant stress for ice deformation and heating except for basal sliding (_b) in which case the effective normal stress is important and the frictional heating occurs at the basal interface. Vertical velocities are computed from the horizontal velocity profile (as given in Appendix I), and both vertical and horizontal advection are computed explicitly via finite differences from surrounding points. Space is not adequate here to present a complete coverage of the model but further details are given in the Appendices. The work is an advance from the earlier three-dimensional modelling of Jenssen and others (1985). Budd and Jenssen (1987), and Radok and others (1987), in that flow lines are no longer explicitly required (but can be derived) since the model reaches a solution by iteration over the whole domain. A steady-state solution may be independently obtained by running the model for a long time (500 ka) with constant boundary conditions. Non-steady-state conditions can be derived from an initial state by computing forward in time with some prescribed varying boundary conditions.

Examples of the use of the model have been to run steady-state solutions for the present regime with different geothermal gradients, different past accumulation regimes, different ice-flow parameters, and different sea-level or bedrock conditions. Only a few examples are given here to examine the model's simulation of the present regime. By comparing both the model-computed balance velocity and dynamics velocity with observations, it should be possible to assess the degree to which steady state currently holds and the degree to which the computed velocities match reality. To compute dynamics velocities for the present regime, we first compute temperature distributions from prescribed computed balance velocities. Similarly, the temperature distribution has also been computed from lower accumulation rates in the past, as indicated from the Vostok ice-core results.

Assessment of Output

To begin with. Figure 1 shows the computed balance velocities based on the present regime (involving the surface elevations, accumulation rate and ice thickness). It should be noted that a change of a constant factor to the accumulation rate would simply result in the same constant-factor scale change in the balance velocities. In this regard, lower accumulation during the ice age would correspond to proportionally lower velocities. The Vostok ice-core data indicate that past accumulation rates have ranged between the present levels and about half those levels. The model results are therefore initially obtained for the present configuration with prescribed velocities first corresponding to the present values and then to half the present values. This allows us to examine the effects of the velocity distribution on the temperature. Similarly, the model can be computed with different prescribed accumulations. The model can also be run with a time series of the input variables changing as prescribed over the entire surface which is an improvement on the flow-line approach of Budd and Young (1983).

Fig. 1. Balance velocities. The distribution of balance velocities is shown as required for steady state with the present accumulation, ice thickness, and elevation contours.

In regard to the basal ice temperatures which have such a strong impact on the computed dynamics velocities, the prescribed geothermal heat flux (ϕ _G) is of great importance. By examining the model results for a range of different values of this gradient, it appears that a reasonable central value to use is about ϕ _G = 5.17 × 10⁻² W m⁻² [1.23 μcal cm⁻² s⁻¹]. This value converts to a basal temperature gradient (depending on the temperature) of roughly γ_G = 2.25°C/10² m (see Table I). A geothermal gradient of this amount has also been found to provide a close fit to the observed deep temperature profile at Vostok (Jenssen and Budd, in press). Figure 2 shows the pattern of

Table I Geothermal heat flux and temperature gradients in ice

Ice conductivity, at 0°C = 2.24 W m⁻¹ °C⁻¹.

Geothermal heat flux unit 1 μcal cm⁻² s⁻¹ = 4.187 × 10⁻² W m⁻².

computed basal temperatures derived from the prescribed geothermal fluxes ϕ _G = 4.62, 5.17, 5.75 × 10⁻²W m⁻², based on the steady state under the present regime. Similar calculations have been carried out with the accumulation rate half its present value or the velocities half the balance velocities or both of these changes together, as shown in Figure 2(d).

The dynamics velocity _d is computed from the internal deformation (V _i) plus the sliding velocity (V _s)

The internal velocity is computed from the integration of the shear strain-rate (

) through the ice related to the shear stress (τ_xy) and temperature (θ). The flow relation used is of the form

where the temperature-dependent factor (_θ) is given by the table determined for tertiary flow as reported by Budd and Jacka (1987, in press). The numerical values are similar to those given by Paterson and Budd (1982) but increased by about a factor of 3 for tertiary flow rates.

The sliding velocity (_s) is taken as the general form derived empirically by Budd and others (1979), with an additional factor for the effect of temperature:

where is the effective normal stress (above buoyancy) relative to pressure-melting point, and are simple monotonie functions, θ is the temperature, and ν ≈ 0.1°C⁻¹. The particular form used here is given in Appendix II, Equation (A₂-4).

The total dynamics velocity is shown in Figure 3 and the sliding component in Figure 4. The sliding velocity is generally small for temperatures below pressure melting. The sliding velocity also increases strongly towards the coast as the normal stress decreases and thus allows high-speed flow in the ice streams even though the shear stress may decrease.

A comparison of the total dynamics velocity with the balance velocity shows general broad similarity but with some notable regions of difference. The largest such region is east of the Filchner Ice Shelf where the dynamics velocities are much below the balance velocities. The discrepancy is more than can be reasonably expected from the uncertainty in the ice-flow parameters or the temperature calculations. It is suggested that either the input data (from the SPRI maps where coverage is sparse) are grossly in error or the ice sheet there is grossly out of balance and may be still building up in response to the post-ice-age accumulation increases.

Fig. 2. Basal temperature. The computed basal temperatures are shown derived for steady state based on the present regime for the following geothermal flux in mW m⁻²: (a) ϕ _G = 46.2, (b) ϕ _G = 51.7, (c) ϕ _G = 57.5. In (d) the corresponding results are for ϕ _G = 51.7 but for the accumulation rate and balance velocity distributions at one-half of their present-day values.

Fig. 3. Dynamics velocity. The velocity distribution computed from the temperature distribution obtained using ϕ _G = 50.4 mW m⁻² and the flow and sliding relations given in Appendix I. There are broad similarities with the balance velocities but notable regions of difference such as up-stream of the Shackleton Range in Dronning Maud Land.

Fig. 4. Sliding velocity. The distribution of sliding velocity computed from the shear and effective normal stresses shows low values inland but high values near the coast particularly in the region of the deep ice streams and outlet glaciers.

Fig. 5. Comparison with observations in East Antarctica. For the region of the IAGP surveys, the observed velocities from Hamley and others are shown in (b). The computed balance velocities are shown in (c) and the dynamics velocities in (d). The model elevations are shown in (a).

Fig. 6. Balance and dynamics velocities versus observed. The balance velocities are plotted against the observed velocities from Figure 5 in (a). The corresponding dynamics velocities are plotted against the observed in (b). Although considerable scatter exists, both show reasonable trends which can be used to assess balance and the effective flow parameters.

Finally, in Figure 5 a comparison is made with the East Antarctic region covered by the IAGP (International Antarctic Glaciology Project). Here, the observed velocities from Budd and Young (1979), Hamley and others (1985), and Hamley (unpublished) in Figure 5a are compared with the balance velocities (Figs 5b and 6a) and the computed dynamics velocities in Figures 5c and 6b. Although there is considerable variation between the results, the general trend of increase of velocity towards the coast, particularly over the deep bedrock regions, is common to all three.

This is not a trivial result because it involves stresses which on some flow lines increase from 0.5 bar to over 1.5 bar then decrease to 0.5 bar while the velocity continues to increase (due to sliding) near the coast.

Since there is inherent error or uncertainty in each of the velocity distributions, statistical regression is useful to compare the observations with both the balance and dynamics velocities. The first conclusion is similar to that expressed by Hamley and others (1985) that the observations in this region do not differ substantially from balance. This involves also some uncertainty in the comparison for the ratio of the average column velocity to surface velocity which is expected (from the modelling) to be generally in the range of 80-90% until the high-sliding zones are reached where it may approach 100%. The second conclusion is that the present dynamics velocity formulation gives a reasonable representation of the observed velocities. The regions where large differences appear between balance and dynamics velocities may well be areas in which the ice sheet is not in balance. On the other hand, in those regions, there may be errors in the data used as input to the model which affect the dynamics velocities. Observed velocities in these regions, as well as improved observations for the input data, would be particularly valuable for clarifying both the state of balance and the dynamics.

In the interim, the comparison between observed and computed dynamics velocities provides some encouragement to use the computed dynamics velocities interactively for the time-varying past conditions through the ice-age changes.

Appendix I

Deformation Velocity and Flow Relation

For a coordinate system with axes horizontal in the line of flow, vertical, and across the flow, it is assumed for the large horizontal scales of the grid that the flow is governed primarily by the horizontal shear stress τ_xz, given by

(A₁-1)

where ρ is the ice density, is the gravitational

Table II Ice-flow rates

at minimum for isotropic ice and stress τ₀ = 1 bar (octahedral)

Note: The temperature used for flow in the model is the temperature difference from pressure-melting point at each depth.

acceleration, α is the surface slope, and is the depth below the surface.

The flow relation used is of the form which relates the shear strain-rate at depth z,

, to the shear stress by

(A₁-2)

where = 3 and _θ is dependent on the temperature and ice properties.

For this context, the dependence on temperature (θ) is given by the data of Budd and Jacka (in press) as shown in Table II, shown for octahedral shear stress τ₀ = 1 bar, and strain-rate

, from laboratory measurements.

The relations between the octahedral and horizontal shear stress and strain-rate are

(A₁-3)

Table II gives the strain-rate at minimum for isotropic ice tests. A further factor ₁ is used with the flow law to allow for anisotropy, tertiary flow, other factors which may affect the flow of ice, or to match with observed surface velocities.

The velocity depth gradient d/d is then given by

(A₁-4)

The horizontal velocity at depth is given by

(A₁-5)

where _s is the surface velocity minus the basal sliding and _θ is determined at each depth using Table II from the temperature profile from the heat conduction equation as given in Appendix II.

The Vertical Velocity Profile

Following Budd and others (1970), the usual vertical coordinate, , is replaced by the relative coordinate ξ = (E – z)/Z, where Z is the ice thickness and is the ice-surface elevation. The continuity equation for an ice element is then

(a)

where

is the horizontal velocity, and (= dξ/dt) is the relative vertical velocity.

The continuity equation for the column as a whole is

(b)

where

is the average horizontal velocity in the column, and is the surface-accumulation rate. It is assumed that

(c)

where ψ is a function of ξ (or ) which may be computed from the flow law used. Use of Equations (a) and (c), with substitution into Equation (b) yields, after some manipulation,

Replacing ξ by (1 - ζ) and integrating from bedrock (ζ = 0) to any general ζ gives, since ζ = 0 at the ice base,

(d)

where

(e)

Appendix II

Heat-Conduction Equation Including Internal Deformation Heating and Advection

The basic thermodynamic equation used in the numerical model is

(A₂-1)

where θ is the temperature (°C), is the ice thickness (m),

is the horizontal velocity (m a⁻¹), is the heat generation from deformation as given below, κ is the thermal diffusivity, κ′ is essentially the derivative of κ with respect to temperature (m² a⁻¹°C⁻¹), ξ is a relative vertical coordinate with 0 at the surface and 1 at the base, and is the relative vertical velocity at ξ (w = dξ/dt).

The quantities Z and A are not functions of depth but other quantities are. In addition, κ, κ′, q, and are functions of temperature so that Equation (A₂-l) is highly non-linear.

The horizontal velocity u is computed from the flow law as given in Appendix I, involving the temperature and stress at any point. The average velocity through the column is u .

The heating q is considered in two parts: q _i which is internal deformational heating, and q _b which is basal frictional heating from sliding.

(A₂-2)

where is the thermal conductivity (kg m a⁻³ °C⁻¹) and

are computed from the ice-flow law as given in Appendix I.

For the sliding component,

(A₂-3)

where _b is the basal sliding velocity. Note is a function of the base temperature.

The computation of basal sliding is based on the empirical sliding studies of Budd and others (1979) and the comparison of dynamics and balance velocities in West Antarctica by Budd and others (1984).

The current form used is given by

(A₂-4)

where _* is the ice thickness above buoyancy (i.e. that required for floating) and _{2 3}, and ₄ are constant, ₂ = 1.5 × 10⁸ m³ year¹ bar¹, ₃ = 50 m², ₄ = 0.0035 m⁻¹, θ² is the temperature relative to pressure-melting point at ν = 0.1 °C⁻¹.

For cases in which the model is run with prescribed (surface or balance) velocities, the total velocity is normalized to the prescribed velocities but the fractions of sliding and internal deformation are determined by the above relations.