2.1. A Stokes variational inequality

The finite element method can be used to compute approximate solutions of the Stokes equations with the contact boundary conditions (2.7a–c). As a first step, we must find a weak formulation by multiplying (2.1) by suitable test functions and integrating by parts. Weak formulations of partial differential equations are formulated in terms of weak derivatives, and the solutions are generally sought in Sobolev spaces (Adams & Fournier Reference Adams and Fournier2003). In order to keep the notation simpler, we assume in this subsection and in § 2.2 that we enforce the Neumann boundary condition (2.9) instead of (2.8). In this case, a suitable space of test functions is the vector-valued Sobolev space $V = \boldsymbol {W}^{1,\kappa }(\varOmega )$ for the velocity, and the Lebesgue space $Q = L^{\kappa '}(\varOmega )$ for the pressure. Here, $\kappa = 1 + 1/n$, where $n$ is the coefficient in Glen's law, and $\kappa ' = \kappa /(\kappa -1)$ is the Hölder conjugate of $\kappa$.

Due to the contact boundary conditions (2.7a–c), the associated weak formulation of the Stokes equations is a variational inequality. Given the convex set of admissible velocity fields

(2.11)\begin{equation} K = \left\lbrace \boldsymbol{v}\in V : \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{n} \leq 0 \ \text{on}\ \varGamma_a \right\rbrace, \end{equation}

the weak formulation of (2.1) is: find $(\boldsymbol {u},p)\in K\times Q$ such that

(2.12a)$$\begin{gather} a(\boldsymbol{u},\boldsymbol{v} - \boldsymbol{u}) - b(p, \boldsymbol{v} - \boldsymbol{u}) \geq f(\boldsymbol{v} - \boldsymbol{u}) \quad \forall\ \boldsymbol{v}\in K, \end{gather}$$

(2.12b)$$\begin{gather}b(q,\boldsymbol{u}) = 0 \quad \forall\ q\in Q. \end{gather}$$

Here, the operators $a : V\times V\to \mathbb {R}$ and $b:Q\times V\to \mathbb {R}$ are defined by

(2.13a,b)\begin{equation} a(\boldsymbol{u},\boldsymbol{v}) = \int_\varOmega 2 \eta\,(|\varepsilon(\boldsymbol{u})|)\,\varepsilon(\boldsymbol{u}) : \varepsilon(\boldsymbol{v})\,\mathrm{d} x \quad \text{and} \quad b(q,\boldsymbol{v}) = \int_\varOmega q (\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v})\,\mathrm{d} x, \end{equation}

and the linear functional $f:V\to \mathbb {R}$ is defined by

(2.14)\begin{equation} f(\boldsymbol{v}) = \int_{\varGamma_t} \tau_b (\boldsymbol{v}\boldsymbol{\cdot}(1,0)^{{\text{T}}})\,\mathrm{d} s - \int_{\varGamma_t} p_i (\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{n})\,\mathrm{d} s - \int_{\varGamma_a\cup\varGamma_d} p_w (\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{n})\,\mathrm{d} s. \end{equation}

As we assume that the water pressure $p_w$ is constant along the bedrock, the incompressibility of the flow velocity can be used to rewrite the right-hand side as

(2.15)\begin{equation} f(\boldsymbol{v}) = \int_{\varGamma_t} \tau_b (\boldsymbol{v}\boldsymbol{\cdot}(1,0)^{{\text{T}}})\,\mathrm{d} s - \int_{\varGamma_t} N (\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{n})\,\mathrm{d} s, \end{equation}

where $N = p_i - p_w$ is the effective pressure. Consequently, the solution depends only on $N$, not on the particular values of $p_i$ and $p_w$.

The variational inequality (2.12) is equivalent to the Stokes equations (2.1) in the sense that if the velocity $\boldsymbol {u}$ is at least twice differentiable and the pressure $p$ is differentiable, then integration by parts can be performed in order to arrive at (2.1) with the boundary conditions introduced in the previous section. Of the three conditions in (2.7a–c), only the kinematic condition $\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {n} \leq 0$ is enforced explicitly in the definition of $K$; the remaining conditions are enforced implicitly in (2.12).

