Introduction

The character of ice-shelf modelling is distinctive in three important ways. First, problems studied with ice-shell models tend to emphasize flow dynamics and comparison with observed flow fields as model-performance measures. This emphasis is motivated by the availability of field observations and by the fact that a well-accepted statement of the governing equations, boundary conditions and constitutive equations is yet to be established. Secondly, ice-shelf models must contend with the ocean and its effect on the basal melting rate. This effect has a direct role in the time-dependent mass balance part of ice-shelf modelling. It also has an indirect influence because the temperature field of ice shelves is strongly influenced by basal melting. In many ways, the difficulty of measuring, or predicting from first principles, a basal melting rate is comparable to that of measuring or predicting a basal sliding velocity in a grounded ice sheet. Thirdly, ice-shelf models must contend with two technical obstacles. The stress-balance equation which governs ice-shelf flow is non-local, i.e. is elliptic in the two horizontal dimensions, and the mass-balance equation is not diffusive, i.e. is hyperbolic, and is thus capable of propagating shocks. These obstacles put stringent demands on the numerical methods applied in ice-shelf models and greatly complicate efforts to couple ice-shelf models with models of inland ice sheets (e.g. Huybrechts, 1990).

The participants of the 1994 EISMINT model-intercomparison workshop field in Bremerhaven decided that ice-shelf model validation and intercomparison could be simplified if a standard test could be developed which separately considered the solution of the elliptic stress balance equations apart from all other considerations, such as the mass balance. Here we present a standard ice-shelf model test which accomplishes this goal. For realism, and to assess ice-shelf model performance under conditions similar to actual research applications, we have based the test on the Ross Ice Shelf, Antarctica. Inter-comparison of five independently developed ice-shelf models presented here illustrates the effects of spatial resolution (e.g. variable resolution in finite-element methods vs uniform resolution in finite-difference methods) and of the form of model equations and boundary conditions (e.g. treatment of the ice-front boundary condition).

Notation

Typical Model Equations

To illustrate the testing package, we compared five ice-shelf models (referred to by their city of origin and by the initials of their creator: Bremerhavenl (PH), Bremerhaven2 (JD), Chicago1 (DRM), Chicago2 (CLH) and Grenoble (VR and CR)) using the Ross Ice Shelf modelling data base. All models except one (Bremerhaven2) are based on the following simplified, or reduced, stress-equilibrium equations.

(1)

and

(2)

where variables are defined in the notation list. The equations on which Bremerhaven2 is based differ from the above equations in that the ice thickness h does not appear on either side of the equations and horizontal viscosity gradients (∂v/∂x and ∂v/∂y) are not accounted for. The principal simplification embodied by all models is that horizontal flow is depth-independent (e.g. Sanderson and Doake, 1979). The validity of this assumption is contingent on the condition that [H]/[L] ≪ 1 where [H] and [L] are characteristic vertical and horizontal distance scales respectively (e.g. Morland, 1987; Morland and Zainuddin, 1987).

Boundary conditions needed to solve Equations (1) and (2) consist of kinematic (velocity specified) conditions at ice-shelf/inland ice boundaries and dynamic (pressure of sea water specified) conditions at the seaward ice front. After Morland (1987; see also Morland and Zainuddin (1987)), the boundary condition at the seaward ice front is idealized by applying a depth-integrated sea-water pressure force along a contour in the horizontal plane that follows the real ice front. In essence, the ice-front boundary condition disregards three-dimensional effects expected within the neighborhood of the ice front.

(3)

where n is the outward-pointing unit normal to the ice-front contour in the horizontal plane.

Equations (1) and (2) use an effective viscosity v to embody Glen’s flow law. The definition of v involves a temperature dependent rate constant and a flow-law exponent n:

(4)

The flow-law rate constant represents a depth-averaged value of the flow-law rate constant which can vary considerably over the Ross Ice Shelf as a result of temperature and density variation (e.g. Thomas and MacAyeal, 1982; MacAyeal and Thomas, 1986).

Ross Ice Shelf Database

