Thermo-mechanically coupled ice-sheet response — cold, polythermal, temperate

Abstract

Classical mixture concepts are the appropriate vehicle for describing the dynamics of ice masses containing some water. We review and derive, respectively, the theoretical formulations of cold, polythermal and temperate ice masses, emphasize the peculiarities of the model equations and point to difficulties that were encountered with the proposed models. The focus is both on the adequate physical motivation of the models and the consistency of their mathematical representation. The paper also has a tutorial character.

As usual, cold ice is treated as a single-component incompressible heat-conducting viscous fluid, while two different models are presented for temperate ice. When it arises in a polythermal ice mass, the water content is small and a simple diffusive model for the moisture content suffices. This diffusive model is further simplified by taking its appropriate limit, when the moisture diffusivity tends to zero. Temperate ice in a wholly temperate — Alpine — glacier is treated as a two-phase flow problem, i.e. the momentum-balance laws of both constituents ice and water are properly accounted for. Such Darcy-type models are suggested because the water arises in a greater proportion; so its dynamic role can no longer be ignored.

The constituent ice is treated as an incompressible non-linearly viscous isotropic body with constitutive properties similar to those of cold ice. The interstitial water is a density-preserving ideal or perfect fluid. The two interact with an interaction force that is proportional to the “porosity” and the seepage velocity. Internal melting that arises will lead to a generalization of the familiar Darcy law.

When water is present, the boundary and transition conditions across internal singular surfaces take special, more complicated forms and involve statements on drainage to the base. These conditions are also discussed in detail.

List of Symbols

1. Introduction

In computations of glacier and ice-sheet flow, glaciologists generally employ simple continuum models with nonlinear, temperature-dependent rheological properties. This leads to the well-recognized fact that, in cold ice, the temperature distribution affects the ice velocity. What is not so easily recognized, however, are the unique complications which arise when ice temperature in these models reaches the melting point. The thermo-mechanical processes which determine the distribution of wet, temperate ice are not easily simplified. Thus, ice-sheet and glacier models typically resort to simple numerical tricks as means to enforce the upper bound of ice temperature. While this pragmatic approach can have its place in studies of short-term ice-sheet behaviour, it is woefully inappropriate in long-term climate studies where the growth and decay of large ice masses depend fundamentally on the thermodynamical behaviour of temperate ice.

In this paper, I shall pull together the experience gained from the derivation and use of models which explicitly treat the thermo-mechanical processes of cold, polythermal and temperate ice masses. In part, my goal is to review the theoretical work carried out over the past 10 years which explicitly addresses the roles of water in basal sliding and of partial melting in ice deformation. This review will describe various competing theories and will discuss their strengths and weaknesses. Specifically, the following section discusses cold ice, how it is treated and where amendments are needed. Then, we turn our attention to polythermal ice and outline the various theoretical formulations by which it is modeled. Beyond this review, I shall present an entirely new description of temperate ice which pulls together the various disparate elements of the previous theories.

Essential to the discussion of cold, polythermal and temperate ice masses is the appreciation of the many roles water plays in ice dynamics. In fact, there exist numerous studies on the significance of water in the basal-sliding mechanism. Sliding with and without cavity formation has been treated with expertise, both physically and mathematically; soft-and hard-bed concepts have been introduced, and sediment-bed-layer instabilities and formation of linked cavity systems, respectively, have been suggested as possible causes of the “catastrophic” surge-type glacier advances; glacier hydrology is recognized as playing an important role, etc. However, a true, interactive dynamic coupling of the motion and deformation of the ice and the seepage-type flow of the water through it, and influencing it, is generally not treated.

This description makes it plausible that the role played by the water in ice may be described by more than a single model. When the water content is small, which is the case in the temperate part of a polythermal ice mass, then a considerable amount of it diffuses along the grain boundaries through the ice. Thus, an adequate model in this case may be simply a diffusive model, in which the balance of moisture mass is governed by its production and diffusion. The water mass within the ice and its motion is too small to play a significant dynamic role. Such a diffusive model has been derived by Fowler and Larson (1978), was corrected by Hutter (1981), and left untouched until Heinz Blatter resurrected it and made an attempt to compute cold-temperate-transition surfaces in glaciers of the Canadian Arctic (Hutter and others, 1988; Blatter, 1991; Blatter and Hutter, 1991).

On the other hand, when the water content is larger such that the pore space, moulins, crevasses, etc. are filled up to a certain level with water, then a more detailed description of the dynamics played by the water is necessary. Fowler (1984) derived such a theory in which moisture flux was no longer treated by Fickian diffusion but rather as a Darcy-type dispersive mechanism, and the role played by the salt impurities was incorporated. Fowler derived his model for polythermal ice and did not stress its applicability to wholly temperate ice masses; that was done by Hutter and Engelhardt ( 1988a, b). However, these models are still flawed somewhat and in particular do not fully treat the boundary and transition conditions.

2. Cold Ice Sheets

The theory of cold ice has been understood both physically and mathematically for over 20 years. Nevertheless, cold ice masses have been the center of interest in the last decade as a result of the need for a systematic derivation of simplified equations. This systematic derivation is neccessary to deduce rationally from “first” principles the glaciologist’s first commandment: “Basal shear stress equals ice density times gravity times overburden depth times surface slope.“ Equally important is the rational derivation of the equations used to describe the growth and decline of cold ice sheets through the various ice ages when the climatic forcing is prescribed. To achieve these principal goals, an asymptotic analysis was required to determine the effects of appropriate depth, length, stress and velocity scales. This scale analysis, put forward in various steps by Fowler, Hutter and Morland was not immediately accepted, as can, for example, be seen from the title of one of my papers in 1981 involving it: The effect of longitudinal strain on the shear stress of an ice sheet. In defense of using stretched coordinates (Hutter, 1981). (I obviously had to overcome the resistance of a stubborn referee. The second part of the title was only added in the revised version of the paper.) The ultimate product of this analysis is the reduced lowest order theory (Hutter, 1983; Morland, 1984; Hutter and others, 1986), referred to as the shallow-ice approximation. Its equations are today used by a number of (younger) scientists in attempts to quantify the role of the large ice sheets through the various climates (Hutter and others, 1986; Hindmarsh and others, 1987, 1989; Herterich, 1988, 1990; Hindmarsh and Hutter, 1988; Calov, 1989, 1990; Huybrechts, 1990a, b; Huybrechts and Oerlemans, 1990; Letréguilly and others, 1990a, b,c). These equations and analogous ones for polythermal and temperate ice will certainly play their role in future analyses of general circulation models when atmosphere-ocean-cryosphere couplings will have overcome the numerical difficulties they are presently still exposed to. The scale analysis, however, also offers the possibility of improving the model equations (in perhaps other situations) by employing a systematic perturbation procedure of the scaled equations.

2.1. Equations

In cold ice, the temperature of any particle is everywhere below the melting point and changes according to the heat that is advected and conducted. Constitutively, cold ice in a glacier or ice sheet is postulated to be a viscous heat conducting, incompressible non-Newtonian fluid. Note that this postulate prevents a priori the modeling of stress-induced anisotropies. The local balance laws of mass, momentum and energy, and the constitutive relations, are in this case (for definitions and notations, see the list of symbols):

(2.1)

Here, accelerations in the momentum equations have been disregarded, the internal energy has been assumed to be a function of temperature only, and Fourier’s law of heat conduction has been used. The last term in the energy equation is the dissipation, tr (tD), and Equation (2.1)₆ is the non-linear stress-stretching relationship with a temperature-dependent rate factor A(T) (usually of Arrhenius type) and a creep-response function

, (2.2)

with appropriate values for a_j (j = 1,..., N) (note we have set t _II=x). We refrain from adhering to the power law for two reasons: first, a finite viscosity law (a ₀≠0) is essential at low stresses* and also has mathematical advantages when stress is expressed in terms of stretching and, secondly, only minimal simplifications emerge when the power law is imposed. Furthermore, the rate factor A(T), though it is written as a function of absolute temperature, is actually a function of the homologous temperature, i.e. the difference between the absolute temperature and the melting temperature T _M.

Equation (2.1) and (2.2) are valid in the ice-sheet domain D. Along its boundaries one has, on the free surface •D_f defined by F_f(x,t) = 0:

(2.3)

