2.1 Glen–Stokes model

The Stokes model consists of the momentum conservation equation when acceleration terms are ignored, together with the incompressibility condition:

where $\unicode[STIX]{x1D748}$ is the Cauchy stress tensor, $\boldsymbol{g}=(0,0,-g)$ , $g$ is the gravitational constant and $\boldsymbol{u}=(u_{x},u_{y},u_{z})$ is the velocity field. Call $\unicode[STIX]{x1D70F}$ the deviatoric stress tensor defined by

where $I$ is the identity tensor and $P$ is the pressure field, and by the requirement that $\text{tr}(\unicode[STIX]{x1D70F})=0$ so that $P=-(1/3)\text{tr}(\unicode[STIX]{x1D748})$ . Glen’s flow law (Glen Reference Glen1953), which describes the mechanical behaviour of ice, consists of the following nonlinear relation:

where $\unicode[STIX]{x1D63F}(\boldsymbol{u})$ denotes the strain rate tensor defined by

$\unicode[STIX]{x1D707}$ is the viscosity defined by

where $|Y|:=\sqrt{(Y:Y)}$ denotes the norm associated with the scalar product $(\;:\;)$ , $B>0$ and $p\in (1,2]$ are two parameters. Observe that equation (2.2) implies that $\unicode[STIX]{x1D63F}(\boldsymbol{u})$ has zero trace.

It is common in glaciology to express the coefficient $B$ and exponent $p$ in terms of different parameters. Define Glen’s exponent $n=1/(p-1)$ and then define a new rate constant $A=(\sqrt{2})^{n-1}B^{-n}$ . (Note that the common convention in the literature is to instead set $A=B^{-n}$ (e.g. Blatter Reference Blatter1995), but our convention reduces clutter in the equations.) In reality, $A$ and $B$ are not constant since they depend on the ice temperature (Glen Reference Glen1953; Paterson Reference Paterson1994), but, for the sake of simplicity, this paper assumes that the ice is isothermal.

2.2 Boundary conditions

The boundary conditions that supplement (2.1), (2.2) are the following. No force applies on the ice–air interface,

where $\boldsymbol{n}$ is an outer normal vector along $\unicode[STIX]{x1D6E4}_{s}$ . Along the lower surface interface, the no-slip condition is

and the nonlinear friction condition reads (Hutter Reference Hutter1983; Schoof & Hewitt Reference Schoof and Hewitt2013)

for $k\in \{x,y\}$ , where $s\in (1,p]$ , $C=C(x,y)>0$ and $\boldsymbol{n}$ is the outward normal unit vector to $\unicode[STIX]{x1D6E4}_{m}$ . Relation (2.10) relates the basal shear stress $[(I-\boldsymbol{n}\boldsymbol{n}^{\text{T}})\unicode[STIX]{x1D70F}]\boldsymbol{\cdot }\boldsymbol{n}$ to the sliding velocity $(I-\boldsymbol{n}\boldsymbol{n}^{\text{T}})\boldsymbol{\cdot }\boldsymbol{u}$ , both of them being projected onto the tangential plane.

2.3 Functional space and variational formulation

We now reformulate the Stokes problem into variational and minimization problems (Glowinski & Le Tallec Reference Glowinski and Le Tallec1989; Jouvet & Rappaz Reference Jouvet and Rappaz2011). For this, we consider the following divergence-free velocity space incorporating Dirichlet conditions where present (Girault & Raviart Reference Girault and Raviart1986):

see Adams & Fournier (Reference Adams and Fournier2003) for an introduction to the Sobolev space $W^{1,p}$ . In addition, we introduce several norms: $\Vert \cdot \Vert _{L^{p}}$ and $|\cdot |_{W^{1,p}}$ denote the $L^{p}$ norm and the $W^{1,p}$ seminorm, defined by $\Vert \boldsymbol{v}\Vert _{L^{p}}^{p}:=\int _{V}|\boldsymbol{v}|^{p}\,\text{d}V$ and $|\boldsymbol{v}|_{W^{1,p}}=\Vert \unicode[STIX]{x1D735}\boldsymbol{v}\Vert _{L^{p}}$ respectively. Furthermore, $\Vert \cdot \Vert _{W^{1,p}}$ denotes the $W^{1,p}$ norm defined by $\Vert \boldsymbol{v}\Vert _{W^{1,p}}=\Vert \boldsymbol{v}\Vert _{L^{p}}+|\boldsymbol{v}|_{W^{1,p}}$ .

