Numerical Models of Ice-Shelf Flow: Ideal/Real

Abstract

We present model calculations that describe the flow of ice shelves of different configurations. We consider “ideal” models with well-defined boundary conditions and simple geometry in order to explore the response of an ice shelf to transient variations in ice-stream input. Gradually increasing the complexity of these simple models allows a better understanding of ice-shelf behavior without the complications that would arise in considering natural ice shelves. We find that the dissipation of ice-thickness variations caused by ice-stream transience is strongly influenced by ice rheology. The presence of an ice rise significantly alters the velocity field of the adjacent ice, when changes in ice-stream input occur. With models of “real” ice shelves, we demonstrate the ability of numerical models to test successfully working hypotheses on ice-shelf thickness distributions. Ice velocities, obtained by diagnostic models of Filchner–Ronne Ice Shelf that use different ice-thickness distributions, are compared with measured ice velocities. This comparison demonstrates that the model employing regions of thin ice in the central part of the ice shelf yields velocities significantly different from the field data. We therefore conclude that zones of thin ice on Filchner–Ronne ice Shelf are unlikely. This conclusion has recently been confirmed by field measurements.

1. Introduction

The need to understand the dynamics and mass balance of the large ice shelves of Antarctica (Ross, Filchner–Ronne, Amery ice shelves) arises because they:

drain large sectors of the Antarctic inland ice, thus reflect the overall behavior of the ice sheet as a whole, and as such serve as sensitive indicators of large-scale and long-term changes that might primarily be induced by global climatic changes.

However, identification of climatic signals in the flow of an ice shelf are complicated by other changes in boundary conditions, which might completely mask the climatically induced alterations in the state of the ice shelf. Ice-stream transience, i.e. the temporal variation in the amount of ice discharge through ice streams into an ice shelf such as observed on the Ross Ice Shelf (Jezek, 1984; Shabtaie and Bentley, 1987), represent a prominent example for such short-term changes in boundary conditions (cf. Lange and MacAyeal, 1988a; MacAyeal and Barcilon, 1988; MacAyeal and Lange, 1988). Thus, besides assessing the role of long-term variations on ice-shelf flow, there clearly exists the need to separate these effects from those caused by short-term changes in boundary conditions.

Because of their complexity, adequate representation of natural ice shelves requires the use of numerical models. In order for models to produce realistic answers though, appropriate boundary conditions have to be known or estimated. However, sufficient boundary conditions are seldom available when features of the size of the major ice shelves are studied. Thus, being forced to use numerical models as the only practical tool to describe large-scale ice-shelf flow, the absence of sufficient boundary conditions might render the modelling results unreliable if not useless.

Given this dilemma, we decided to take a two-fold approach in modelling ice-shelf flow. As a first approach, we construct models of relatively simple, well-defined, “ideal” ice shelves. Simulating natural ice-shelf flow usually leads to complications that are a result of non-linear feed-backs between different elements of the model, which usually cannot be traced back to a specific cause. Using our simple models, we can gradually increase their complexity, while subjecting them to identical variations in specific boundary conditions. Therefore, at each step in this procedure, we are able to assess the exact response of our model ice shelves. This enables us to get an insight into the nature and magnitude of changes that might occur in, for example, ice-velocity or ice-thickness fields in response to changing boundary conditions.

As the first problem, we addressed ice-stream transience in simple one or two ice-stream models. Having acquired a better understanding as to what likely changes in the major properties of an ice shelf will be, when subjected to alternating ice-stream inflow, we will have more confidence in simulating similar processes in models of natural ice shelves.

As a second approach, whenever sufficient boundary conditions are available, we will use models of “real” ice shelves to address specific questions of ice-shelf flow. The ice-thickness distribution of the central part of the Filchner–Ronne Ice Shelf is one of the problems we have tackled with the aid of “real” models. Here, the basic question is whether extended zones of reduced ice thickness, as seen in the maps of Robin and others (1983) and Crabtree and Doake (1986), are real features or based on misinterpretations of radar measurements, as suggested by Thyssen (1988).

While some of our model studies have already been published (Lange and MacAyeal, 1986; 1988a, b; MacAyeal and Barcilon, 1988; MacAyeal and Lange, 1988), the primary aim of this paper is to demonstrate the potential of our two-fold approach in simulating ice-shelf flow. Rather than giving a complete description of either modelling technique, we present selected results of “ideal” and “real” models.

