2.1 Cold ice

The standard model to describe ice sheet flow is that of an incompressible power-law viscous fluid, governed by equations of mass and force balance,

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\tau}_{ij}}{\partial x_{j}}-\frac{\partial p}{\partial x_{i}}=-{\it\rho}g_{i}, & \displaystyle\end{eqnarray}$$

(2.3a-c ) $$\begin{eqnarray}\displaystyle {\it\tau}_{ij}=A^{-1/n}D^{1/n-1}\unicode[STIX]{x1D60B}_{ij},\quad \unicode[STIX]{x1D60B}_{ij}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial u_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right),\quad D=\sqrt{\unicode[STIX]{x1D60B}_{ij}\unicode[STIX]{x1D60B}_{ij}/2}, & & \displaystyle\end{eqnarray}$$

where repeated indices are summed over. Here, $\boldsymbol{u}=(u_{1},u_{2},u_{3})$ is the velocity vector, $\unicode[STIX]{x1D60B}_{ij}$ is the strain rate tensor, $D$ is its second invariant, ${\it\tau}_{ij}$ is the deviatoric stress tensor, $p$ is the pressure, ${\it\rho}$ is the density of ice and $\boldsymbol{g}=(g_{1},g_{2},g_{3})$ is the acceleration due to gravity. $A$ and $n$ are coefficients in the flow law. Temperature enters the problem because the coefficient $A$ depends strongly on temperature, usually with an assumed Arrhenius dependence. Thus, the model is coupled to an energy equation,

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}c\left(\frac{\partial T}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}T\right)-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left(K\boldsymbol{{\rm\nabla}}T\right)=a, & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle a={\it\tau}_{ij}\unicode[STIX]{x1D60B}_{ij}, & \displaystyle\end{eqnarray}$$

in which $c$ is the specific heat capacity, $K$ is the thermal conductivity and $a$ is the viscous dissipation.

2.2 Temperate ice

When ice reaches the melting point $T_{m}$ , there are three potentially important changes required to this model. Firstly, the temperate ice is porous, and with a non-zero porosity ${\it\phi}$ (volume fraction), the flow-law coefficient may depend on porosity as well as temperature $A=A(T,{\it\phi})$ (Duval Reference Duval1977; Lliboutry & Duval Reference Lliboutry and Duval1985). Secondly, the energy balance must be altered to account for the latent heat required to melt or refreeze the interstitial water (and hence determine the evolution of the porosity). Thirdly, the porous ice matrix may be permeable, and the ability for water to drain means that the ice–water mixture is no longer necessarily incompressible.

To address these changes, we consider the general mixture form of the mass, momentum and energy conservation equations for permeable temperate ice, broadly working by analogy with models for melt extraction from the Earth’s mantle (e.g. Hewitt & Fowler Reference Hewitt and Fowler2008). The ice is assumed to have porosity ${\it\phi}$ , and we denote the internal melting rate (the rate of mass conversion from ice to water) by $m$ . In terms of the average pore-water velocity $\boldsymbol{u}_{w}$ we also define the relative moisture flux (or Darcy transport velocity)

(2.6) $$\begin{eqnarray}\displaystyle \boldsymbol{j}\equiv {\it\phi}(\boldsymbol{u}_{w}-\boldsymbol{u}), & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}$ is still the velocity of the ice; $\boldsymbol{j}$ is commonly referred to as the ‘diffusive’ moisture flux, but we refrain from that terminology here since it is not (in our formulation) strictly diffusive. Conservation of ice and water mass are then expressed as

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial (1-{\it\phi})}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }[(1-{\it\phi})\boldsymbol{u}]=-\frac{m}{{\it\rho}}, & \displaystyle\end{eqnarray}$$

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\phi}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }[{\it\phi}\boldsymbol{u}+\boldsymbol{j}]=\frac{m}{{\it\rho}_{w}}, & \displaystyle\end{eqnarray}$$

where ${\it\rho}_{w}$ is the density of water, and combining these gives an expression for overall mass conservation,

(2.9) $$\begin{eqnarray}\displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{j}-\frac{{\it\rho}_{w}-{\it\rho}}{{\it\rho}_{w}{\it\rho}}m. & & \displaystyle\end{eqnarray}$$

This expression holds in cold ice too, where $\boldsymbol{j}=0$ and $m=0$ so it reduces to the incompressibility constraint (2.1).

Accounting for the latent heat of melting $L$ , the overall energy balance is expressed as

(2.10) $$\begin{eqnarray}\displaystyle & & \displaystyle mL+\left[{\it\rho}c(1-{\it\phi})+{\it\rho}_{w}c_{w}{\it\phi}\right]\frac{\partial T}{\partial t}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\left[{\it\rho}c(1-{\it\phi})\boldsymbol{u}+{\it\rho}_{w}c_{w}({\it\phi}\boldsymbol{u}+\boldsymbol{j})\right]\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}T-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left(K\boldsymbol{{\rm\nabla}}T\right)=a,\end{eqnarray}$$

