On the numerical computation of temperature in an ice sheet

Abstract

Quadrature solutions for temperature fields in the central regions of ice sheets where flow is by internal deformation are presented and compared with numerical computations and the European Ice Sheet Modelling Initiative (EISMINT) Benchmark standards. The solution appears to be more accurate than the EISMINT Benchmark standards for thermomechanically uncoupled temperature distributions. The ability of finite-difference and pseudo-spectral methods to reproduce this solution is considered.

Errors in the EISMINT Benchmark are larger than expected. The possibility that they could arise from inaccurate evaluation of the vertical velocity is considered, and found to be unlikely. Formulae for computing the vertical velocity in Jenssen's σ coordinate are complicated and may lead to programming errors. A simple form of the heat-transport equation in the σ coordinate for flow by sliding and internal deformation is derived arid presented, and it is shown how this form is particularly suitable for finite-element computations. It is also shown that only two quadratures are necessary to compute the vertical and horizontal heat transport and the ice flux.

1. Introduction

The recent European Ice Sheet Modelling Initiative (EISMINT) Benchmark series of ice-sheet modelling intercomparisons has provided an invaluable set of tests for comparing and validating ice-sheet models in ranges of parameter space where no analytical solutions exist (Huybrechts and others, 1996). While the analytical solution for temperature has long been known (Robin, 1955) in regions where horizontal flow can be neglected in the heat equation and where vertical velocity is proportional to elevation as in a plug flow, it has not been generally appreciated that accurate solutions only requiring a quadrature exist for flow by internal deformation, for both linear and non-linear rheologies, although much of the relevant theory already exists (Zotikov, 1986). As with the Robin solution, one needs to solve a second-order boundary value problem, and a first integral is readily obtained analytically. The solution can then be obtained by a quadrature, which plays the same role as the error function in the Robin solution.

If the ice-sheet thickness and upper surface vertical velocity are known, then this quadrature solution can be applied, which has the consequence that the EISMINT Benchmark Level 1 temperature calculations, which are not coupled to the ice-sheet rheology, can be compared with this solution in areas where horizontal flows are negligible (e.g. the central areas of ice sheets). Results from these areas have been reported (Huybrechts and others, 1996) and are compared here with the analytical solution. This brings us to the main point of this paper: the mean EISMINT Benchmark results are not particularly accurate, and if one were to seek to validate one's model by obtaining the mean EISMINT Benchmark value, this would be possible only by using rather inaccurate finite-difference (FD) approximations. The solution presented in this paper is more accurate and it is also shown that second-order schemes with ten FD grid intervals can reproduce this solution to the expected 1% accuracy, as compared with the 5% accuracy of the mean value of the EISMINT Benchmark experiment.