2. Methods, Numerical Techniques, and Boundary Conditions

2.1. General

We will only give a short outline of basic methods and techniques used in our model calculations; details can be found in MacAyeal (1987), MacAyeal and Thomas (1982, 1986), and MacAyeal and Lange (1988).

Two model classes are distinguished: diagnostic models that yield “snap-shot” solutions of ice-shelf flow under prescribed boundary conditions and prognostic models, which treat the temporal evolution of an ice shelf. In both classes, we solve the momentum-balance equation

(1)

using a proper constitutive relation (Paterson, 1981)

(2)

where Τ and

are the stress and strain-rate tensors, ρ is ice density, g is the gravitational acceleration, e_z is the unit vector orthogonal to the geoid, τ is the second invariant of the stress tensor, Β the depth-integrated stiffness constant, and n is a constant (taken as 3 in most cases; see below). Mass balance is expressed as

(3)

where H is ice thickness, t is time, u is ice velocity, Ȧ is the surface accumulation or ablation rate, and Ṁ is the bottom-melting or freezing rate (∇_H is the horizontal divergence operator). Ice velocities u are assumed to be independent of depth (no shear) and thermodynamics are not treated in the prognostic models. Instead, a constant stiffness parameter Β is used. The nominal value of Β = 1.6 × 10⁸ Pa s^1/3 is derived from the temperature profile at J9 on the Ross Ice Shelf, probably a representative value for Antarctic ice shelves. In order to study the effects of a uniform warming of the ice shelf, we used a value of 0.8 × 10⁸ Pa s^1/3 in some of our models. While most models are treated with non-linear rheology (i.e. n > 1), we also simulated Newtonian rheology (i.e. n = 1) in some models.

In the diagnostic models, we compute temperature profiles based on the analytical solution of Crary (1961) for a ten-layer ice shelf at each grid cell. The temperature profile is used to derive a weighted stiffness parameter Β for each cell.

Ice input in all models is via ice streams and constant snow accumulation, but bottom freezing is ignored. Iceberg calving (not explicitly modelled here) is assumed to balance exactly the ice flow towards the edge.

Fig. 1. Configurations of the one ice-stream (a) and two ice-stream (b) models. In (b), the entire numerical domain, including a 10 km boundary, which is required in our calculations, is shown and will be displayed in subsequent figures. Note the non-symmetric position of the ice-stream entrances in (b) and the position of the ice rise (size = 20 km × 20 km).

2.2. “Ideal” models

2.2.1. Model configurations and boundary conditions

Figure 1 gives the configuration of the “ideal” one and two ice-stream models (for further details, see Table I). The ice-shelf domain is represented by rectangular grid elements of constant sizes for each model (Table I). As for boundary conditions, ice is treated as stagnant along the boundary between ice shelf and surrounding inland ice (no-slip condition), and ice-stream input is varied between limits, which are regarded as minimum and maximum cases, i.e. between 0 and 40 km³/a, respectively (this corresponds to ice velocities at the ice-stream entrance of 0 and 1000 m/a, respectively). This will result in ice-thickness and ice-velocity variations of the ice shelf that are also bounded. The magnitude of the thickness variations defines what we call the ice-thickness anomaly envelope (see below).

In the two ice-stream models, we include a central ice rise (Fig. 1). Ice is treated as being stagnant over the area of the ice rise.

Table 1 Configuration and boundary conditions for “Ideal” models

*Ice-stream discharges and inflow velocities are given for the minimum and maximum cases, respectively; for the two ice-stream model, the values apply to both ice streams simultaneously.

†Ȧ and Ṁ are the accumulation rate and the bottom-melting rate of Equation (3), respectively (in ice equivalent).

‡n and Β are the constants in the constitutive relation (Equation (2)).

2.2.2. Ice-thickness anomaly envelopes

Ice-thickness anomalies reflect the instantaneous ice-thickness deviation from the long-term mean. Since the ice-thickness variations will be produced by ice-stream discharge fluctuations (see above), ice-thickness anomalies will fall within an envelope determined by the maximum and minimum thickness attainable from ice-stream forcing alone. Whenever particular details of ice-stream transience are unknown, the spatial distribution of this envelope provides the best available measure of ice-shelf response. Thus, regions sensitive solely to atmospheric and oceanic forcing can be identified.

We define the absolute and relative ice-thickness anomaly envelopes Η_a and H_r, respectively;

