Computational Ice-Divide Analysis of a Cold Plane Ice Sheet Under Steady Conditions

Abstract

The dimensionless form of the field equations and boundary conditions governing plane flow of a grounded cold ice sheet emerge from balance statements of mass, momentum, and energy. They constitute an amended version of a reduced model of ice-sheet flow, due to Morland (1984) and Hutter (1983), and circumvent the restrictions imposed by the reduced model, namely the neglect of the longitudinal stretching effects. The amended version permits satisfaction of mass balance at the ice divide for arbitrary basal sliding conditions and gives a better reproduction of the local flow features.

Under very mild simplifying assumptions, namely that horizontal thermal conduction can be ignored close to the divide, we present a numerical analysis of the ice divide which has second-order accuracy. This analysis permits determination of the temperature profile, velocity, and stress distributions in a symmetric ice divide, provided that the ice-divide height, the local behavior of the accumulation and surface-temperature functions, and the geothermal heat flow are prescribed.

1. Introduction

Ice divides are areas where important surrogate information about the climate or glacier history may be obtained. For example, they are the best areas for obtaining cores for isotope analysis (Raymond, 1983), and their flow regime can have a profound influence on erratic trajectories (Boulton and others, 1985), especially if one includes saddle points under the class of divides. However, ice divides present certain mathematical difficulties which (i) prevent the reduced model of Morland (1984) and Hutter (1983) from describing their evolution and characteristics, and (ii) require high-order precision in their numerical analysis when reliable local solutions of the stress, velocity, and temperature distributions are being sought.

Steady-state conditions of the reduced model were numerically studied by Yakowitz and others (1986) and Hutter and others (1986, 1987). These solutions were constructed with the aid of a marching procedure by starting at the divide, selecting a divide height, performing an ice-divide analysis, and proceeding column by column up to the margin (at which the height must vanish) was reached. The pre-selccted ice-divide height was varied until the integrated accumulation function vanished for the solution to satisfy total mass balance. But, unfortunately, a large amount of basal slip was required for the predictions to be accurate. The reason for the inadequacy of the model was the neglect of longitudinal stretching effects, which we know to be significant in the vicinity of ice divides (Raymond, 1983).

Our new improved model equations incorporate the longitudinal stretching effects. They are derived from the exact continuum equations by ignoring the horizontal thermal diffusion in the heat equation (this term is known to be small). With this model, strong basal friction is now permissible, but a finite-difference solution scheme must possess at least second-order accuracy to admit an acceptable local solution. Even higher-order accuracy is needed if a complete ice-sheet solution is to be found.

We will demonstrate this analysis for plane-flow conditions and a locally symmetric ice divide. The analysis of a dome of a locally axisymmetric ice sheet is analogous and only requires minor changes. Basic equations will only be briefly stated and not derived, as they can be found, for instance, in Hutter (1983).

2. Governing Equations

The equations of the model prior to any approximations are as follows.

2.1. Field equations

When restricted to plane flow, the dimensionless velocities U, W, stresses σ_x, σ_z, τ, and the temperature Τ are governed by the equations (see Fig. 1):

(1)

Fig. 1. Symmetrie plane ice sheet indicating the coordinates, the Tree surface, and the bottom profiles, including a sketch of the finite-difference approximation close to the divide.

Table 1a Scales used to non-dimensionalize the field equations and boundary conditions. Numerical values are typical for the situation. A range of values is discussed in the text

Table 1b Physical parameters

Table 1c Dimensionless parameters based on the above scales

(In the reduced model equations v = ε² and θ = Ψ/s.)

Table 2a MATERIAL RESPONSE FUNCTIONS FOR POL YCR YST ALLINE ICE

Table 2b Accumulation-rate and surface-temperature functions (examples)

Table 2c. Basal Sliding Law and Geothermal Heat Flow

Fig. 2. Typical shape of the temperature distribution at the divide (Z = 0 corresponding to the bottom, Ζ = Η_D to the top).