There exist many different approaches for solving a variational inequality with the finite element method (Glowinski, Lions & Trémoliéres Reference Glowinski, Lions and Trémoliéres1981). For example, in the penalty method, the constraint $\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {n} \leq 0$ in the definition of $K$ is enforced via the addition of a nonlinear penalty term in (2.12a), and the variational inequality is transformed into a variational equality over $\boldsymbol {W}^{1,\kappa }(\varOmega )$. This is the approach used in Stubblefield et al. (Reference Stubblefield, Spiegelman and Creyts2021). We find that what works best in combination with the advection equation (2.10) is a mixed method with Lagrange multipliers that enforces an average-wise contact condition. We remark that such a condition could also be enforced via a penalty term, but in this case the contact boundary conditions would not be enforced exactly.

The contact boundary conditions (2.7a–c) can be enforced explicitly in (2.12) with a Lagrange multiplier that essentially represents the sum of the normal stress and the water pressure, $\sigma _{nn} + p_w$, along the attached region $\varGamma _a$. In order to introduce this reformulation of (2.12), we need to define a suitable set of admissible Lagrange multipliers. We define the range of the normal trace operator onto $\varGamma _a$ by

(2.16)\begin{equation} \varSigma = \left\lbrace \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{n} |_{\varGamma_a} : \boldsymbol{v}\in V \right\rbrace. \end{equation}

The space of admissible Lagrange multipliers is then the convex cone of elements in the dual space $\varSigma '$ satisfying a weak equivalent of the dynamic contact boundary condition $\sigma _{nn} + p_w \leq 0$ in (2.7a–c):

(2.17)\begin{equation} \varLambda = \left\lbrace \mu\in\varSigma' : \mu(\xi) \geq 0 \ \forall \xi\in\varSigma\ \text{such that}\ \xi\leq 0\ \text{on}\ \varGamma_a \right\rbrace. \end{equation}

The mixed formulation of the variational inequality (2.12) may be written as: find $(\boldsymbol {u},p,\lambda )\in V\times Q\times \varSigma '$ such that

(2.18a)$$\begin{gather} a(\boldsymbol{u},\boldsymbol{v}) - b(p, \boldsymbol{v} ) - \lambda(\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{n}|_{\varGamma_a}) = f(\boldsymbol{v}) \quad \forall\ \boldsymbol{v}\in V, \end{gather}$$

(2.18b)$$\begin{gather}b(q,\boldsymbol{u}) = 0 \quad \forall\ q\in Q, \end{gather}$$

(2.18c)$$\begin{gather}\mu(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{n}|_{\varGamma_a}) \geq 0 \quad \forall \mu\in\varLambda,\quad \lambda\in\varLambda \quad \text{and} \quad \lambda(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{n}|_{\varGamma_a}) = 0. \end{gather}$$

Equation (2.18a) can be obtained by multiplying (2.1a) by a test function $\boldsymbol {v}$, integrating by parts, and setting

(2.19)\begin{equation} \lambda = \sigma_{nn} + p_w. \end{equation}

Formulation (2.18) offers several advantages for finite element approximation in comparison to (2.12). First, there is no longer any need to build an approximation of the convex set $K$ because we seek the velocity in the space $V$. Second, the conditions (2.18c), which are a weak analogue of the contact boundary conditions (2.7a–c), can be rewritten as a nonlinear equation at the discrete level; see (2.22c) below. This means that the contact boundary conditions can be enforced exactly by solving the finite element system analogous to (2.18). As we explain shortly, this leads to a robust and stable numerical method for solving the complete subglacial cavitation system.

Under certain conditions on the geometry of the domain and the linear functional $f$, the variational inequality (2.12) has a unique solution. A proof can be found in de Diego, Farrell & Hewitt (Reference de Diego, Farrell and Hewitt2021), together with an analysis of the finite element approximation of the mixed system (2.18).

2.2. A numerical method for subglacial cavitation