where $c_{w}$ is the specific heat capacity of the water and the viscous dissipation $a$ should strictly include additional terms due to the energy dissipated by the relative water flow, although these can be shown to be small so are ignored (see e.g. Fowler Reference Fowler1984). Again, this equation reduces to (2.4) in cold ice, where $m$ , ${\it\phi}$ and $\boldsymbol{j}$ are all zero. For the temperate ice, since the temperature is by definition constrained to be at the melting point $T_{m}$ , this equation serves to determine the melting rate $m$ and hence the source terms in (2.8) and (2.9). In particular, if $T_{m}$ is taken as constant, we simply have

(2.11) $$\begin{eqnarray}\displaystyle m=a/L. & & \displaystyle\end{eqnarray}$$

If the effect of pressure and impurities on the melting point are accounted for, these lead to corrections to this relationship.

In order to solve (2.8), the volume flux of interstitial water must be specified. We adopt Darcy’s law with permeability $k_{0}{\it\phi}^{{\it\alpha}}$ to write

(2.12) $$\begin{eqnarray}\displaystyle \boldsymbol{j}=k_{0}{\it\phi}^{{\it\alpha}}\left({\it\rho}_{w}\boldsymbol{g}-\boldsymbol{{\rm\nabla}}p_{w}\right), & & \displaystyle\end{eqnarray}$$

where $k_{0}$ and ${\it\alpha}$ are positive constants. In our analysis later, we assume that ${\it\alpha}\geqslant 2$ ; such values can be justified from the small porosity limit of a Carman–Kozeny relationship, with the precise choice of ${\it\alpha}$ depending on the partitioning of water between veins (triple junctions at crystal boundaries) and nodes at which such veins meet (Nye & Frank Reference Nye and Frank1973).

Although it is tempting to assume that the pore pressure is equal to the ice pressure, $p_{w}=p$ , this is too restrictive an assumption. Instead we follow the theory of viscously compacting media (McKenzie Reference McKenzie1984) and define the compaction pressure (or effective pressure; Terzaghi Reference Terzaghi1923; Biot Reference Biot1941), $p_{e}$ , according to the relationship

(2.13) $$\begin{eqnarray}\displaystyle p_{e}\equiv p-p_{w} & & \displaystyle\end{eqnarray}$$

and relate it to the compaction rate according to

(2.14) $$\begin{eqnarray}\displaystyle p_{e}=-\frac{\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}}{{\it\zeta}}, & & \displaystyle\end{eqnarray}$$

where ${\it\zeta}={\it\eta}/{\it\phi}$ is the bulk viscosity and ${\it\eta}=(1/2)A^{-1/n}D^{1/n-1}$ is the effective viscosity. This relationship is equivalent to the statement

(2.15) $$\begin{eqnarray}\displaystyle \frac{\partial {\it\phi}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}=\frac{m}{{\it\rho}}-\frac{(1-{\it\phi}){\it\phi}}{{\it\eta}}p_{e}, & & \displaystyle\end{eqnarray}$$

which expresses pore-scale mass conservation on the assumption that creep closure of the pores occurs at a rate proportional to the pressure difference $p_{e}$ divided by the effective viscosity (cf. the standard equations for melt opening and creep closure of a Röthlisberger channel or a borehole (Nye Reference Nye1953; Fowler Reference Fowler1984)). Although $p_{e}$ is the analogue of the ‘effective pressure’ $N$ , commonly used to describe subglacial drainage conditions, we adopt the notation $p_{e}$ so as to avoid any confusion (i.e.  $N$ is reserved for use at the ice sheet or glacier bed, whereas $p_{e}$ is defined within the ice).

At the microscopic scales relevant to compaction, surface energy effects may also become relevant (e.g. Bercovici, Ricard & Schubert Reference Bercovici, Ricard and Schubert2001). At triple grain junctions (where connected liquid veins are likely to exist), surface tension will reduce the compaction rate for a given compaction pressure. A somewhat more general version of (2.14)–(2.15) would be

(2.16a,b ) $$\begin{eqnarray}\displaystyle p_{e}-p_{s}({\it\phi})=-\frac{\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}}{{\it\zeta}},\quad \frac{\partial {\it\phi}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}=\frac{m}{{\it\rho}}-\frac{(1-{\it\phi}){\it\phi}}{{\it\eta}}\left[p_{e}-p_{s}({\it\phi})\right], & & \displaystyle\end{eqnarray}$$

