Polythermal ice sheets and glaciers contain both cold ice and temperate ice. We present two new models to describe the temperature and water content of such ice masses, accounting for the possibility of gravity- and pressure-driven water drainage according to Darcy's law. Both models are based on the principle of energy conservation; one additionally invokes the theory of viscous compaction to calculate pore water pressure, and the other involves a modification of existing enthalpy gradient methods to include gravity-driven drainage. The models self-consistently predict the evolution of temperature in cold ice and of water content in temperate ice. Numerical solutions are described, and a number of illustrative test problems are presented, allowing comparison with existing methods. The suggested models are simple enough to be incorporated in existing ice-sheet models with little modification.

This paper concerns the modelling of temperature and water content in polythermal ice masses. These are ice sheets and glaciers that contain both cold ice, below its pressure melting temperature, and temperate ice, which is a two-phase mixture of ice and water at the melting temperature. Such models are an important component of ice-sheet simulations because the rheology of ice and hence its ability to flow depend strongly on temperature and water content. The specific focus of this paper is on models that allow for gravity-driven drainage of meltwater in temperate ice.

There is already an extensive literature on the subject of polythermal ice
masses

Ice temperature is governed by an advection–diffusion equation, with dissipation due to viscous creep acting as a heat source. If the temperature reaches the melting point, further dissipative heating leads to internal melting, and the resulting temperate ice can be described using mixture theory. Combining both options leads to a free-boundary problem for temperature and water content (the two are often combined, as enthalpy).

A key ingredient of this theory is a constitutive law for water transport in
the temperate ice. There is little experimental constraint on this.

Two previously suggested models for the water transport are to make use of
either Darcy's law or a diffusive flux

In a companion paper

Enthalpy gradient models have become widely used in the literature

In Sect.

We assume that ice velocity satisfies the usual Stokes equations,

Here

It is a basic assumption of this study that these equations provide an
adequate description of ice deformation even when the temperate “ice” is
actually a two-phase mixture of ice and pore water (the additional
possibility of partial saturation, with air as well as water in the pore
space, is ignored). As discussed by

Conservation of energy is expressed as

Here