From now on, we assume that $\unicode[STIX]{x1D6E4}_{0}$ has a non-zero measure, so that Poincaré inequality (Brezis Reference Brezis1999; Adams & Fournier Reference Adams and Fournier2003) ensures the equivalence of the norms $\Vert \cdot \Vert _{W^{1,p}}$ and $|\cdot |_{W^{1,p}}$ on the functional space ${\mathcal{X}}$ while Korn’s inequality (e.g. Quarteroni & Valli Reference Quarteroni and Valli1994; Ern & Guermond Reference Ern and Guermond2004) implies the equivalence of the norms $\Vert \cdot \Vert _{W^{1,p}}$ and the one defined by $\Vert \unicode[STIX]{x1D63F}(\cdot )\Vert _{L^{p}}$ on the space ${\mathcal{X}}$ . (One can get rid of this assumption. However, one then has to replace ${\mathcal{X}}$ by its quotient set induced by rigid body motions (Schoof Reference Schoof2010), otherwise Poincaré inequality is no longer valid. Doing so makes the analysis substantially more complex.) The Poincaré and Korn inequalities can be summarized by the following lemma, whose proof is found in the above references.

We now formulate the weak form of the Glen–Stokes model. After using (2.3), multiplying (2.1) by a test function $\boldsymbol{v}\in {\mathcal{X}}$ , integrating over $V$ , and using Gauß’s theorem twice, we obtain

First, by (2.4) and (2.6), the first term of (2.14) is

Second, using (2.7a,b ), (2.8), (2.9), (2.10) and $\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{n}=0$ , the second term of (2.14) becomes

after noting that $\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{n}=0$ implies $([(\boldsymbol{n}\boldsymbol{n}^{\text{T}})\unicode[STIX]{x1D70F}]\boldsymbol{\cdot }\boldsymbol{n})\boldsymbol{\cdot }\boldsymbol{v}=0$ and $[(\boldsymbol{n}\boldsymbol{n}^{\text{T}})\boldsymbol{\cdot }\boldsymbol{u}]\boldsymbol{\cdot }\boldsymbol{v}=0$ . Finally, since $\boldsymbol{v}$ satisfies the constraints of the functional space ${\mathcal{X}}$ (see (2.11)), the third and fourth terms of (2.14) are cancelled out such that the pressure field $P$ vanishes from the variational formulation, which simplifies to

The above variational formulation (without pressure) is usually called the ‘reduced’ formulation. In contrast, there exists a ‘mixed’ formulation, in which the pressure $P$ plays the role of a Lagrange multilayer for the constraint $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0$ ; see, e.g., Girault & Raviart (Reference Girault and Raviart1986). In this paper, we work only with the ‘reduced’ formulation.

2.4 Minimization formulation and well-posedness

The problem defined by (2.17) can be associated with the following minimization problem (Jouvet & Rappaz Reference Jouvet and Rappaz2011):

where

Indeed, it is easy to verify that (2.17) corresponds to $\langle {\mathcal{J}}^{\prime }(\boldsymbol{u}),\boldsymbol{v}\rangle =0$ , where $\langle {\mathcal{J}}^{\prime }(\boldsymbol{u}),\boldsymbol{v}\rangle$ is the Gâteaux derivative of ${\mathcal{J}}$ at $\boldsymbol{u}$ in the direction $\boldsymbol{v}$ . One can prove the strong continuity and the convexity of ${\mathcal{J}}$ in ${\mathcal{X}}$ directly from (2.19); see lemma 2.4 and lemma 2.5 in Jouvet & Rappaz (Reference Jouvet and Rappaz2011). On the other hand, one needs the Poincaré and Korn inequalities (2.12)–(2.13) in order to prove the strict convexity of ${\mathcal{J}}$ and the coercivity property,

where $c_{1},c_{2}$ are positive constants; see lemma 2.5 and lemma 2.6 in Jouvet & Rappaz (Reference Jouvet and Rappaz2011). Following arguments of convex analysis (Girault & Raviart Reference Girault and Raviart1986), this implies the existence and the uniqueness of a minimizer of ${\mathcal{J}}$ in ${\mathcal{X}}$ (Jouvet & Rappaz Reference Jouvet and Rappaz2011), and, equivalently, the well-posedness of the variational problem (2.17). Let us note that one can show the existence of a unique pressure field $P$ such that $(\boldsymbol{u},P)$ satisfies the ‘mixed’ formulation associated with (2.17), supposing an inf–sup condition for the functional spaces of velocity and pressure; see Jouvet & Rappaz (Reference Jouvet and Rappaz2011).

