Large-Scale Numerical Modelling of the Antarctic Ice Sheet

Abstract

A large-scale dynamic numerical model of the Antarctic ice sheet has been developed to study its present state of ice flow and mass balance as well as its response to long-term changes of climate or sea-level.

The flow of ice over a two-dimensional grid is determined from the ice thickness, the basal shear stress, the bedrock depth, and ice flow parameters derived from velocities of existing ice sheets. The change in ice thickness with time is governed by the continuity equation involving the ice flux divergence and the ice accumulation or ablation. At the ice sheet seaward boundary, a floating criterion and floating ice thinning rate apply. Bedrock depression with a time-delayed response dependent on the history of the ice load is also included.

A 61 × 61 point grid with 100 km spacing has been used to represent the ice-sheet surface, bedrock, and accumulation rate. The model has been used to simul a te the growth of the present ice sheet and i ts reaction to changes of sea-level, bedrock depression, accumulation rate, ice flow parameters, and the iceshelf thinning rate.

Preliminary results suggest that the present ice sheet is not in equilibrium but rather is still adjusting to changes of these parameters.

1. Objectives

The work presented here represents one step in a program to develop successively improved numerical models of the Antarctic ice sheet. The main objective at this stage is to develop a reasonable model of the whole Antarctic ice sheet which can be used as an interactive component within a global climate model to study long-term climatic change. Attention has been focused on the adequate representation of the large-scale growth and decay responses of the ice sheet to environmental changes.

A similar model for the North American region was developed by Budd and Smith (1981) to examine the orbital radiation theory of the ice ages. In that region, it was found that the southern limit of the ice sheet on land, and the areal extent, were largely governed by the ablation regimes. For the Antarctic ice sheet, net surface ablation is negligible within the total budget, and the bulk of the ice sheet generally has a positive surface net budget extending to the seaward edge. The main limits to the extent of the ice sheet appear to include the increasing depth below sea-level of the bedrock, causing the ice to float, strain thinning of floating ice shelves, basal riielting, and calving of icebergs.

Successful modelling of the ice sheet therefore requires an adequate representation of these coastal processes in so far as they affect the limit to the extent of the ice sheet.

Since, in the past, the climate has been generally colder rather than warmer, changing surface ablation rates could not be expected to have been a major factor in the changing extent of the Antarctic ice sheet. The most important factors seem to have been changing sea-level, changing bedrock depression (with delayed response), and, possibly, changing precipitation.

The objective so far has been to develop a model of the ice sheet which responds realistically to changing sea-level or precipitation and simulates the response of the bedrock to varying loads. Such a model will be useful for studying the role of Antarct ica in the changes of climate, ice volume, and sealevel over the last 120 ka. In addition, it is expected that a considerable amount can be learned from the model about the present state of the Antarctic ice sheet, the response times for changes, possible past regimes, and expected future trends.

2. Methodology

The model used in this study is similar to that used for the North American ice sheet model by Budd and Smith (1981). The bedrock elevation distribution with respect to sea-level is specified over a twodimensional grid b₀(x,y) for some initial time t0. The initial accumulation distribution A₀(x,y) is also specified over the grid and may change with time depending on other factors such as the elevation E. An initial ice thickness distribution Z₀ may or may not be specified.

for the north american ice sheet study of budd and smith (1981), a grid scale of 200 km was used and, although some important features (such as individual valleys) were smoothed out, the gross features were reasonably represented. for the current antarctic study, a 61 × 61 grid with a spacing of 100 km was adopted which captured the essential features of the antarctic region reaching out to beyond the continental shelf. again, although major outlet valley glaciers were smoothed out, the large-scale features were well portrayed.

On the large scale, the basal shear stress τ_b, is considered to be the primary driving force on the ice, and (with other stresses neglected) is given by the down-slope gravitational

(1)

stress where Z is the ice thickness, ρ is the ice density, g is the gravitational acceleration, and ᾱ is the largescale average surface slope. This slope is calculated from the x and y components (α_x,α_y) by

(2)

The surface velocity V is based on an empirical relation given by Budd and Smith (1981):

(3)

where n ~ 4 and k₁ ~ 0.03 bar⁻ⁿ a⁻¹, but can be considered to be dependent mainly on the temperature of the basal ice, and set accordingly.

For the North American ice sheet, Budd and Smith (1981) found the ice-sheet dimensions depended only slightly on the value of k₁. Hence, in this study, the k₁ value was kept constant for each calculation and sensitivity studies were carried out by using different values of k in separate calculations. A basal sliding option for the model as given by Budd and Smith (1981) is also available but comparative results including this option are not given here.