Following standard practice, it is convenient to combine
Eqs. (

The energy and moisture Eqs. (

Enthalpy is transported by advection and conduction in the cold regions, and by advection and water transport in the temperate regions (as well as conduction due to variations of the pressure-dependent melting point).

The jump condition (Stefan condition) that naturally arises from
Eq. (

In the following subsections, we describe three specific models for the water
flux

This model assumes that water transport occurs according to Darcy's law, driven by gravity and pressure gradients. We write

If we approximate the ice pressure as hydrostatic,

Substituting this into Eq. (

A common approach is to suppose that the water transport is diffusive, so

This model does not allow for significant drainage and can result in the
build-up of large amounts of water in the ice. A cap on porosity and removal
of additional water are often implemented by adding a term

Motivated by the desire to capture gravity-driven drainage of the meltwater,
this model takes the gravitational component of Darcy's law from
Eq. (

This model has the advantage of allowing water to drain from the ice, as in
the compaction pressure model, but without the need to solve an additional
equation for

The model in Eq. (

This section explains the numerical method used to solve the energy model
consisting of Eqs. (

We employ a first-order finite volume discretization with a simple operator
splitting method. Temperature, porosity, and effective pressure are
discretized on a regular grid, while ice velocity and water flux are defined
on a staggered grid. Using a harmonic average to define the permeability on
the staggered grid means that the condition

At each time step, Eq. (

For clarity, we summarize the time discretization here, with superscripts
referring to the time step (ignore the hats on some of the variables, which
are included for later reference):

Values of parameters used for example solutions. The
value of

Boundary conditions are discussed by

This algorithm has been tested against steady-state solutions calculated
using a shooting method as well as asymptotic methods

Since ice sheets and glaciers typically have a shallow aspect ratio, it is
common to approximate conduction terms using the fact that vertical (

In Sect.

The model can accommodate a pressure-dependent melting point

Parallel-sided slab setup. Arrows indicate the bed-parallel velocity, and shading indicates the region of temperate ice (to be found).

A simple one-dimensional test problem for these models is the case of a
parallel-sided slab of ice, as seen in Fig.

For a genuinely parallel moving slab the normal velocity

The top surface

One of the reasons for exploring this simple setup is that semi-analytical
solutions are possible, allowing the algorithms to be tested. Constructing
these solutions involves straightforward algebraic manipulations and solution
of a transcendental equation for the position of the cold-temperate
transition (details may be found, for example, in

Steady-state solutions for the parallel-sided slab,
showing temperature and porosity as a function of height. Normal velocities
are

Figure

The standard enthalpy gradient model closely approximates the no-water-flux
solution for both downward- and upward-moving ice, the diffusive term
naturally smoothing out the discontinuity in porosity at the freezing
interface (provided

For the case of downward-moving ice (panel a), the location of the
cold-temperate boundary is almost identical between all the models. This is
reassuring, since the location of the interface in this case should be
determined by the twin conditions

For the case

For the case of upward-moving ice (panel c), the location of the
cold-temperate boundary in the two models with gravity-driven drainage is now
different from the no-water-flux solution. This is because the drainage
results in lower porosities and thus less advected latent heat transport
than in the no-water-flux case; the jump condition Eq. (

Steady-state solutions for the parallel-sided slab
with the compaction pressure model, showing porosity and effective pressure
as a function of height. Normal velocities are

Temperate glacier setup, showing boundary
conditions for the compaction pressure model. Solid lines show ice
streamlines and red shading shows the magnitude of viscous dissipation. Here

Temperate glacier solutions for

The effective pressure in the compaction pressure model is shown in
Fig.

Overall, we see that inclusion of gravity-driven drainage with either model reduces the porosity, as expected, and alters the location of the cold-temperate interfaces due to changes in the advected latent heat. These basic features of the models carry over into more complex problems in more dimensions.

We illustrate the effect of water transport further by considering an
idealized mountain glacier as seen in Fig.

The effective viscosity and dissipation rate are given by
Eq. (

Both the upper surface

It is worth pointing out that in a more complete model an energy balance and
firn compaction model at the upper surface would prescribe both

Ice cap setup. Solid lines show ice streamlines,
red shading shows the magnitude of viscous dissipation, and blue shading
shows the region of temperate ice (to be found). Here

Figure

We see that the compaction pressure model (panel d) and the modified enthalpy
gradient model (panel c) give very similar results. The meltwater flux at the
bed is of course larger than for the standard enthalpy gradient model. The
flux is similar in the two gravity-driven drainage models to what reaches the
bed if the drainage function is included in the standard model (panel b), but
it does have slight differences; notably, it is transported further down
glacier, because the water is being advected as well as draining vertically,
rather than draining instantaneously to the bed. Of course, the actual amount
that drains is sensitive to the drainage function

For the final example we take a symmetric ice-sheet profile as seen in
Fig.

Ice cap solutions for

The upper surface

Figure

The predicted region of temperate ice is similar between the models, since it
is largely controlled by the location of dissipative heating. However, there
are some differences: notably the downstream end of the temperate region
extends closer to the ice surface in the standard enthalpy gradient method
without drainage, because of the large advected latent heat transport, which
is all released at the freezing interface. This is essentially the same
behaviour as seen in the upward-moving parallel slab in
Fig.

The distribution of porosity is very similar between the modified enthalpy gradient method and the compaction pressure model. It increases towards the bed because of the accumulated drainage of water from above, whereas the porosity distribution in the standard model with drainage is almost uniform because it is evacuated directly to the bed.

Ice cap solutions as in Fig.

We have presented methods for modelling the water content of temperate ice that make use of Darcy's law and thus allow water to drain under gravity. The models provide an alternative to commonly used enthalpy gradient methods, and we have shown that they can be solved in a similar way using an enthalpy-based formulation. The gravity-driven water fluxes are coupled to latent heat transport in temperate ice, so they can result in changes to the location of cold-temperate interfaces from what is predicted by other models. The new models correctly handle jump conditions at both melting and freezing cold-temperate interfaces.

We envisage that the compaction pressure method could quite easily be
implemented in existing ice-sheet models, requiring many fewer alterations
than a full two-phase formulation. Alternatively, the
modified enthalpy gradient method
has been shown to give very similar results in many cases. That model assumes
that the pore pressure is equal to the ice pressure, and the buoyancy

Although the compaction pressure model requires the addition of a new variable, the effective pressure, the computational cost of this new variable is relatively small. In fact, the compaction pressure model was often cheaper to solve than the modified enthalpy gradient model using our algorithm, the latter generally requiring smaller time steps to be taken. The advantage of including the effective pressure is that it allows for a connection with the subglacial drainage system and for the possibility that the pore water pressure becomes close to hydrostatic (in which case gravity-driven drainage does not occur as rapidly as one might otherwise think).

The compaction pressure model is essentially a stripped-down version of the
more detailed theory presented by

The permeability of the temperate ice evidently plays an important role in
this theory. It is unlikely that a particular value for

The length scale over which the effective pressure gradients play a role is
seen from Eq. (

The nature of the solutions in the case of a small compaction length is
explored by

This situation of upward advection into a temperate region is unlikely to be common, if indeed it ever occurs, but one can imagine that such behaviour might occur in an overdeepening, where de-pressurization results in freeze-on of subglacial water. The margins of ice streams are another place where this may be relevant if lateral advection brings cold ice underneath ice that has been heated in the shear margin.

We end by commenting on straightforward extensions to this model. One
possibility is to include the term

A similar amendment that would be straightforward in a coupled model is to
include the deviatoric stress components in Darcy's law (

The fact that deviatoric stresses in the ice could lead to differential
compaction of different orientations of veins

The authors declare that they have no conflict of interest.

Ian J. Hewitt is supported by a Marie Curie FP7 Career Integration Grant within the European Union's 7th Framework Programme. Christian Schoof acknowledges the support of NSERC grants 357193-13 and 446042-13 and a Killam Faculty Research Fellowship at UBC. Edited by: O. Gagliardini Reviewed by: A. Aschwanden and one anonymous referee