where the term $p_{s}({\it\phi})$ is a curvature- and hence porosity-dependent stress that serves in this case to reduce compaction rates; it is analogous to the term ${\it\sigma}\text{d}{\it\alpha}/\text{d}{\it\phi}$ in Bercovici, Ricard & Schubert (Reference Bercovici, Ricard and Schubert2001). In part motivated by the fact that there is no obvious constitutive relation for $p_{s}$ , beyond the requirement that $p_{s}$ should decrease with ${\it\phi}$ (and presumably be positive), we persist with our original, simpler equations (2.14)–(2.15). However, we recognize that the reduced model we develop below could be generalized to retain $p_{s}$ , and that this may become necessary: in practice, $p_{s}$ should scale as ${\it\sigma}/R$ , where ${\it\sigma}$ is surface free energy and $R$ is the radius of curvature, and will become comparable to $p_{e}$ when $R\sim {\it\sigma}/p_{e}\sim 5\times 10^{-5}$  m for ${\it\sigma}=3\times 10^{-2}$  J m $^{-2}$ and $p_{e}\sim 2\times 10^{3}$  Pa as estimated below. Such radii of curvature are not impossible at triple grain boundaries, especially where porosities are small (see e.g. Nye & Mae Reference Nye and Mae1973).

Using (2.14) to eliminate $p_{w}$ , Darcy’s law (2.12) becomes

(2.17) $$\begin{eqnarray}\displaystyle \boldsymbol{j}=k_{0}{\it\phi}^{{\it\alpha}}\left(({\it\rho}_{w}-{\it\rho})\boldsymbol{g}+\boldsymbol{{\rm\nabla}}p_{e}+\boldsymbol{{\rm\nabla}}p_{r}\right). & & \displaystyle\end{eqnarray}$$

Here we have substituted $\boldsymbol{{\rm\nabla}}p={\it\rho}\boldsymbol{g}+\boldsymbol{{\rm\nabla}}p_{r}$ , thereby defining the reduced (i.e. non-hydrostatic) ice pressure $p_{r}$ .

Allowing for compaction of the ice matrix strictly speaking not only changes the divergence-free nature of the velocity field, but also the form of the remaining Stokes equations (2.2)–(2.3). The standard definition of the deviatoric stress requires a vanishing trace ${\it\tau}_{ii}$ , and we additionally have to account for the stress supported by pore water. Conforming with these requirements, we put

(2.18) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial x_{j}}[(1-{\it\phi}){\it\tau}_{ij}]-\frac{\partial }{\partial x_{i}}[(1-{\it\phi})p+{\it\phi}p_{w}]=-(1-{\it\phi}){\it\rho}g_{i}-{\it\phi}{\it\rho}_{w}g_{i}, & \displaystyle\end{eqnarray}$$

(2.19) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\tau}_{ij}=AD^{1/n-1}\left(\unicode[STIX]{x1D60B}_{ij}-{\textstyle \frac{1}{3}}D_{kk}{\it\delta}_{ij}\right), & \displaystyle\end{eqnarray}$$

where $A=A(T,{\it\phi})$ . Note that in the limit of ${\it\phi}=0$ , we have $D_{kk}\equiv \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0$ , and (2.18)–(2.19) again reduce to their incompressible cold-ice versions, (2.2)–(2.3).

2.3 Cold–temperate boundaries

At any boundary between cold and temperate ice it is necessary to apply jump conditions to ensure that mass, momentum and energy are conserved. These conditions are derived by balancing the fluxes and stresses at the boundary. In addition, we assume that the ice velocity must be continuous, since there would otherwise be non-integrable stress singularities (just as temperature must be continuous or there would be a heat flux singularity). Denoting the normal to the boundary by $\boldsymbol{n}$ and its migration velocity by $\boldsymbol{v}$ (the tangential component of $\boldsymbol{v}$ is defined to be the same as that of the ice velocity $\boldsymbol{u}$ ), and denoting the temperate side of the boundary by $^{+}$ and the cold side by $^{-}$ , the conditions are expressed as

(2.20a,b ) $$\begin{eqnarray}\displaystyle \boldsymbol{u}^{+}=\boldsymbol{u}^{-},\quad T^{+}=T^{-}=T_{m}, & & \displaystyle\end{eqnarray}$$

(2.21) $$\begin{eqnarray}\displaystyle & \displaystyle [{\it\rho}(1-{\it\phi})(\boldsymbol{u}-\boldsymbol{v})]^{+}\boldsymbol{\cdot }\boldsymbol{n}+[{\it\rho}_{w}{\it\phi}(\boldsymbol{u}-\boldsymbol{v})+{\it\rho}_{w}\,\boldsymbol{j}]^{+}\boldsymbol{\cdot }\boldsymbol{n}=[{\it\rho}(\boldsymbol{u}-\boldsymbol{v})]^{-}\boldsymbol{\cdot }\boldsymbol{n}, & \displaystyle\end{eqnarray}$$