In this subsection, we present a numerical algorithm for solving the complete subglacial cavitation problem that results from coupling the Stokes equations (2.1) to the time-dependent advection equation (2.10). We write the discrete counterparts to these equations and explain how they are coupled by deforming the domain according to a contact criterion. The resulting algorithm is summarised in Algorithm 1. The solver is implemented in Firedrake (Rathgeber et al. Reference Rathgeber, Ham, Mitchell, Lange, Luporini, McRae, Bercea, Markall and Kelly2016), using the version available at zenodo/Firedrake-20211103.0 (2021). The code for the viscous contact solver presented in Algorithm 1 is openly available, as described in the data availability statement below.

Subglacial cavitation is a time-dependent problem, so its discretisation requires a partition of a given time interval $[0,T]$ into intervals of length ${\rm \Delta} t$. At each time step $t_k = k {\rm \Delta} t$ for $k = 0,1,\ldots, N_T = T/{\rm \Delta} t$, we must consider a discrete cavity roof $\theta _h^{k}$ and domain $\varOmega ^{k}$ given by

(2.20)\begin{equation} \varOmega^{k} = \lbrace (x,y)\in\mathbb{R}^{2} : 0\leq x \leq L,\ \theta_h^{k}(x) \leq y \leq H \rbrace. \end{equation}

Additionally, the lower boundary is the union of the attached and detached regions $\varGamma _a^{k}$ and $\varGamma _d^{k}$, respectively. These are also time-dependent and are determined at each time step by the contact criterion (2.29) below.

Algorithm 1:

Solution procedure for the subglacial cavity problem

2.2.1. Finite element approximation of the Stokes variational inequality

On each domain $\varOmega ^{k}$, we seek a finite element approximation to (2.18). To this end, we consider a triangulation of $\varOmega ^{k}$ and the finite element spaces $V_h$, $Q_h$ and $\varSigma _h$ given by piecewise continuous quadratic vector fields, piecewise constant scalar fields and piecewise constant scalar functions on $\varGamma _a^{k}$, respectively. This velocity–pressure pair is known to be inf-sup stable in two dimensions (Boffi, Brezzi & Fortin Reference Boffi, Brezzi and Fortin2013). Furthermore, as we explain in § 2.2.2, working with a piecewise constant Lagrange multiplier enables the use of a very simple upwinding scheme for solving the advection equation (2.10) in a manner consistent with the discrete contact boundary conditions.

We write $V_h = \operatorname {span}{\{\boldsymbol {v}_i\}}_{i=1}^{N_v}$, $Q_h = \operatorname {span}{\{q_j\}}_{j=1}^{N_q}$ and $\varSigma _h = \operatorname {span}{\{\mu _k\}}_{k=1}^{N_\mu }$, where $N_v = \dim {V_h}$, $N_q = \dim {Q_h}$ and $N_\mu = \dim {\varSigma _h}$. For the functions $(\boldsymbol {u}_h, p_h, \lambda _h) \in V_h\times Q_h\times \varSigma _h$, we denote the respective degrees of freedom in $\mathbb {R}^{N_v}$, $\mathbb {R}^{N_q}$ and $\mathbb {R}^{N_\mu }$ by $\mathbf{u}$, $\mathbf{p}$ and $\boldsymbol{\lambda }$. In order to write an algebraic counterpart of (2.18c), we need the operator

(2.21)\begin{equation} \boldsymbol{\gamma_n} : \mathbb{R}^{N_v} \to \mathbb{R}^{N_\mu} \end{equation}

that returns the average normal components of a vector $\boldsymbol {v}_h\in V_h$ along the edges on $\varGamma _a^{k}$. The discrete counterpart to (2.18) is

(2.22a)$$\begin{gather} \boldsymbol{A}(\mathbf{u}) - \boldsymbol{B} \mathbf{p} - \boldsymbol{D}\boldsymbol{\lambda} = \boldsymbol{f}, \end{gather}$$

(2.22b)$$\begin{gather}-\boldsymbol{B}^{{\text{T}}}\mathbf{u} = 0, \end{gather}$$

(2.22c)$$\begin{gather}\boldsymbol{\lambda} + \boldsymbol{C}(\boldsymbol{\lambda},\mathbf{u}) = 0 . \end{gather}$$