where a(Τ) is a temperature-dependent rate function, w(τ_II/s²) is a creep response function, and

(2)

is the second stress-deviator invariant (see Hutter and others, 1986). s, v, and ε are dimensionless constants, which are independent here but were given by s = εθ^−1/2 and v = ε² in the reduced model. We still may assume that ε and v are smaller than unity and that θ is 0(1) or larger, but an inter-relationship is not presently suggested. In addition, α and β are given constants; they are called an energy-dissipation and a thermal diffusion number, respectively. In terms of physical constants and typical scales, these constants are given in Table I. Numerical values have a considerable spread.

In Equations (1), statements (a) and (b) are a horizontal force balance, (c) is the incompressibility condition, while (d) and (e) are constitutive relations of stress relating strain-rate in a non-linear fashion to the stress deviator. The non-linearity is expressed by the temperature-dependent rate factor a(T) and the creep response function w(τ_II/s²). Finally, (f) is the heat equation that incorporates advection, thermal diffusion, and strain heating. As shown by Morland (1984), all these terms are significant in thermo-mcchanical coupled ice-sheet dynamics.

Rate factors a(T) between 263 and 273 Κ satisfy a polynomial relation a(T) = P₃(T). Likewise, w(τ_II/s²) is often a power law, but numerous laboratory and field studies (see Hutter, 1983, chapter 2) have indicated that polynomial expressions fit available data better. We shall adopt a “finite viscosity” expression to avoid mathematical singularities at the free surface (Johnson and McMeeking, 1984) but leave both a(T) and w(τ_II/s²) arbitrary to adjust for local constitutive behavior. A choice is given in Table II.

The constitutive relations of Equations (Id and e) do not incorporate any stress-induced anisotropics (i.e. directional dependencies). To our knowledge, such anisotropy effects have so far not been incorporated in any glaciological context. However, inhomogeneitics (i.e. positional dependencies) can easily be incorporated in the rate factors as well as the creep-response function by making a and ω also position-dependent. So, our formulation permits incorporation of the differences of the softening effect of the dust particles in Pleistocene ice by an internal variable (see Hutter and Vulliet, 1985) or simply by making a and w position-dependent.

2.2. Boundary conditions

The field Equations (1) and (2) are complemented by boundary conditions at the free surface

(3)

and at the base

(4)

In Equations (3), (a) is a kinematic statement expressing the evolution of the surface profile; (b) and (c) are surface-stress conditions relating σ_x, σ_z, and τ to the atmospheric pressure but ignoring the small wind-stress term. Relation (d) is the Dirichlet condition of temperature, i.e. we prescribe the surface temperature through T_s(X,H,t).At the base, a sliding condition is applied. This is formally implemented in Equations (4) by prescribing the horizontal velocity component (Equation (e)) and requesting tangency of the sliding velocity along the base Ζ = F_B(X), expressed in Equation (f). Statement (g) equates the heat flow from the basal surface into the ice to the dimensionless frictional heat plus the dimensionless geothermal heat flow

which is positive as a flow towards the ice. Orders of magnitude for

are 0.1–1.0. The term involving γ represents generation of frictional heat along the base; γ is a dimensionless quantity with the order of magnitude ≤ 10²; it is also defined in Table I.

We emphasize that the boundary condition in Equations (3)-(4) are exact statements valid for any value of ε. The climatology input is provided by p_atm (usually = 0) and the two forcing functions A_cc(.) and T_s(.), and the interaction with the substratum is described by the geothermal heat flow

and the sliding law U_F(.).

The accumulation rate has generally a strong dependence on H. Denoting by H_E the snow-equilibrium height, we have A_cc < 0 if Η < H_E and A_cc > 0 if Η > H_E with A_cc = 0 at Η = H_E. A polynomial dependence on the variable (H - H_E) may be appropriate, but an X- or t-dependence can also be incorporated by choosing H_E = H_E(X,t).

