2.1. The equations used

The core problem is that of integrating the three-dimensional heat-conduction equations:

Here q (units K s⁻¹) is the heating within the ice due to friction, T (K) is the absolute temperature of the ice, p (kg m⁻³) is the density, c (J kg⁻¹ deg⁻¹) is the specific heat capacity, k (W m⁻¹ deg⁻¹) is the thermal conductivity, and ν is the three-dimensional velocity.

where u, v and w are the components of velocity and i, j and k are unit vectors along x,y and z.

The coordinates for this equation are the common rectangular (x,y, z, t) system with z directed downwards from an origin at sea-level. Since the upper and lower boundaries of the ice are not horizontal—inasmuch as the surface is assumed to have any arbitrary shape, as does the bedrock beneath, and melting is a function of the computations—there will be parts of these boundaries which cross the axes, leading to difficulties in the computational stage when derivatives are replaced by finite differences. Furthermore, the changes in horizontal and vertical disposition of the ice (as the ice mass changes shape in time) will mean that some particular points of the finite-difference grid will sometimes be in air, sometimes in water, sometimes in ice. Thus boundary errors may be quite large unless great care is taken, generally by decreasing the grid size (and therefore the time step), and using highly accurate finite differences near the edges of the ice.

On the other hand, this problem may be overcome quite simply by a change in the vertical coordinate. The equations which result are considerably more complex than those above, but since the number of points in the grid does not have to be increased, nor the time step decreased, nor special finite differences applied near the boundaries, the computation is considerably quicker and more accurate.

A new coordinate, ζ, is defined:

where M is the total melt beneath the ice, D is the depth of the ice, F is the bedrock height and E the ice surface elevation (all measured above the arbitrary, fixed reference level: sea-level). It is clear that all the ice is always contained between the levels ζ = 0 (the surface) and ζ = 1 (the lower boundary between ice and water or ice and bedrock), and that as the ice grows or shrinks, the vertical coordinate automatically expands or contracts.

This “mapping” of the (x, y, z, t) system into a non-orthogonal set of axes (x, y, ζ, t) is quite analogous to the σ-coordinate system developed by Reference PhillipsPhillips (1957) for numerical weather forecasting. Both new systems share the relative disadvantage that simple derivatives become transformed into much more complex forms. Techniques whereby these transformations may be developed will be omitted here, but may be found for example in Reference HaltinerHaltiner (1971). It may be shown that:

and

where , and the subscripts refer to the coordinate systems. No subscript indicates independence of either set of axes.

The increased complexity of the terms of Equation (1) in the ζ-system involves only a moderate increase in programming and computation time, which is far outweighed by the very decided advantage of automatic removal of some of the more difficult boundary problems.

The following equations, needed to complete the model, will not be discussed in detail; their interpretation can be found in Reference Budd, Budd, Jenssen and RadokBudd and others (1971) and Budd and Jenssen ([1975]).

Heat input into the ice mass both from below and from internal viscous dissipation is given, for basal heating only, by

where g is the acceleration due to gravity, and Γ is the geothermal heat flux. The second term represents the contribution made to the base gradient by frictional warming, in which ν _b and τ_b are the basal horizontal velocity and the base stress.

The heating by frictional dissipation inside the ice mass is

where ∂ν/∂ζ has a number of analytic forms or is determined directly from the observed and experimental data. In this latter case the lower boundary condition is best chosen to be not Equation (4a) but rather

Thus heating in the general case results from viscous dissipation and from the input to the base due to the geothermal heat flux.

When phase changes occur, Equations (4a) and (4c) are replaced by

as a lower boundary condition. The pressure melting temperature T_pm may be taken (see, for example, Reference DorseyDorsey, 1940) to be a simple function of depth: T _pm = 273 K — (0.007 K/m) D.

It may be shown quite readily (see Reference Budd, Budd, Jenssen and RadokBudd and others, 1971) that the amount of melt is simply proportional to the difference between the computed temperature gradient [∂T/∂ζ]_e (using Equation (5)) and the basal heat input. Thus

Here L is the latent heat of fusion, and [∂T/∂ζ]_b is given by either Equation (4c) for generalized flow, or is given by Equation (4a) for columnar motion. The quantity on the right-hand side of Equation (6) is a positive quantity for melting, and negative for freezing.

