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 Reference Fowler and LarsonFowler and Larson (1978), was corrected by Reference HutterHutter (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 (Reference Hutter, Blatter and FunkHutter and others, 1988; Reference BlatterBlatter, 1991; Reference Blatter and HutterBlatter 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. Reference FowlerFowler (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 Reference Hutter and EngelhardtHutter and Engelhardt ( 1988a, Reference Hutter and Engelhardtb). 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 (Reference HutterHutter, 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 (Reference HutterHutter, 1983; Reference MorlandMorland, 1984; Reference Hutter, Yakowitz and SzidarovszkyHutter 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 (Reference Hutter, Yakowitz and SzidarovszkyHutter and others, 1986; Reference Hindmarsh, Morland, Boulton and HutterHindmarsh and others, 1987, Reference Hindmarsh, Boulton and Hutter1989; Reference HerterichHerterich, 1988, Reference Herterich1990; Reference Hindmarsh and HutterHindmarsh and Hutter, 1988; Reference CalovCalov, 1989, Reference Calov1990; Reference HuybrechtsHuybrechts, 1990a, Reference Huybrechtsb; Reference Huybrechts and OerlemansHuybrechts and Oerlemans, 1990; Reference Letréguilly, Reeh and HuybrechtsLetréguilly and others, 1990a, Reference Letréguilly, Huybrechts and Reehb,Reference Letréguilly, Reeh and Huybrechtsc). 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):
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)6 is the non-linear stress-stretching relationship with a temperature-dependent rate factor A(T) (usually of Arrhenius type) and a creep-response function
with appropriate values for aj (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≠0) is essential at low stressesFootnote * 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 •Df defined by Ff(x,t) = 0:
and along the base •Df,defined by 
:†
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)2 requires the free surface to be stress-free, while Equation (2.4)2 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, Reference WeertmannWeertman, 1957, Reference Weertmann1961, Reference Weertmann1964, Reference Weertmann1967, Reference Weertmann1971, Reference Weertmann1979; Reference LliboutryLliboutry, 1964, 1968, Reference Lliboutry1975, Reference Lliboutry1979, Reference Lliboutry1987; Reference NyeNye, 1969, Reference Nye1970; Reference KambKamb, 1970; Reference MorlandMorland, 1976a,Reference Morlandb; Reference FowlerFowler 1981a, Reference Fowlerb, Reference Fowler1986, Reference Fowler1987; Reference ShreveShreve, 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)3, 4 is of the Neumann type if the basal temperature does not reach melting, or of the Dirichlet type if the melting temperature
is reached. Here T
0(= 0.01°C) is the melting temperature at p
0(= 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)2 with the boundary conditions in Equations (2.3)2 and Equation (2.4)2 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)4 subject to the velocity distribution
. -
(3) Update the temperature field by solving Equation (2.1)3 subject to the boundary conditions in Equation (2.3)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).
)
.
. (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(t0). The temperature may simply be prescribed (because for Τ an initial value problem is solved), v(x,t0) follows from the solution of the momentum equation at given T(x,t0). 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 (Reference HutterHutter, 1983; Reference MorlandMorland, 1984, Reference Hutter, Yakowitz and SzidarovszkyHutter 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. Reference HerterichHerterich, 1988, or Reference Budd and JenssenBudd 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
and along the base, defined by z=h f(x,y,t):
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)4 ), 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
and the velocities become
with
in which (u, v) b are the horizontal components of the basal velocity vector evaluated in Equation (2.10)1 as prescribed in Equation (2.8)1 ,ᐁ denotes the two-dimensional horizontal gradient operator. Cs and Cg 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
and for the polynomial law in Equation (2.2),
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 (Reference HutterHutter, 1980, Reference Hutter1981, Reference Hutter1983) 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 (Reference HutterHutter, 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 (Reference HerterichHerterich, 1988; Reference Dahl-JensenDahl-Jensen, 1989a, Reference Dahl-Jensenb). 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)
otherwise, in which horizontal thermal diffusion has been ignored. Thus, Equation (2.12) is parabolic in x and y, suggesting forward-marching procedures (Reference Hutter, Yakowitz and SzidarovszkyHutter 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)1 :
q is the volume flux and can easily be determined from Equation (2.10)1 :
where
For a power-law fluid, this becomes
and for the polynomial law in Equation (2.2)
where each of the QjS 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
in which Δ2 is the horizontal (two-dimensional) Laplacian operator and {•} is the curly bracketed term in Equation (2.14). As long as Cs and Q are bounded, which we shall now assume, the term {•} will vanish at hf = hb , making Equation (2.18) singular at the margin. For the no-slip condition, Cs = 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
So with
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
α=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 = xm is approached, unless the basal drag coefficient C ∞(h f–h b)–1 is close to the margin. Reference Morland and JohnsonMorland 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 Reference Hindmarsh, Boulton and HutterHindmarsh (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 Reference RaymondRaymond (1983) and Reference Hutter, Yakowitz and SzidarovszkyHutter 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 Reference Szidarovszky, Hutter and YakowitzSzidarovszky 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 (Reference RaymondRaymond, 1983). Reference Szidarovszky, Hutter and YakowitzSzidarovszky 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 (Reference Hutter, Yakowitz and SzidarovszkyHutter and others, 1986; Reference Hindmarsh and HutterHindmarsh and Hutter, 1988; Reference Hindmarsh, Boulton and HutterHindmarsh 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 (Reference Elliot and OckendonElliot 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,
where Τ — TM is the homologous temperature and β and γ are constants. Reference Hindmarsh, Boulton and HutterHindmarsh and others (1989) performed computations for
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 (Reference Hutter and EngelhardtHutter and Engelhardt, 1988a, Reference Hutter and Engelhardtb). 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 (Reference MüllerMü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 (Reference Fowler and LarsonFowler and Larson, 1978; Reference HutterHutter, 1982; Reference FowlerFowler, 1984; Reference Hutter and EngelhardtHutter and Engelhardt, 1988a, Reference Hutter and Engelhardtb; Reference Hutter, Blatter and FunkHutter and others, 1988; Reference BlatterBlatter, 1991; Reference Blatter and HutterBlatter 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 (Reference Blatter and HutterBlatter, 1991; Reference BlatterBlatter 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:
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 Reference LliboutryLliboutry (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
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)5 ). The moisture production M as given in Equation (3.1)6 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 downFootnote *
then Equation (3.1)6 would follow from Equation (3.1)2 and Equation (3.3) by ignoring moisture flux j and heat flux q. Reference Fowler and LarsonFowler and Larson (1978) used Equations (3.3) and (3.1)2 , and requested moisture flux and heat flux in temperate ice to vanish and then wrote Equation (3.1)6 . Reference HutterHutter (1982) ignored Equation (3.3), included a non trivial moisture flux and postulated Equation (3.1)6 . A model not involving a closure condition for moisture production, Equation (3.1)6 , 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, Reference HutterHutter’s (1982) model is used.
If relations in Equation (3.1)5,6 are substituted in Equation (3.1)2 , the mass-balance equation for the moisture content reads
an equation which is formally the same as the energy Equation (2.1)3 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 Reference Blatter and HutterBlatter (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 V1 is the ice velocity, then along a non-deforming bed one has
where V1 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
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
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)
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 waterFootnote *
Here, Equation (3.6) has been used, and in the last line the contribution from the kinetic energy has been ignored. With
Equation (3.9) finally assumes the form
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)6 . 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
but this equation is only meaningful when the ice above the temperate base is temperate.
The remaining boundary condition is the sliding law
as prescribed already in Equation (2.4)2 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 (Reference HobbsHobbs, 1974; Reference LliboutryLliboutry, 1987b, Reference Lliboutry1993; Reference Alts and HutterAlts and Hutter, 1988a, Reference Alts and Hutterb,Reference Alts and Hutterc, Reference Alts and Hutter1989). 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 (Reference Stauffer and OeschgerStauffer and Oeschger, 1979).
In classical continuum thermodynamics, a phase change surface is defined as a surface at which temperature and tangential velocity are continuous (Reference HutterHutter, 1983; Reference MüllerMü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:
Furthermore, the jump conditions of moisture content, momentum, energy and entropy must hold (Reference HutterHutter, 1982, Reference Hutter1983):
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)1
, 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)2
is only approximate because it ignores the contribution 
which is very small. Equation (3.15)3
is the entropy balance when 
has been used (see Reference HutterHutter, 1982; Reference Hutter, Blatter and FunkHutter and others, 1988; Reference BlatterBlatter, 1991). The energy balance is not needed, because it reduces to the same statement as the entropy balance. Indeed, from
this claim is immediate. Combining Equations (3.15)1 , and (3.15)3 , finally yields
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 Reference Fowler and LarsonFowler 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. Reference HutterHutter (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 Τ < TM , 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
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.
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.
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,
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
the boundary-layer solution on the temperate side of the CTS, see Figure 2, takes the form
Furthermore,
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. Reference BlatterBlatter (1991) and Reference Blatter and HutterBlatter 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)2
). 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.
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
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
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 Reference Hutter, Blatter and FunkHutter and others (1988), Reference Blatter and HutterBlatter and Hutter (1991), and Reference BlatterBlatter (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 (Reference Blatter and HutterBlatter 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 (Reference FowlerFowler, 1984; Reference Hutter and EngelhardtHutter and Engelhardt, 1988a, Reference Hutter and Engelhardtb) has so far largely been ignored.
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 (Reference RöthlisbergerRöthlisberger, 1972; Reference ShreveShreve, 1972; Reference WeertmannWeertman, 1972; and many others) with basal sliding mechanisms in parts based upon linked cavity configurations of the basal water conduit system (Reference WalderWalder, 1986; Reference KambKamb, 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 (Reference Vulliet and HutterVulliet 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 saltsFootnote * 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
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
in which M α is the constituent mass production and ma is the constituent momentum production (= specific interaction force), which satisfy the conservation constraints
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.
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
this shows that their sum is made up of the pore pressure and the constitutive parts of the constituent stresses, tw and tI (this sum is not identical to the mixture stress, but it agrees with it, if the diffusive momentum fluxes are ignored (see Reference Hutter and EngelhardtHutter and Engelhardt, 1988a, Reference Hutter and Engelhardtb). 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
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 tI and tw, 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 tw = 0. Alternatively, the partial stress tensor tI is postulated as
Accordingly, the materially dependent stress tensor tI 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 DI need not be traceless, despite the incompressibility assumption,Footnote 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
This is a materially objective vector quantity, i.e. vF transforms as a vector under (non-inertial) observer transformations, even though vw and vI 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
that is objective. We therefore generalize Darcy’s law by postulating that
with a constant of proportionality μ/k*.μ is a viscosity
and k* is the so-called coefficient of absolute permeability (m2). Combination of Equations (4.8) and (4.9) yields the generalized Darcy law
It involves a contribution of the volumetric melting rate. In earlier formulations (Reference FowlerFowler, 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,
Here, use has been made that tw = 0 and. Alternatively, Equation (4.5)4 (with Equation (4.10) substituted for the interaction force) can be solved for the seepage velocity
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
and may then write Equation (4.12) in the form
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−1) 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
Here, all simplifying assumptions have been incorporated, and it is understood that the expressions for vF, 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), pw (1), v1 (3), vw (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)6 must then be replaced by
, an equation which guarantees the tracelessness of, (4.16).
-
(2) It is usually justified to assume that
. In this case the seepage velocity is approximately given by vF = nvw, and the generalized Darcy law may be written as. (4.17)Substituting this expression into Equation (4.15)1 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.
.
. In this case the seepage velocity is approximately given by vF = nvw, and the generalized Darcy law may be written as
Let us now turn to the energy equation for the mixture as a whole; it may be stated as
in which the “material” derivative d(·)/dt is defined as the time rate of change for an observer following the barycentric velocity
∈, 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
Here, we have set T
1=T
w=T
M, and have assumed that 
being the latent heat of fusion and C1
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
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
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
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 vw — vI, with an additional dependence on ρ
I and/or n. Thus,
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,Footnote but ice may support a deviatoric viscous component. So, we assume
where, at a first approximation, A is assumed constant. Moreover, the constitutive part of σ 1 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
and constitute 15 equations for the 16 fields 
. The equation that is still missing is, of course, the balance of energy,
in which heat conduction is ignored as before. A derivation along the lines of the previous deductions yields
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 Reference Hutter and EngelhardtHutter and Engelhardt (1988a, Reference Hutter and Engelhardtb).
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.Footnote 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, VI.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 (Reference LliboutryLliboutry, 1968, Reference Lliboutry1975, Reference Lliboutry1979b). We note, finally that Equation (4-30b) automatically satisfies the tangency condition vI · 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
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,
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
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)~−1
, 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.
which arises in the kinematic 
by the factor (1 — n)~−1
, which would be expected.
, then 
, as it must be.
consists only of surface run-off, then 
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
(A second condition is not needed because integration of Equation (4.27)4 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, Reference Hutter and EngelhardtHutter 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 = hph(x,y,t), whose evolution equation may be written as
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
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
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
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
As far as the jump conditions of linear momenta are concerned, they read
for the ice and water constituents, respectively. For ice this yields
and for water
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 implyFootnote
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
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.
Acknowledgements
This paper was presented at the conference held to mark the occasion of Professor L. Lliboutry’s retirement (23 November 1990). It was submitted to the Laboratoire de Glaciologie, Grenoble, on 20 December 1990 for publication and now appears, more than 2 years later, with few amendments and alterations. I thank R. Hindmarsh and D. MacAyeal for their constructive review and H. Blatter and R. Grève for inspirative discussions. The revised paper was drafted while I was holding a Visiting Professorship at the Geographical Institute, ΕΤΗ Zürich.
I thank Mrs Danner and Mr Schneider for their diligent production of the typed manuscript.
The accuracy of references in the text and in this list is the responsibility of the author, to whom queries should be addressed.
The accuracy of references in the text and in this list is the responsibility of the author, to whom queries should be addressed.
= constant, 
= constant, and does not imply that vI or vW are solenoidal.
is also imposed in 
is to survive.