and along the base •D_f,defined by :†

(2.4)

The first of these are kinematic statements and express a local mass balance and the tangency of the ice velocity to the basal surface. Equation (2.3)₂ requires the free surface to be stress-free, while Equation (2.4)₂ connects the basal sliding velocity to the shear traction with a coefficient that may depend on the absolute value of this shear stress and upon the pressure acting perpendicular to the surface. The functional form of this sliding coefficient has been the center of activity through the last three decades of glaciological research, especially in the presence of water (see, for instance, Weertman, 1957, 1961, 1964, 1967, 1971, 1979; Lliboutry, 1964, 1968, 1975, 1979, 1987; Nye, 1969, 1970; Kamb, 1970; Morland, 1976a,b; Fowler 1981a, b, 1986, 1987; Shreve, 1984); to the numerical modeler these works are only marginally useful as he needs an explicit functional relation for C. The thermal condition in EquationS (2.4)₃, ₄ is of the Neumann type if the basal temperature does not reach melting, or of the Dirichlet type if the melting temperature

(2.5)

is reached. Here T ₀(= 0.01°C) is the melting temperature at p ₀(= 630 Pa) and c is the Clausius-Clapeyron constant. The climalological input enters through the accumulation-rate function (also called mass or snow balance) and the surface temperature T _f( x ,t).

2.2. Solution strategy for cold ice masses

The above-stated initial boundary-value problem comprises a complicated free boundary-value problem for the domain geometry {D(t),•D _f(t)} and the temperature and velocity fields. To solve this problem, either approximations must be pursued or numerical integration techniques employed. As for the latter, an iterative solution technique is advantageous and now used by nearly all numerical modelers. The idea is as follows:

• (1) For a fixed time t, assume that the domain D(t) and the temperature distribution T(x,t) within it are known. We shall explain below how this step is initiated.

• (2) Use the momentum Equation (with the prescribed temperature distribution) Equation (2.1)₂ with the boundary conditions in Equations (2.3)₂ and Equation (2.4)₂ to obtain the velocity field v in D(t). This is a non-linear elliptic boundary-value problem with mixed boundary conditions. Since it is not self-adjoint, care is taken in the selection of the solution algorithm.

• (2) Determine the new domain D(t+Δt) by integrating in time the kinematic Equation (2.3)₄ subject to the velocity distribution .

• (3) Update the temperature field by solving Equation (2.1)₃ subject to the boundary conditions in Equation (2.3)₃ and (2.4)_3,4 . This problem, parabolic in time and elliptic in space, yields . (Since this step uses v(x, ί), a further iteration of the velocity and temperature distributions may be needed at fixed D(t+Δt). )