Continuity of mass for the local element, for the entire ice flowing through a column above the bedrock, and for the amount of melt water is governed by the equation

(note that here w = dz/dt and not dζ/dt) and

where m is the instantaneous amount of melt, m = ∂M/∂t and A is the surface rate of accretion (or ablation) of matter.

The total amount of melt at any location above the bedrock is simply:

where m is found from Equation (6).

Finally, as described by Budd and Jenssen ([1975]) the horizontal velocity at any point within the ice mass is found by integrating the shear strain-rate ∂ν/∂ζ using a table of empirical values of that quantity as a function of temperature and stress. With the stress defined by:

the table look-up is easily effected and leads to:

Here ν _b is the velocity at the lower ice boundary, and, for no slip, is zero.

With the horizontal velocity known, the vertical component of motion is found using the equation of continuity for an element of the ice (Equation (7)) and integrating vertically:

where is the vertical velocity at the base of the ice and w _a is the (x, y, z, t) vertical velocity at ζ = a. Once again it must be remembered that w = dz/dt and not dζ/dt.

2.2. The flow law used

As mentioned above, the horizontal velocity is determined by integration of the shear strain-rate, which is in turn a function of temperature and stress. The relation used is an empirical one, based on laboratory and field measurements, and presented in graphical form by Reference BuddBudd (1969) in his figure 2.2. For the model, the curves of that diagram were transformed into tabular form (Table I).

With stress computed by Equation (10) and temperature computed through Equation (1), the Table allows strain-rate to be found for any point. Integration along a vertical in the mass (Equation (11)) will give the horizontal velocity ν which is broken into its x and y components.

2.3. Temperature conduction of bedrock

Over long periods of time, the temperature profile in the underlying bedrock will vary with time in what may in all probability be a complex interaction with the temperature profile of the ice mass above. Since this pilot model was to be used to test its feasibility, and was to be used only for short periods of time (less than 5 000 years), heat conduction in the bedrock was ignored. There was, then, the further assumption that this effect was a second-order one. On the practical side, inclusion of points in the bedrock would increase either (or both) the time of the computation or the number of points in the grid. Since both time and space of the computer used were being stretched simply by the basic model, the refinement of bedrock conduction was not added.

3.1. The basic procedure

The use of the heat-conduction and continuity equations on a computer requires additional theory, amendment and information. The exact derivatives must be replaced by finite differences (so that a three-dimensional mesh must be chosen to cover the area of interest); the time step of the integration procedure must be consistent with the mesh intervals in order to avoid computational instability; and suitable boundary and initial conditions must be chosen. In addition, and with the object of saving computer time, the integration for this particular model does not proceed at the same pace for all grid points: some regions will have many integrative steps, with a small time interval, whilst others will have few steps, but with a large time interval.

The skeleton of the computational procedure is as follows:

• i. Assuming that temperatures everywhere within the grid are known (see below), the velocity distribution (both horizontal and vertical) is computed.

• ii. With a value for the time step involved in the integration correctly chosen so as to satisfy the stability criterion at the given location (see below), the temperature changes within the ice are found and hence the new temperature distribution is determined. The heat-conduction equation used for this step is a modified and simplified form of Equation (1) in the ζ-coordinate system.

• iii. The continuity equations for depth and melt are then satisfied, and any base temperatures in need of correction are adjusted.

• iv. The ancillary computation of trajectory locations is made.

• v. Steps (i) to (iv) are repeated until the integration time is some arbitrary amount.

• vi. The ancillary computation of dielectric absorption, both per unit depth and throughout the entire ice depth, is made.

• vii. Relevant parameters are then printed out. These include the two-dimensional fields: surface and basal temperature distribution; minimum temperatures at all verticals within the ice, and the depth at which these low values occur; the surface and basal computed temperature gradients; distribution of melt water; the mean temperature of the ice in any vertical; the effective “mean absorption temperature”, that is the temperature a column would have to have, if, when isothermal, the dielectric absorption were to be the same as that computed at that location; the surface horizontal streamlines ; and the vertical velocity distribution. See Reference Budd, Budd, Jenssen and RadokBudd and others (1971) for a further discussion of these quantities.

  Another important output statistic is the mass budget for the entire ice mass—and this includes accumulation and ablation areas, loss due to melting, and loss of ice at the horizontal edges of the ice due to “calving”. Local mass budgets for regions of particular interest are also made. Floating ice, attached to the main mass, is followed until calving occurs.