Two other sources of error are investigated: the numerical integration required to compute the vertical velocity, and programming error. Firstly, it seems that the error in the EISMINT Benchmark Level 1 is too large to be explicable by integration error associated with the computation of the vertical velocity. Secondly, when directly transformed into the terrain-following coordinate (Jenssen's (1977) σ coordinate) usually used in glaciology, the heat-transport equation becomes extremely complicated in the sense that it contains a large number of terms. However, many of these terms can be shown to cancel, and a much simpler form of the heat equation thereby derived. This form, a development of work by MacAyeal (1997) and Hulbe (1998), lends itself naturally to consistent definition of flux and vertical heat transport, and also means that the computation of the heat transport and ice continuity requires only two quadratures.

2. Theory

The steady-state diffusion advection equation for an ice sheet near its centre, where horizontal advection and dissipation can be neglected, is given by

(1)

where β= K/ Ha,H is the thickness of the ice, a is the accumulation rate, K is the thermal diifusivity of ice, θis the temperature, ζ — z/H, z is the elevation above the bed, and w is the vertical velocity w normalised by the accumulation rate such thatw w(ζ= 1) = -1 (see Equation (A23)). This relation has a first integral (Zotikov, 1986, section 4.1)

(2)

(3)

where superscript b stands for evaluation at the base. We can thus write the formal solution

(4)

In the thermally uncoupled case, w is simply a function of ζand the computation of this solution is therefore just a quadrature, i.e. a numerical integration invoking ζonly and not θ. The upper surface boundary condition can be inserted by evaluating the quadrature at the upper surface ζ= 1 and eliminating θ^b.For plug How, where the horizontal velocity is independent of ζ,w( ζ) = - ζ and we readily obtain

(5)

which is the Robin solution.

For flow by internal deformation, we proceed as follows. We need to demonstrate that we have sufficient information from the EISMINT Benchmark ice-sheet geometries and forcing functions to solve the temperature equation where horizontal velocities arc negligible. A crucial point to understand is that small horizontal velocities do not imply small horizontal gradients of horizontal velocity, which must be large enough to balance the accumulation.

It is shown in the Appendix (Equation (A22)) that the vertical velocity w near the centre of an ice sheet is given by

(6)

where

(7)

is the integral horizontal flux of ice for the interval 0 to ζand u is the horizontal velocity. We term this the partial flux, to distinguish it from the quantity q^s( ζ) = generally known in glaciology as the flux. The horizontal divergence operator ▽_H is taken along surfaces of constant ζ. Implicitly, the velocity, partial flux and flux are all functions of a horizontal vector coordinate ρ, over which the divergence is taken, but we do not show this explicitly in the following derivation, where we are interested in the vertical variation of u and q.

Since

(8)

We can also write

(9)

We now need to evaluate the partial flux ▽_H·q (ζ), and in particular we need to be able to show that if we know ▽_H·q ^s we can thereby compute ▽_H·q (ζ).This is possible under the shallow-ice approximation (Hutter, 1983). Firstly, we need to show how to compute u near the ice-sheet centre in order to compute its horizontal divergence, which is not small, even though u is small. For thermally uncoupled flow we use the Glen rheology together with the shallow-ice approximation to obtain the shear relation

(10)

Where

(11)

A is a rate factor, ρ is the density of ice, g is the acceleration clue to gravity and we are using the reduced constitutive relationship

(12)

(13)

With these assumptions we have, in common with the EISMINT Benchmark Level 1, specifically excluded the divide of an ice sheet, where there is a different mechanical zone (Raymond, 1983).

Now, using Equation (10) and the definition of the partial flux (Equation (7)) we can compute

(14)

which means that at the upper surface ζ= 1,

(15)

We can use Equations (6) and (14) to write

(16)

and since from Equation (15) we know that

(17)

we can substitute this into Equation (16) to arrive at

(18)

This functional dependence has been derived previously (Raymond, 1983). As a check it can be seen easily that this gives the correct answer w( ζ) = - ζ for plug-like shear flow (n →∞). These steps have eliminated our need to know the rate factor; given a geometry obtained from the EISMINT Benchmark experiments, the accumulation rate and the geothermal heat gradient, we have all the information needed to solve the heat Equation (1).

In order to compute the solution for the temperature from the format solution (Equation (4)), we need to integrate w( ζ), and obtain

(19)

This can be substituted into the formal solution (Equation (4)) and the vertical velocity computed by a quadrature.

In the applications below, Riemann sums (a second-order method) have been used. We divide the interval [0, ζ_m] into m — 1 grids, not necessarily of equal size, with nodal points numbered i =(1, m). Then, the quadrature is defined by

(20)

(21)

Since we are dealing with summations, numerical problems arising from finite word lengths are less acute than in difference solutions, meaning that fine grids, leading to very accurate evaluations, are in principle possible. In contrast, in difference solutions of the original equations, overly fine grids lead to the subtraction of very similar values of the temperature, and a concomitant loss of accuracy as a result of finite word lengths. In the applications below we specify m = 1001,which implies a relative error of 10^-6 for a second-order method.

3. The Eismint Benchmark Temperature Calculation

The properties of the advection-diffusion equation for internal deformation are well known and we simply present some examples of the solutions. We take the EISMINT Benchmark Level 1 parameters n = 3, a = 0.3 m a^-1 and consider the thickness computed by the fixed-margin type-I "3d" models (Huybrechts and others, 1996, table 2), where the mean thickness is H = 3419.9 m. The geothermal heat flux is 0.042 W m ^-2 and the thermal conductivity of ice is 2.1 W m ^-1 K^-1, which gives a basal temperature gradient of 0.02 K m^-1. The upper surface temperature is -34.15 K. Figure 1 presents solutions for nϵ {1,3,5,10,30, ∞}; the last case corresponds to the case where the shear layer is infinitesimally thick, indistinguishable from a plug flow, and the temperature is thus described by the Robin solution.

Fig. 1. Illustrating analytical solutions to the ice-sheet temperature equation for various n indicated. "Robin" refers to the plug-flow solution due to Robin (1955).

The EISMINT Benchmark Level 1 temperatures are presented in what is termed a "homologous temperature" θ _E, which in Robin (1955) is defined as

(22)

where γ corresponds to the EISMINT Benchmark β and γH represents the melting-point depression caused by an overburden pressure corresponding to ice thickness H. In those experiments, γ was set to be 8.7 x 10 ^-4 Km^-1,which is the appropriate value for air-saturated ice (Paterson, 1994, p.212). The homologous temperature is more commonly defined as

(23)

where θ_M is the melting point at standard temperature and pressure, and the definition used in the EISMINT Benchmark Level 1 is only approximately consistent with this definition when θ_H ≈ 1 i.e. near the melting point. Indeed, the EISMINT Benchmark definition can lead to negative absolute temperature, and if one wants to avoid this problem and also define the homologous temperature in K, a consistent definition is

(24)

It is not clear that all the modellers were using the same definition of the homologous temperature, which at a temperature of around 264 K would lead to an error of 0.11 K. The results of these various corrections applied to the analytical solution are shown in Table 1. The analytical solution is accurate to at least two decimal places. One should compare the column headed EISMINT Benchmark to the analytical solution headed corrected analytic 1 as these should have had the same corrections applied.

Table 1. Basal temperatures for EISMINT Benchmark Level 1. Results are expressed in K - 273.15, which is roughly equivalent to °C

There is clearly some deviation between EISMINT Benchmark Level 1 and the analytical solutions. The accuracy achievable is a function of the grid interval used in the vertical, something not reported in the EISMINT Benchmark Level 1 experiments. The accuracy achievable is illustrated here using two FD methods and one pseudo-spectral method (Fornberg, 1996) which solve Equation (1). The FD methods have (i) a first-order upwinded difference Operator and (ii) a second-order central difference operator for the advective term, and second-order operators for the diffusive terms. The pseudo-spectral method uses the weight-generat-ing algorithm described by Fornberg and uses "Chebyshev clustering", i.e. a grid density specified by a Chebyshev function. It has spectral accuracy. Eleven points were used in each case.

The results are shown in Figure 2. As expected, the pseudo-spectral method shows very high accuracy, while the second-order FD method shows fair accuracy and gives a good estimate of the basal temperature. The first-order method is not of high accuracy. Closer inspection of the results shows that errors in the first-order FD solution are about one-tenth of the temperature range (23 K), while for the second-order solution they are about one-hundredth; these errors are consistent with error estimates for the FD methods and the grid spacing of Δ ζ = 0.1. Note that the temperature solution does not require node clustering at the base, as the curvature and higher derivatives are low: Taylor-expansion error analysis for FD methods shows that errors are proportional to the higher derivatives. Such clustering might be needed for accurate integrations of the velocity and flux near the base, when there is thermomechanical coupling. Errors for the pseudo-spectral method are <0.01 K. Inaccuracies in the EISMINT Benchmark results, where decreased basal temperatures compared with the analytical solution are reported, cannot be due to the use of first-order methods, as these lead to warmer basal temperatures (Fig. 2); this is predictable, as first-order solutions increase numerical diffusion, which decreases the significance of advection and increases basal temperatures.

Fig. 2. Illustrating solutions to the ice-sheet temperature equation and the analytical solution. FD O(Δζ) and FD O[(Δζ)²] are first- and second-order FD methods. PS is a pseudo-spectral method and "Analytical" is the method described in this paper. Marked points correspond to solution nodes; the FD methods are on evenly spaced grids, while pseudo-spectral nodes exhibit "Chebyshev clustering". Pseudo-spectral and analytical solutions are coincident on the plot and differ by <0.01 K.

4. Computation of Vertical Velocity

It is likely that some of the reported error in the EISMINT Benchmark experiments arises from the numerical computation of the vertical velocity, which in the present paper has been computed analytically. Numerical models generally use some possibly rather disguised form of the integrated continuity relationship

(25)

(26)

to determine the vertical velocity, using numerical integration to approximate the integral defining q.

As was recognised by Hutter and others (1986), there is never any need to perform an expensive double quadrature as might be suggested by Equation (26), as integration by-parts can be used to turn the expression for q into a single integral. For completeness, we treat the case where A depends upon ζ, and define a new quantity

(27)

Then, from Equation (14), we see that

(28)

where, in an obvious notation, we have integrated by parts using

(29)

For the present case, where A is constant, we obtain

(30)

and we infer

(31)

or

(32)

To investigate quadrature error, we approximate these integrals by the Riemann sums, with the same discretisation used in Expression (20). We obtain

(33)

and

(34)

(35)

For ten grids (i.e. 11 points) the two numerical methods of evaluation (Equations (33) and (34)) gave identical answers, and the maximum absolute difference in the velocities between the numerical integration and the analytical solution at any point was <1 %. Since the grid spacing is 0.1, this error magnitude is consistent with the fact that Riemann sums are second-order accurate. Maximum relative error, where the difference between the two velocities was normalised by the velocity, was around 5%, but this occurred at the first gridpoint from the base, where vertical velocity is small (-0.023) and advection therefore does not contribute signifi-Cantly to heat transport. To put it another way, we might expect the relative error in points near the base to be higher, as a relatively coarser grid has been used in their computation (one gridpoint in this case). However, the small velocity near the base mitigates this error. From these results, it seems unlikely that error in the numerical integration of the vertical velocity led to the reported errors in the EISMINT Benchmark results.

Are there any aspects of the problem which might tend to lead to programming error? Straightforward application of Jenssen's (1977) "σ" transform (i.e. using a normalised vertical coordinate) of the heat equation with a computation of the vertical velocity based on the continuity Equation (A4) and the numerical integrations in Equation (33) or (34) leads to quite complicated formulae for the vertical heat transport (see Equation (A15)), and it is conceivable that this has resulted in programming errors in previous work. In fact, combination of continuity equation and the σ transform allows cancellation of several terms; this has already been noted for plug flow by MacAyeal (1997) and Hulbe (1998). We present here a relatively simple form of the advection equation in normalised coordinates which generalises the MacAyeal-Hulbe formula to include flow by internal deformation.

The heat-transport equation in normalised coordinates can be written

(36)

where q^s is the total flux through the thickness (i.e. q^s(ρ) = q(ρ,ζ = 1)· Its derivation is quite lengthy and is given in the Appendix. The MacAyeal-Hulbe formula is obtained for the special case of plug flow, since for such a flow regime q = ζq^s because horizontal velocity is independent of vertical elevation. Equation (36) has advantages for thermomechanically coupled flows, since the general effect of increased basal temperatures is to make flows more pluglike, meaning that the dominant term in the vertical heat transport is usually aζ. Errors in Computing the fluxes by numerical integration should therefore be of less significance in computing the vertical heat transport. An alternative form of Equation (36) is

(37)

Equations (36) and (37) also lead to simple design of finite-element schemes for computing (he vertical heat transport (see Hulbe (1998) for a discussion of the difficulties). With the present formulation, flux divergence terms in the vertical heat transport can be multiplied by a test function ϕ and Green's theorem can then be applied to obtain weak forms like

(38)

where the areal integrals are over the element, n represents the outward-pointing normal to the boundary of the element, and the second term on the righthand side is the integral around the boundary Γ of the element. This expression obviates the need for elements with special properties in the horizontal in order to allow computation of the vertical heat transport. The other terms in the heat equation can be treated by the usual methods (sec MacAyeal, 1997; Hulbe, 1998).

An important point about Equation (36) is that it naturally leads to consistent computations of the total flux (used in the equation describing the evolution of the ice-sheet profile) and the vertical heat transport through the computation of the partial flux. Finally, we show that only two quadratures are necessary in the computation of the flux divergence for the ice-sheet profile and the heat-transport equation. From Expressions (10), (27) and (28) we can see that

(39)

(40)

(41)

(42)

We approximate these integrals by the Riemann sums, with the same discretisation used in Expression (20).The integral in Equation (40) can be computed by the Riemann sum

(43)

while the second term on the righthand side of Equation (42) can be computed by the Riemann sum

(44)

In these expressions A_i__1/2 refers to A evaluated at the point ζ_i-1/2. Obviously one performs the discretisation such that the quadratures for successive ζ_m can be computed recursively from the preceding value, e.g.

(45)

etc.

5. Discussion and Conclusions

A quadrature solution for ice-sheet temperature equation for flow by internal deformation, corresponding to Robin's (1955) solution for plug flow, has been found. This solution is of much higher accuracy than the EISM1NT Benchmark Level 1 series. Good accuracy of numerical solutions can be obtained with central difference estimates of the temperature gradient, and very high accuracy can be obtained by using pseudo-spectral methods. Overly fine discretisations of the vertical transport equation using higher-order schemes can lead to numerical artefacts (specifically, "wiggles") but this does not seem to be a problem for coarse grids.

E1SMINT Benchmark 1 temperature results do not seem to be as good as they should be; the error is larger than one would expect. It is unlikely that numerical errors resulting from integrating the vertical velocity are the cause of this. It is remarked upon that the complicated expressions resulting from evaluation of the vertical velocity and substitution of this into the σ transform can lead lo programming errors, and a simpler formula for heat transport in a σ coordinate is presented. This formula shows that only two quadratures are needed to compute the ice-sheet evolution and the heat transport.

There is independent evidence that the computed ice-sheet geometries in the EISMINT Benchmark Level 1 were accurate, as codes which gave typical results in these experiments also gave good results when compared against a map-plane analytical solution for the case n = 1 (Hindmarsh and Payne, 1996).