The change in ice thickness with time t at a point (x,y) is given by the continuity equation integrated through the ice thickness

(4)

where A is the net budget or accumulation rate in units of ice depth per unit area per time, V̄ is the average velocity through the ice thickness, and where the basal melt rate for grounded ice has been neglected.

A simple formulation of the variation of accumulation with time as the elevation changes based on the present accumulation rate distribution A₀ is prescribed as follows: for each 1 km increase in elevation E above the present surface E_p the accumulation rate is reduced by a factor of 2:

(5)

It is notable that this elevation relation also provides a useful approximation to the present accumulation distribution over the grounded part of the ice sheet starting from a uniform base of ~0.6 Mg m⁻².

The simple bedrock depression formulation given by Budd and Smith (1981) has also been used here, viz.,

(6)

where b^*t is the bed depression at time t, and Z_i is the ice thickness above that required for grounding the ice in sea water at i thousand years before t. This formula gives the equilibrium depression of one quarter of the ice load and has a delayed response of ~5 ka. If it were supposed that the present ice sheet is in isostatic equilibrium, then the initial bed for zero ice thickness could be taken as the present bed raised up by one quarter of the present ice load. Although the present bedrock may not be in equilibrium, this assumption has been used within the present study.

3. Input data

Bedrock data has been obtained from the Soviet 1:15 000 000 Antarctic bathymetric map of 1974, which shows both ice surface and bedrock elevation contours. With the aid of a digitizer, these contours were read into a computer and the data then interpolated onto the nodes of a square (x,y) grid representing an area of 6000 × 6000 km², extending to beyond the continental shelf. In order to concentrate on the large-scale features rather than small-scale variations, the bed was smoothed with a nine point filter. The resultant bed is shown in Figure 1(a). A similar process was carried out for the present ice-sheet surface. These distributions were then used to compute the present ice thickness distribution, the equilibrium bedrock depression, and, hence, the raised (relaxed) bed.

Fig.1. Input data within the 61 × 61 point 100 km spacing grid around the South Pole (+). The border markings are at 1 000 km spacing. (a) Antarctic bedrock. The smoothed bedrock distribution used in the model is shown with contour intervals at 500 m ranging from -1 500 m to +1 500 m. Tha regions above sea-level are shown shaded. The interior regions below 1 000 m are shown stippled. The edge of the present ice sheet is indicated by a dashed line. (b) Antarctic accumulation rate. Contours are shown for the gridded net accumulation rate data used in the model as the basic distribution for the present ice sheet. Closed contours of 0.6 m a⁻¹ (of ice) occur near the Antarctic Peninsula. The remaining contours decrease inland as follows: 0.5, 0.4, 0.3, 0.2, 0. 1, 0.05, 0.025 m a^{− 1}.

For the present accumulation rates, the map of Kotlyakov and others (1974) has been used with additions and modifications from Allison (1979), Budd and Young (1979), Clausen and others (1979), Raynaud and others (1979), Young (1979), Morgan and Jacka (1981).

These data were also digitized and interpolated onto the grid giving the distribution shown in Figure 1(b).

The two data sets form the basis for a series of experiments with the model aimed at simulating the present Antarctic ice sheet and its response to changing boundary conditions.

4. Seaward boundary

It has been found that an important factor controlling the extent of the ice sheet is the boundary condition applying at the transition zone between grounded and floating ice. Within the present model, floating is prescribed when the depth below sea-level reaches some specified fraction, say, λ of the ice thickness. For existing floating ice shelves, a value of λ between 0.9 and 0.84 is appropriate. For the results presented here, the value of λ is set at 0.9.

Once the ice becomes afloat, the basal shear stress falls to zero so that Equation (3) no longer controls the flow. For free-floating ice shelves, Weertman's (1957) ice-shelf creep formulation can be expected to apply, for which the vertical creep rate ε_z increases with ice thickness Z according to

(7)

where n is the power flow law index for ice creep.

For ice shelves bounded at their sides, it has been found that the strain-rate thinning actually decreases with increasing ice thickness away from the ice front, cf. Budd (1966), Budd and others (1967), Thomas 1973[a], 1973[b]. In this case, the ice flow is largely governed by the downward surface slope towards the ice front and the transverse strain-rate across the ice shelf, cf. Sanderson (1979), Budd and others (1982).