2.5 Partial elimination of the vertical velocity

We could continue working with the reduced weak formulation (2.17). However, it is advantageous to eliminate the vertical component of the velocity, $u_{z}$ , from the formulation whenever this is possible. A closer look at (2.17) shows that this is already possible for the two last terms.

From now on, we assume that

where $\boldsymbol{x}=(x,y)$ denotes the horizontal coordinates, $z$ denotes the vertical coordinate, and $\text{}\underline{s}(x,y)$ and $\overline{s}(x,y)$ are the elevations of the lower and upper ice surfaces; see figure 1. The ice thickness is denoted by $h$ and satisfies $h=\overline{s}-\text{}\underline{s}$ .

In fact, $\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}=0$ and $\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{n}=0$ on $\unicode[STIX]{x1D6E4}_{m}$ imply

where $\unicode[STIX]{x1D735}_{\boldsymbol{x}}$ is the gradient with respect to the horizontal coordinates $\boldsymbol{x}=(x,y)$ and $\boldsymbol{n}^{\text{T}}=((\unicode[STIX]{x1D735}_{\boldsymbol{x}}\text{}\underline{s})^{\text{T}},-1)$ . This defines a scalar product $(\boldsymbol{a},\boldsymbol{b})_{M}:=(\unicode[STIX]{x1D648}\boldsymbol{a})\boldsymbol{\cdot }\boldsymbol{b}$ and its associated norm $|\boldsymbol{a}|_{M}:=\sqrt{(\boldsymbol{a},\boldsymbol{ a})_{M}}$ , where the $3\times 3$ matrix $\unicode[STIX]{x1D648}$ is defined by

so that

Since the last row of $\unicode[STIX]{x1D648}$ in (2.23) is filled with zeroes, (2.24) shows that the second term in (2.17) does not actually depend on the last component $u_{z}$ .

On the other hand, using the Gauß theorem we have, for any $\boldsymbol{v}\in {\mathcal{X}}$ (see (2.11)),

As a consequence, the right-hand side in (2.17) can be formulated without the vertical component of the velocity:

In summary, the variational formulation (2.17) can now be rewritten as

while the functional to be minimized is

It must be stressed that only the first term still depends on the vertical velocity in both formulations (2.27) and (2.28).

2.6 Shallow ice flow models

In this section, we infer the classical shallow ice flow models (Blatter, SIA and SSA) by applying simplifications to the Glen–Stokes model at the variational level.

Let us define $\unicode[STIX]{x1D63F}_{0}(\boldsymbol{u})$ , $\unicode[STIX]{x1D63F}_{1}(\boldsymbol{u})$ and $\unicode[STIX]{x1D63F}_{2}(\boldsymbol{u})$ by splitting the strain rate tensor as follows:

where

where $\unicode[STIX]{x1D63F}_{0}$ and $\unicode[STIX]{x1D63F}_{1}$ depend only on the horizontal velocity components $(u_{x},u_{y})$ while $\unicode[STIX]{x1D63F}_{2}$ depends on the horizontal derivatives of the vertical component. We also define $\unicode[STIX]{x1D63F}_{01}(\boldsymbol{u}):=\unicode[STIX]{x1D63F}_{0}(\boldsymbol{u})+\unicode[STIX]{x1D63F}_{1}(\boldsymbol{u})$ , and $\unicode[STIX]{x1D63F}_{02}(\boldsymbol{u}):=\unicode[STIX]{x1D63F}_{0}(\boldsymbol{u})+\unicode[STIX]{x1D63F}_{2}(\boldsymbol{u})$ .

A hierarchy of ice sheet models is associated with an expansion in the small parameter $\unicode[STIX]{x1D716}=[h]/[\boldsymbol{x}]$ , the aspect ratio of the ice geometry $V$ . Here, $[h]$ and $[\boldsymbol{x}]$ denote its typical height and length. It is easy to verify that in the splitting (2.29), $\unicode[STIX]{x1D63F}_{0}$ scales with $1/\unicode[STIX]{x1D716}$ , $\unicode[STIX]{x1D63F}_{1}$ scales with $1$ (i.e. it is dimensionless) and $\unicode[STIX]{x1D63F}_{2}$ scales with $\unicode[STIX]{x1D716}$ . From this observation, one can obtain submodels by neglecting the second-order component $\unicode[STIX]{x1D63F}_{2}$ or/and the first-order component $\unicode[STIX]{x1D63F}_{1}$ in (2.29). Figure 2, and the paragraphs below, describe this hierarchy by identifying replacements for the full tensor $\unicode[STIX]{x1D63F}(\boldsymbol{u})$ . Looking forward, the estimators in § 3 are in one-to-one correspondence with the arrows in figure 2.