This procedure requires initiation of the velocity and temperature fields in D(t₀). The temperature may simply be prescribed (because for Τ an initial value problem is solved), v(x,t₀) follows from the solution of the momentum equation at given T(x,t₀). Often one assumes an initial temperature distribution which corresponds to a steady temperature field corresponding to this velocity distribution. It is no more physical than any other initial temperature distribution. To my knowledge, this general thermo-mechanically coupled Stokes-flow problem with freely evolving boundary (and not imposing the shallow-ice approximation) has for the first time been solved but not sufficiently exploited by Shilun Qin and K. Hutter (paper in preparation) by using the finite-element method. To be sure, there are commercially available codes for three-dimensional time dependent Stokes (low or heat conduction, and these can be coupled to a kinematic equation for an evolving surface; however, these are generally not optimally designed for the kinematic and rheological peculiarities of the problem at hand. Qin has been able to symmetrize the emerging matrix equations in his FE-solution without introducing the adjoint operator equation to the above system, which speeds up computations; for details see paper in preparation by Shilun Qin and K. Hutter. Since required CPU times are generally high, optimization of the computations in this regard is essential.

Finally, computational routines of the above equations for arbitrary geometry (without the imposition of the shallow-ice approximation) are needed in three-dimensional flow problems of cold glaciers. Furthermore, it is known that the shallow-ice approximation is invalid at the ice margin, at ice-sheet-ice-shelf transition zones (near grounding lines) and at ice divides (domes). However, it is yet to be seen whether such models can compete with the simpler finite-difference models for the shallow-ice equations in a simulation of ice-sheet formation and decay through ice ages.

2.3. Shallow-ice approximation

The dimensional analysis leading to the shallow-ice approximation has been performed several times in the past (Hutter, 1983; Morland, 1984, Hutter and others, 1986). We shall not repeat it here; instead, we give a heuristic derivation of the reduced equations. This is what glaciologists do anyway (e.g. Herterich, 1988, or Budd and Jenssen, 1989). However, I wish to stress here why applied mathematicians do emphasize such scaling arguments: they naturally suggest an entire hierarchy of approximations and allow identification of the ranges of validity for each of these. We shall shortly come back to this point.

Basic to the shallow-ice approximation is the recognition that ice sheets are long and wide but shallow so that variations of any field quantity in the horizontal direction are small in comparison with variations in the vertical direction:

,

for Some fields g.

This assumption can then be used to prove that

.

Thus, there is negligible horizontal shearing and the normal stresses are isotropic and reduce to

.

In other words, there are no materially dependent normal stresses in the shallow-ice approximation; normal stress effects are ignored and can be accounted for only in a higher-order model approximation. With the above assumptions, the reduced field equations are

(2.6)

(2.7)

and along the base, defined by z=h _f(x,y,t):

(2.8)

These equations comprise the shallow-ice approximation and are substantially simpler than the original non-linear Stokes-flow equations. Evidently, the pressure distribution is hydrostatic (Equation (2.6)₄ ), which makes the drastic simplifications possible. (It is at this point where it is probably most easily recognized that longitudinal stress or stretching effects are ignored.) Upon formal integration the stresses are

, (2.9)

and the velocities become

with

, (2.10)

in which (u, v) _b are the horizontal components of the basal velocity vector evaluated in Equation (2.10)₁ as prescribed in Equation (2.8)₁ ,ᐁ denotes the two-dimensional horizontal gradient operator. C_s and C_g are positive scalars so that (u,v) at any point in an ice sheet is anti-parallel to ∇h _f. Flowlines of the horizontal velocity components must be the orthogonal trajectories of the lines of equal height of the free surface and this must be so for any point on a vertical line. Whenever observations indicate a different direction, the shallow-ice approximation cannot be the appropriate model. Obviously, horizontal strain-rate effects must come into play under these circumstances.

For a prescribed domain D with boundary and a given temperature field, the formulas Equations (2.9) and (2.10) are simple, and computations involve a quadrature only. For a power-law fluid

(2.11a)

and for the polynomial law in Equation (2.2),

(2.11b)

Let us pause for some remarks. It follows from Equation (2.9) that to lowest order a varying rate factor does not cause the glaciologists’ first commandment to be violated. This was never, apparently, an obvious question to the precursors of the shallow-ice approximation. Moreover, when the basal topography is less smooth to the extent that the flow within a boundary layer close to the base no longer follows exactly the direction of steepest descent of the free surface, or when typical lengths of topographic undulations are of the order of the ice thickness, then such an ad hoc motivation of the shallow-ice approximation must be replaced by the systematic treatment (Hutter, 1980, 1981, 1983) that allows derivation of the next-order corrections. In a plane-flow situation, I have done this and derived higher-order model equations for the time dependent surface elevation of an ice slope (Hutter, 1980). It was that calculation which showed that the perturbation solutions constructed with this systematic scheme breaks down when a power-law rheology is used. A finite-viscosity law offers remedy, and retaining the power law requires matched asymptotic expansions to construct the correct solution.

It is only recently that such higher corrections are accounted for in model calculations; however, systematic derivations are not generally attempted (Herterich, 1988; Dahl-Jensen, 1989a, b). This is certainly a promising avenue for future research.

The remaining real numerical problem is the determination of the temperature field: in D(x,y,z,t)

(2.12)

otherwise, in which horizontal thermal diffusion has been ignored. Thus, Equation (2.12) is parabolic in x and y, suggesting forward-marching procedures (Hutter and others (1986) have used this in a steady-state two-dimensional ice-sheet study). Computationally, however, diffusion stabilizes the numerical schemes so that it is advantageous to use a diffusive term in the expression k div grad T=k ΔT despite the expected smallness of the contribution .

There remains the derivation of the domain-evolution equation from a vertical integration of Equation (2.6)₁ :

(2.13)

q is the volume flux and can easily be determined from Equation (2.10)₁ :

(2.14)

where

. (2.15)

For a power-law fluid, this becomes

(2.16)

and for the polynomial law in Equation (2.2)

(2.17)

where each of the Q_jS is given by a corresponding formula Equation (2.16). Combining Equations (2.13)–(2.17) yields a nonlinear diffusion equation for ftf which, however, is in general singular at the margin. To see this singularity, substitute Equation (2.14) into Equation (2.13) to obtain

(2.18)

in which Δ₂ is the horizontal (two-dimensional) Laplacian operator and {•} is the curly bracketed term in Equation (2.14). As long as C_s and Q are bounded, which we shall now assume, the term {•} will vanish at h_f = h_b , making Equation (2.18) singular at the margin. For the no-slip condition, C_s = 0 and q = Q∇hf. The strength of the singularity follows from Equation (2.9), if we request the basal shear stresses to be bounded. In a local coordinate system in which χ is in the direction of steepest descent of the surface and y perpendicular to it, we have

. (2.19)

So with

(2.20)

boundedness of σxz(x _m) implies , with . A similar local analysis for the term in Equation (2.18) using Equations (2.14), (2.15) and (2.16) shows that

, if the basal sliding coefficient is bounded, C < ∞.

α=1, if C becomes singular as (h _f–h _b)^–1, as x→x _m.

In the shallow-ice approximation, the surface slope has a square-root singularity when the margin x = x_m is approached, unless the basal drag coefficient C ∞(h _f–h _b)^–1 is close to the margin. Morland and Johnson (1980) proposed such a sliding law and thus made their shallow-ice solutions uniformly valid but, of course, such a law is less clearly based on known physics (Fowler, 1990). The singularity (when C< ∞) is evidently not due to any rheological peculiarity (as, for example, the infinite viscosity at small stretchings of the power law) but to the use of the shallow approximation. It is a real physical effect that must be coped with.

Fowler (1990) showed for two-dimensional isothermal flow that in the vicinity of the margin the full Stokes-fiow problem must be solved; and at that margin, slopes are large but not infinite. This tells the numerical modeler, who usually smears over this singularity, that mesh sizes close to the margin should be selected sufficiently small, as otherwise the discretized solution will be inaccurate.

A different kind of singularity occurs at a divide. Fowler (1990) quoted Hindmarsh (1989) for proving that the reduced model does have infinite curvature of the surface at the divide if a power-law fluid is used with n > 1. The problem had already been recognized by Raymond (1983) and Hutter and others (1986) in a steady-state plane ice-sheet analysis. These solutions were constructed with the aid of a marching procedure column by column up to the margin. However, as stated by Szidarovszky and others (1989) “... 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 the longitudinal stretching effects, which we know to be significant in the vicinity of ice divides (Raymond, 1983). Szidarovszky and others (1989) constructed a solution including longitudinal stretching effects (and thus avoiding infinite curvature when Glen’s flow law is used), and Fowler (1990) sketched a scale analysis. The numerical modeler can learn from this that computations should preferably not be started at a divide and that constitutive models incorporating finite viscosity at zero stretching is advantageous. Furthermore, it is advisable to select small grid sizes close to ice divides.

I also ought to mention the limitations of the cold ice-sheet model that are based on pure physical grounds. The foremost limitation is the restriction that the ice sheet muit be cold. This restriction is enforced computationally by setting T≣T _M whenever Τ should rise above the melting point. As long as melting temperatures are only reached at basal points, the numerical implementation is consistent with the model. On theoretical grounds, however, finite regions of temperate ice near the ground are possible, and are likely to occur close to the margin where strain heating is large (Hutter and others, 1986; Hindmarsh and Hutter, 1988; Hindmarsh and others, 1989). In the latter situation, the enthalpy density of the ice is increased by increasing the water content. This may be redistributed by movement of the interstitial water. Then two-phase flow will arise; models for it will be described in sections 3 and 4.

A simpler remedy is to adopt the “enthalpy” method (Elliot and Ockendon, 1982) where the effect of the latent heat is represented by an increased specific heat over a narrow temperature range. The idea is to let the specific heat capacity tend to infinity as the temperature approaches the melting point. This is equivalent to letting the water drain from the ice without any adjustment to mass balance and provides an infinite sink for heat entering the system.

The formal specification of the enthalpy method is as follows. Let £ represent the internal energy of the ice. Then we choose the specific heat capacity in the form, for example,

(2.21)

where Τ — T_M is the homologous temperature and β and γ are constants. Hindmarsh and others (1989) performed computations for

(2.22)

Notice that C(T)→∞ when the melting temperature is approached. Computationally, the effect of Equation (2.21) is that the melting temperature is never reached; the model formally never breaks down and the latent heat that would be consumed by the melting processes is approximately taken into account.

3. Polythermal Ice Masses

3.1. General Remarks

The concepts on which the governing equations of polythermal or temperate ice masses are based are the classical mixture concepts; in fact, they can be used in glaciology in many more situations in addition to those mentioned here (Hutter and Engelhardt, 1988a, b). Basic to mixture theories is that the interpénétration of the various constituents is conceptually idealized by assuming that each point in a body is simultaneously occupied by all constituents, and that for each of these constituents balance relations of mass, momentum, energy and entropy hold. The only difference between this point of view and that of a one-component body is that constituent mass, momentum and energy are not conserved; only their sum must be conserved. As is the case with one-component continua, the respective equations do not suffice for the description of the field variables involved, and so some variables must be expressed in terms of others by materially dependent constitutive quantities, e.g. the partial stress tensors of each constituent can be postulated as a functional of the constituent deformation fields. Even in cases in which these equations can be written down explicitly, they are in general very complicated; so simplifications are needed.

To this end, mixture concepts were developed in which only some — the most important ones — of all balance laws are used. For instance, when water percolates through cold firn, one must assume that firn and water have different temperatures, implying that two energy balances are needed for the firn and the water separately; however, when water is present in temperate ice, the assumption of equal temperatures may be quite reasonable; a single energy balance for the mixture as a whole may be sufficient. Similarly, when studying impurities enclosed in ice, their momentum can be lumped into that of ice. It is evident this principle leads to a hierarchy of mixture models each having its own complexity and being applicable under restricted physical situations only. The most important models in glaciological applications have the following structure (Müller, 1985):

• Class I. The balance laws of mass of all constituents are used, but momentum and energy-balance statements are only formulated for the mixture as a whole. Often energy considerations are ignored.

  These models are typical for the description of the diffusive motion of any particulate substance contained as a contaminant or impurity in another substance. Equations have advective and diffusive structure.

• Class II. The interpenetrating constituents have comparable mass fractions and their velocities are sufficiently distinct from each other. Here, the balance laws of mass and momenta are formulated for all the constituents, but only an energy balance for the mixture as a whole is used. In other words, the various constituents have sufficiently distinct momenta or velocities but a single temperature suffices for all. Often, again, energy considerations may be unimportant.

  These models go beyond the classical diffusive models (which can only describe dilution and not concentration). The interaction forces between the constituents are important and make the separation of the constituents possible. In ground-water flow, Darcy’s law expresses this momentum exchange mechanism: the force exerted by the water on to the soil and, alternatively, that exerted by the soil on the water.

• Class III. The next level is then the full-mixture thermodynamics with all constituent balance laws of mass, momenta and energies.

  The creep deformation of cold firn under the influence of percolating surface meltwater would have to be formulated by a theory of this complexity.

We will describe temperate ice in a polythermal glacier as belonging to the first class of mixture models. The motion of the water through the ice is then basically a diffusive process. In wholly temperate ice masses, the concept is one of the second class. Water moves through ice much like ground water through soil, and so a law similar to Darcy’s law must be derived to describe its motion. Detailed descriptions have been given in earlier publications (Fowler and Larson, 1978; Hutter, 1982; Fowler, 1984; Hutter and Engelhardt, 1988a, b; Hutter and others, 1988; Blatter, 1991; Blatter and Hutter, 1991), and so we only outline the highlights here.

3.2. Polythermal Ice Masses

Polythermal ice is assumed to consist of several disjoint regions of ice masses, each exhibiting physically distinct behavior. Some parts are cold, others are temperate. Together, the disjoint sets form the polythermal ice mass. We shall define the ice as cold at a particular point, if the temperature of the particle at that point is below the melting temperature, defined by Equation (2.5) and if there is no moisture present. The definition of the term “temperate” is difficult, in particular when salts are present. We shall ignore them because impurities that are soluble in water do not, in general, arise in greater proportion (exceptions are marine ice sheets and glaciers), and may then define ice as temperate whenever the local temperature reaches the melting point, which in this case depends on pressure alone. Under this restricted assumption, cold and temperate ice regions are separated by a surface, the cold-temperate transition surface across which the water content may jump from zero to a finite value. Observations provide support for such a simplified view point (Blatter, 1991; Blatter and Hutter, 1991).

For many polythermal glaciers and ice masses, the cold parts are the regions close to the free surface and the temperate regions touch the bed. Their existence is then due to both sufficient geothermal heat and strain heating (i.e. dissipative heat due to viscous flow). The amount of water that is generated this way is likely to be small so that crevasses and other openings are free of water, but the ice is “wet”. This means that the water collects in and moves along the grain boundaries rather than in larger cracks, and it may also be trapped in bubbly inclusions. Thus its momentum will not considerably differ from that of ice: in short, a diffusive model is appropriate.

A mathematical description of polythermal ice therefore consists of the field equations in the cold and temperate sub-regions, respectively, and boundary conditions on the free surface, the base and the cold temperate transition surface.

Field Equations

Field equations for cold ice and boundary conditions along those parts of the free surface and the bed where the ice is cold have been given in section 2, Equation (2.1)–(2.5). They will not be repeated.

Temperate ice is postulated to be a binary mixture consisting of ice and water, in which the concentration (by mass) of water is very much smaller than that of ice (perhaps approximately 1–3%). Two balance laws of mass for the constituents ice and water (or equivalently, and used here, for water and the ice-water mixture) and only one balance law of momentum for the mixture comprises the adequate mixture model (of Class I). If we further ignore volume changes under phase changes, an assumption that is very well satisfied at such small concentrations, the incompressibility assumption can also be imposed. This suggests that the following field equations and definitions hold: Mass balance for mixture: div v = 0, Mass balance for moisture content w: , Momentum balance for mixture:

, (3.1)

Stress-stretching relation: , Fick’s law of moisture flux: , Moisture production: Barycentric velocity: .

Here, w=ρW/ρ is the moisture content; the ratio between the water density per unit volume of mixture and the mixture density, ρ = pw + p1. Furthermore, L is the latent heat of fusion and v is a small moisture diffusivity. A is a moisture-dependent rate factor that is only poorly known. According to Lliboutry (1976), it increases with increasing w, but its dependence on w is not as strong as that of A(T) on temperature. As of today, because of lack of better knowledge, one still assumes A is independent of w.

The moisture flux can be shown to have the form

, (3.2)

but this expression is only needed for transformation purposes, since we have postulated a constitutive relation for j (Fick’s law as in Equation (3.1)₅ ). The moisture production M as given in Equation (3.1)₆ is essentially a constitutive relationship that has been prescribed in that form because the energy balance has been ignored. If the energy balance for ice plus water is written down*

, (3.3)

then Equation (3.1)₆ would follow from Equation (3.1)₂ and Equation (3.3) by ignoring moisture flux j and heat flux q. Fowler and Larson (1978) used Equations (3.3) and (3.1)₂ , and requested moisture flux and heat flux in temperate ice to vanish and then wrote Equation (3.1)₆ . Hutter (1982) ignored Equation (3.3), included a non trivial moisture flux and postulated Equation (3.1)₆ . A model not involving a closure condition for moisture production, Equation (3.1)₆ , but instead proposing the energy law as in Equation (3.3) is free from such ad hoc assumptions; it has not been proposed before. Ensuing developments are independent of this choice, but I might mention that having a non-vanishing moisture flux is crucial. I shall shortly return to this point. In what follows, Hutter’s (1982) model is used.

If relations in Equation (3.1)_5,6 are substituted in Equation (3.1)₂ , the mass-balance equation for the moisture content reads

, (3.4)

an equation which is formally the same as the energy Equation (2.1)₃ for cold ice. It is of the advection diffusion reaction type. In particular, because of the small (but non vanishing) diffusivity ν, it is parabolic in time and space but becomes hyperbolic if diffusion is ignored. This difference is important when boundary conditions are formulated. The diffusion equation requires boundary conditions along the entire closed boundary of the domain where the ice is temperate. The hyperbolic equation allows them only along the upstream positions. This has far-reaching consequences. We proceed for the moment with the diffusive case.

Boundary Conditions

We assume that the free boundary is entirely cold; its boundary conditions are then given by Equation (2.3). (The more general case has been discussed in detail by Blatter (1991).) Alternatively, the basal boundary is assumed to consist of a cold part and a temperate part . Along the former, the boundary conditions are stated in Equation (2.4); they include the case for which the base just reaches the melting point.

Along the temperate bed , above which a nonempty domain with temperate ice exists, the basal surface melting rate cannot be neglected in the formulation of the moisture-flux boundary condition. If V₁ is the ice velocity, then along a non-deforming bed one has

(3.5)

where V₁ is the ice velocity. Equation (3.5) says that all the water is instantly drawn to the bed, an assumption that is likely to be realistic, given the extremely small amounts of meltwater at the base. With the definition of the barycentric velocity and with the aid of Equation (3.2), it is then readily shown that Equation (3.5) has the alternative form

. (3.6)

This shows that the barycentric velocity deviates from being tangential to the base by two contributions, the melting rate and a weighted moisture flow through the base. The boundary condition describing moisture flux at the base follows from the jump conditions associated with the balance laws in Equation (3.1). That of the moisture content is

(3.7)

where and ƒ⁺ is the value of ƒ immediately below the base, whereas f⁻ is the corresponding value of ƒ just above it. Thus, Equation (3.7) may be written as (we omit the negative signs on the ice side of the basal surface)

. (3.8)

Use was also made of Equation (3.6), and the basal drainage function has been introduced. It is the occurrence of this function that makes the theory of polythermal glaciers practically difficult, as it must be known if a boundary-value problem is to be solved. One could instead prescribe the moisture content at the base, but this seems to be as difficult as prescribing .

The moisture-flux boundary condition Equation (3.8) contains the basal melting rate for which a governing equation must be established. This equation is the energy-jump condition corresponding to the global balance law of energy for ice and water*

. (3.9)

Here, Equation (3.6) has been used, and in the last line the contribution from the kinetic energy has been ignored. With

(3.10)

Equation (3.9) finally assumes the form

. (3.11)

For cold ice which reaches the melting temperature at the base, w = 0 and j = 0, for whichEquation (3.11) reduces to Equation (2.4)₆ . Alternatively, when the ice above the base is temperate, the conductive heat may be ignored (q = 0). In any case, Equations Equation (3.8) and Equation (3.11) together define the moisture-flux boundary condition along the temperate base. For non-vanishing w, the two conditions can be combined and then yield the equation

, (3.12)

but this equation is only meaningful when the ice above the temperate base is temperate.

The remaining boundary condition is the sliding law

(3.13)

as prescribed already in Equation (2.4)₂ for the special case that v.n = 0. This sliding law identically satisfies the tangency condition.

Cold-Temperate Transition Surface (CTS)

Contrary to intuition, the transition from cold to temperate ice is not smooth, but takes place at a singular surface where the moisture content and the temperature gradient may suffer a finite jump. For this discontinuity to exist, the salt and impurity content of the ice must be sufficiently small. Even if this were not the case, the transition layer is so thin that it could be treated as a discontinuity anyway. Support for a smooth transition would also come from thermostatics of isolated lens-type water inclusions in an ice matrix. Because of surface energy and surface tension, such water inclusions can exist at temperatures down to —10°C; however, only for extremely small lenses (Hobbs, 1974; Lliboutry, 1987b, 1993; Alts and Hutter, 1988a, b,c, 1989). Most water collects in veins, whose local mean curvatures are smaller than lenses. Thus, surface-tension induced supercooling of the water is correspondingly less for veins than for lenses. So, the assumption of an abrupt change of the temperature gradient at the transition from cold to temperate ice is probably reasonable. Corroboration to this effect is also provided by some measured temperature profiles in boreholes from Greenland, which suggest a discontinuity in the vertical temperature gradient at the suspected location of the CTS (Stauffer and Oeschger, 1979).

In classical continuum thermodynamics, a phase change surface is defined as a surface at which temperature and tangential velocity are continuous (Hutter, 1983; Müller, 1985). We now adopt and extend this definition and define the CTS as a singular surface at which the (melting) temperature and tangential barycentric velocity are continuous:

. (3.14)

Furthermore, the jump conditions of moisture content, momentum, energy and entropy must hold (Hutter, 1982, 1983):

(3.15)

In these equations, is the melting-freezing rate (> 0 for melting), and we have taken the ()⁻ sign to signify the cold side. According to Equation (3.15)₁ , the jump in moisture flow is related to the jump in moisture content, and the amount of melting of ice at the CTS. Equation (3.15)₂ is only approximate because it ignores the contribution which is very small. Equation (3.15)₃ is the entropy balance when has been used (see Hutter, 1982; Hutter and others, 1988; Blatter, 1991). The energy balance is not needed, because it reduces to the same statement as the entropy balance. Indeed, from

(3.16)

this claim is immediate. Combining Equations (3.15)₁ , and (3.15)₃ , finally yields

; (3.17)

therefore, the kink in the temperature profile is related to that in the moisture-content profile. Ignoring the moisture flux ab initio (j = 0) makes the temperature gradient continuous at the CTS, and, via Equations Equation (3.15)_1,3 , implies there. (We must exclude the other possibility , as it would correspond to the thermodynamic equilibrium.) Since w⁻ = 0, this would imply that the moisture content would also have to vanish on the temperate side.

On the other hand, by ignoring a Clausius—Clapeyron adjustment of the melting temperature, Equations (3.15)_1,3 imply that at a melting CTS the moisture content must vanish on either side, , implying that also dw⁺/dy = 0 (on the temperate side of the CTS.) and dT/dy = 0 (on the cold side of the CTS)*. For a freezing CTS, a non-vanishing moisture jump is possible. This implies that the build-up of temperate ice zones is very smooth and therefore probably very hard, while the freezing of temperate patches is accompanied with moisture jumps at the CTS.

On Reduced Models

With Equation (3.17), we have reached the point where discussion of reduced models is becoming meaningful. Before we begin, we note that reduced models are needed, because a full diffusive model forces field glaciologists to determine the basal drainage function, or the basal moisture content, which is virtually impossible. In Fowler and Larson’s (1978) model, diffusion is ignored ab initio, i.e. V = 0 implying j = 0, and so the temperature profile has a continuous derivative throughout, and the moisture content is continuous at the CTS. This prevents a simple Stefan problem ↑ from being solvable. Hutter (1982) has shown that the motion of the CTS through a block of ice initially consisting of two parts, one being temperate with known uniform moisture content and the other being cold at Τ < T_M , cannot be determined from such a reduced model. In spite of this, observations indicate the temperature profile to have a kink at the CTS, making the Fowler-Larson model inappropriate. Another restriction of the Fowler-Larson theory (pointed out already by Fowler and Larson themselves) follows from the hyperbolicity of the moisture-evolution Equation (3.4), which, with vanishing moisture diffusivity, takes the form

. (3.18)

Its characteristics are precisely the streamlines. According to Equation (3.18), the moisture content will grow along any streamline, and so the CTS can never cross a streamline twice, because the moisture content would have to vanish in this theory at the intersection points, which is impossible if the moisture content has to grow from one intersection point to the other. Configurations of CTSs as shown in Figure 1 are therefore not possible according to this theory.

Fig. 1. A section through an ice sheet with three possible cold-temperate transition surfaces and streamlines that cross these surfaces more than once. Case 3 is thought to arise close to the margins of a glacier or ice sheet ( schematic).

A different limit theory, which combines the practical advantages of not having to prescribe a basal boundary condition on moisture flux or content, yet saving the property of allowing a jump in moisture content at the CTS is obtained if the diffusive model is taken at the limit as the diffusivity becomes infinitely small everywhere except in a very thin boundary layer close to the CTS. In such a distinguished limit, the moisture mass-balance equation agrees with Equation (3.18), is hyperbolic and thus only needs prescribed boundary conditions at the upstream boundaries. Along the CTS, where w must vanish on the cold side, but can have a finite-positive value on the temperate side, Equation (3.15) survive to their full length. Here, an expression for the jump must be found.

To this end, consider Figure 2, which shows a close-up of the boundary layer at the CTS. Accordingly, and the asymptotic value of w at the outer edge of the boundary layer is w _∞. In a model which ignores the boundary layer, it will be the quantity that is recognized as the moisture jump across the CTS. An approximate solution of the moisture balance equation motivates the form , where a is a phenomenological parameter. If this is so, the moisture jump in Equation (3.15) must be replaced by , in order to account properly for the local jump relations.

Fig. 2. Boundary-layer profile of the moisture content w close to the CTS.

To derive the approximate form of the mass balance for the moisture content in Equations Equation (3.1)_2,4 in the boundary layer close to the CTS, we conjecture that the local time derivatives of w, the advection of moisture parallel to the CTS, the moisture diffusion parallel to the CTS, the moisture production can all be ignored in comparison to moisture diffusion and moisture convection perpendicular to the CTS. Hence,

, (3.19)

in which y is the coordinate perpendicular to the CTS and is the barycentric velocity component in the same direction into the cold ice, which is assumed not to vary across the boundary layer. With the boundary conditions

, (3.20)

the boundary-layer solution on the temperate side of the CTS, see Figure 2, takes the form

. (3.21)

Furthermore,

, (3.22)

which may be non-zero also when v→0.

This calculation motivates the choice , where a is a phenomenological parameter that may be estimated from simple model calculations of the full problem. Blatter (1991) and Blatter and Hutter (1991) used a = 2 on the basis of a similar boundary-layer computation as the one presented above (but an unrealistic boundary condition replacing Equation (3.20)₂ ). a > 1 means that Fickian diffusion will drive moisture back to the CTS while a < 1 corresponds to the case driving it away from it. Both cases seem to be possible, suggesting that a ought to depend on other field quantities defined on the CTS. I will give some further remarks on this below.

This theory is now used as follows: in all equations the diffusivity ν is set to zero, and the moisture jump arising in Equation (3.15) is replaced by , e.g.

(3.15a)

In this way, the reduced theory computes the far field while properly accounting for the CTS boundary layer.

This model therefore bears the advantage of not having to prescribe a boundary condition on moisture flux or moisture content at the temperate basal surface, yet to preserve the jump properties at the CTS. In this theory, the CTS may cross any streamline as many times as demanded by the remaining equations; and, in particular, isolated temperate patches may exist in a cold environment and vice versa.

The crucial part of this distinguished limit model is the selection of the coefficient a. It certainly depends on the dynamic state of the CTS and may thus be postulated in the form

(3.23)

Further dependencies on are, in principle, also possible; however, the dynamics of these variables are described by the Equation (3.15) and a further inclusion in Equation (3.23) is not felt necessary. We also anticipate that under freezing conditions water will be driven towards the CTS, whereas under melting conditions this diffusion is away from the CTS; thus, , where ƒ is now a new functional of the form of (Equation 3.32). The simplest expression would be

, (3.24)

with appropriately selected coefficients m and p. Whether Equation (3.24) is reasonable must follow from a comparison of solutions of the full boundary value problem v∞0) with solutions of corresponding boundary value problems using the limit theory (v→0) with an explicit choice of a like Equation (3.24).