• viii. The process recycles back to step (i) until the full computation time has been reached.

The computational problems involved in this procedure are discussed in detail in a separate report (to be published by D. Jenssen). The main problems concern the stability of the finite-difference algorithms and the parameterization of the spread of the ice sheet over the oceans.

3.2. Computational stability

For the simple one-dimensional case, it can be shown that in order to preserve stability and to damp errors, preventing exponential growth, the time and space increments must satisfy the following relation:

The non-linearity, and three-dimensional nature, of Equation (1) make the determination of its stability criterion very much more difficult. However, if it is realized that the time step depends on the smallest grid spacing, which is the vertical, then Inequality (13) must be a good approximation to the true stability criterion. (The horizontal spacing is measured in tens of kilometres, the vertical in tens of metres.) The non-linearity of the heat-conduction equation will probably reduce the maximum time step slightly, if at all. Pragmatically, Inequality (13) has been used in the preliminary three-dimensional model calculations to determine the maximum time step which can safely be used, and no evidence of computational instability has been found. With Inequality (13) violated, however, forward time integration did in fact very quickly lead to the meaningless results typical of computational instability.

Now examine the stability criterion (13) in more detail. Since the number of points in the vertical is the same for all points of the grid, Δζ is constant, and thus the time step, Δt, is roughly proportional to the square of the ice thickness. It would appear that to avoid computational instability, the choice of the time step is governed by the thinnest ice which will be encountered in the grid which has been arbitrarily set at 5o m. If, say 50% of the ice is of reasonable thickness, then a large time step (100 years or so) may be used. For the thinnest ice, which probably will not occupy more than a maximum of 10% of the grid, the time step must be less than 1 year—and may even drop to 0.1 years. Thus if the same value of Δt is used over the entire mesh, gross inefficiency in computation time will result: the heat conduction equation for the points of thick ice will be integrated perhaps a thousand times when in fact one integration with a value of Δt a thousand times larger would have produced the same temperatures.

To avoid this inefficiency, different values of the time step were used in regions of different ice thickness, so that if the largest time step permissible anywhere in the grid is Δt _max, then at any arbitrary point the time step Δt is such that At. Δt _max = NΔt, where N is an integer. This means that while the integration proceeds at different paces in different spots, every Δt _max all these independent integrations “mesh” together. Of course, none of the Δt used violate the stability criterion, so that not only is computational stability preserved, but also the computation of temperatures proceeds everywhere at an optimum speed.

3.3. Changes in horizontal extent of the grid

During the course of the computation, as the ice grows or shrinks horizontally, or as it flows over the ocean, the type of material at a grid point may change, so that the computations performed there will change in frequency and nature, and “shocks” at the boundary may be set up. The coarser the grid, the more severe will be these shocks, unless special “tracking” procedures (equivalent to having a very fine grid along the boundaries of the more coarse one) are utilized.

For any complete study, the ice flow past the grid boundaries must be treated from time step to time step, much in the same manner as the glacier snout was followed computationally in the model of Budd and Jenssen ([1975]). For the present preliminary three-dimensional computations, however, no such treatment has been made. The discussion of the results obtained (see below) will show the effect of neglecting this.

3.4. Ice shelves and calving icebergs

As mentioned earlier, the minimum ice depth considered has been about 5o m. Any thinner ice results in so small a time step that computation time becomes prohibitively large—especially when it is noted that total integration times of the order of hundreds of thousands of years may be required.

This restriction can of course be varied according to the size of the ice sheet, and with the times over which an integration is to be performed. It has the advantage of allowing calving to be treated in a very simple manner. Whenever the depth becomes smaller than the cut-off value chosen, all ice at that point is removed from the grid (provided that the point is in the ocean and that at least one of its eight closest neighbours is also in the ocean, and has a depth less than the minimum figure). The ice removed in this manner is thought of as a detached iceberg, and is no longer treated.