Geophysical and glaciological observations of the Ross Ice Shelf have been made since the early part of this century (e.g. Bentley, 1984), We compiled a data base of these observations primarily from measurements made during the Ross Ice Shelf Geophysical and Glaciological Survey, 1973–78 (RIGGS) and the Siple Coast Project 1983–90 (see J. Glaciol., 39(133)). These observations include ice thickness (Bentley and others, 1979; MacAyeal and others, 1987), sea-bed depth (Bentley and Jezek, 1982), horizontal surface velocity U (Thomas and MacAyeal, 1982; Thomas and others, 1984; MacAyeal and others, 1987), annual average surface (10 m) temperature (Thomas and others, 1984), accumulation rate (Thomas and others, 1984), and density/depth and temperature/depth profiles at selected locations (Kirchner and others, 1979; Thomas and MacAyeal, 1982). Ice-stream and glacier input velocity boundary conditions were compiled from MacAyeal and others (1987) and MacAyeal and Thomas (1986). Ice-shelf modellers may access this. Ross Ice Shelf data base at the World Data Center for Glaciology A, Boulder, Colorado, U.S.A.

Numerical Domain

Two standard numerical discretizations of a 465 000 km² area of the Ross Ice Shelf are provided in the test data base (Fig. 1). A regularly spaced 147 × 111 array of grid points is designed for finite-difference models. An irregular mesh of triangular elements is designed for finite-element models. The spatial resolution of the finite-difference grid is 6.822 km. The spatial resolution of the finite-element mesh varies from 0.5 km near boundaries to 20 km in central areas (Fig. 1). The finite-difference grid and the finite-element mesh were constructed to possess approximately the same number of grid points or nodes. (For the finite-difference grid, there are 16 317 grid points, of which 11 067 are active. For the finite-element mesh, there are 16 231 nodes and 30 639 elements.)

Fig. 1. Map of Ross Ice Shelf (upper panel) and finite-difference and finite-element discretizations of the Ross Ice Shelf near the Byrd Glacier outlet (lower panels). Coordinates are grid latitude and longitude as described by Bentley and others (1979). Region A (shaded) is the Ross Ice Shelf. Region B is the artificial 1 m thick ice-shelf region added to the finite-difference domain to create a rectilinear ice-front continu (note that it cuts off a part of the ice shelf near Ross Island). Label C denotes the ice-front contour in the finite-element domain. The part of the open ice front west of Ross Island (the McMurdo Ice Shelf) is treated as a closed boundary in the finite-difference model. Labels D denote locations of outlet glaciers and ice streams where kinematic boundary conditions are non-zero. Grid resolution is 6.822 km for the finite-difference domain. Mesh resolution for the finite-element domain varies from 0.5 to 20 km.

Two options for treating the ice-front boundary are provided to accommodate finite-difference models which require simple, rectilinear geometry at boundaries where dynamic conditions are imposed. In both options, the ice shelf is artificially extended beyond the mid 1970s ice-front position (sec Fig. 1). In the Grenoble model, the ice thickness in this artificial region is reduced to 1 m, while in the two Bremerhaven models the ice thickness is extrapolated from the existing ice shelf.

Test results presented below show that artificial extension of the ice shelf for finite-difference models generates little degradation of model performance within the confines of the real ice shelf. However, along the actual ice-front contour (separating the real and artificial parts of the ice shelf), finite-difference models cannot strictly enforce zero tangential stress. Thus, in circumstances where ice-front phenomena are of interest, finite-element models may be preferable because of their capacity to handle curved geometry.

Model Performance Assessment

Test results are assessed through comparison of the derived surface velocity and observations. At the simplest level, a single scalar measure of velocity misfit, X ², is computed according to the formula

(5)

where u is the model velocity, subscript i identities the RIGGS or Siple Coast Project station where a surface velocity is observed. The mask variable δ _i ; is 0 or 1 depending on whether station i is within the numerical domain and whether there are valid velocity observations available for this station. The weighting variable σ _|U| is the assumed uncertainty of the observed velocity, which varies from ±5 to ±30 m a⁻¹ among the RIGGS stations (Thomas and others, 1984). For the model intercomparison presented below, we adopt a uniform value of σ _|U| = 30 m a⁻¹. The sum defining X ² in Equation (5) is not divided by N = 210, the number of stations where model and data can be compared, because we wish to identify X ² with the χ ² variable used in statistical analysis of model/data mismatch (e.g. Menke, 1989, p. 32).