Model Computations

The full model equations of polythermal ice masses that are outlined above have never been solved numerically. The construction of solutions is bound to be difficult because the evolution of two moving surfaces must be determined along with the velocity and temperature fields in the cold-ice region, and the velocity and moisture content fields in the temperate-ice region. Approximations involve a scale analysis and will essentially lead to equations appropriate for the shallow-ice approximation. Such calculations for the case v∞0 were done by Hutter and others (1988), Blatter and Hutter (1991), and Blatter (1991). The reduced model equations were derived for two-dimensional flow and then solved for a strictly parallel-sided slab consisting of a cold layer on top of a temperate bottom layer. Later computations for a plane flow model of Laika Glacier, Canada, involved fixed domain mappings of the governing equations prior to their finite-difference representation. With the discretized equations, a CTS was successfully constructed and provento exist in Laika Glacier as a remnant of the Little Ice Age (Blatter and Hutter, 1991).

Physically, the computations have shown that the proposed model provides promising indications to predict adequately polythermal structures in ice sheets and glaciers. In fact, the Laika Glacier results give reason to assume that its polythermal structure is the manifestation of transient processes. Numerically, it has been demonstrated that the proposed equations can be solved under idealized two-dimensional but nevertheless realistic flow situations. I am optimistic that, before long, ice-sheet dynamics will be analysed through the ice ages on the basis of a polythermal model. Conceptually, the long and rather painful development of the final model equations with its distinguished limit v→0 is a demonstration of the fruitful interplay between practical possibilities and theoretical desires. Field glaciologists feel uncomfortable with the determination of the moisture diffusivity and prescription of a basal boundary condition on moisture flux, and theoreticians cannot do without a moisture jump at the CTS. The reduced model equations satisfy both in an ideal fashion. Its test is urgent.

