We consider two important and widely studied problems in glaciology that involve contact. The first is that of the grounding line of an ice sheet flowing from the continent and into the sea, where the ice floats and loses contact with the bedrock. The stability of the position of this grounding line was first questioned by Weertman in 1974  and since then numerous analyses have attempted to prove or disprove the possibility of an instability, see e.g. [28, 38, 24]. The second problem is that of subglacial cavitation, where the ice detaches from the bedrock along the lee side of an obstacle. Subglacial cavitation occurs at a much smaller scale, along the interface between the ice and the bedrock, and is usually formulated as a boundary layer problem. Lliboutry first proposed the possibility of cavities forming between the ice and the bedrock in 1968 . Since then, subglacial cavitation has been recognised as a fundamental mechanism in glacial sliding, attracting the attention of both theoretical [17, 37] and experimental  studies.
A precise understanding of both grounding line dynamics and subglacial cavitation is of great relevance for predicting future sea level rise and therefore comprehending large scale climate dynamics [36, 20]. The two contact problems described above are modelled by coupling a Stokes problem for the ice flow with a time-dependent advection equation for the free surface. At each instant in time, a Stokes problem must be solved with contact boundary conditions that allow the detachment of the ice from the bed. These contact conditions transform the instantaneous Stokes problem into a variational inequality.
Numerous finite element simulations of these equations have been carried out, see e.g. [19, 13, 16, 27]. However, to the best of our knowledge, no formal analysis of this problem and its approximation exists in the mathematical literature. Moreover, we believe that the discretisations used in these computations can be improved upon, by exploiting the structure of the variational inequality. Although the Stokes variational inequality is superficially similar to the elastic contact problem, which has been widely studied [34, 26], the Stokes problem includes three substantial difficulties that must be addressed carefully: the presence of rigid body modes in the space of admissible velocities, the nonlinear rheological law used to model ice as a viscous fluid, and the nonlinearity of the boundary condition used for the sliding law when modelling the grounding line problem.
In this work we analyse the instantaneous Stokes variational inequality and its approximation and focus on the first two difficulties due to rigid body modes and the nonlinear rheology. The presence of rigid body modes renders this problem semicoercive. Although semicoercive variational inequalities have been studied in the past [25, 40, 2], existing analyses use purely indirect arguments which give very limited information on the effects on the discretisation of the finite element spaces used. Here, we present a novel approach based on the use of metric projections onto closed convex cones to obtain constructive proofs for the discretisation errors. On the other hand, the nonlinear rheology complicates the estimation of errors for the discrete problem. Here, we use the techniques from [5, 30] to establish a convergence analysis.
We propose a mixed formulation of the Stokes variational inequality where a Lagrange multiplier is used to enforce the contact conditions. This formulation permits a structure-preserving discretisation that explicitly enforces a discrete version of the contact conditions, up to rounding errors. This allows for a precise distinction between regions where the ice detaches from the bed and those where it remains attached. This precision is extremely useful when coupling the Stokes variational inequality with the time-dependent advection equation for the free surface.
1.1 Outline of the paper
In Section 2, the Stokes variational inequality and its mixed formulation are presented. We prove a Korn-type inequality involving a metric projection onto a cone of rigid modes that will be used throughout the analysis, and we demonstrate that the mixed formulation is well posed. In Section 3, we analyse a family of finite element approximations of the mixed problem and present error estimates in terms of best approximation results for the velocity, pressure and Lagrange multiplier. Finally, in Section 4, a concrete finite element scheme involving quadratic elements for the velocity and piecewise constant elements for the pressure and the Lagrange multiplier is introduced. We then present error estimates for this scheme and provide numerical results for two test cases. For the first test, we solve a problem with a manufactured solution to calculate convergence rates and compare these with our estimates. For the second test, we compute the evolution of a subglacial cavity to exhibit the benefits of using a mixed formulation in applications of glaciological interest.
Given two normed vector spacesand and a bounded linear operator , the dual of is denoted by and the dual operator to by . The range of is denoted by and its kernel by . The norm in is denoted by and the pairing between elements in the primal and dual spaces by for and . We will work with the Lebesgue and Sobolev spaces , where and , defined as the set of functions with weak derivatives up to order which are -integrable. When we write . The space of polynomials of degree over a simplex (interval, triangle, tetrahedron) is denoted by . The space of continuous functions over a domain is given by . Vector functions and vector function spaces will be denoted with bold symbols, e.g. and .
2 Formulation of the problem
In this section we introduce the semicoercive variational inequality that arises in glaciology and present its formulation as a mixed problem with a Lagrange multiplier. An auxiliary analytical result concerning a metric projection is then proved. Finally, we analyse the existence and uniqueness of solutions for the mixed problem.
2.1 A model for ice flow
We denote by a bounded, connected and polygonal domain which represents the glacier. The assumptions of the domain being two-dimensional and polygonal are made in order to simplify the analysis, but we expect the essential results presented here to extend to three dimensional domains with smooth enough boundaries. Ice is generally modelled as a viscous incompressible flow whose motion is described by the Stokes equation : equationparentequation
In the equations above, represents the ice velocity, the pressure and
is a prescribed body force due to gravitational forces. The tensoris the symmetric part of the velocity gradient, that is,
The coefficient is the effective viscosity of ice, which relates the stress and strain rates. A power law, usually called Glen’s law , is the most common choice of rheological law for ice:
Here, represents the Frobenius norm of a matrix: for with components we have . The field is a prescribed function for which . The parameter is constant and is usually set to ; for we recover the standard linear Stokes flow. From now on, we simply write
where is in and satisfies a.e. on for some . Moreover, is in for . This expression for reveals the -Stokes nature of the problem when considered as a variational problem in the setting of Sobolev spaces.
2.2 Boundary conditions
For a given velocity and pressure field, we define the stress tensor by
where is the identity tensor field. Let denote the unit outward-pointing normal vector to the boundary . We define the normal and tangential stresses at the boundary as
The boundary is partitioned into three disjoint open sets , and . The subset represents the part of the boundary in contact with the atmosphere (and the ocean or a water-filled cavity in the case of a marine ice sheet and a subglacial cavity, respectively). Here we enforce
where represents a prescribed surface traction force. On we enforce the contact conditions which allow the ice to detach from but not penetrate the bedrock. In particular, detachment can occur if the normal stress equals the subglacial water pressure, which is defined everywhere along a thin lubrication layer in between the ice and the bedrock. In order to write these in a simplified form, we assume that the water pressure along the bed is constant and we measure stresses relative to that water pressure. We also assume that the ice can slide freely and impose no tangential stresses. Then, the boundary conditions at are given by equationparentequation
Finally, represents a portion of the boundary on which the ice is frozen and hence we prescribe no slip boundary conditions:
As explained in Section 1, one of the challenges we wish to address with this work is the case when rigid body modes are present in the space of admissible velocities. In these cases, the problem becomes semicoercive instead of coercive and the structure of the problem changes significantly. This occurs whenever is empty. For this reason, we assume that can be empty and require the subsets and to have positive measures.
2.3 The mixed formulation
We now present the mixed formulation whose analysis and approximation is the focus of this work. To do so, we first write (1) with boundary conditions (4)-(6) as a variational inequality. Then, we introduce the mixed formulation by defining a Lagrange multiplier which enforces a constraint that arises due to the contact boundary conditions (5a). In Appendix A we specify and prove the sense in which these different formulations are equivalent.
We denote by the normal trace operator onto . This operator is built by extending to the operator on , defined on smooth functions. The closed convex subset of is then defined by
We also introduce the operators and defined by
Moreover, the action of the applied body and surface forces on the domain is expressed via the function , defined as
In the mixed formulation, the constraint on is enforced via a Lagrange multiplier. We denote the range of by and equip this space with the norm. We assume the geometry of and to be sufficiently regular for this space to be a Banach space, see [34, Section 5], [26, Chapter III] and [1, Chapter 7] for discussions on normal traces and trace spaces. The Lagrange multiplier is sought in the convex cone of multipliers
The equivalent mixed formulation of (9) is: find such that equationparentequation
2.4 A metric projection onto the cone of rigid body modes
In this section we present a projection operator onto the cone of rigid body modes inside and prove a Korn-type inequality involving this operator. This result, which we believe to be novel, allows us to prove the well-posedness of the continuous and discrete problems and to obtain estimates for the velocity error in the -norm.
We define the space of rigid body modes inside by
We also introduce the subspace of rigid body modes inside , which we denote by . Note that and for all . In fact, it can be shown that , see [34, Lemma 6.1], and therefore .
The fact that coincides with the kernel of complicates the construction of error estimates in the -norm. Our solution is to make use of the metric projection onto the closed convex cone , which we shall denote by . This metric projection assigns to each function a rigid body mode for which
In Appendix B we explain that metric projections are well-defined on uniformly convex Banach spaces and that the range of is the closed convex cone , see (58) for the definition of a polar cone. Since the Sobolev space is uniformly convex for , the operator has a closed convex range. This property is exploited below to prove Theorem 1.
We first prove a preliminary Korn inequality on closed convex sets, Lemma 1. For this preliminary result, we need the following generalised Korn inequality from : for a bounded and Lipschitz domain there is a constant such that
Let be a bounded and Lipschitz domain and a closed convex subset which satisfies . Then, there is a constant such that
By using (12), it suffices to show that
Assume by contradiction that (13) does not hold. Then, there is a sequence in such that and as . By (12), we see that is bounded in and therefore we may extract a subsequence, also denoted , which converges weakly to a in and strongly in . Since is closed and convex, it is also weakly closed by [7, Theorem 3.7], so . Moreover, by the lower semicontinuity of , we also have that and . However, by construction , a contradiction.
Let denote the metric projection onto the closed convex cone . There is a such that
2.5 Well posedness of the mixed formulation
Questions on the existence and uniqueness of solutions of the mixed system (10) can be answered by studying an equivalent minimisation problem. This equivalence depends on the so-called inf-sup property holding for the operators and . Let
These inf-sup conditions can be stated as
Condition (15) is proved in [32, Lemma 3.2.7] and (16) follows from the inverse mapping theorem because is surjective onto the closed Banach space . We also define the space of divergence-free functions and the convex set as
Then, (10) is equivalent to the minimisation of the functional
over , see Appendix A.
If , which in our setting occurs when , then the Korn inequality in [9, Lemma 3] implies that is a norm on equivalent to . Since
it follows that the problem is coercive. However, whenever the set of rigid body modes is not the zero set, we then say that the operator is semicoercive with respect to the seminorm because the bound (18) does not hold for . Below, in Theorem 2, we show that a consequence of the semicoercivity of is that (10) will have a unique solution only when the following compatibility condition holds:
Condition (19) will not only allow us to establish the well-posedness of (10), but will also be required for proving velocity error estimates in the -norm. This is possible because the map from is a continuous map defined over a compact set. Therefore, whenever (19) holds, the inequality
follows with the constant defined as
The importance of the compatibility condition (19) is well-known in the study of semicoercive variational inequalities, see [31, 40, 34] in the context of general variational inequalities and [39, 9] in a glaciological setting.
If , then or tend to infinity; as a result, and the functional is coercive. Uniqueness of the solution follows from the strong convexity of the functional over , see the proof for [9, Theorem 1].
Regarding the necessity of (19) for the existence and uniqueness of solutions to (10) when , assume by contradiction that is a unique solution and that there is a rigid mode for which . By testing with in (10a), we find that
and therefore we must have and . It is then straightforward to see that is also a solution to (10), contradicting our initial assumption.
We can bound the norm of the pressure by using (22), the equality
The compatibility condition (19) becomes redundant whenever the portion of the boundary with no slip boundary conditions has a positive measure, i.e. . On the other hand, (19) cannot be satisfied if there exist non zero rigid body modes tangential to . In this case, if solves (10), then for any , with on , the triple will also solve (10).
3 Abstract discretisation
In this section we propose an abstract discretisation of the mixed system (10) built in terms of a collection of finite dimensional spaces satisfying certain key properties. We can then introduce a discrete system analogous to (10) and investigate the conditions under which we have a unique solution. Then, we prove Lemmas 2, 3, and 4, which establish upper bounds for the errors of the discrete solutions.
3.1 The discrete mixed formulation
For each parameter , let , and be finite dimensional subspaces. We assume that to avoid the need of introducing discrete compatibility conditions. We define the discrete cone
and the discrete convex set
An immediate consequence of the definitions of and is that but unless . Additionally, we have that thanks to the assumption .
The discrete analogue of the variational inequality (9) is: find such that
This discrete variational inequality can be written as a mixed problem by introducing a Lagrange multiplier. This results in the discrete mixed formulation that is the counterpart of (10): find such that equationparentequation
An advantage of using a mixed formulation at the discrete level is that we explicitly enforce a discrete version of the contact conditions (5a). Just as in (11), it is possible to show that the conditions and (24c) are equivalent to
In order to state a minimisation problem equivalent to (24), we must introduce the subspace of of discretely divergence-free functions and the discrete convex set :
Then, the discrete mixed problem (24) is equivalent to the minimisation over of the functional defined in (17), provided that two discrete inf-sup conditions hold. For , these discrete conditions can be stated as
When the conditions (26) and (27) hold, then (23), (24) and the minimisation of over are equivalent problems. The proofs for such equivalences require the same arguments as the proofs presented in Appendix A. If we additionally assume to be coercive over , then admits a unique minimiser over and the discrete mixed formulation is well-posed. For these reasons, the discrete inf-sup conditions guarantee a unique solution for (24) and set constraints on the choice of spaces , and used when approximating solutions of (10).
Analogously to the continuous case, the coercivity of over hinges on the compatibility condition (19).
Assume that (26) and (27) hold. Whenever , the discrete mixed problem (23) has a unique solution. On the other hand, if , then (24) admits a unique solution if and only if the compatibility condition (19) holds. Additionally, if (19) holds when applicable, any discrete solution to (24) is uniformly bounded from above, that is,
where the constant depends on .
The converse statement and the bound (28) again follow by taking the same steps as in the continuous case. These steps can be taken due to our definitions of and and the assumption .
3.2 Upper bounds for the velocity error
An important tool presented in [5, 30] for establishing error estimates for non-Newtonian flows is the use of the function111This operator is defined differently in [5, 30]. In these references, the term is denoted by . defined by
This operator is closely related to the operator . The inequalities
hold uniformly for all for two constants . The following variation of Young’s inequality,
is valid for any and , with the constant depending on . Additionally, the seminorm is connected to via the inequality
which holds for any . See [30, Lemmas 2.3, 2.4] for proofs of (30) and (32), and [5, Lemma 2.7] for (31). It is also worth mentioning that the quasi-norm presented in  can be written in terms of , see [5, Remark 2.6].
We next present what we believe is a novel approach for bounding the velocity error in the -norm which uses the ideas presented in Section 2.4. By Theorem 1 and (32), the velocity error can be decomposed into two components as
where the constant depends on . For the first term on the right of (33), which represents the rigid component of the error, we present the following result:
Set and test equation (10a) with . Reordering, we find that
The constant depends on the norms of the solution and these are bounded from above by the norm of , see (21).
As mentioned in the introduction, previous analyses of finite element approximations of semicoercive variational inequalities either only consider the error in a seminorm  or use indirect arguments to prove the convergence of the approximate solution in the complete norm [42, 25, 40, 2, 8]. In these cases, arguments by contradiction involving a sequence of triangulations are used. In Lemma 2, on the other hand, we provide a fully constructive proof for bounding the rigid component of the velocity error from above. This result is a key ingredient in obtaining the error estimates for the finite element scheme presented in the next section.
If the pair is divergence free in the sense that for all implies that , then the term in inequality (35) can be removed.
3.3 Upper bounds for the pressure and Lagrange multiplier errors
We finalise the analysis of the abstract discretisation by bounding the errors for the pressure and the Lagrange multiplier from above.
From the inf-sup condition (26) it follows that