Figure 2.

The hierarchy of the ice flow models.

The Blatter model (Blatter Reference Blatter1995) is obtained by neglecting $\unicode[STIX]{x1D63F}_{2}(\boldsymbol{u})$ , or equivalently by replacing $\unicode[STIX]{x1D63F}(\boldsymbol{u})$ by $\unicode[STIX]{x1D63F}_{01}(\boldsymbol{u})$ in the formulation (2.27) or (2.28). As a consequence, the variational formulation is

where $\boldsymbol{u}_{1}$ denotes the Blatter solution. By doing so, the variational equality no longer involves any vertical component of the velocity. As a consequence, the Blatter model solves the horizontal components of the velocity only, and $u_{z}$ can be computed a posteriori from the horizontal components by integrating the incompressibility condition,

for $k\in \{x,y\}$ .

In solving (2.33), one can remove the third component from the definition of the functional space ${\mathcal{X}}$ together with the incompressibility condition. Let $\hat{{\mathcal{X}}}$ be the space ${\mathcal{X}}$ but with the vertical velocity component removed. One can then show that any solution $\boldsymbol{u}_{1}\in {\mathcal{X}}$ of (2.33) minimizes (2.28), with $\unicode[STIX]{x1D63F}$ being replaced by $\unicode[STIX]{x1D63F}_{01}$ , and vice versa with the solution of (2.33) from $\hat{{\mathcal{X}}}$ being extended by solving (2.34). (For simplicity, we keep the notation ${\mathcal{X}}$ to denote $\hat{{\mathcal{X}}}$ in what follows.) In addition, the functional conserves all of its properties such that the existence and uniqueness of a minimizer are automatically ensured; see Colinge & Rappaz (Reference Colinge and Rappaz1999) and Schoof (Reference Schoof2010). The traditional Blatter model derivation involves the ‘hydrostatic approximation’ assumption. It is used to get rid of the pressure and to transfer the gravitational forces on the right-hand side of the momentum conservation equations from the vertical to the horizontal components. In the current approach, the ‘hydrostatic approximation’ assumption is not needed since the pressure $P$ is eliminated by restricting the functional space for the velocity field to divergence-free fields (§ 2.3), while this transfer of the gravitational force components results from the application of the Gauß theorem (§ 2.5), as already mentioned in Dukowicz (Reference Dukowicz2012). Thus, working with the weak form of the Stokes equation instead of the strong original form makes the derivation simpler. Although one can rewrite (2.33) in a simplified subspace of ${\mathcal{X}}$ , we deliberately keep the formulation intact in order to obtain a unified formulation for all submodels.

Going one step further, the SIA (Hutter Reference Hutter1983) is obtained by neglecting $\unicode[STIX]{x1D63F}_{1}$ , or equivalently by replacing $\unicode[STIX]{x1D63F}_{01}$ by $\unicode[STIX]{x1D63F}_{0}$ in (2.33). As a consequence, the variational formulation for the SIA is

where $\boldsymbol{u}_{0}$ denotes the SIA solution. As in the Blatter model, (2.35) does not involve $u_{z}$ . Furthermore, the first term in (2.35) now contains only vertical derivatives, but no more horizontal derivatives. The removal of these horizontal derivatives renders the overall problem substantially easier to solve. In particular, the SIA solution can be written analytically, as independent vertical integrals for the horizontal components of velocity; see, e.g., Greve & Blatter (Reference Greve and Blatter2009). Because of the loss of the $\unicode[STIX]{x1D63F}_{1}$ term in the variational inequality, one cannot state the existence and the uniqueness of a solution as before, since Korn’s inequality no longer holds with $\unicode[STIX]{x1D63F}_{0}$ instead of $\unicode[STIX]{x1D63F}$ . However, since the associated functional contains no horizontal derivatives, one can permute ‘ $\text{min}$ ’ and horizontal integration, and then solve the problem pointwise in $\boldsymbol{x}$ . Call $\unicode[STIX]{x1D6FA}$ and $\unicode[STIX]{x1D6FA}_{k}$ the projections of $V$ and its boundaries $\unicode[STIX]{x1D6E4}_{k}$ on the horizontal plane $(Oxy)$ for $k\in \{s,0,m\}$ . Given $\boldsymbol{x}\in \unicode[STIX]{x1D6FA}$ , the pointwise minimization problem