4. Temperate Glaciers

We now consider ice masses that are wholly temperate through all or most of the seasonal cycle and may only have a cold surface layer in the firn zone during winter time (which will conceptually be ignored here). Summer surface melt rates are large in these cases; surface water is collected in surface channels (and also percolates through the firn) and quickly reaches the deeper parts of the glacier through moulins and crevasses. In a typical warm summer day such a glacier is saturated with water in its deeper part, and comprises in its upper part of ice and metamorphosed snow through the pore volume of which the water, fed from the surface, flows in an unsaturated fashion (Fig. 3). This source of water keeps the deeper, saturated, part of the glaciers alive. The domains are separated by the phreatic surface, which is the free surface of the “ground-water fluid” in the saturated part at depth. Since the surface melt rate follows the daily radiative input, the phreatic surface performs large oscillations with a dominant period of 24 h, and may at night occasionally reach the base. The large flucutations in interstitial (water) pressure cause corresponding daily variations in ice velocity through the dependence of the sliding law on water pressure. In short, the system “glacier and its water” directly responds to the availability of surface water; corresponding time-scales are days, perhaps weeks. Of course, this is all common glaciological knowledge, but its mathematical modeling (Fowler, 1984; Hutter and Engelhardt, 1988a, b) has so far largely been ignored.