A thorough treatment of floating ice shelves in the model would therefore require horizontal stress integrations, both across and along the ice shelf. For the present context, we are primarily interested in the extent and shape of the grounded ice, and so the floating ice is only of concern here in so far as it is involved in the spreading and coalescence of the grounded ice. A second feature of the floating ice that is relevant to the extent of the ice sheet is the basal melt rate. Here again, the complex mass and heat exchanges between ice and ocean water have been omitted and, instead, the strain- and melt-rate thinning have been included in a prescribed thinning rate for floating ice. For most of the model runs discussed here, the standard thinning rate used is 1 × 10⁻² a⁻¹. Two examples of comparable runs with values of 2 × 10⁻³ and 2 × 10_-2 a⁻¹ respectively are also considered. Table I

Table I. Observed ice-shelf strain-rates

shows observed strainrates and ice thicknesses for existing ice shelves which are generally free-floating and exhibit the expected increasing strain-rate with ice thickness. The values in parentheses for the Amery Ice Shelf show the decrease of strain-rate inland for a bounded ice shelf. The values of the measured strain-rates largely span the range of 2 × 10⁻³ to 2 × 10⁻² a⁻¹ studied with the model here.

A theoretical discussion of the importance of iceshelf thinning rate on the extent of ice sheets with bedrock depression has been given by Weertman (1974).

5. Numerical experiments

Although the model described above is relatively simple, it can still be used to investigate a number of interesting questions concerning the Antarctic ice sheet, such as the following. What is the present state of the ice sheet and what is the equilibrium to which it may be tending? How does the ice sheet react to changes in such factors as: sea-level, accumulation rate, bedrock depression, the ice flow law parameters? How long does the ice sheet take to react to external changes and to approach a new equilibrium? What is the expected pattern of growth of the ice sheet from an initial bare rock state? What is the cause of the existing "non-equilibrium" state of the ice sheet? How does the past history of the Antarctic ice sheet fit into the picture of global changes of climate and sea-level?

In order to investigate these types of questions, a large number of computer model runs (-100) have been made with the model by setting some initial conditions (bedrock, accmulation rate, etc.) and letting the model run for a long period of simulated time (e.g. 100 ka). The time step generally used in these simulations is 100 a. Only a few of the results of this work can be reported here so those most relevant to the above questions have been selected for illustration in Table II

Table II. Summary of results of selected model runs^*

For example, no.4 was started with the present relaxed (i.e. adjusted for zero load) bedrock distribution and with no initial ice. The initial accumulation rate was taken as the present distribution which then decreased when the modelled ice elevation grew above the present day elevation according to Equation (5). Sea-level was set at the present level and the computation carried out for 150 ka of simulated time, by which stage a steady state was approximately achieved (cf. Figs. 2, 3, and 4). The final volume was then ~36 × 10⁶ km³, compared to ~23 × 10⁶ km³ for the present ice-sheet volume, as obtained from the digitization here which also agrees with the volume for the present grounded ice given by Bardin and Suyetova (1967)

Fig.2. Growth of the Antarctic ice sheet. The growth of the modelled ice sheet to steady state for run no.4 is illustrated by the surface elevation distribution at the following times: (a) 5 ka, (b) 10 ka, (c) 40 ka, and (d) 150 ka. Contour intervals are 500 m and range from 0 to 4 000 m. The most seaward interval represents floating ice. The highest elevations are shaded

Fig.3. Profiles during growth to equilibrium. The ice surface and bedrock elevation profiles along a line following the meridian 90°-270° are shown for the model sequence of Figure 2. The times are indicated in 10³ a.

Fig.4. Ice volume curves. The changes in ice volume during growth to equilibrium are shown for various model calculations with run numbers as given in Table II. In some cases the run time was extended to 200 ka to establish a clearer equilibrium.

The features of the other runs are indicated in the table. For example, no.3 was started with both the existing ice sheet and bedrock. In no.5, sealevel was maintained at 150 m below present sea-level, whereas in no.6 the bedrock was started 150 m lower than the relaxed bedrock distribution. This was equivalent to using a sea-level 150 m higher than present.

6. Discussion of Results

The growth of the modelled ice sheet volumes towards equilibrium is illustrated in Figure 4. The approach to equilibrium is approximately logarithmic with a time constant of ~25-50 ka. It can be seen that the final steady state is not reached for those cases starting with zero ice until about 150 ka. For those cases starting with the present accumulation and bedrock, the ice sheet grows larger than the present ice sheet, particularly in West Antarctica. Final steady-state cross-sections for different cases are shown in Figure 5.

Fig.5. Equilibrium profiles. The final steady-state surface and bedrock profiles for various run' (indicated by the run numbers from Table I) are shown for the meridians 90^°-270^* (in (a)) and 0^* -180^* (in (b)), compared with the profiles for the present ice sheet (P).

