On the finite element approximation of a semicoercive Stokes variational inequality arising in glaciology

by   Gonzalo G. de Diego, et al.
University of Oxford

Stokes variational inequalities arise in the formulation of glaciological problems involving contact. Two important examples of such problems are that of the grounding line of a marine ice sheet and the evolution of a subglacial cavity. In general, rigid modes are present in the velocity space, rendering the variational inequality semicoercive. In this work, we consider a mixed formulation of this variational inequality involving a Lagrange multiplier and provide an analysis of its finite element approximation. Error estimates in the presence of rigid modes are obtained by means of a novel technique involving metric projections onto closed convex cones. Numerical results are reported to validate the error estimates and demonstrate the advantages of using a mixed formulation in a glaciological application.



There are no comments yet.


page 1

page 2

page 3

page 4


Mixed methods for the velocity-pressure-pseudostress formulation of the Stokes eigenvalue problem

In two and three dimensional domains, we analyze mixed finite element me...

Numerical approximation of viscous contact problems applied to glacial sliding

Viscous contact problems describe the time evolution of fluid flows in c...

Superconvergent mixed finite element methods for the pseudostress-velocity formulation of the Oseen equation

We present a priori and superconvergence error estimates of mixed finite...

Mixed Variational Finite Elements for Implicit, General-Purpose Simulation of Deformables

We propose and explore a new, general-purpose method for the implicit ti...

The Hodge Laplacian on Axisymmetric Domains

We study the mixed formulation of the abstract Hodge Laplacian on axisym...

A mixed parameter formulation with applications to linear viscoelasticity

In this work we propose and analyze an abstract parameter dependent mode...

Crack propagation simulation without crack tracking algorithm: embedded discontinuity formulation with incompatible modes

We show that for the simulation of crack propagation in quasi-brittle, t...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

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 [43] 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 [35]. Since then, subglacial cavitation has been recognised as a fundamental mechanism in glacial sliding, attracting the attention of both theoretical [17, 37] and experimental [44] 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.

1.2 Notation

Given two normed vector spaces

and 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 [18]: equationparentequation


In the equations above, represents the ice velocity, the pressure and

is a prescribed body force due to gravitational forces. The tensor

is 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 [22], 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.

In order to build a weak formulation of (1) and (4)-(6), we must first define suitable function spaces in which to seek the velocity and the pressure. Let and

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 order for (8) to make sense, we require and . Then, (1) with boundary conditions (4)-(6) can be reformulated as the variational inequality: find such that


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


The Lagrange multiplier essentially coincides with on . Indeed, if the solution to (10) is sufficiently smooth for integration by parts to hold, we arrive at on . Moreover, the conditions and (10c) are equivalent to


which is a weak representation of the contact boundary conditions (5a).

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 [41]: for a bounded and Lipschitz domain there is a constant such that

Lemma 1

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.

The main result of this section, Theorem 1, follows from Lemma 1 together with the fact that the range of is a closed, convex set whose intersection with is .

Theorem 1

Let denote the metric projection onto the closed convex cone . There is a such that


By Lemma 10, the range of is the closed convex set , defined as in (58). Moreover, we clearly have that

Therefore, we may apply Lemma 1 to prove the theorem.

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.

Theorem 2

If , then a solution to (10) exists and is unique. If , then there is a unique solution to (10) if and only if the compatibility condition (19) holds. Moreover, provided (19) holds when applicable, any solution of (10) is bounded from above, i.e.


where the constant depends on .


By Lemmas 7 and 8, the mixed problem (10) is equivalent to the minimisation of , defined in (17). A minimiser will exist if is coercive, see [15, Theorem 2, Section 8.2]. By Theorem 1 and (19),

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.

Finally, we can show that a solution to (10) is bounded from above by first setting in (10) and noting that

We have that , so it follows from Theorem 1 and (18) that


We can bound the norm of the pressure by using (22), the equality

and the inf-sup condition (15). The same follows for the norm of using (16). The last term to bound is the norm of . From (10a) we deduce that

We complete the proof by using (19), Theorem 1 and the bounds for the norms of and .


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).

Theorem 3

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 .


This proof proceeds in exactly the same way as the proof for Theorem 2. If (19) holds, we can obtain the bound

for any by using Theorem 1 and (20). This shows that a solution exists. The proof for the uniqueness is the same as in the continuous case.

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 [4] 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:

Lemma 2

Assume that the compatibility condition (19) holds. Let be the solution to (10) and to (24). Then, there is a constant depending on such that



Set and test equation (10a) with . Reordering, we find that

Since the compatibility conditions hold, we may apply (20) on the term to the left. Moreover, by adding and subtracting , applying Hölder’s inequality and Theorem 1, we establish the inequality

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 [31] 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.

The second term on the right of (33) can be bounded from above with inequality (18) and by using the properties of .

Lemma 3

Let the triples and be solutions to (10) and (24), respectively. Then there is a which depends on such that



From (10a) and (24a), we deduce that

holds for any . Using (30) and (31), we have

for an arbitrary . Additionally, by using Young’s inequality, we deduce that

for any . Then, via (32), and by setting and sufficiently small, inequality (35) is established.


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.

Lemma 4

Assume that the discrete inf-sup conditions (26) and (27) hold. Let be the solution of (10) and of (24). Then, there is a constant depending on and such that


for all and .


Since and are subsets of and respectively, we can obtain the following equality from (10a) and (24a):


The inf-sup condition (26) for the pressure space holds over the space of vector fields with a normal component vanishing on . Therefore, for , and from equation (38) we derive

From the inf-sup condition (26) it follows that