Fig. 3. Sketch of a wholly temperate glacier with the ice saturated by water in its deeper part, its unsaturated ice at the upper part and the phreatic surface separating the two.

The reason is that, commonly, the flow of water through glaciers is treated as a pipe flow through a system of intraglacial channels (Röthlisberger, 1972; Shreve, 1972; Weertman, 1972; and many others) with basal sliding mechanisms in parts based upon linked cavity configurations of the basal water conduit system (Walder, 1986; Kamb, 1987). These latter situations are certainly another possible approach to the physics of glaciers in the presence of ample water, and they are more appropriate in explaining glacier-surge mechanisms. However, they are not ideally suited to cope with a true interaction of ice and water. I take the view that it is permissible to distribute all the intraglacial channels, moulins and cracks smoothly over the space that is occupied by the ice and to model their dynamic effects by a Darcy-type permeability. This implies that the differential lengths of such a theory are at least of the order of a mean distance of the cracks.

Temperate ice will therefore be assumed to be a porous continuum consisting of ice and interstitial water (saturated), or ice, water and air (unsaturated). The continuity assumption of this porous medium may somewhat stretch the circumstances arising in reality; it is presently the only compromise by which the dynamics of temperate glaciers may become mathematically tractable.

Because the amount of water moving through the ice is relatively large and the velocity of the interstitial water exceeds that of the ice, considerations of balance of momentum for the water are important. We treat the saturated and unsaturated regions separately and discuss boundary and transition conditions last. The model I shall present is very close to earlier formulations of the description of the deformation of natural slopes of soil (Vulliet and Hutter, 1988). In fact, what is new in the temperate-ice situation is the occurrence of melting within the ice mass.

4.1. Temperate Ice Saturated With Water

We assume a binary mixture concept of Class II for ice and water, ignore the role played by the salts* and deal with a single temperature for the constituents water and ice. The pore space along grain boundaries, cracks and the larger intraglacier openings, such as moulins and crevasses, is lumped into the single variable “porosity”, n. The bulk (true) densities, , assumed to be constant (incompressibility) and the constituent densities in the mixture, are then related to each other by

. (4.1)

Alternative variables, common in the soil mechanics literature, are the mass fraction, or concentration of water (moisture content), w and the void ratio, e= water volume per solid volume, the transformation from one to the other being defined as indicated in Table 1.. Balances of mass and momenta read

(4.2)

in which M ^α is the constituent mass production and m^a is the constituent momentum production (= specific interaction force), which satisfy the conservation constraints

. (4.3)

A scale analysis can be performed with Equations Equation (4.2) which indicates that inertial terms can be ignored. This assumption is certainly well satisfied for the constituent ice but may, however, be questionable under certain circumstances for the water. This is particularly so, when the intraglacial water flows mainly through a few conduits; but, in that case, the continuity assumption is questionable anyway; I will ignore this case here mainly because I have no clue myself how to treat it.

Table 1. Transformations between porosity n, void ratio e and mosture content w

Because of the incompressibility assumption constitutive relations cannot be postulated for β ^I and β ^W, but only for their differences with an interstitial pressure pw, known as pore pressure. Distributing pw among the constituents according to their volume fraction, one obtains for the water and ice stresses

; (4.4)