(4)

(5)

where

(6)

Here, H₊(x,y) and H₋(x,y) are the ice thicknesses at each grid point, corresponding to maximum and minimum ice-stream influx, respectively. The ice thicknesses used in these calculations are obtained with the aid of a prognostic model, by starting from an arbitrary configuration, with a prescribed ice-stream influx, after it has reached equilibrium (usually after about 1500 a; see MacAyeal and Lange, 1988). The ice-thickness anomaly envelopes are evaluated for each grid point and mapped for the ice shelf under consideration.

2.3. “Real” models

We test the hypotheses on the ice-thickness distribution of the Filchner–Ronne Ice Shelf. In model 1, we use ice thicknesses given by Robin and others (1983). In model 2, we arbitrarily thicken the thin central area of the Robin and others’ map by adjusting the thicknesses in the thin zone to those of the adjacent thicker lobes of the ice shelf (see Lange and MacAyeal (1986) for further details and their figure 1 for details of the model configuration). The ice shelf is represented by rectangular and triangular elements of constant size (length = 1.6 × 10³ m).

Ice input occurs through the ice streams entering Filchner-Ronne Ice Shelf with velocities of 400 m/a (Lange and MacAyeal, 1986) and by constant snow accumulation (0.15 m/a ice equivalent). Other boundary conditions are similar to those used in the “ideal” models.

Since ice thicknesses determine the major driving force of ice-shelf flow, different velocity fields emerge from the two ice-thickness distributions. The corresponding ice velocities are then compared with existing field data. We also compare the resulting strain-rate fields of the two models with a map of surface features (e.g. crevasses; Swithinbank and others, 1988). The ice thicknesses, together with these features, will help in selecting the most likely thickness distribution of Filchner–Ronne Ice Shelf.

3. Results

3.1. “Ideal” models

We will present only results of models which involve ice-stream transience. However, the temporal variation will only be discussed in terms of ice-thickness anomaly envelopes for cases of maximum versus minimum ice-stream input (see above; for results on models with time-varying ice-stream input, see MacAyeal and Lange, 1988; paper in preparation by M.A. Lange and D.R. MacAyeal).

Fig. 2. Absolute (lower panels, in m) and relative (upper panels) ice-thickness anomaly envelopes for the one ice-stream models. 1 and 2 are results of models with non-linear rheology, and 3 represents results of a model with Newtonian rheology.

Figure 2 gives absolute and relative ice-thickness anomaly envelopes for models 11–3 of Table I. Because of symmetry, we use the upper half of the diagram for relative anomaly envelopes and the lower part for absolute anomaly envelopes, respectively. As can be seen, ice rheology strongly influences the response of an ice shelf to varying ice-stream input. While in model I1 a substantial fraction of the ice shelf shows significant changes in absolute ice thicknesses (as seen in the absolute anomaly envelopes), in model I2, in which a reduced stiffness parameter is used, the region affected by these variations is more restricted to the first 100 km down-stream of the ice-stream entrance in model I2. Hence, the effects of varying ice-stream input are removed more readily in the ductile ice than in model I1. However, because of overall thinner ice in model I2, which is a result of the reduced stiffness parameter, relative anomaly envelopes are remarkably similar in both models. Thus, in both cases, maximum discharge through the ice stream will result in a thickening by a factor of 1.67 for a substantial part of the respective ice shelf.

These results also imply that the effects of ice-stream transience on absolute ice thicknesses will be less significant than changes in the overall condition of the ice shelf. Thus, for these models, the climatic signal dominates any signals caused by ice-stream transience.

An ice shelf obeying a Newtonian flow law (cf. Doake and Wolff, 1985) reacts drastically differently to ice-stream transience than ice shelves complying with a power flow law with n > 1. Constant viscosities lead to enhanced changes in ice thicknesses. These changes are significantly different from those of models I1 and I2, both in magnitude and spatial distribution of the ice-thickness anomaly envelopes. Increasing the ice-stream input to its maximum value leads to three times thicker ice for a substantial part of the ice shelf (i.e. for regions with H_r > 1.0). Maximum changes now occur further down-stream from the ice-stream entrance, approximately in the mid-section of the ice shelf. Thus, while ice-thickness anomalies are dissipated fairly close to the ice-stream inlet in models I1 and I2, they maintain a high value for most of the ice shelf when linear rheology is assumed. This is probably caused by the high strain-rates of model I3, which result from higher friction along the margins of the ice shelf. Thus, in order to maintain the balance between ice-stream inflow and calving rates at the ice-shelf front, greater ice thicknesses are required in its mid-section.