(2.22) $$\begin{eqnarray}\displaystyle & \displaystyle [(1-{\it\phi}){\it\tau}_{ij}n_{j}-(1-{\it\phi})p{\it\delta}_{ij}-{\it\phi}p_{w}{\it\delta}_{ij}]^{+}n_{j}=({\it\tau}_{ij}-p{\it\delta}_{ij})^{-}n_{j}, & \displaystyle\end{eqnarray}$$

(2.23) $$\begin{eqnarray}\displaystyle & \displaystyle \left[{\it\rho}_{w}L{\it\phi}(\boldsymbol{u}-\boldsymbol{v})+{\it\rho}_{w}L\boldsymbol{j}-K\boldsymbol{{\rm\nabla}}T\right]^{+}\boldsymbol{\cdot }\boldsymbol{n}=-K\boldsymbol{{\rm\nabla}}T^{-}\boldsymbol{\cdot }\boldsymbol{n}. & \displaystyle\end{eqnarray}$$

The last condition is a Stefan condition: the flow of latent heat into the boundary from within the temperate subdomain gives rise to a potential discontinuity in conductive heat flux across it.

2.4 Non-dimensionalization and approximation

A complete model for the temperate ice is given by (2.8)–(2.10), (2.15)–(2.17), together with (2.18)–(2.19). For cold ice, (2.1)–(2.4) replace these equations (in fact, the temperate ice model reduces to (2.1)–(2.4) at zero porosity), while the conditions (2.20)–(2.23) apply at any cold–temperate boundary.

We now make two particular approximations that provide a simpler but still comprehensive model for the temperate ice. In making these approximations it is helpful to non-dimensionalize the model, but before losing the dimensions, we highlight the approximations to be made. The first is a Boussinesq approximation; that is, the density difference will be ignored except in the buoyancy term ( $({\it\rho}_{w}-{\it\rho})$ in (2.17)). Secondly, we will assume that the relative moisture flux is significantly less than the ice flux, so that the porosities remain small and the ice matrix is still approximately incompressible, $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}\approx 0$ in (2.9). In addition, we will assume below that $T_{m}$ is constant. Although this restriction could easily be relaxed to include pressure dependence, it helps to clarify the discussion.

To non-dimensionalize the variables, we write

(2.24a-d ) $$\begin{eqnarray}\displaystyle x=[x]x^{\ast },\quad \boldsymbol{u}=[u]\boldsymbol{u}^{\ast },\quad t=([x]/[u])t^{\ast },\quad a=[a]a^{\ast }, & & \displaystyle\end{eqnarray}$$

(2.25a-e ) $$\begin{eqnarray}\displaystyle T=T_{ref}+[T]T^{\ast },\quad {\it\phi}={\it\varepsilon}{\it\phi}^{\ast },\quad \boldsymbol{j}=[j]\boldsymbol{j}^{\ast },\quad p_{e}=[p_{e}]p_{e}^{\ast },\quad p_{r}=[p_{r}]p_{r}^{\ast }, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $T_{ref}$ is a reference temperature included mostly for future convenience for efforts to try to generalize the model to a pressure-dependent melting point; otherwise, we can simply set $T_{ref}=T_{m}$ . We choose natural scales for the temperature, porosity, the compaction pressure and the reduced pressure as

(2.26a-e ) $$\begin{eqnarray}\displaystyle [T]=\frac{[a][x]^{2}}{K},\quad {\it\varepsilon}=\frac{{\it\rho}c[T]}{{\it\rho}_{w}L},\quad [j]=\frac{[a][x]}{{\it\rho}_{w}L},\quad [p_{e}]=\frac{[{\it\eta}][a]}{{\it\varepsilon}{\it\rho}_{w}L},\quad [p_{r}]=\frac{[{\it\eta}][u]}{[x]}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $[{\it\eta}]$ is a typical value for the effective viscosity ( $[{\it\eta}]=(1/2)A^{-1/n}([u]/[x])^{1/n-1}$ ), and values for the typical dissipation rate $[a]$ and for $[u]$ and $[x]$ are assumed to be given by the geometry, driving stresses and sliding behaviour relevant to the ice flow at hand.

Table 1.

Typical values of parameters and scales for the variables. $[{\it\phi}]={\it\varepsilon}$ is a ‘scale’ for the porosities, which are by definition dimensionless. There is considerable uncertainty over the permeability constants $k_{0}$ and ${\it\alpha}$ , and consequently ${\it\kappa}$ .

As an example, we estimate values for these scalings in table 1, using a typical driving stress of approximately $10^{5}~\text{Pa}$ and strain rate $[D]\approx 0.1~\text{m}~\text{y}^{-1}$ to infer order-of-magnitude estimates of the dissipation rate $[a]$ and viscosity $[{\it\eta}]=(1/2)A^{-1/n}[D]^{1/n-1}$ . A vertical length scale of $100~\text{m}$ and velocity scale of $1~\text{m}~\text{y}^{-1}$ are used in anticipation of our one-dimensional test cases in § 4.