over the set of functions

admits a unique solution since the functional in (2.36) is continuous, strictly convex and coercive in the functional space (2.37); see Schoof (Reference Schoof2010).

Finally, the SSA (Morland Reference Morland, Veen and Oerlemans1987) is obtained from the Blatter model by replacing the space of solutions ${\mathcal{X}}$ by its subspace of vertically constant horizontal component functions

As a consequence, the variational formulation for the SSA is

where $\overline{\boldsymbol{u}}_{1}$ denotes the SSA solution. Again, (2.39) does not involve $u_{z}$ . Note that, if $\boldsymbol{v}\in \overline{{\mathcal{X}}}$ , then $\unicode[STIX]{x1D63F}_{01}(\boldsymbol{v})$ simplifies into $\unicode[STIX]{x1D63F}_{1}(\boldsymbol{v})$ in (2.39). However, we keep this form intact in the following. Since the space $\overline{{\mathcal{X}}}$ is a closed and convex subset of ${\mathcal{X}}$ , ${\mathcal{J}}$ admits a unique minimum in $\overline{{\mathcal{X}}}$ , $\overline{\boldsymbol{u}}_{1}$ the unique solution of (2.39); see also Schoof (Reference Schoof2006).

2.7 Stability with respect to data

The next lemma shows that all of the solutions are bounded by functions of the model data.

3.1 Preliminary lemmas

Before starting, we need to recall two lemmas that will be useful later. The first lemma consists of inequalities (see Lemma 2.1 in Barrett & Liu (Reference Barrett and Liu1993), slightly modified with the triangle inequality $|\unicode[STIX]{x1D709}|+|\unicode[STIX]{x1D702}|\leqslant 2(|\unicode[STIX]{x1D709}|+|\unicode[STIX]{x1D709}-\unicode[STIX]{x1D702}|)\leqslant 4(|\unicode[STIX]{x1D709}|+|\unicode[STIX]{x1D702}|)$ ).

The second lemma introduces the concept of solution-dependent quasinorm (see Barrett & Liu (Reference Barrett and Liu1993)).

The quasinorms (3.7) and (3.8) will be used to deal with the nonlinearities of the problem (when $p\neq 2$ and $s\neq 2$ ) to obtain the estimates (3.1)–(3.4). In contrast, the derivation gets simpler in the case $p=2$ and $s=2$ (Newtonian fluid and linear sliding law) since one can work with the standard $W^{1,2}$ -norm directly instead of any quasinorm so that the inequalities (3.5), (3.6) and (3.9) are actually not needed.

3.2 Proof of estimate (3.1) – Blatter–SSA

As seen earlier, the SSA solution is obtained by changing the space of ${\mathcal{X}}$ to an approximation subspace $\overline{{\mathcal{X}}}$ , and not by changing the variational formulation (as needed in order to obtain the Blatter or SIA solution). Thus, standard techniques can be applied to derive a posteriori error estimates in order to establish (3.1). Before doing so, we need some approximation properties of $\overline{{\mathcal{X}}}$ in ${\mathcal{X}}$ , as stated in the following lemma; see the proof in Jouvet & Picasso (Reference Jouvet and Picasso0000).

By definition of the quasinorm (3.8) and inequality (3.6), we have

Since $\boldsymbol{u}_{1}$ solves (2.33), then taking $\boldsymbol{v}=\boldsymbol{u}_{1}-\overline{\boldsymbol{u}}_{1}$ , we deduce

Since $\overline{\boldsymbol{u}}_{1}$ solves (2.39), where $\overline{\boldsymbol{v}}\in \overline{{\mathcal{X}}}$ is arbitrary, we obtain

which becomes after using Gauß’s theorem and the symmetry of $\unicode[STIX]{x1D63F}_{01}$

By applying Hölder’s inequality, we obtain