Here, we have introduced the matrices $\boldsymbol {B}\in \mathbb {R}^{N_v\times N_q}$ and $\boldsymbol {D}\in \mathbb {R}^{N_v\times N_\mu }$, the vector $\boldsymbol {f}\in \mathbb {R}^{N_v}$, and the nonlinear operators $\boldsymbol {A}:\mathbb {R}^{N_v}\to \mathbb {R}^{N_v}$ and $\boldsymbol {C}:\mathbb {R}^{N_\mu }\times \mathbb {R}^{N_v}\to \mathbb {R}^{N_\mu }$. The matrices are given by the elements $\boldsymbol {B}_{ij} = b(\boldsymbol {v}_i,q_j)$ and $\boldsymbol {D}_{ij} = \int _{\varGamma ^{k}_a} \mu _j (\boldsymbol {v}_i \boldsymbol {\cdot } \boldsymbol {n})\,\mathrm {d} s$, and the vector $\boldsymbol {f}$ is given by

(2.23)\begin{equation} \boldsymbol{f}_i ={-} \int_{\varGamma_t} p_i(\boldsymbol{v}_i\boldsymbol{\cdot}\boldsymbol{n})\,\mathrm{d} s - \int_{\varGamma_a\cup\varGamma_d} p_w(\boldsymbol{v}_i\boldsymbol{\cdot}\boldsymbol{n})\,\mathrm{d} s + \int_{\varGamma_t} \tau_b (\boldsymbol{v}_i\boldsymbol{\cdot}(1,0)^{{\text{T}}})\,\mathrm{d} s. \end{equation}

The nonlinear operator $\boldsymbol {A}$ is given by

(2.24)\begin{equation} [\boldsymbol{A}(\mathbf{u})]_i = \int_\varOmega 2 \eta_\epsilon(|\varepsilon(\boldsymbol{u}_h)|)\, \varepsilon(\boldsymbol{u}_h) :\varepsilon(\boldsymbol{v}_i)\,\mathrm{d} x, \end{equation}

where $\eta _\epsilon$ is a regularised form of (2.2) that prevents infinite viscosity at zero strain rates, defined by

(2.25)\begin{equation} \eta_\epsilon(|\varepsilon(\boldsymbol{u})|) = \tfrac{1}{2}A^{{-}1/n} \left( \tfrac{1}{2} |\varepsilon(\boldsymbol{u})|^{2} + \epsilon^{2} \right)^{{(1-n)}/{2n}}, \end{equation}

for a regularisation parameter $\epsilon > 0$ that we set to $\epsilon = 10^{-2}$. The nonlinear function $\boldsymbol {C}$ is

(2.26)\begin{equation} \boldsymbol{C}(\boldsymbol{\lambda},\mathbf{u}) = \max{\left\lbrace 0, -\boldsymbol{\lambda}+ c (\boldsymbol{\gamma_n}\mathbf{u}) \right\rbrace} \end{equation}

for an arbitrary $c>0$ that we set to $c = 1$. We find the computations carried out in this paper to be insensitive to the value of $c > 0$ when $c\in [10^{-6},10^{6}]$. In (2.26), the $\max$ operation is understood to be carried out for each of the elements in the vector $-\boldsymbol{\lambda } + c (\boldsymbol {\gamma _n}\mathbf{u})\in \mathbb {R}^{N_\mu }$. The use of (2.26) is a common way of expressing contact conditions; a particular advantage is that the nonlinear system (2.22) can be solved with a semi-smooth Newton method (Hintermüller, Ito & Kunisch Reference Hintermüller, Ito and Kunisch2002). Moreover, (2.22c) is equivalent to enforcing

(2.27a–c)\begin{equation} \boldsymbol{\gamma_n} \mathbf{u} \leq 0,\quad \boldsymbol{\lambda} \leq 0 \quad \text{and} \quad (\boldsymbol{\gamma_n} \mathbf{u})\boldsymbol{\cdot} \boldsymbol{\lambda} = 0, \end{equation}