The resulting dimensionless parameters are

(2.27a-e ) $$\begin{eqnarray}\displaystyle r=\frac{{\it\rho}}{{\it\rho}_{w}},\quad {\it\varepsilon}=\frac{{\it\rho}c[T]}{{\it\rho}_{w}L},\quad \mathit{Pe}\,=\frac{{\it\rho}c[u][x]}{K},\quad {\it\kappa}=\frac{k_{0}{\it\varepsilon}^{{\it\alpha}}({\it\rho}_{w}-{\it\rho})g}{[j]},\quad {\it\delta}=\frac{[p_{e}]}{({\it\rho}_{w}-{\it\rho})g[x]}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

which represent the density ratio, the porosity scale (also in this case an inverse Stefan number), the Péclet number (the ratio of advection to conduction), a permeability parameter ${\it\kappa}$ (a dimensionless measure of both permeability and gravity) and a compaction parameter ${\it\delta}$ , which measures how much compaction pressure gradients contribute to relative moisture flux.

We drop the asterisk decorations immediately with the understanding that all variables are dimensionless. The scaled versions of force balance (2.18) and constitutive law (2.19) (recall that we can reduce these to (2.2)–(2.3) in the cold domain by setting ${\it\phi}=D_{kk}=0$ ) are

(2.28) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial x_{j}}[(1-{\it\varepsilon}{\it\phi}){\it\tau}_{ij}]+\frac{{\it\varepsilon}}{\mathit{Pe}\,}\frac{\partial }{\partial x_{i}}({\it\phi}p_{e})-\frac{\partial p_{r}}{\partial x_{i}}=\frac{{\it\varepsilon}}{\mathit{Pe}\,{\it\delta}}{\it\phi}g_{i}, & \displaystyle\end{eqnarray}$$

(2.29) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\tau}_{ij}=A(T,{\it\phi})D^{1/n-1}\left(\unicode[STIX]{x1D60B}_{ij}-{\textstyle \frac{1}{3}}D_{kk}{\it\delta}_{ij}\right), & \displaystyle\end{eqnarray}$$

where $A$ is now dimensionless.

In the cold ice domain, mass balance and conservation of energy are

(2.30) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(2.31) $$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Pe}\,\left(\frac{\partial T}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}T\right)-{\rm\nabla}^{2}T=a. & \displaystyle\end{eqnarray}$$

In the temperate domain,

(2.32) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=-\frac{{\it\varepsilon}}{\mathit{Pe}\,}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{j}-\frac{{\it\varepsilon}(1-r)}{r\mathit{Pe}\,}a, & \displaystyle\end{eqnarray}$$

(2.33) $$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Pe}\,\left(\frac{\partial {\it\phi}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}+{\it\phi}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}\right)+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{j}=a, & \displaystyle\end{eqnarray}$$

(2.34) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{j}=\frac{{\it\phi}p_{e}}{{\it\eta}}-\frac{1-r}{r\mathit{Pe}\,}a, & \displaystyle\end{eqnarray}$$

(2.35) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{j}={\it\kappa}{\it\phi}^{{\it\alpha}}\left(\boldsymbol{g}+{\it\delta}\boldsymbol{{\rm\nabla}}p_{e}+\mathit{Pe}\,{\it\delta}\boldsymbol{{\rm\nabla}}p_{r}\right), & \displaystyle\end{eqnarray}$$

${\it\eta}=(1/2)AD^{1/n-1}$ being a dimensionless viscosity. These represent respectively, overall mass conservation, moisture conservation, the compaction relation and Darcy’s law. The non-dimensional gravity is $\boldsymbol{g}=(0,0,-1)$ , assuming the axes are oriented with the $z$ -direction pointing vertically upwards. The jump conditions at the boundaries between cold ( $^{-}$ ) and temperate ( $^{+}$ ) domains become

(2.36) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{j}^{+}\boldsymbol{\cdot }\boldsymbol{n}=-(1-r)\mathit{Pe}\,[{\it\phi}(\boldsymbol{u}-\boldsymbol{v})]^{+}\boldsymbol{\cdot }\boldsymbol{n}, & \displaystyle\end{eqnarray}$$

(2.37) $$\begin{eqnarray}\displaystyle & \displaystyle \left[(1-{\it\varepsilon}{\it\phi}){\it\tau}_{ij}+\frac{{\it\varepsilon}}{\mathit{Pe}\,}{\it\phi}p_{e}{\it\delta}_{ij}-p_{r}{\it\delta}_{ij}\right]^{+}n_{j}=({\it\tau}_{ij}-p{\it\delta}_{ij})^{-}n_{j}, & \displaystyle\end{eqnarray}$$