The surface temperature T_s is the temperature measured 10 m below the surface. It also has a strong dependence on height, which in a first approximation can be assumed to be linear.

The geothermal heat flow

is often assumed to be constant (corresponding to 3°C/100m). However, it may equally be assumed to be a function of position and time. Computationally, it is also advantageous to prescribe the geothermal heat through a “virtual” horizontal plane,

(5)

rather than

. Values of Q^geoth(X, F_B(X)) are generally larger in the central region than towards the margin. This follows simply from the insulating effect that is exhibited by the ice sheet.

The sliding law is probably the most complicated of all. It must be selected with caution because both physical and mathematical arguments will determine its form (Hutter, 1983). For a finite slope margin (wedge type), one must have

as the margin position is approached. In other words, the sliding velocity is linear in the shear traction and the overburden pressure.

The explicit forms of all these functional relations used in this paper are given in Table II.

We now give a heuristic proof that Equations (1) can circumvent the restriction that basal sliding must be large. To this end, recall that, from Hutter and others (1986), the curvature of the surface profile at the ice divide, computed according to the reduced equations is given by

(6)

where

and Η_D = H(0), F_B = F_B(0), σ_X(0) = σ_Z(0), and thus, τ_II(0) = 0 since τ(0) = 0. With τ_II = 0 at the divide, w = w(0) which vanishes for a power flow law implying GLIDE = 0. This necessarily requires SLIDE ≠ 0. We have chosen w(0) ≠ 0 and thus would make H″ dependent upon a “finite viscosity at zero stress”, a quantity which is poorly known. In addition, GLIDE is very small which again requires that SLIDE is of order unity. The point is that the shallow-ice approximation yields doubtful results close to the divide. For the improved model, longitudinal stretching is significant so σ_X(0) ≠ σ_Z(0) cannot vanish even with τ(0) = 0.

Balancing both sides of Equation (1d) then requires

so that Equations (lc and d) with

At a symmetric ice divide, Equations (la) and (le) are identically satisfied, and we may write in accordance with mass balance

and then deduce from above

This relation defines both ∂U/∂x and

at the divide, if T and Η are prescribed. With

, one has τ_II ≠ 0, w(τ_II) ≠ 0, and thus GLIDE ≠ 0. In fact, we shall demonstrate that w ≤ 0(1) under usual circumstances and so SLIDE may vanish without making

infinite.

Henceforth, steady-state conditions will be considered and time derivatives in Equations (If) and (3a) will be omitted.

2.3. Numerical solution strategy

Suppose we ignore the terms

in the energy Equation (If). The soundness of this approximation can be tested a posteriori. Let the finite-difference approximation of the derivative of f with respect to X be denoted by δf. (We shall define δf below.) Then, Equations (1) and (2) may be written as

(7)

with

(8)

Note that Equation (7a) is a re-arrangement of Equation (Id) and Equation (7b) is a combination of Equation (7a) and the definition of τ_II, in Equations (2).

The boundary conditions are as follows:

(9)

and at Ζ = H(X) :

(10)

Observe that for any fixed value of x Equations (7) and boundary conditions (9) and (10) define a two-point boundary value problem (TPBVP) for the unknown functions τ, σ_z, W, U, T, and T′, provided that the approximating derivatives δσ_x, δτ, δU, δW, δT, and δF_B are known.

In order to take advantage of this property, the solution strategy is as follows. For any given value of the divide height Η_D, we shall proceed in three steps:

Step 1. All unknowns will be determined at the divide and in the column next to the divide. (That is, for x = 0 and for x = Δx.)

Step 2. Using a marching procedure, the unknown functions will be determined for all x = kΔx(k Δ 2), such that H(kΔx) > 0. The margin position x_M(H_D) is reached when for the first time H((k + 1)Δx) ≤ 0.

Step 3. The divide height H_D is determined by the condition that the integrated squared accumulation must vanish. That is, the value of H_D is obtained by solving the equation (T_a is called the target function)

(11)