a discrete equivalent of (2.7a–c). By solving (2.27a–c), we see that on each edge on $\varGamma _a^{k}$, we have that the edge wishes either to detach ($\boldsymbol {\gamma _n} \mathbf{u} < 0$) or to remain attached to the bed ($\boldsymbol{\lambda } \leq 0$ and $\boldsymbol {\gamma _n} \mathbf{u} = 0$). At the continuous level, one cannot write (2.18c) as in (2.22c) because in general the Lagrange multiplier $\lambda$ is only an element of the dual space $\varSigma '$ and may not even be a function on $\varGamma _a$ (Stadler Reference Stadler2007).

2.2.2. The discrete advection equation and a contact criterion

In order to write a discretisation of (2.10), we introduce some further notation. We denote the points along the lower boundary by $(x_i,\theta _i)$, with the index $i$ increasing from left to right as in figure 2. The edge between $(x_{i-1},\theta _{i-1})$ and $(x_i,\theta _i)$ is denoted by $e_i$. We write $u^{k}_{n,i}$ for the value of $\boldsymbol {\gamma _n}\mathbf{u}$ associated with $e_i$ at time $t_k$ (by the definition of $\boldsymbol {\gamma _n}$, this value corresponds to the average value of $\boldsymbol {u}_h\boldsymbol {\cdot }\boldsymbol {n}$ along $e_i$).

Figure 2.

Notation and ordering for nodes $\{(x_i,\theta _i)\}$ and edges $\{e_i\}$ along the attached region $\varGamma _a^{k}$.

A numerical scheme that solves the advection equation (2.10a) should satisfy two conditions. First, it must be stable, as described in standard textbooks (LeVeque Reference LeVeque2007) (in Appendix A.2 we present results with and without stabilisation of (2.10a)). Additionally, the scheme should be compatible with the discrete contact boundary conditions (2.27a–c), in the sense that if for two edges $e_i$ and $e_{i+1}$ we have that $u_{n,i}^{k} = 0$ and $u_{n,i+1}^{k} = 0$, then we should have that $\theta ^{k+1}_{i} = \theta ^{k}_{i}$; that is, the node in between the edges remains unchanged. This means that the discrete counterpart of $\boldsymbol {u}\boldsymbol {\cdot } (\partial \theta /\partial x, -1)^{{\text{T}} }$ should be defined in terms of $\boldsymbol {\gamma _n} \mathbf{u}$. If not, then we generally find that all edges will detach within a few time steps due to approximation errors. Here, we propose the following explicit scheme:

(2.28)\begin{equation} \frac{\theta_i^{k+1} - \theta_i^{k}}{{\rm \Delta} t} = \left(\left(\frac{\theta^{k}_i - \theta^{k}_{i-1}}{x_i - x_{i-1}}\right)^{2} + 1 \right)^{1/2} u^{k}_{n,i}, \end{equation}

where the velocity at each node is taken from the edge located immediately upstream. This results in an upwinding scheme with the property that penetration cannot occur along $\varGamma _a$, because the average values of $\boldsymbol {u}_h\boldsymbol {\cdot }\boldsymbol {n}$ along the edges are non-positive; see (2.27a–c). However, a node on a previously detached edge can penetrate the bedrock. To avoid this and enforce (2.10b), we simply set $\theta ^{k+1}_i = b(x_i)$ whenever $\theta ^{k+1}_i < b(x_i)$.

Given a cavity roof $\theta _h^{k}$ at $t_k$, we must set the attached and detached regions $\varGamma _a^{k}$ and $\varGamma _d^{k}$ to solve the Stokes variational inequality at $t_k$. This is carried out by assigning each edge $e_i$ with endpoints $(x_{i-1},\theta ^{k}_{i-1})$ and $(x_i,\theta ^{k}_i)$ to either $\varGamma _a^{k}$ or $\varGamma _d^{k}$. We do this by looking at the node downstream of the edge, in accordance with the scheme for the advection equation:

(2.29)\begin{equation} e_i\ \text{belongs to}\ \varGamma_d^{k}\ \text{if and only if}\ \theta^{k}_i - b(x_i)> \mathrm{tol}, \end{equation}

where the tolerance $\mathrm {tol}$ is set to $10^{-9}$. Consequently, given a cavity roof $\theta ^{k}_h$ at time $t_k$, the attached and detached regions $\varGamma ^{k}_a$ and $\varGamma ^{k}_d$ are fully determined.