(2.38) $$\begin{eqnarray}\displaystyle & \displaystyle r\mathit{Pe}\,[{\it\phi}(\boldsymbol{u}-\boldsymbol{v})]^{+}\boldsymbol{\cdot }\boldsymbol{n}=-\boldsymbol{{\rm\nabla}}T^{-}\boldsymbol{\cdot }\boldsymbol{n}, & \displaystyle\end{eqnarray}$$

with $\boldsymbol{u}$ and $T$ still continuous, so $T=T_{m}$ . The first of these jump conditions indicates that if a finite quantity of water is carried into the boundary (and therefore freezes there and expands), mass conservation requires that there must be a compensating moisture flux away from the boundary. The second is a statement of force balance while the third is the Stefan condition, which relates freezing at the boundary to the rate of heat transfer away from it.

To complete our approximation, we now take the limits ${\it\varepsilon}\rightarrow 0$ , and $r\rightarrow 1$ . ${\it\varepsilon}\ll 1$ corresponds to the assumption of a low moisture content, which ensures that the flow of the ice matrix remains divergence free. This follows from (2.32), assuming $\mathit{Pe}\,$ does not simultaneously go to zero; so treating the ice as incompressible in this way makes sense provided the ice moves predominantly due to background flow of the ice sheet, rather than through compaction alone (the latter would be an unusual circumstance). Setting $r\ll 1$ implies that the effects of contraction and expansion upon melting and freezing are ignored in so far as they induce moisture movement. The effect of the density difference causing pore water to drain under gravity is still included in (2.35); this is the Boussinesq approximation.

The limit of small ${\it\delta}$ is also appropriate, as table 1 indicates that ${\it\varepsilon}$ and ${\it\delta}$ are typically of similar size. We choose to retain the term ${\it\delta}\boldsymbol{{\rm\nabla}}p_{e}$ in the definition of $\boldsymbol{j}$ for now, as it leads to the formation of boundary layers, as discussed in § 4 and appendix D. These boundary layers involve large gradients in $p_{e}$ , but not gradients in $p_{r}$ , as the latter is prescribed from the ice flow solution. We can therefore legitimately omit the term $\mathit{Pe}\,{\it\delta}\boldsymbol{{\rm\nabla}}p_{r}$ in (2.35), provided the Péclet number is not too large. Specifically, we shall assume that $1\gg \mathit{Pe}\,{\it\delta}\gg {\it\epsilon}$ , which allows us also to omit the right-hand side of (2.28).

2.5 Summary

Here we summarize the reduced model. At leading order the equations for ice flow, to be solved in both cold and temperate domains, are

(2.39) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(2.40) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\tau}_{ij}}{\partial x_{j}}-\frac{\partial p_{r}}{\partial x_{i}}=0, & \displaystyle\end{eqnarray}$$

(2.41a-d ) $$\begin{eqnarray}\displaystyle {\it\tau}_{ij}=2{\it\eta}\unicode[STIX]{x1D60B}_{ij},\quad {\it\eta}=\frac{1}{2}A(T,{\it\phi})^{-1/n}D^{1/n-1},\quad \unicode[STIX]{x1D60B}_{ij}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial u_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right),\quad D=\sqrt{\unicode[STIX]{x1D60B}_{ij}\unicode[STIX]{x1D60B}_{ij}/2}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The energy equation for cold ice is

(2.42a-c ) $$\begin{eqnarray}\displaystyle \mathit{Pe}\,\left(\frac{\partial T}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}T\right)-{\rm\nabla}^{2}T=a,\quad {\it\phi}=0,\quad T\leqslant T_{m}, & & \displaystyle\end{eqnarray}$$

where $a={\it\tau}_{ij}\unicode[STIX]{x1D60B}_{ij}$ and the equivalent equations for temperate ice are

(2.43a-c ) $$\begin{eqnarray}\displaystyle \mathit{Pe}\,\left(\frac{\partial {\it\phi}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}\right)+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{j}=a,\quad T=T_{m},\quad {\it\phi}>0, & & \displaystyle\end{eqnarray}$$

(2.44) $$\begin{eqnarray}\displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{j}=\frac{{\it\phi}p_{e}}{{\it\eta}}, & & \displaystyle\end{eqnarray}$$

(2.45) $$\begin{eqnarray}\displaystyle \boldsymbol{j}={\it\kappa}{\it\phi}^{{\it\alpha}}\left(\boldsymbol{g}+{\it\delta}\boldsymbol{{\rm\nabla}}p_{e}\right). & & \displaystyle\end{eqnarray}$$

The conditions to apply at the cold–temperate boundaries are

(2.46) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{j}^{+}\boldsymbol{\cdot }\boldsymbol{n}=0, & \displaystyle\end{eqnarray}$$