If the ice depth is less than 5o m but the point is over the bedrock, or if all its neighbouring points have ice thicker than this minimum figure, the point in question merely represents a local deep depression, and cannot possibly calve away from the main ice mass. Such special points are tagged, and either arbitrary temperature distributions assigned to them (of a steady-state nature) or they are given temperatures which are weighted averages of the surrounding ice temperatures. Thus velocities may be computed for these anomalous points, and continuity requirements may be satisfied as well. In this manner ice may grow over the bedrock, or extend, without calving, over the oceans to form ice shelves.

For the ice over the land, the temperature and velocity are not computed by the main equations of the model until the thickness exceeds 50 m. Over the sea, ice velocity and temperature are always weighted averages of the surroundings. This simplistic approach to ice-shelf dynamics and thermodynamics is one feature of the model which must be treated more realistically in any larger-scale study.

4.1. Grid, time steps, initial conditions

The model outlined above has been tested by means of a coarse representation on the Greenland ice sheet. For this pilot study an IBM 7044 machine was used. Its restricted capacity imposed the choice of a 12 X 12 X 10 (vertical) point grid; the east-west grid spacing was 100 km, the north-south spacing 200 km, while the vertical spacing varied from 5 to 300 m.

The boundary of this grid is shown in Figure 1. Also shown in this figure is the land extent of the grid: any points between the interior and exterior boundaries were considered to be oceanic. This is obviously not true for those points midway on the east coast of Greenland, but the main purpose of this preliminary study was a thorough testing of the model and of the program based on it. The uneven land boundary on eastern Greenland was deliberately chosen to examine the effects of various backward finite-difference operators on the calculations. Note the almost isolated point seven rows from the northern boundary: a similarly isolated point exists at the southernmost land point of the Greenland grid.

Fig. 1. Boundaries of the grid, and maximum permissible time steps. The heavy outer line gives the horizontal extent of the grid (12 X 12 points), while the thin straight line gives the boundary of the (assumed) land area within the model. Circles indicate mesh points with an assumed initial ice thickness less than 1 m. Dashed curved lines give isochrones (in years) of the maximum permitted time steps which can be used in the forward integration of the model. Larger time steps will lead to computational instability.

Fig. 2. The model accumulation distribution. Contours of accumulation (Mg m⁻³) from Reference MockMock (1967). Accumulation over oceanic areas is assumed to be zero.