this shows that their sum is made up of the pore pressure and the constitutive parts of the constituent stresses, t^w and t^I (this sum is not identical to the mixture stress, but it agrees with it, if the diffusive momentum fluxes are ignored (see Hutter and Engelhardt, 1988a, b). It can also be easily demonstrated that the neglect of the momenta in the momentum-balance laws implies that these diffusive momentum fluxes can be ignored; so the two assumptions are consistent with each other.

With Equations (4.1), (4.3), (4.4) and upon ignoring the time rates of change of momenta, the balance laws of mass and momenta take the form

(4.5)

These will be regarded as field equations for n, pw, v _I and v _w, respectively. Thus, materially dependent statements must be established for the stresses t^I and t^w, the volumetric melt rate M and the interaction force ρβ, where ρ is the density of the mixture.

As for the stresses, we assume the pore fluid to be inviscid, implying t^w = 0. Alternatively, the partial stress tensor t^I is postulated as

(4.6)

Accordingly, the materially dependent stress tensor t^I is decomposed into an isotropic pressure term p ^Il and a deviatoric contribution . This deviatoric contribution is related in a non-linear fashion to the deviatoric part of the stretching tensor . A is a moisture-dependent rate factor and ƒ is the creep-response function, assumed to depend upon the second invariant of . Since the stretching tensor D^I need not be traceless, despite the incompressibility assumption, we have also postulated a (non-linear) viscous relationship between and p ^I; the pressure-dependent parameter ζ is a bulk viscosity.

Due to lack of detailed knowledge, one must still postulate the stress-stretching relation in the form known for polythermal ice. Thus A will be assumed constant (independent of ω) and ƒ may be prescribed in the form of Equation (2.2), while the bulk viscosity ζ is assumed to vanish identically. This corresponds to the Stokes hypothesis and owing to Equation (4.6) necessarily requires p ^I = 0, or div ν^I = 0.

Establishing a constitutive relation for the interaction force amounts to formulating Darcy’s law. To this end, we introduce the so-called seepage velocity defined by

. (4.7)

This is a materially objective vector quantity, i.e. v^F transforms as a vector under (non-inertial) observer transformations, even though v^w and v^I do not. In Darcy’s law, this seepage velocity is related to the interaction force pß; however, when mass production of water does not vanish (M≠0), ρβ does not transform as an objective vector under such changes of observer frames. It is the combination

(4.8)

that is objective. We therefore generalize Darcy’s law by postulating that

(4.9)

with a constant of proportionality μ/k*.μ is a viscosity and k* is the so-called coefficient of absolute permeability (m²). Combination of Equations (4.8) and (4.9) yields the generalized Darcy law

. (4.10)

It involves a contribution of the volumetric melting rate. In earlier formulations (Fowler, 1984), this contribution is missing.

In the soil mechanics literature, one does not work with Equation (4.5) but rather a set of equations that can be derived from them. First, by adding the two force-balance statements in Equation (4.5)_3,4 one obtains a force balance for the mixture as a whole,

. (4.11)

Here, use has been made that t^w = 0 and. Alternatively, Equation (4.5)₄ (with Equation (4.10) substituted for the interaction force) can be solved for the seepage velocity

. (4.12)

In ensuing developments, the term involving grad η will be ignored; in other words, variations in porosity are assumed to be an order of magnitude smaller than those of the pore pressure. We also introduce the gravitational potential ϕ and the piezometric head h, according to

(4.13)

and may then write Equation (4.12) in the form

. (4.14)

This is the reduced momentum-balance law for the constituent water closest in form to Darcy’s law. In fact, k is Darcy’s permeability coefficient with dimension (m s⁻¹) and Equation (4.14) is identical to Darcy’s law when M = 0.

The field equations describing the motion of water and ice in the saturated region of a temperate glacier are now given by

(4.15)

Here, all simplifying assumptions have been incorporated, and it is understood that the expressions for v^F, w and h are substituted. The functions A, ƒ and ζ are assumed to be explicitly known. With this understanding, the set in Equation (4.15) forms 14 independent equations, the counting being indicated in parentheses on the right of each line in Equations Equation (4.15). Unknown field variables are n (1), p_w (1), v¹ (3), v^w (3), (5), p ^I (1) and Μ (1), a total of 15 fields, and so one additional statement is needed (to describe the evolution of the melting rate M). This will be the energy equation, derived below.

Two additional simplifications might to advantage be imposed in a first computational attempt to solve the system Equation (4.15):

• (1) If the Stokes hypothesis is made, then p ^I = 0, so the viscous ice pressure drops out of Equation (4.15), reducing the number of unknown fields to 13. Equation (4.15)₆ must then be replaced by

  

  , (4.16)
  , an equation which guarantees the tracelessness of .

• (2) It is usually justified to assume that . In this case the seepage velocity is approximately given by v^F = nv^w, and the generalized Darcy law may be written as

  

  . (4.17)

  Substituting this expression into Equation (4.15)₁ yields

  

  , (4.18)

  which is a diffusion equation with a diffusivity that depends on the melting rate M. With M = 0, this equation is also known in soil mechanics.

Let us now turn to the energy equation for the mixture as a whole; it may be stated as

, (4.19)

in which the “material” derivative d(·)/dt is defined as the time rate of change for an observer following the barycentric velocity

. (4.20)

∈, q and Φ are the internal energy, heat flux and dissipation function, all of the mixture as a whole. We shall ignore the conductive heat flow (q≃0) because, owing to the Clausius-Clapeyron equation, the melting temperature varies very weakly with pressure: so cannot assume appreciable values. If we also ignore the kinetic energy of the diffusive motion, one has

(4.21)

Here, we have set T ₁=T _w=T _M, and have assumed that being the latent heat of fusion and C₁ the specific heat of ice at constant volume.

Next, consider the dissipation function Φ = tr(tD) where t is the mixture stress and D the stretching of the barycentric velocity v. With (in which the diffusive momentum flux has been ignored), and with Equation (4.20), one may readily show that

. (4.22)

This expression could be put into a different form by substituting Equation (4.15)_5,6 ; however, no further physical insight is gained thereby. According to Equation (4.22), the dissipation function consists of three terms: (i) the work by the pressure done on the dilatational part of the barycentric motion, (ii) that of the ice stresses on the stretching of the water motion, and (iii) the work done by the ice stresses on the stretching of the ice motion. It is not clear whether any of these contributions is negligibly small but it is likely that the first contribution is an order of magnitude smaller than the other two. We also expect the second term to be smaller than the third. Pilot calculations will have to clear this point.

With Equations (4.21) and (4.22), the energy equation takes its final form

. (4.23)

This equation is the missing link and serves as “closure condition” for Equation (4.15).

4.2. Temperate, Non-Saturated Ice Region

Above the phreatic surface there is insufficient water available that the pore volume would be completely filled. It is not clear to me how this region is best treated. So, I give a first, simplistic approach.

I shall regard the ice-water-(air) mixture as a binary mixture of ice and water in which both constituents are treated in a dynamic fashion, but acceleration terms are ignored. This yields

. (4.24)

Both ρ _w and ρ _I are unknown, but ρ _I can be related to the porosity, which is also unknown. Constitutive relations must be established for the stresses, the interaction force and the melting rate. These relations must take into account that percolation through a non-saturated matrix is associated with some notion of turbulence. So, we should require that , will quadratically depend on v^w — v^I, with an additional dependence on ρ _I and/or n. Thus,

, (4.25)

where C is a dimensionless drag coefficient, which may carry further dependencies.

For the stresses it may be convenient to introduce a pressure p, and to suppose that the water phase cannot sustain shear stresses, but ice may support a deviatoric viscous component. So, we assume

(4.26)

where, at a first approximation, A is assumed constant. Moreover, the constitutive part of σ ¹ is automatically deviatoric; this corresponds to neglecting the isotropic viscous stress (Stokes hypothesis) and requires the velocity field of the ice to be solenoidal (see discussion in connection with formula Equation (4.16)). With the representations in Equations (4.25) and (4.26), the field Equation (4.24)take the form

(4.27)

and constitute 15 equations for the 16 fields . The equation that is still missing is, of course, the balance of energy,

, (4.28)

in which heat conduction is ignored as before. A derivation along the lines of the previous deductions yields

. (4.29)

This does not contain any contribution from the interaction force which was chosen in the form of Equation (4.25) to account for the action of turbulence. This fact appears to be a weakness of the present model; viscous stresses in the constituent water may have to be incorporated to have turbulence in the water contribute to the melting rate.

4.3 Boundary and Transition Conditions

These must be formulated at the base, the free surface and the phreatic surface, respectively. We shall only outline the highlights and refer for details to Hutter and Engelhardt (1988a, b).

Basal Surface

The bed will be assumed to consist of hard, solid rock that may be partly fractured; however, the flow of water from the saturated ice into the solid rock will be assumed slow on the time-scales of daily fluctuations of surface melt rates. Such an assumption is, perhaps, realistic for alpine glaciers on a rock bed. It permits us to request the basal water discharge to vanish.

On this basis, the boundary conditions that one wishes to impose at the base are as follows:

• (1) Tangency of the water velocity to the rock bed,

  

  . (4.30a)

• (2) Tangency of the ice velocity to the rock bed, V_I.n = 0. However, this condition is only approximate, as some ice will melt at the base, and the corresponding amount of water production is so small in comparison to the surface-melt production that it can safely be ignored.

  As a consequence, the sliding law is chosen in conformity with this tangency condition

  

  (4.30b)

  in which C(·) is a drag coefficient and the dot signifies unspecified further dependencies. For instance, it is likely that C depends on ||τ^*|| and the pore pressure p ^w at the base. In fact, the determination of this latter dependence is one of Lliboutry’s major achievements (Lliboutry, 1968, 1975, 1979b). We note, finally that Equation (4-30b) automatically satisfies the tangency condition v_I · n = 0.

• (3) A complete set of boundary conditions would also necessitate the formulation of an energetic boundary condition involving geothermal heat and basal melt rate; however, on time-scales of decades and less the latter is dynamically of no relevance and so one can dispense with it.

It should be mentioned that, in the above, we impose a boundary condition on basal discharge, a condition which field glaciologists are not comfortable with. In a way, the situation is less critical here than it is for polythermal ice, because the abundance of water flow from the surface into the ice makes it possible to ignore a possible loss of water to the rock bed, because this loss is likely to be felt only on time-scales that are large in comparison to those dictated by surface-melt rates. Thus, field glaciologists feel probably more comfortable here with the imposition of the boundary condition of vanishing basal water discharge than with a corresponding moisture-boundary condition at a temperate base of polythermal ice. This is fortunate, because formulation of some sort of boundary condition on the water mass loss is mathematically necessary. On the other hand, a proper treatment would have to incorporate the water-saturated basal rock or soil beneath the ice and analyse the seepage flow through this porous matrix along with that of the ice. It would requin estimates of the coefficient of permeability of the bed, and a possible surface of impermeability would have to be imposed farther down. Through such an extension of the model, the uncertainties at the base may be “shifted farther away” from the ice and possibly be less perturbing but, obviously, the model will become even more complicated.

Free Surface

Let us consider next the free surface of the temperate ice unsaturated with water. With z=h _f(x,y,t) and

.

its kinematic equation is

; (4.31)

however, in view of the very short times (days, weeks, seldom months) on which these processes take place, one may ignore the geometric variation of the top surface and simply assume that z = hf(x,y,t) is evaluated for a fixed time. Nevertheless, (or a) is the surface net accumulation (melt) rate. To find an interpretation of , consider the jump condition of mass for the mixture as a whole at the free surface,

(4.32)

in which ()⁺ denotes the atmospheric side and ()⁻ the ice side, and ν is the barycentric velocity, ; then is the volume flux of ice from the atmosphere to the ice mass and is the corresponding mass flow. Physically, this flow consists of the accumulation of snow and the loss of mass at the surface due to evaporation and/ or sublimation as well as run-off of surface melt water. It must, on the ice side, be equal to the snow mass deposited and the mass of melted ice at the surface minus the surface water infiltrating the unsaturated ice mass. This interpretation also follows from Equation (4.32); indeed, by substituting the expression of the barycentric velocity into Equation (4.32) one obtains

(4.33)

where

(1) = accumulation-evaporation-sublimation-surface run-off

(2) = snow addition-surface melting

(3) = infiltration.

Here, is the accumulation function on the ice side, which is the sum over the added ice volume by solid precipitation and the melting rate. Three special cases are immediate from Equation (4.33):

• (1) When no melting takes place and therefore no water is present, then ρ_w=0 and

  

  . (4.34a)

  The accumulation-rate function which arises in the kinematic Equation (4.31), is larger than by the factor (1 — n)~⁻¹ , which would be expected.

• (2) When there is no net accumulation, , then Equation (4.33) implies

  

  , (4.34b)

  i.e. the infiltration rate is prescribed by the melting rate; in particular , as it must be.

• (3) When consists only of surface run-off, then Equation (4.33) says

  

  (4.34C)

  which is also obvious.

From an observational point of view, two of the three terms in Equation (4.33) must be measured; this is the fact that makes the application of this theory difficult and somewhat questionable. It is very difficult if not impossible to separate the total amount of surface water from that which runs off superficially.

The boundary condition of momentum is formulated as the momentum-jump condition of the mixture as a whole. Ignoring the jump of the convective momentum flux in comparison to the jump in the traction vector, this condition reads

. (4.35)

(A second condition is not needed because integration of Equation (4.27)₄ yields an expression for which all variables at the surface are already fixed.)

In the case when there is no accumulation and only ablation is taking place, Hutter and Engelhardt (1988a) derived from the entropy-jump condition a relation between the jump of the normal heat flow and the melting rate Oj.. We shall not go into any details here, because prescribing the atmospheric heat flow is conceptually just the same as prescribing The condition does not free us from a separate prescription of the surface run-off.

Phreatic Surface

Consider next the phreatic surfacez = h_ph(x,y,t), whose evolution equation may be written as

, (4.36)

in which

We shall adopt the convention that the saturated ice is on the (-) side of the phreatic surface and the non-saturated ice is on its (+) side; is then the accretion-rate function and represents the nourishment of the “ground-water table” from above. As long as no water percolates through the non-saturated upper ice region, one has and the phreatic surface is material with respect to the constituent water.

In general, however, is a dynamic quantity which must be related to the other basic fields such as the constituent densities and velocities on either side of the phreatic surface. Thus, additional relations are needed. These follow from the jump conditions of mass across the phreatic surface for the mixture as a whole and for the constituents. We begin with the jump condition of mass for the mixture as a whole, viz.

,

and transform its lefthand side (LHS) and righthand side (RHS) as follows

By equating these relations, we deduce the following relation

(4.37)

This is a first equation involving the accretion rate and the jumps of the densities, the porosities and the normal velocities.

Two further relations are obtained from the jump conditions of the individual constituents. Considering the water component first, we may write

(4.38)

It follows from this that, if either the jump of water density or the accretion rate vanish, then there is no jump in the normal component of the water velocity. Analogously, for the ice constituent

,

or after some re-arrangement

. (4.39)

It follows from Equations (4.37), (4.38) and (4.39) that if , then vw. n is continuous across the phreatic surface as is (1-n)(vI — vw). n; the water density pw and n may still suffer a jump across the phreatic surface in this case. On the other hand, if the porosity is continuous, , then Equations (4.37)–(4.39) imply that VI. n is continuous and nothing more. Both, vw. n and vI. n are continuous if .

This analysis suggests that the following jump conditions (of mass) may be reasonable to impose at the phreatic surface

(4.40)

As far as the jump conditions of linear momenta are concerned, they read

for the ice and water constituents, respectively. For ice this yields

(4.41)

and for water

(4.42)

This equation tells us that is parallel to n and thus . Ignoring in Equations (4.41) and (4.42) the momentum jumps of the ice and water, these imply

, (4.43)

in which continuity of the porosity n has also been assumed.

We close with a brief remark on energy. Its field equation has the form of Equation (4.28) both in the saturated and the non-saturated regions. In this equation, diffusive energy flux is ignored, so the energy-jump condition reduces to

. (4.44)

The internal energy in this approximation must be continuous across the phreatic surface.

This completes the formulation of the transition conditions at the phreatic surface.

To summarize the results of this section: we have presented a mathematical model of a wholly temperate ice mass (glacier) that consists of two disjoint sub-regions (1) ice, of which the pore volume (cracks, crevasses, grain boundaries) is saturated with water and (2) ice whose pore space is only partly filled with water. In this region, capillary effects are ignored.

For the region where the ice is saturated with water, the governing field equations comprise Equations (4.15) and (4.23) for the 15 fields, n (porosity), pw (pore-water pressure), VI and Vw (ice and water velocities), (deviatoric ice stress), p ^I (ice pressure) and M (melting rate). They are obtained from two constituent mass-and momentum-balance laws, constitutive relations for the ice stress and the ice-water interaction force (Darcy’s law), and an energy balance for the mixture as a whole (Equation (4.23)).

The temperate, non-saturated region is described by physically similar Equations (4.27) and (4.28) now for 16 fields. The difference is that the peculiar densities for ice and water are now independent. Boundary conditions comprise kinematic and dynamic statements at the free surface (Equations (4.31) and (4.35)), at the base (Equation (4-30a)) and at the phreatic surface (Equations (4.36), (4.43) and (4.46).

5. Concluding Remarks

In the foregoing analysis, I have presented detailed analyses of cold, polythermal and temperate ice masses, emphasizing thereby conceptual rigor and applicability of the deduced theoretical model.

For cold ice masses, the theoretical model is farthest advanced of all. The governing equations are understood, their numerical solution has been pushed ahead by a few computational enthusiasts and the thermo-mechanically coupled response of academic and realistic ice sheets to climatological forcings can be analyzed through inter-glacial cycles.

Calculations on the response of Antarctica and Greenland as cold ice sheets as well as other ones on Arctic glaciers indicate that some, if not many, of our terrestrial ice masses possess temperate patches in an otherwise cold environment. Ice masses of this sort are called polythermal. For these, several variants of diffusive theories were presented. It was made clear that deduction of the adequate theory had to satisfy the theoretician’s requirement of being consistent with a Stefan-type condition at the cold-temperate transition surface, plus the practitioner’s wishes of being minimal to the extent that no prescription of a boundary condition was mathematically required that could not be measured. Such a theoretical formulation is in our hands. First computational results have been obtained with this formulation. Further, full-scale analyses are in sight and the theory promises to be the next generation of ice-sheet modeling in future climate dynamics.

Least advanced of all is a theoretical model for temperate ice masses. Formulations for these may at present perhaps be overly sophisticated for field and applied glaciologists — most of these seem to be involved with questions of essentially local nature — however, I think the time is ripe to start a serious attempt to deduce a theoretical model for the dynamical interplay of the ice and the interstitial water whose water table performs large oscillations on a daily time-scale. The theory presented in section 4 is such a proposal and pilot computations ought to be performed with it in order to check its mathematical adequacy but equally so also to demonstrate its usefulness.