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):

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.

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

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

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 (Reference Johnson and McMeekingJohnson 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

and at the base

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,

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

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

with

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:

and at Ζ = H(X) :

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)

For example, the secant method (Reference Yakowitz and SzidarovszkyYakowitz 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.1. Divide analysis

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

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

with boundary conditions

is obtained. In the above

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

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

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

where

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

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

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),

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

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

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 τ

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

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)

and therefore, since Δ is known at the divide,

Furthermore, we know that

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

or, after re-arrangement,

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

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

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

Thus, all functions at level Ζ are determined, since

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.