Using the approximation inequalities of $\overline{{\mathcal{X}}}$ in ${\mathcal{X}}$ , namely (3.11) and (3.12) with $\boldsymbol{v}=\boldsymbol{u}_{1}-\overline{\boldsymbol{u}}_{1}$ , we obtain

In addition, (2.40) and (3.9) lead to

The final result (3.1) follows from (3.19) and (3.20).

3.3 Proof of estimate (3.3) – Blatter–SIA

Taking the difference between (2.33) and (2.35), together with the test function $\boldsymbol{v}=\boldsymbol{u}_{1}-\boldsymbol{u}_{0}$ , we obtain

By adding and subtracting $\int _{V}B|\unicode[STIX]{x1D63F}_{01}(\boldsymbol{u}_{0})|^{p-2}(\unicode[STIX]{x1D63F}_{01}(\boldsymbol{u}_{0}):\unicode[STIX]{x1D63F}_{01}(\boldsymbol{u}_{1}-\boldsymbol{u}_{0}))\,\text{d}V$ , it follows (using $\unicode[STIX]{x1D63F}_{0}:\unicode[STIX]{x1D63F}_{0}=\unicode[STIX]{x1D63F}_{0}:\unicode[STIX]{x1D63F}_{01}$ ) that

From (3.5), (3.6), the definition (3.8), and using $\unicode[STIX]{x1D63F}_{01}(\boldsymbol{u}_{0})=\unicode[STIX]{x1D63F}_{0}(\boldsymbol{u}_{0})+\unicode[STIX]{x1D63F}_{1}(\boldsymbol{u}_{0})$ , we have

Using the Hölder inequality, we obtain

From (2.40) and (3.9), we have

We conclude (3.3) from (3.24) and (3.25).

3.4 Proof of estimate (3.4) – Stokes–Blatter

Taking the difference between (2.33) and (2.27), together with the test function $\boldsymbol{v}=\boldsymbol{u}_{1}-\boldsymbol{u}$ , we obtain

By adding and subtracting $\int _{V}B|\unicode[STIX]{x1D63F}(\boldsymbol{u}_{1})|^{p-2}(\unicode[STIX]{x1D63F}(\boldsymbol{u}_{1}):\unicode[STIX]{x1D63F}(\boldsymbol{u}_{1}-\boldsymbol{u}))\,\text{d}V$ , and using $\unicode[STIX]{x1D63F}=\unicode[STIX]{x1D63F}_{01}+\unicode[STIX]{x1D63F}_{2}$ , we get

where

From (3.6) and the definition (3.7), we obtain

Next, $A_{1}$ and $A_{2}$ are treated separately. In $A_{1}$ , one can simplify $(\unicode[STIX]{x1D63F}(\boldsymbol{u}_{1}):\unicode[STIX]{x1D63F}_{2}(\boldsymbol{u}_{1}-\boldsymbol{u}))$ into $(\unicode[STIX]{x1D63F}_{02}(\boldsymbol{u}_{1}):\unicode[STIX]{x1D63F}_{2}(\boldsymbol{u}_{1}-\boldsymbol{u}))$ since $\unicode[STIX]{x1D63F}_{1}:\unicode[STIX]{x1D63F}_{2}=0$ . We now integrate by parts on $\unicode[STIX]{x1D6FA}$ (the projection of $V$ on the horizontal plane $(Oxy)$ ), so that

and the boundary term $[(B|\unicode[STIX]{x1D63F}(\boldsymbol{u}_{1})|^{p-2}(\unicode[STIX]{x1D63F}_{02}(\boldsymbol{u}_{1}))\boldsymbol{\cdot }\boldsymbol{n}]\boldsymbol{\cdot }(\boldsymbol{u}_{1}-\boldsymbol{u})$ vanishes on $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}\times [\text{}\underline{s},\overline{s}]$ . It should be noted that only the vertical component $(\boldsymbol{u}_{1}-\boldsymbol{u})$ is involved in (3.31). Using the Hölder inequality,

is obtained. Using the Hölder inequality, it can be shown that

Since we assumed no basal sliding, $\boldsymbol{v}(\boldsymbol{\cdot },\text{}\underline{s}(\boldsymbol{\cdot }))=\mathbf{0}$ , so that

using the incompressibility. Combination of (3.31) and (3.34) leads to

On the other hand, from (3.5) and the Hölder inequality, we obtain

On combining (3.30), (3.35) and (3.36),

is obtained. From (2.40) and (3.9), we have

The result of (3.37) and (3.38) is (3.4).