Fig. 3. Model surface elevation and bedrock. Budd and others (unpublished) after Reference BaderBader (1961) and Holtzscherer and Bauer ([[956]). Continuous lines are surface-elevation contours, broken lines are for bedrock: underscored contours are for surface elevation. Units are metres (above sea-level).

Fig. 4. Initial and + 1 000 year ice thicknesses. The dashed line gives the initial assumed ice thickness in metres (consistent with the data of Figure 3), and the full line gives contours of thickness obtaining after a forward time integration of 1 000 years. Underscored contours are for + 1 000 years.

Fig. 5. Model surface temperatures. Contours of the initial surface ice temperature (in °C), from Budd and others (unpublished). This pattern will change during the integration procedure as the ice rises or falls.

Fig. 6. Steady-state base ice temperatures (from the three-dimensional model). The temperatures (ºC) obtaining at the ice/bedrock boundary computed from the model when there is no change at any grid point of temperature with time. See text for further details. Hatched area indicates region of melting.

Another feature of Figure 1 is that those grid points which initially had no ice (thickness less than 50 m) are shown as open circles. Figure 1 also shows isochrones of the initial time step (which changed with the ice thickness during the integration). In order to speed up the computation, points with less than 50 m of ice and with local time-step values less than 1% of the maximum permissible for the grid were assigned steady-state temperatures or the average of the surrounding temperatures, and treated only in the continuity equation, until both the ice thickness exceeded 50 m and the local time step had become larger than 1% of the permissible grid maximum.

The surface accumulation, based on Reference MockMock (1967) is shown in Figure 2. Oceanic values were assumed to be zero in the pilot calculation; this produced strong accumulation gradients along the western and south-eastern coasts and had important effects on the calculations.

Figure 3 shows the bedrock (broken lines) and surface elevation contours (solid lines) used in this study. These are based on data from Reference BaderBader (1961) and Holtzscherer and Bauer ([1956]). The corresponding ice thickness distribution is given in Figure 4 as the dashed lines: the solid lines are computational results, and are discussed below.

Surface temperature is a boundary condition; values of this parameter are given in Figure 5, based on data in Reference Mock and WeeksMock and Weeks (1965). However the full temperature field must initially be prescribed, and this distribution should be consistent with the model equations, in order to avoid spurious rates of temperature change which may amplify during the time integration (especially in view of the velocity/temperature feedback) to produce eventually small but significant variations in the temperature and motion fields.

To prevent this from happening, the initial temperatures were computed with the one-dimensional steady-state model of Reference Budd, Budd, Jenssen and RadokBudd and others (1971), and used as a first guess to the true temperatures by inputting them to the three-dimensional model. The surface temperature was forced to be constant in time, and, as well, no change in the shape, extent, or thickness of the ice field was allowed. Mass outflow from the edge of the ice was permitted, however. With these restrictions, the temperature distribution adjusted so that at any point decreased in magnitude and approached zero, that is reached a steady state. Since the velocity distribution was still being calculated (and used as a source of frictional heating), it remained consistent with the temperature pattern.

This computational approach to the steady state is rather slow: for the present case, some 120 000 years of simulated time were required to have a rate of change in temperature everywhere of less than 0.01 deg/100 years.

It is not feasible to show graphically the entire three-dimensional temperature field, but a sensitive parameter which has been extensively used in previous studies is the basal temperature. This is shown in Figure 6, and, for comparison, the initial basal temperatures of the one-dimensional steady-state model are given in Figure 7. Considerable differences in the extent of basal melting are indicated between the (assumed) initial and (computed) steady states; these are due almost entirely to the coarse nature of the vertical grid which could not cope adequately with the rapidly changing temperature gradients, and hence heat fluxes, near the base of the ice. Even so the actual temperature differences are only of the order of a few degrees and readily reducible (though not completely removable), by increasing the vertical resolution.

The velocity distribution corresponding to the computed temperatures is shown in Figure 8; broken lines give the initial value and full lines those developing subsequently. For comparison Figure 9 gives the velocities deduced from assuming a balanced mass flow over Greenland. Both patterns show a region of very high velocity near the west coast with quite low velocities over the central and part of the eastern portion of Greenland. The computer model is unable to match velocities along the east coast because the ice there is so thin that it is excluded from the velocity computations.

Fig. 7. Steady-state base ice temperatures (from a one-dimensional model). Temperatures (°C) obtaining at the ice/bedrock boundary computed from a simple model which ignores all horizontal temperature effects, and which has a fine vertical resolution. These temperatures are used as a first guess to those of Figure 6. Hatched area indicating melt region. See text for further details.

Fig. 8. Initial (steady-state) and + 1 000 year surface velocities. Isotachs (m/year) of surface velocity initially are computed by running the model in its steady-state mode. These isotachs are the dashed lines, and they are consistent with the temperature distribution of Figure 6. Velocities after 1 000 years are shown as full lines, and these are consistent with the data of Figure 11. Underscored values are for + 1 000 year contours.

Fig. 9. Balance velocities. Isotachs (m/year) computed by assuming a mass balance over the Greenland continent with mass transport being determined by the accumulation distribution, and ice thickness remaining constant at its present values.

Fig. 10. One-dimensional transient-state-model surface velocities. Using the temperatures computed from a one-dimensional model, these surface velocities are found in the manner described in the text. Units are m/year.

Apart from this difference there remain enough discrepancies between Figures 8 and 9 to suggest that the ice velocity distribution in the Greenland ice sheet is not consistent with the assumption of balanced flow. This is further underlined by the results of a transient one-dimensional temperature and velocity calculation along the lines described by Reference Budd, Budd, Jenssen and RadokBudd and others (1971). The resulting velocities are reproduced in Figure 10 and (disregarding magnitudes which involve a simple adjustable parameter) the general pattern is very close to that of the computed steady-state velocities in Figure 8.

4.2. Results of a 1 000 year calculation

With the grid initial inputs and time steps of Section 4.1, the integration of the basic equations (Section 2) allows temperatures, velocities, ice shape, and ice extent to be computed at any future time.

Fig. 11. Base temperatures after 1 000 years. The basal temperatures (°C) obtained by a forward time integration of 1 000 years. These are consistent with the ice thickness data of Figure 4, and the isotach pattern of Figure 8. Hatched area indicates basal temperatures greater than the pressure melting point.

The full lines in Figure 4 show the ice thicknesses which have developed after 1 000 years of unrestricted flow. The 3 000 m contour has enlarged, and the maximum thickness of the ice has shifted slightly towards the west, where all ice thicknesses have increased markedly. Growth of the ice to the east has not been nearly so pronounced, due to the thinner and colder ice there, which will lead to much smaller velocities—see also Figure 8. A very rapid build-up of ice has occurred to the west at about lat. 68° where a small pocket of depths greater than 3 000 m can be observed. This points out a deficiency in the present model associated with the very coarse grid used (100 km west to east and 200 km north to south): in the integration of the continuity equation, the large distances over which the finite differences are taken tend to produce a smoothing effect, and to reduce the mass flow at any particular point. Near the edges of Greenland, in spite of allowing mass transfer to the oceans, this will produce the bank-up of material of Figure 4.

The base temperatures (Fig. 11) show an interesting phenomenon over the west of Greenland; due to the increased ice thicknesses and higher surface elevations, the consequent decrease in surface temperatures has produced a marked decrease in base temperature, and contracted basal melting to a much narrower region. In fact on the west coast the temperatures have fallen to below —5°C. The velocities at 1 000 years (full lines in Figure 8) show a general increase over most of Greenland and a somewhat more complex pattern with the formation of a region of high speeds to the south-west. The areas of greatest velocity near the coast also correspond to the growth of oceanic ice shelves.

5. Discussion

Since Figure 11 shows that the base of the ice has cooled near the western extremities, the velocity increases there must be due to the greater ice thickness with the consequent increased surface slopes.

On the other hand the rapid decrease of velocity between these maxima and the edge of the grid arises directly from the colder ice temperatures. This slowing down of the ice flow will also contribute to the banking-up of the ice near the coast.

We have therefore the following sequence of events for this particular case: due largely, it is thought, to the coarse grid, but also to the accumulation distribution and possibly the flow law used, the ice grew rapidly near the west coast; the surface slope increased and so therefore did the velocity; at the coast, however, the increased surface elevation led to colder temperatures and lower velocities—the ice flow was impeded; the ice depth increased and the surface slope increased; and so forth. Thus the velocities near the coast became larger rapidly, whilst those at the coast did not increase so quickly, and mass accumulated at the western edge of Greenland.

Restated, in computational terms for the inadequate form of the model here considered:

• i. The coarse vertical spacing leads to incorrect lower temperatures, with the implication of large heat inputs: this, in turn, will lead to a spuriously high velocity distribution.

• ii. The coarse horizontal spacing leads to smoothing out of features inland with a consequent lowering of velocity there.

• iii. The coarse horizontal spacing near the coast exaggerates the surface slope, and hence the velocity, creating a large mass flow. This is prevented from leaving the grid because

• iv. the coarse horizontal spacing at the coast, coupled with the inferior method of allowing for floating ice shelves cannot transport mass away from the central land mass.

• v. The coarse horizontal spacing also contributes to the mass barrier on the west because it implies a large horizontal gradient of accumulation near the coast.

These deficiencies of the model interact with each other in an unstable feed-back mechanism so that mass builds up excessively on the western coast ever more quickly. Shortly after 1 000 years, in fact, the surface slopes computed within the model became so large that the concomitant base stresses exceeded the maximum allowed for in the table, and the computation was terminated.

On the credit side it may be claimed that the equations of the model are realistic since they lead to three-dimensional distributions of temperature and velocity consistent with those obtained from two other, independent, one-dimensional models (the steady and transient state models of Reference Budd, Budd, Jenssen and RadokBudd and others (1971)). Differences between these three models are merely those to be expected when a basic model is successively made more realistic and hence more complex.

For the first time, velocity and temperature have been made mutually interactive for a three-dimensional ice mass, and their distributions at any time are consistent in regard to temperature, motion, continuity, ice thickness, and mass flow. Moreover the model permits the generation of floating ice to be followed, as well as mass loss by calving.