Once we have computed a new cavity roof $\theta ^{k+1}_h$, the Stokes problem is then solved over a new domain $\varOmega ^{k+1}$. This means that the mesh is deformed according to $\theta ^{k+1}_h$. Since we use vertically extruded meshes for the computations in this paper, a simple algorithm is used that moves each node vertically in a linear manner with respect to the original position of the nodes at $t = 0$. If we denote the mesh nodes over $x_i$ by $y_{ij}$, and let $y = H$ be the upper boundary of the domain, then we compute the new vertical positions of the nodes via

(2.30)\begin{equation} y_{ij}^{k+1} = \theta^{k+1}_i + \frac{H - \theta^{k+1}_i}{H - \theta^{0}_i} ( y^{0}_{ij} - \theta^{0}_i ) \end{equation}

for given $y^{0}_{ij}$, $\theta ^{0}_i$ and $\theta ^{k+1}_i$.

3.1. Steady subglacial cavities

The steady state of (2.1) and (2.10) with the boundary conditions $u = u_i$ and $\sigma _{nn} = -N$ on the top boundary $\varGamma _t$ can be found with Algorithm 1 by evolving the cavity from an initial state until the norm of the discrete derivative $(\theta ^{k+1} - \theta ^{k})/{\rm \Delta} t$ is below a prescribed threshold $10^{-4}$. In this subsection, we find the steady cavity for a Newtonian flow ($n = 1$ in Glen's law (2.2)) over a sinusoidal bed of small amplitude $r = 0.01$ given by $b(x) = rL \cos {(2{\rm \pi} x/L)}$. The non-dimensional effective pressure is set to $N = 0.3$, and the non-dimensional horizontal velocity to $u_i = 1$. We first consider a linear rheology and a bedrock of small amplitude in order to compare our results with the analytical solution of the linearised cavitation problem considered in Fowler (Reference Fowler1986) and Schoof (Reference Schoof2005). A brief description of this method is included in Appendix B.

We use five different meshes with $n_e$ vertices distributed uniformly along the lower boundary. In table 1, we present the non-dimensional basal shear stress $\tau _b$ and sliding speed $u_b$ along the cavities. These values are calculated with

(3.1a,b)\begin{equation} \tau_b ={-} \frac{1}{L} \int_{\varGamma_a\cup\varGamma_d}(\sigma_{nn} + p_w)n_x\,\mathrm{d} s \quad \text{and}\quad u_b = \frac{1}{L}\int_{\varGamma_a\cup\varGamma_d}u\,\mathrm{d} s, \end{equation}

where $\boldsymbol {n} = (n_x,n_y)$. The expression for $\tau _b$ can be derived from the expression for the force exerted by the ice on the bed as in Schoof (Reference Schoof2005). Recall that the Lagrange multiplier $\lambda$ represents the normal stress $\sigma _{nn} + p_w$ along $\varGamma _a$ (see (2.19)); therefore $\tau _b$ can be calculated from $\lambda$ via

(3.2)\begin{equation} \tau_b ={-} \frac{1}{L} \int_{\varGamma_a}\lambda n_x\, \mathrm{d} s. \end{equation}

As we show in the next subsection, the values of $\tau _b$ and $u_b$ can be used to construct a sliding law.

Table 1.

Information about the different meshes used to compute the steady cavity states together with calculations of the basal shear stress $\tau _b$, the basal sliding speed $u_b$, and the detachment and reattachment points of the cavity.

The formula for $u_b$ presented in (3.1b) might seem strange if the subglacial cavity domain is interpreted as a boundary layer between an ice sheet and the bedrock. This would suggest that we take the sliding speed to be the average value of the horizontal velocity along the top boundary $\varGamma _t$. However, our computations indicate that if the height of the domain $H$ is sufficiently large, then we can expect the shear stress to approach a constant value and the horizontal velocity to vary with $y^{n}$ as $y$ approaches $H$, where $n$ here is the exponent in Glen's law (2.2). Therefore, the horizontal velocity along the top boundary depends strongly on the height of the domain. For this reason, and following Gagliardini et al. (Reference Gagliardini, Cohen, Raback and Zwinger2007), we use (3.1b) to calculate $u_b$. In this case, we find that $u_b$ is independent of $H$ for sufficiently large values of $H$. In particular, throughout this paper we set $H=L$. In agreement with Gagliardini et al. (Reference Gagliardini, Cohen, Raback and Zwinger2007), we find this value of $H$ to be sufficiently large.

In figure 3, we present the steady cavity shape and normal stresses $\sigma _{nn}$ along the attached region $\varGamma _a$ for the different meshes. We can see from these figures that the cavity shape is computed accurately even with the coarsest mesh. Additionally, we also present the stress distribution obtained from the linearised theory, which is determined uniquely for an effective pressure $N$ and a sliding speed $u_b$. The result from the linearised theory plotted in figure 3 is computed with the value of $u_b$ calculated with the most refined mesh.

Figure 3.

(a) Steady cavity shapes and (b) normal stresses along the attached region computed with different mesh sizes. The stress distribution obtained from the linearised theory is shown in black. Here, the amplitude of the obstacle is set to $r = 0.01$, the scaled effective pressure to $N = 0.3$, and the scaled horizontal velocity at the top to $u_i = 1$.

The plot for the normal stress distribution demonstrates that the contact conditions (2.7a–c) are satisfied exactly at the discrete level for all of the meshes, because $\sigma _{nn} + p_w \leq 0$. This plot also exhibits the singularity of the normal stresses at the reattachment point. This singularity complicates the approximation of the normal stresses along the attached region and can lead to very inaccurate computations of the sliding law in largely cavitated states. However, figure 3 also indicates that with increasing mesh refinement, the solver appears to converge towards the linearised solution.

At the discrete level, the Lagrange multiplier $\lambda$, which we use to calculate $\tau _b$, is piecewise constant on each edge along the lower boundary. In figure 3, these values, the degrees of freedom of $\lambda$, are plotted at the midpoints of each edge. Due to the contact criterion (2.29), the edge immediately upstream of the first reattached node is treated as part of $\varGamma _a$ (this is required in order to allow for the possibility of subsequent detachment there), which explains why there are non-zero values of $\sigma _{nn} + p_w$ left of the reattachment node in figure 3. This is particularly visible for the coarsest mesh with $n_e = 16$.

3.2. Computation of the linear and nonlinear steady sliding law

We next perform calculations similar to those of § 3.1 but for varying effective pressure $N$ and power-law exponent $n$. This allows us to map out a steady sliding law for ice sliding over a hard bed with cavitation as in Gagliardini et al. (Reference Gagliardini, Cohen, Raback and Zwinger2007). The steady sliding law can be mapped out by varying only $N$ and keeping $u_i$ fixed, since dimensional analysis of the steady problem shows that the scaled basal shear stress $\tau _b/N$ depends only on the ratio $u_b/(ALN^{n})$, and not independently on $u_b$ or $N$ (Fowler Reference Fowler1986) (the same will not be true of the unsteady problem in § 4).

Several previous studies have suggested what form the sliding law should take, both with and without cavitation (Kamb Reference Kamb1970; Fowler Reference Fowler1981, Reference Fowler1986; Gudmundsson Reference Gudmundsson1997; Schoof Reference Schoof2005; Gagliardini et al. Reference Gagliardini, Cohen, Raback and Zwinger2007). The law proposed in Gudmundsson (Reference Gudmundsson1997) for the uncavitated case can be written as

(3.3)\begin{equation} \left(\frac{\tau_b}{rN}\right)^{n} = \alpha(n)\,\frac{r}{AL}\,\frac{u_b}{N^{n}}, \end{equation}

where $r = a/L$, and $\alpha (n)$ is a function depending on $n$. The function $\alpha (n)$ is related to the parameter $c_0$ (which also depends on $n$) considered in Gudmundsson (Reference Gudmundsson1997) by

(3.4)\begin{equation} \alpha(n) = \frac{(2{\rm \pi})^{n+2}}{2c_0}. \end{equation}

For a Newtonian flow, the complex analysis method presented in Fowler (Reference Fowler1986) and Schoof (Reference Schoof2005) yields an exact solution to the linearised problem. In particular, for high effective pressures, no cavitation occurs and a linear sliding law as in (3.3) with $c_0 = 1$ is found. For effective pressures lower than a critical value $8{\rm \pi} ^{2} r \eta u_b / L$, cavitation occurs and the sliding law becomes nonlinear, varying with $N$ as well as $u_b$.

We compute the sliding law over a sinusoidal bed of different amplitudes $r$ and for $n=1$, 3 and 5, and we plot the results in figure 4 using the scaling suggested by (3.3). The mesh has 192 vertices along the lower boundary, and 7296 cells. The locations of the cavity endpoints are also plotted in figure 4, along with the solution to the linearised problem calculated with the method from Fowler (Reference Fowler1986) and Schoof (Reference Schoof2005). For each $n$, the parameter $\alpha (n)$ is computed by calculating the slope of the curve near the origin (where there is no cavitation) for the lowest value of $r$. The corresponding values of the parameter $c_0$ can then be calculated from (3.4); these values can be found in table 2 together with those obtained in Gudmundsson (Reference Gudmundsson1997) and Gagliardini et al. (Reference Gagliardini, Cohen, Raback and Zwinger2007). We see that the results obtained in these works are broadly similar to ours and that the value $c_0 = 1$ obtained in the linearised theory is approached in all cases when $n = 1$.

Figure 4.

(a) Computed sliding law for steady glacial sliding with cavitation. (b) Cavity endpoints. The parameter $n$ in Glen's law (2.2) is set to $n = 1$ (green), 3 (blue) and 5 (red). For each $n$, we compute the sliding law for beds of amplitude $r = 0.01$ (circles, light), 0.04 (squares) and 0.08 (triangles, dark). For these computations, we set the scaled velocity at the top of the domain to $u_i = 1$ and the scaled viscosity parameter to $A = 0.5$. We use a mesh with 192 cells along the lower boundary. The parameter $\alpha (n)$ is computed from the slope of the curve near the origin for the lowest $r$.

Table 2.

Value of the parameter $c_0$ associated with the sliding laws. This parameter is computed with the slope of these curves near the origin.

The computed sliding laws with cavitation in figure 4 are multivalued for $\tau _b/(rN)$ as expected (Fowler Reference Fowler1986; Schoof Reference Schoof2005). This aspect of the law justifies the use of the Dirichlet boundary condition $u = u_i$ instead of the alternative Neumann boundary condition. We find that if we use the Neumann boundary condition $\boldsymbol {\sigma }_{nt} = \tau _b(1,0)^{{\text{T}} }$ and initiate the cavity from a fully attached state, then the solver always evolves to the steady state associated with the upsloping region of the curve (see also § 4). From figure 4, we also deduce that the validity of the sliding law (3.3) along the linear segment of the curves (where either little or no cavitation has occurred) decreases with increasing values of $r$ and $n$. For example, when $n = 5$, one can observe that the linear segment of the sliding law for $r = 0.08$ clearly does not collapse onto the corresponding linear segment for $r = 0.01$. As soon as the cavity size increases and the sliding laws cease to be linear, the aspects of these curves largely differ for different values of $n$.

In figure 4, we use a scaling different to the one used in Gagliardini et al. (Reference Gagliardini, Cohen, Raback and Zwinger2007). In Gagliardini et al. (Reference Gagliardini, Cohen, Raback and Zwinger2007), the computed maximum value reached by $\tau _b/N$ is included in the scaling for the sliding law. In this way, the maximum value reached by the scaled sliding law equals 1 by design. However, we preferred the scaling based on (3.3) because it contains fewer terms that are unknown a priori. It is also worth mentioning that for different values of $n$, the curves in figure 4 do not collapse into a single curve when plotted with the scaling from Gagliardini et al. (Reference Gagliardini, Cohen, Raback and Zwinger2007).

For the linear case with $n = 1$, the numerical results computed with the finite element solver highly resemble those obtained with the linearised solution. For $r=0.08$, a slight difference from the linearised solution can be seen near the peak of the sliding law. This difference is probably a consequence of nonlinear effects that are accentuated with increasing amplitudes of bedrock roughness.