Other measures of model performance can be devised by the test-package user to emphasize model performance around specific physical features, such as ice rises, shear margins and the ice front, or to emphasize differences that occur in model intercomparison exercises. These measures are likely to take the form of contour, vector and tensor diagrams that can be constructed as diagnostics of the stress balance equations governing ice-shelf flow.

Table 1. Model performance index, X² (non-dimensional)

Intercomparison of Five Ice-Shelf Models

We intercompare five ice-shelf models introduced above. Three models, Grenoble, Bremerhaven1 and Bremerhaven2, are based on the finite-difference method. Two models, Chicago1 and Chicago2, are based on the finite-element method. Grenoble is currently under development for future use in large-scale ice-sheet modelling (Rommelaere and Ritz, 1996). Bremerhaven1 is associated with large-scale coupled dynamic/thermodynamic models of Antarctica (Huybrechts, 1990). Bremerhaven2 has been applied to the Ronne–Filchner and Ekström Ice Shelves (Determann, 1991). Chicago1 is under development for studying ice-shelf/ice-stream interaction in Hudson Strait and the Labrador Sea off the North American continent during the glacial period. Chicago2 has been used to simulate the effects of recent iceberg calving off the Larsen Ice Shelf, Antarctica. None of these models has been specifically set up to study Ross Ice Shelf flow.

Fig. 2. Comparison between model-derived velocity magnitude (using Chicago1 and a spatially uni form estimate of the depth-averaged flow-law rate constant) and observed velocity magnitude for 156 of the 210 RIGGS and Siple Coast Project stations contained within the finite-element model domain (upper left). Misfit would be improved if the depth averaged flow-law rate constant were spatially variable. Comparison between model-derived velocity magnitude for Grenoble and Chicago1 (lower left), Grenoble and Bremerhaven1 (upper right) and Bremerhaven2 and Bremerhaven1 (lower right) at the 156 station.

In this intercomparison we use a spatially uniform value of (1.9 × 10⁸ ). We adopt this simplification, because spatially varying values estimated by Thomas and MacAyeal (1982) have been found to yield a poor fit between model and observation. This inadequacy is due to the assumptions Thomas and MacAyeal made about the basal-melting rate and to their simplified, one-dimensional heat-flow model used to determine internal ice-shell temperatures. Despite the flow-law rate constant used in the test, the results of all five models are sufficiently close to the observed velocity that an intercomparison of the models appears justified.

The X ² values associated with the five model runs are listed in Table 1. Figure 2 presents comparison plots of | u _i | vs | U _i |, i = 1, …, 210, for the output of Chicago1. As explained above, agreement between observed and modelled velocity was not an intended outcome of the intercomparison test. Pairwise intercomparison of the |u _i |, i = 1, …, 210, obtained by three of the models (Chicago2 is sufficiently close to Chicago1 that its results need not be specially highlighted) are also shown in Figure 2. Models agree exactly when their respective values of |u _i |, i = 1, …, 210, plot on a 45° line in Figure 2. According to Figure 2, results of Grenoble and Chicago1 agree for points at the high end of the velocity range; but, for points in the low-velocity range, there is considerable scatter (on the order of ±100 m a⁻¹). Results of Grenoble and Bremerhaven1 agree well in the low velocity range and are scattered off the 45° line in the high-velocity range. Bremerhaven1 and Bremerhaven2 show considerable scatter off the 45° line in all velocity ranges.

Figure 3 displays contour maps of |u| associated with Chicago1, Grenoble, Bremerhaven1 and Bremerhaven2. All models produced a velocity pattern that is similar to that observed (Thomas and MacAyeal, 1982). The flow maximum is located at the ice front and strong shear layers are produced at the sides of ice rises and along coasts. In all models, the maximum velocity achieved at the ice from exceeds the observed flow by about 25%. We expect this misfit to be corrected once better estimates of the flow-law rate constant are available.

Differences between the models (figs 2 and 3) are apparent, First, the effect of model-equation differences between Bremerhaven2 and the other models (eliminating h from both sides of Equations (1) and (2), and disregarding ∂v/∂x and ∂v/∂y terms) is nicely demonstrated by comparing panels G and H of Figure 3 with the other panels. The main effect of the model-equation differences is to uncouple velocity gradients from thickness gradients. This results in a weaker velocity gradient close to the grounding line and a stronger velocity gradient in the relatively flat area downstream of Crary Ice Rise in the Bremerhaven2 results compared to the other model results.