For example, the secant method (Yakowitz and Szidarovszky, 1986) can be used for the numerical solution of this equation. The advantage of using the secant method is the fact that this method does not use derivatives of the target function. Note that the computation of the target function at any value of Η_D requires the completion of steps 1 and 2 for this fixed value of Η_D.

In this paper, we shall not discuss the marching process in general; we shall leave this step to a succeeding paper. But step 1 will be discussed in detail.

We now demonstrate how the ice-divide analysis is performed.

3. Numerical Treatment of the Ice Divide

For simplicity, a symmetric configuration will be considered.

3.1. Divide analysis

At a symmetric ice divide the following symmetry conditions must be fulfilled:

(12)

Substituting these relations into the steady-state versions of the governing field Equations (1) and boundary conditions (2) and (3), the following system of equations

(13)

with boundary conditions

(14)

is obtained. In the above

(15)

and W^(0)(Z) is a prescribed function. Our initial — and later improved — ice-divide analysis will therefore be based on an approximate representation of the vertical velocity profile, the first estimate being based upon the linear relation

(16)

Note also that, since w is increasing in Δ² and Δ is a multiplier, there is a unique solution Δ of Equation (13e) for any values of H(0) and T. For the suggested numerical procedure this is important. The procedure is now as follows:

Equations (13c, d, and e) with the definitions in Equation (15) and the boundary conditions in Equations (14a and b) yield a TPBVP for T′ and T. In the solution of this TPBVP, Δ is obtained by solving the non-linear Equation (13e) using thereby the estimate in Equation (16); τ_II then follows from Equation (15a,b). Finally, we may recall that U = τ = 0 at the divide.

Note that the functions σ_X and σ_Z are not determined; only their difference, function Δ, is obtained. As we shall see later, in the column next to the divide a simultaneous computation of all variables with those at the divide is possible which makes all the variables available.

3.2. Second-order derivative approximation

We list here without proof a first X-derivative approximation of second order, i.e. the error is 0(ΔX²). Let X = 0 be the point of symmetry. Consider the derivative ∂f/∂x Then to second order,

where

(17)

This formula will be used repeatedly in what follows. We also need the backward second derivative in Ζ

where

(18)

In the above, even functions satisfy the property f(X) = f(-X) while for odd functions one has f(X) = -f(-X).

3.3 Column next to the divide

Recall that at the divide we know the functions W, T, T′, Δ, and U = τ = 0.

Assume now that X = ΔX. Then, up to second order, the even functions are the same at the divide and in the column next to the divide. Thus

(19)

The other functions U, τ, σ_x, and σ_II at the divide and next to the divide can be obtained simultaneously in several steps.

Step 1

From Equation (7c) and relations (17c and d), we may deduce the approximations

(20)

Thus, U(ΔX,Z) is now known for all Z.

Step 2

By using the boundary conditions in Equations (3) at Ζ = H(0), the values of σ_x(ΔΧ,Η(0)), σ_z(ΔΧ,Η(0)), τ(ΔΧ,Η(0)) can be determined with an estimated value of H′(ΔX). (H′ (ΔX) evaluated from Equation (3a) would be zero if second-order accuracy is used.) Then Equations (3b and c) and the definition of Δ imply for X = ΔX and Ζ = H(0),

(21)

all of which can be computed.

Step 3

We iterate on Equations (21) by constructing improved values of H′. From a Taylor-series expansion of H′ we have

(22)

From a Taylor-series expansion in X of the steady version of Equation (3a), it follows that

(23)

The derivative terms of A_cc(.) can be computed to any desired order of numerical accuracy, because this function is prescribed.

is given by Equation (16), but

is still unknown; it can be determined in the following way.

By differentiating Equation (le) with respect to X and imposing the symmetry relation τ

(24)

so that with the aid of Equations (17c) and (18) we obtain for Ζ = H(0)

(25)

With Equation (25), an improved estimate of