These results demonstrate that ice-thickness changes due to varying ice-stream input may provide a clue as to the

Fig. 3. Relative ice-thickness anomaly maps for the two ice-stream models. Both models use non-linear constitutive relations. Model 2 uses a reduced stiffness parameter B. The asymmetry in the shapes of the isolines is caused by the non-symmetric positions of the ice-stream entrances.

dominant rheology of the ice shelf under consideration. In contrast to the results for models with non-linear rheology, the results for models with linear rheology also imply that ice-stream transience might indeed mask other changes, i.e. of large-scale oceanic and/or atmospheric nature.

Figure 3 gives relative ice-thickness anomalies for models II1 and II2 (because of the non-symmetric distribution of the ice-stream inlets relative to the axis of the ice shelf, the entire map of ice-thickness anomalies is shown). Reduction of the stiffness parameter by one-half (corresponding to an overall warming of 12°C for the whole ice shelf) leads to similar relative ice-thickness anomalies for both models. However, the absolute anomaly envelopes are different for both cases, similar to the results for the one ice-stream models (cf. Fig. 2).

As expected, ice velocities change significantly in response to ice-stream discharges. This is most clearly seen directly at the entrance of the ice streams but also throughout the whole ice shelf. Reducing the stiffness parameter leads to smaller absolute velocity changes compared to the nominal case. However, the ice-stream signal dominates any changes caused by the altered rheology (not shown here; see paper in preparation by M.A. Lange and D.R. MacAyeal). Thus, based on observations of ice velocities, ice-stream transience will result in changes that might exceed those seen when overall boundary conditions (such as ice temperatures) are altered.

Fig. 4. Ice velocities of the nominal two ice-stream model. Panel (a) displays the velocities for the minimum-discharge case; panel (b) displays the velocities for the maximum-discharge case.

Fig. 5. Predicted (numbered isolines) and measured (numbered arrows) ice velocities (in m/a) of Filchner–Ronne Ice Shelf. Panel (a) displays results of model 1, where reduced thicknesses in the central part of the ice shelf (approximately between 57° and 50°W and between 72° and 83°S) were used. Model 2 (panel (b)) uses an ice-thickness distribution, where the thin central zone is artificially filled by meteoric ice.

The presence of an ice rise has relatively modest effects on the ice-thickness anomaly envelopes. However, if we analyze the velocity fields (note that we present absolute ice velocities and not anomaly envelopes as in the case of ice thicknesses), the ice rise is clearly seen, once ice-stream input reaches its maximum. Figure 4a and b gives results for model II1 for minimum and maximum ice-stream discharge, respectively. If we were to remove the ice rise in models II, the velocity distribution of the central part of the ice shelf as seen in Figure 4b would resemble that of the minimum-discharge case (Figure 4a; however, not with respect to absolute values of ice velocity).

3.2. “Real” models

Figure 5a, b gives computed ice-velocity contours for the Filchner–Ronne Ice Shelf based on our two ice-thickness distributions used in the computations, i.e. models 1 and 2, respectively. As can be seen, the velocity distribution is similar in both models. However, the zone of thin ice modelled in run 1 results in significantly smaller velocities in the near-ice-edge region of the ice shelf.

Also shown in Figure 5 are velocity vectors, which were obtained by comparing satellite positions of selected points on the ice shelf at yearly intervals (Kock and Wiegand, 1986). They provide an independent boundary condition against which our predicted ice velocities have to be compared. It is obvious that the velocity distribution of model 2 yields the better match with the measured velocities. Thus, based on this comparison alone, it appears that a zone of thin ice in the center of Filchner–Ronne Ice Shelf is unlikely. This agrees well with our earlier conclusion (Lange and MacAyeal, 1986), which was based on a comparison between predicted velocities and ice-edge velocities that were derived from direct measurements and indirectly from comparisons of geographical positions of the ice-shelf edge in different years (cf. Lange and Kohnen, 1985; Lange, 1987).