A second, and perhaps major, difference between models arises because of differences in boundary conditions specified at the ice front. In models Chicago1 and Chicago2, the finite-element method is ideally suited to treat curved ice-front boundary contours. The normal and tangential components of the boundary condition presented in Equation (3) are thus readily satisfied along the actual ice-front boundary by both models; and this, perhaps, is why the results of the two Chicago models are virtually identical.

Fig. 3. Model velocity magnitude for Chicago1 (panels A and B), Grenoble (panels C and D). Bremerhaven1 (panels E and F) and Bremerhaven2 (panels G and H). Contour interval for panels A, C, E and G is 100 m a⁻¹ and for panels B, D, F and H is 25 m a⁻¹ up to 200 m a⁻¹ and is 200 m a⁻¹ thereafter.

In the finite-difference models, Grenoble and Bremerhaven1, the exact pressure balance is specified in the direction normal to the ice front (Equation (3)). In both models, the tangential velocity component is determined by the dynamic condition ∂v/∂x = −∂u/∂y, where v is the tangential velocity, u is the normal velocity, x is the coordinate directed perpendicular to the ice front and y is the coordinate parallel to the ice front.

The difference between finite-difference models Grenoble, Bremerhaven1 and Bremerhaven2 is the adoption of a 1 m thick artificial ice-shelf extension in the Grenoble model, and an ice-shelf extension with extrapolated thickness in the Bremerhaven1 and Bremerhaven2 models. In the case of the Grenoble model, small thickness in the artificial ice-shelf extension makes the solution less sensitive to the exact prescription of the boundary condition. This explains why the results of the Grenoble model are almost the same as those of the finite-element models (Chicago1 and Chicago2), which are presumed to satisfy the ice-front boundary conditions most readily along an irregular geometry.

The Bremerhaven2 model has a boundary condition for the normal velocity component associated with a freely expanding ice shelf in the normal direction only. Here, the transverse-velocity component has a zero gradient in the normal direction. These boundary conditions are the primary cause of the high ice-shelf velocities in the Bremerhaven2 model compared to all the other models.

One surprise resulting from the model intercomparison is that difference in spatial resolution between finite-element and finite-difference domains does not appear to be significant in determining large scale flow patterns. This suggests that the numerical efficiency associated with finite-difference methods over finite-element methods would be a leading advantage in designing ice-shelf models to treat large-scale flow. However, subtle differences between finite-element and finite-difference models can be seen in several small-scale features where spatial resolution has an influence on the fidelity of the numerical domain to the natural geometry. In panels B and D of Figure 3, comparing results from the Chicago1 and Grenoble models, for example, ice velocity dowstream of the Mulock Glacier outlet in the Chicago1 simulation is dramatically less than that in the Grenoble simulation. This difference is due to the fact that the finite-element mesh reproduces a narrower and extended glacier outlet channel than that reproduced by the finite difference grid. Apparently, the outflow of Mulock Glacier meets greater resistance in the long, narrow finite-element representation than in the more open finite-difference representation of the same ice-shelf geometry. Regardless of which model results are more accurate, the difference displayed by the Mulock Glacier outlet comparison suggests that models based on finite-element and finite-difference methods may have important differences in applications to regional-scale problems of complex geometry.

Conclusion

The test package presented here offers several benefits to potential ice-shelf model developers. First, it offers a chance to catch coding mistakes. Secondly, developers can periodically recheck their code against the observed flow of the Ross Ice Shelf after making improvements in model physics. Thirdly, differences between valions numerical algorithms and their suitability for specific purposes can be assessed using a standard, controlled test problem. In our experience, all three benefits were realized in performing the intercomparison of five ice-shelf models. We learned that finite-difference and finite-element methods perform equally well in large-scale applications typical in ice and climate research. The finite-element method, however, appears to have an advantage in applications where complex, regional scale phenomena are of interest and where the ice-front boundary condition must be applied rigorously. We, additionally, learned that treatment of ice-front boundary conditions is the major factor in determining differences between ice-shelf model results in our standardized test.