and therefore that of H′ (ΔX) is computed as given by Equations (22) and (23). Step 2 then provides improved values for σ_x(ΔΧ,H), σ_z(ΔΧ,Η), τ(ΔΧ,Η).

Step 4

Similarly, at the divide we have from the boundary condition (14c)

(26)

and therefore, since Δ is known at the divide,

(27)

Furthermore, we know that

(28)

Step 5

Starting from the free surface, we show how a marching procedure down the two columns at the divide and next to it can be developed to obtain all these functions σ_x(0,Ζ), σ.(0,Ζ), σ_x(ΔΧ,Ζ), σ_x(ΔΧ,Ζ), and τ(ΔX,Ζ) for all values of Z. For Ζ = H(0), they are all known. Let now Ζ be given, and assume that all the above functions arc known for Ζ + ΔΖ. Then the procedure is as follows: in Equation (7b), the respective derivative approximations of

and δτ(ΔX,Ζ + ΔΖ) are substituted. Together with Equation (19c) this yields

(29)

or, after re-arrangement,

(30)

Note that Equation (29a) is only a first-order derivative approximation, but after multiplying it by ΔZ and re-arranging the terms in order to obtain relation (30), we finally obtain a second-order approximation of σ_z(ΔΧ,Ζ). Relation (29b) is itself a second-order derivative approximation, and so Equation (30b) gives a third-order approximation for σ_z(ΔΧ,Ζ).

Next, since

is even, one has up to second order

Consequently, Equation (30) also holds at X = 0, so that

(31)

With this we deduce, since Δ(0,Ζ) is already known,

(32)

Finally, function τ(ΔΧ,Ζ) can be updated as follows: in Equation (7a), δσ_x and ∂τ/∂Z at (ΔX,Ζ + ΔΖ) can be approximated. Then, as for relations (30a and b), we obtain

(33)

Thus, all functions at level Ζ are determined, since

(34)

By this marching process, all variables at the divide and in the column just next to the divide are computed.

The procedure described above is based on a prescribed form of function W(Z) of the divide, which was given by Equation (16). Note that by starting from any feasible function W(Z), the above process can be repeated easily. For any function W(Z), the quality of the fit of the resulting solution functions can be measured by adding the squares of the relative errors in satisfying Equations (7) and boundary conditions (9) and (10) at X = 0 and X = ΔX. Then, the corrrect choice will be that function W(Z) which minimizes this overall added squared error. We refrain from presenting any details and pass directly to the presentation of numerical results.

4. Numerical Example

The procedure described in sections 3.1–3.3 was tested in a hypothetical numerical example. The actual values of the parameters were selected as

The step sizes in both directions, X and Ζ were selected as

In order to obtain a feeling of how the numerical results depend on the values of the above parameters, we systematically changed the value of one of the above parameters to obtain six cases, which were as follows:

Case 1: Original data

Case 2: H_D = 0.9

Case 3: α = 0.3

Case 4: β = 0.041

Case 5: v = 0.003776

Case 6: s² = 0.000408.

The temperature distributions for these cases are tabulated in Table III, and a typical shape (for original data) of function T(0,Ζ) is drawn in Figure 2. For the

Table 3 Temperature distribution on the divide (values in centigrade are obtained by multiplying the numbers by a factor of 20)

Fig. 3. Distribution of the longitudinal stresses σ_x and σ_z at the divide.

same case, functions σ_x and σ_z are shown in Figure 3, and Figure 4 illustrates function τ_II. These give an indication of the significance of the longitudinal stretching effects.

The accuracy is governed by the accuracy of the shooting-method computation for determining the solution of the two-point boundary-value problem in Equations (13c, d, and e) and (14a, b, and c). The error at Ζ = Η_D was always about 10⁻⁶–10⁻⁷.

Acknowledgement

We thank Professor D.R. MacAyeal for his critical review of an earlier version of this paper.

Fig. 4. Vertical distribution of the second stress-deviator invariant at the divide. (In the reduced model σ would vanish.)