Computed strain-rates of model 1 suggest a regime of compressional flow, with principal axes oriented roughly parallel to the ice front in the region of thin ice. This feature is completely absent in the model 2 strain-rates (see figures 5 and 6 of Lange and MacAyeal, 1986). Compression in the central zone should lead to buckling and/or crevassing in or along the margin of this zone. However, most of the crevassing, as seen in the map of Swithinbank and others (1988), is clearly related to ice-stream flow regimes and does not seem to comply with the strain-rate field of model 1. Thus, this result also suggests that model 1 ice thicknesses are incorrect.

Finally, Engelhardt and Determann (1987) reported a thickness measurement through a bore hole at 78°19′S and 57°08′E on Filchner–Ronne Ice Shelf. This position corresponds to a region where Crabtree and Doake (1986) found ice thicknesses of about 200 m. However, the measured ice thickness is 465 m. This provides direct evidence for the fact that maps of the thin central zone of Filchner–Ronne Ice Shelf (Robin and others, 1983; Crabtree and Doake, 1986) are probably based on erroneous interpretations of the radar records, as suggested by Thyssen (1988). Engelhardt and Determann (1987) proposed that this zone of the ice shelf might be underlain by non-meteoric ice of as yet unspecified origin. They hypothesized that this ice might represent frazil (i.e. saline) ice formed in the water column beneath the ice shelf that gradually accumulates down-stream of Doake Ice Rumples. Lange and MacAyeal (1988b), on the basis of prognostic models of the central part of Filchner–Ronne Ice Shelf, have demonstrated that such a process provides the most likely mechanism for the accumulation of bottom ice beneath the ice shelf.

4. Conclusions

The results presented here demonstrate the potential of numerical models in assessing ice-shelf dynamics. Models of “ideal” configuration and well-defined boundary conditions yield a valuable insight into the response characteristics of an ice-shelf to ice-stream transience. In this study, we particularly describe differences in ice-thickness anomaly envelopes for ice shelves of different rheology but otherwise identical boundary conditions. Our approach enables us to separate ice-thickness changes caused by variation in ice-stream discharge alone from changes that might arise due to the feed-back of different elements in a model of a “real” ice shelf. Such feed-backs might arise from the complicated geometry and spatially inhomogeneous properties of a natural ice shelf.

We have shown that ice rheology is important in determining the magnitude and distribution of ice-thickness anomalies on the ice shelf in response to ice-stream transience. The most significant differences are seen when models with linear rheology are compared with those of models obeying a non-linear flow law. Ice-thickness anomaly envelopes of non-linearly flowing ice shelves will generally be restricted to a region close to the ice-stream outlet. Reducing the stiffness parameter in these models will only reduce the magnitude of the absolute anomaly envelopes but will result in relative anomaly envelopes similar to the nominal model. Hence, it is expected that these ice shelves will be more strongly affected by overall changes in the state of the ice shelf, such as their mean temperature, than by changes in ice-stream input. This is different for the linearly behaving ice shelves. Here, ice-stream transience will lead to ice-thickness changes that might mask effects of variations in global boundary conditions.

The presence of an ice rise will primarily be seen in changes in the ice-velocity field of the ice shelf. While ice velocities will increase throughout the entire ice shelf as a result of enhanced ice-stream input, the velocity field in the vicinity of an ice rise will additionally be disturbed and will be different from that arising with no ice rise present.

Models of “real” ice shelves provide powerful tools in testing working hypotheses on the properties of an ice shelf. Our results demonstrate that the existence of a central zone of thin ice in Filchner–Ronne Ice Shelf is very unlikely. This conclusion is based on a comparison of modelled ice velocities with field data. Our predicted velocities obtained from models where the thicknesses of the central zone are similar to the surrounding ice yield a significantly better agreement with the measured velocities. Based on the predicted strain-rates of our model, in comparison with a map of surface features on Filchner–Ronne Ice Shelf (Swithinbank and others, 1988), we reach the same conclusion. This conclusion was recently confirmed by direct thickness measurement in the central zone of the ice shelf (Engelhardt and Determann, 1987), which resulted in significantly thicker ice than predicted previously (Crabtree and Doake, 1986).

Acknowledgements

This research was supported by the U.S. National Science Foundation (DPP 8509451) and the NATO Scientific Affairs Division (27–854/85). M.A.L. would like to thank the Deutsche Forschungsgemeinschaft, who supported attendance at the International Glaciological Society symposium in Hobart, where this paper was presented. We thank K. Hutter for a constructive review of an earlier version of this paper. This is Alfred-Wegener-lnstitut für Polar- und Meeresforschung contribution 113.