For the standard version (no.4), Figure 2 shows the sequence of surface elevation maps as the ice sheet grows from zero initial thickness to equilibrium. For the run (no.5) where the sea-level has been reduced from the present by 150 m, the ice sheet grows out further, particularly in the large iceshelf basins. This pattern is illustrated by the profiles shown in Figure 5. It is notable that the equilibrium bedrock is also deeper by ~ 200 m compared to the standard case, and is more than 400 m deeper over the Ross Sea basin where the new thick ice has advanced.

In order to test the effect of bedrock elevation changes within the model, run no.6 was carried out with the bedrock uniformly lower by 150 m relative to that of no.4. The result is illustrated in Figure 5(a) and (b) which show that the resultant ice-sheet shape is very much closer to that of the present ice sheet than those runs based on the present conditions. The effect of sea-level and bedrock variations on the extent of ice in the Ross Sea region is particularly dramatic. The low sea-level case results in the grounded ice reaching the edge of the continental shelf, as inferred for the ~18 ka BP maximum by Hughes and others (1981). The volume difference is similar to that given by Lingle and Clark (1979). In contrast, the lower bedrock condition of no.6 resulted in a grounded ice margin position, remarkably similar to that of the present, and with a similarly low West Antarctic ice sheet surface.

Although the varying sea-level studies are still in progress, it seems from the results so far that the present ice-sheet size and shape could have resulted from a previously larger ice sheet and consequently more depressed bedrock associated with the sea-level minimum around ~18 ka BP. Following the sea-level rise, the ice may have retreated to near its present situation, and since then the bedrock could have risen causing the ice sheet to grow and advance again. This pattern seems to support the inferences regarding the present advance of the ice shelf made by Thomas (1976) and Thomas and Bentley (1978) and the rise of the basement rock in the Ross Ice Shelf region by Greischar and Bentley (1980).

In all the cases described so far, the modelled ice sheets at equilibrium are higher than the present ice sheet. The question therefore arises as to how the present ice sheet came to be as low as it is. Possible causes include lower accumulation rates in the past and/or higher flow rates associated with higher basal ice temperature. In order to test these conjectures, and to test the sensitivity of the model to changes of the parameters A and k, the model was run firstly with the initial accumulation A reduced by one half (no.7), and secondly with k increased to a value more appropriate for temperate ice (no.8). In both cases, the equilibrium elevation and ice thickness are reduced, and the approach to equilibrium is particularly slower in the case of the reduced accumulation rate. In Figure 6, the model equilibrium surface elevation maps resulting for several of these different conditions are shown in comparison with that of the present ice sheet.

Fig.6. Equilibrium ice sheets. The surface elevation contours for the modelled steady-state ice sheets are compared with those for the present Antarctic ice sheet (b). The run numbers are from Table II as follows: (a) no.7, low accumulation; (c) no.8, high flow rate; (d) no.10, low ice-shelf thinning rate. The corresponding map for the "standard" run no.4, is given in Figure 2(d). Contours are plotted for elevations starting at zero and increasing at 500 m intervals.

Although much more work is needed to investigate the possible past history sequences of the Antarctic ice sheet, it appears that the present simple model can be used to investigate the factors which could lead to the present shape and size for the Antarctic ice sheet, such as changes in bedrock, sea-level, accumulation, and effective basal-ice temperature.

7. Future Developments

Studies underway with the existing model include the response of the ice sheet at maximum extent to a rise in sea-level, and the response of the ice sheet to sea-level changes resulting from the northern hemisphere ice volume changes as for example derived by Budd and Smith (1981).

In the future, a number of improvements to the model could be made, e.g. increased resolution, explicit modelling of ice shelves, incorporation of longitudinal and transverse strain-rates, introduction of an ice sliding law over the transition zone approaching floating (cf. Budd and others 1979), quantification of basal melt rates for ice shelves, and calculations calculations of vertical profiles of temperature and velocity through the ice sheet.

The calculation of temperatures could be carried out at each time step or after longer periods (as in Budd and others (1971) or Jenssen (1977)). The basal temperatures would then depend on the ice-sheet history as well as on its properties at any given time. At this stage, however, since the past history is not well known, it is considered a more clear-cut strategy to keep the flow parameter constant with a given run until the model's sensitivity to varying k values has been adequately determined.

Perhaps of more importance is that, as higher resolution or more detail is sought, greater use could be made of more recent or more accurate input data.

Even while these new developments are proceeding much can be learned using the present model for the interpretation of the past history of the ice sheet from ice-core data and other palaeo-indicators.