(2.47) $$\begin{eqnarray}\displaystyle & \displaystyle \mathit{Pe}\,[{\it\phi}(\boldsymbol{u}-\boldsymbol{v})]^{+}\boldsymbol{\cdot }\boldsymbol{n}=-\boldsymbol{{\rm\nabla}}T^{-}\boldsymbol{\cdot }\boldsymbol{n}, & \displaystyle\end{eqnarray}$$

(2.48) $$\begin{eqnarray}\displaystyle & \displaystyle T^{-}=T_{m}. & \displaystyle\end{eqnarray}$$

Note that the temperature and porosity equations in the cold and temperate regions can be combined as an enthalpy equation in the usual way. The dimensionless enthalpy is defined as $h=T+{\it\phi}$ , and conversely temperature and porosity are derived from this as

(2.49) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}ll@{}}T=h,{\it\phi}=0, & \text{if}~h<T_{m},\\ T=T_{m},{\it\phi}=h-T_{m}, & \text{if}~h\geqslant T_{m}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Then (2.42) and (2.43) are equivalent to

(2.50a,b ) $$\begin{eqnarray}\displaystyle \mathit{Pe}\,\left(\frac{\partial h}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}h\right)+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{Q}=a,\quad \boldsymbol{Q}=\left\{\begin{array}{@{}ll@{}}-\boldsymbol{{\rm\nabla}}T,\quad & h<T_{m},\\ \boldsymbol{j},\quad & h\geqslant T_{m}.\end{array}\right. & & \displaystyle\end{eqnarray}$$

Enthalpy is transported by advection and conduction in the cold ice and by advection and moisture transport in the temperate ice. The Stefan condition (2.47) is the natural jump condition that arises from treating this equation in weak form (e.g. Ockendon et al. Reference Ockendon, Howison, Lacey and Movchan2003, § 2.5).

3.1 Interior boundaries

The interior boundaries between cold and temperate regions are free boundaries, whose position must be determined as part of the problem. Their location is determined from the condition $T=T_{m}$ , along with the jump condition (2.47). The manner in which these conditions operate depends upon the direction of ice flow relative to the boundary.

In particular, since $T\leqslant T_{m}$ in the cold region and $\boldsymbol{n}$ is defined to point into the temperate region, we always have $\boldsymbol{{\rm\nabla}}T^{-}\boldsymbol{\cdot }\boldsymbol{n}\geqslant 0$ , and (2.47) therefore implies that ${\it\phi}(\boldsymbol{u}-\boldsymbol{v})^{+}\boldsymbol{\cdot }\boldsymbol{n}\leqslant 0$ . If ice flows from a cold to a temperate region, so $(\boldsymbol{u}-\boldsymbol{v})\boldsymbol{\cdot }\boldsymbol{n}>0$ (we refer to this an ‘inflow’ boundary), then this is only possible if ${\it\phi}^{+}=0$ . Consequently (2.47) implies that in fact $\boldsymbol{{\rm\nabla}}T^{-}\boldsymbol{\cdot }\boldsymbol{n}=0$ . Thus, at any inflow boundary, (2.47) actually requires ${\it\phi}=\boldsymbol{{\rm\nabla}}T\boldsymbol{\cdot }\boldsymbol{n}=0$ . The condition ${\it\phi}=0$ seems intuitive, meaning that ice entering the temperate region starts with no water content. From the perspective of the heat transport problem on the cold side, the boundary is both at the melting temperature $T=T_{m}$ and has zero diffusive heat flux $-\boldsymbol{{\rm\nabla}}T\boldsymbol{\cdot }\boldsymbol{n}$ across it.

These two conditions are equivalent to the conditions applied for a parabolic obstacle problem, with $T\leqslant T_{m}$ defining the obstacle; this is the approach to defining the cold–temperate boundary implicit in the models of Zwinger et al. (Reference Zwinger, Greve, Gagliardini, Shiraiwa and Lyly2007) and Suckale et al. (Reference Suckale, Platt, Perol and Rice2014). With these boundary conditions, the temperature problem decouples from the physics of the temperate region ${\it\Omega}^{+}$ .

However, at an outflow boundary (where $(\boldsymbol{u}-\boldsymbol{v})\boldsymbol{\cdot }\boldsymbol{n}\leqslant 0$ ), unlike an inflow boundary, the boundary conditions on temperature do not decouple from those on the temperate ice side. There can be a discontinuity in ${\it\phi}$ , and the jump condition (2.47) relates this discontinuity to the jump in temperature gradient, in the usual manner of a Stefan condition. This reflects the fact that freezing can occur at an outflow boundary. From the point of view of determining the location of the boundary, the porosity in the temperate region needs to be solved for primarily because it features in this condition. The obstacle problem formulation mentioned above is only appropriate when there are no outflow boundaries.

The jump condition for mass conservation (2.46) provides the necessary boundary condition for the moisture transport problem (2.44)–(2.45) which, for a given instantaneous porosity field ${\it\phi}$ , constitutes an elliptic (Helmholtz-like) problem for $p_{e}$ , explicitly

(3.1) $$\begin{eqnarray}\displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left[{\it\kappa}{\it\phi}^{{\it\alpha}}\left(g+{\it\delta}\boldsymbol{{\rm\nabla}}p_{e}\right)\right]-\frac{{\it\phi}p_{e}}{{\it\eta}}=0. & & \displaystyle\end{eqnarray}$$

The equivalent conditions for this problem when the temperate domain extends to the exterior boundary will be discussed below. However, it is worth pointing out that the elliptic nature of this problem means that the divergence term $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{j}$ in (2.43) acts as a non-local function of the porosity field. The remaining hyperbolic part of that equation suggests that a boundary condition is only required for ${\it\phi}$ on the ‘inflow’ parts of the boundary, which is precisely where the above discussion around (2.47) revealed that ${\it\phi}=0$ .

3.2 Exterior boundaries

The conditions to apply at external boundaries of the domain depend on whether the ice adjacent to the boundary is temperate or cold. In some cases, this distinction may be prescribed as part of the boundary condition (e.g. if the temperature is prescribed to be below the melting point), whilst in others (e.g. if a given heat flux is supplied) it needs to be determined as part of the problem through the imposition of the constraint $T\leqslant T_{m}$ at the boundary, and may change in time.

At cold exterior boundaries, a condition on temperature or temperature gradient (or a relation between them) must be prescribed, in the usual way. Naturally, no conditions are required for the moisture transport problem there.

At temperate exterior boundaries, the appropriate conditions again depend on whether they are ‘inflow’ or ‘outflow’ boundaries, with $(\boldsymbol{u}-\boldsymbol{v})\boldsymbol{\cdot }\boldsymbol{n}$ being positive or negative respectively. Here the normal $\boldsymbol{n}$ points into the temperate region, and $\boldsymbol{v}$ is the velocity of the boundary itself (see figure 1). In the case of an inflow boundary, the porosity ${\it\phi}$ must be prescribed, as well as a condition on either the moisture flux $\boldsymbol{j}\boldsymbol{\cdot }\boldsymbol{n}$ or the compaction pressure $p_{e}$ for the moisture transport problem (2.34)–(2.35). In the case of an outflow boundary, only a condition on $\boldsymbol{j}\boldsymbol{\cdot }\boldsymbol{n}$ or $p_{e}$ is required.

Figure 1.

Temperate and cold domains, ${\it\Omega}^{+}$ and ${\it\Omega}^{-}$ , respectively, separated by interior boundaries ${\it\Gamma}$ with normal $\boldsymbol{n}$ . Thick lines correspond to exterior boundaries.

For example, at the surface of a temperate glacier, the accumulation region corresponds to an inflow boundary, and we may expect to prescribe the porosity and the moisture flux $\boldsymbol{j}\boldsymbol{\cdot }\boldsymbol{n}$ , possibly as the output of a model for the compacting firn layer above. In the ablation region, where the surface corresponds to an outflow boundary, the compaction pressure $p_{e}$ is likely to be zero, with both the pore pressure and ice pressure being atmospheric. (There is in fact a possible complication that arises here: the pore pressure $p_{w}$ cannot drop below atmospheric pressure or, in isolation from the atmosphere, below the triple point pressure of water. This bound may be attained near the glacier surface and lead to a region of partially saturated temperate ice overlying saturated temperate ice, and hence the introduction of another free boundary. The appropriate extension of our model to this case is likely to be similar to an analogous situation in subglacial drainage (Schoof, Hewitt & Werder Reference Schoof, Hewitt and Werder2012), but we do not pursue this further here.)

At the base of a temperate glacier, where basal melting occurs, we have an outflow boundary. A single boundary condition is required. The obvious choice is to equate the englacial water pressure $p_{w}$ with subglacial water pressure. If we define effective pressure at the bed $N$ as normal stress minus subglacial water pressure, then in dimensionless terms

(3.2) $$\begin{eqnarray}\displaystyle p_{e}=N+\mathit{Pe}\,{\it\tau}_{nn}, & & \displaystyle\end{eqnarray}$$

where ${\it\tau}_{nn}$ is the dimensionless deviatoric normal stress at the bed (the difference between $p_{e}$ and $N$ arises because the definition of compaction pressure does not take account of deviatoric stresses or the normal direction to the bed).

If basal freeze-on occurs (for example due to decompression and supercooling of subglacial water), the boundary would be of inflow type, and we additionally expect to apply a condition on the porosity of the frozen-on ice. The porosity condition would presumably depend on precisely how the frozen-on ice is formed; if frazil ice is formed in the subglacial water and accumulates mechanically at the underside of the glacier, a non-zero porosity may result, while direct freeze-on to the base of the glacier may result in a zero porosity boundary condition. Again, we do not investigate the physics involved further here.