Volume 2011 (2011), Article ID 164581, 24 pages
http://dx.doi.org/10.1155/2011/164581
Research Article

## Analysis and Finite Element Approximation of a Nonlinear Stationary Stokes Problem Arising in Glaciology

1Department of Mathematics and Computer Science, Free University of Berlin, 14195 Berlin, Germany
2Chair of Numerical Analysis and Simulation, EPFL, 1015 Lausanne, Switzerland

Received 13 September 2011; Accepted 25 October 2011

Copyright © 2011 Guillaume Jouvet and Jacques Rappaz. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

The aim of this paper is to study a nonlinear stationary Stokes problem with mixed boundary conditions that describes the ice velocity and pressure fields of grounded glaciers under Glen's flow law. Using convex analysis arguments, we prove the existence and the uniqueness of a weak solution. A finite element method is applied with approximation spaces that satisfy the inf-sup condition, and a priori error estimates are established by using a quasinorm technique. Several algorithms (including Newton's method) are proposed to solve the nonlinearity of the Stokes problem and are proved to be convergent. Our results are supported by numerical convergence studies.

#### 1. Introduction

In this paper we consider a model problem that is commonly used by glaciologists to compute the motion of glaciers. Ice is assumed to be an incompressible non-Newtonian fluid governed by Glen's law [1]. Glen's law and the mass momentum equation lead to a nonlinear stationary Stokes problem with a strain-dependent viscosity.

Glacier models based on Glen's law have already been studied by several authors. However, all of them have considered a simplified model, called first-order approximation [2]. This model is obtained by rewriting the Stokes equations into a dimensionless form and by dropping all terms of order , where is the typical aspect ratio of glaciers. This simplification results into a nonlinear elliptic problem for the horizontal velocity field, the vertical component, and the pressure field being determined a posteriori. Colinge and Rappaz first demonstrated the well-posedness of this problem and proved the convergence of the finite element approximation with piecewise linear continuous functions in [3]. Inspired by the work of Baranger/Najib [4] and Barrett/Liu [5] on non-Newtonian problems, a priori and a posteriori error estimates were obtained later in [68].

Unlike the first-order approximation, the original Stokes model which is considered in this paper is a saddle point problem for the velocity and the pressure fields. We prove the existence and the uniqueness of a weak solution using an equivalent minimisation problem and an inf-sup stability condition. Next, we establish a priori estimates for a finite element approximation using a quasinorm technique [5]. Eventually, we investigate several successive approximation algorithms to solve the system nonlinearity. In particular, we upgrade by using Newton's method the fixed point algorithm, given in [3] and proved to be convergent in [6, 9].

Boundary conditions describe the basal sliding phenomena that can significantly influence the glacier ice flows. In [3, 68], the first order approximation model was coupled to a Dirichlet condition. However, this approach requires the basal velocity distribution which is unknown. To overcome this difficulty, several sliding laws—including a Coulomb-type law—were considered in [10]. In our model, we use a sliding law that results in a nonlinear Dirichlet-Robin boundary condition.

This paper is organised as follows: the physical model is presented in Section 2. We prove the well-posedness of the weak problem in Section 3. In Section 4, we apply a finite element method and establish a priori error estimates. Successive approximation algorithms to solve the system nonlinearity are proposed and proved to be convergent in Section 5. In Section 6, convergence studies are performed to support the results of Sections 4 and 5.

#### 2. The Model

Let us suppose that ice occupies the domain , with or 3. Ice can be considered as an incompressible non-Newtonian fluid with negligible inertial effects [11]. It follows that the velocity and the pressure of ice solve the stationary nonlinear Stokes problem in : where denotes the rate of strain tensor, the viscosity of ice, and the gravity force. Here above, the viscosity depends on and is defined by the regularised Glen's flow law [11]. More precisely, for a given velocity field , the viscosity satisfies the following nonlinear equation: where , is a positive parameter, is Glen's exponent, and is a small regularization parameter which prevents infinite viscosity for zero strain ( in the original Glen's law [1]). When , then the viscosity is constant and (2.1) correspond to the classical linear Stokes problem related to a Newtonian fluid. In the framework of glaciology, is often taken equal to 3; see [12].

Let us set the boundary conditions for the system of (2.1). Three mechanical circumstances may occur at the boundary of a glacier: (i) no force applies on the ice-air interface; (ii) ice slides on the bedrock-ice interface; (iii) ice is stuck to the bedrock-ice interface. The boundary of is thus split into three parts: , , and , referring to circumstances (i), (ii), and (iii), respectively. We assume throughout that is bounded, its boundaries and , are and . We consider the free surface condition: where is the unit outward normal vector along the boundary of the domain . We apply the nonlinear sliding condition [10, 13, 14]: where are the orthogonal vectors tangent to the boundary , that is, when and when . Here above, is the sliding coefficient that is given by where is the Euclidean norm of , is Glen's exponent, is a positive parameter, and is a small parameter which prevents infinite for zero velocity. The no-sliding condition writes

Note that the conditions applied on boundaries , , and are Neumann, Robin-Dirichlet, and Dirichlet conditions, respectively. When (Newtonian flow) and , the problem (2.1) with boundary conditions (2.3), (2.6) has already been widely studied; see, for instance, [1517].

#### 3. Existence and Uniqueness

In this section, we prove that there exists a unique weak solution to problem (2.1) with mixed boundary conditions (2.3), (2.4), and (2.6). Pressure is first eliminated from the system by restricting the velocity space to divergence-free fields. Afterwards, the reduced problem is transformed into a minimisation problem. Following [3, 8], its well-posedness is proved by using convex analysis arguments. The existence and the uniqueness of the pressure field are ensured by an inf-sup condition. We now state in the next lemma several properties of the function that will often be used in Sections 3, 4 and 5.

Lemma 3.1. For all , there exists a unique satisfying (2.2). The function is and decreasing. There exist such that:

Proof. The properties of and inequalities (3.1) and (3.2) can be easily deduced from Lemmas  1 and 2 of [7]. Inequality (3.3) is obtained by differentiating (2.2) with respect to . Inequalities (3.4) and (3.5) result from inequality (3.1), Lemma  2.1 in [5] and inequality . Details are given in [12].

Let us notice that property (3.1) was introduced by Barrett and Liu (see [5]) in order to obtain a priori error estimates of a similar problem to the one treated in this paper. Define the Banach spaces: where are conjugate exponents and is Glen's exponent. By using (3.2), we have for all . Then, if , we have . By using the trace inequality for all , see [19 page 197], we obtain . Similarly, we can show . Owing to Hölder's inequality, the mixed formulation of problem (2.1) with boundary conditions (2.3), (2.4), and (2.6) that consists of finding such that is meaningful.

Remark 3.2. If , pressure in (3.8) is defined up to a constant. In that case, is replaced by . Moreover, if , then and is constant and if , then the (linear) problem (3.8) is well posed; see, for instance, [1517].

The next lemma states the equivalence of norms and on space .

Lemma 3.3 (Korn's inequality). If and if , then there exists a constant such that for all such that on .

Proof. We apply Corollary  4.1 in [18] ( being the identity matrix) and Lemma  3.1 page 40 in [16].

We consider the divergence-free velocity space:

In , the pressure field vanishes of the variational formulation (3.8). The reduced formulation consists then of finding such that To transform problem (3.11) into a minimisation problem, we introduce the functional where The functional is Gâteaux differentiable, and its first derivative , at point , in direction , is given by Clearly, any minimiser of in satisfies (3.11). We now establish several lemmas that allow us to prove the existence and the uniqueness of this minimiser in Theorem 3.8. We show the continuity of in Lemma 3.5, the strict convexity of in Lemma 3.6, and the coercivity (in the sense of (3.18)) of in Lemma 3.7. The continuity of requires the following result (Lemma  4 in [3])

Lemma 3.4. Let be a measurable set of and , then one has the following inequality:

Lemma 3.5. The functional is -continuous.

Proof. By using (3.2), (2.5), and , we have, for all These two inequalities together with Lemma 3.4 imply the -continuity of .

Lemma 3.6. The functional is strictly convex on .

Proof. Clearly, and . From (3.3), we have if , and then is strictly convex. Since is an increasing function, is strictly convex. In the same way, we can show that is strictly convex by using (2.5). Let satisfying and . From Korn's inequality (Lemma 3.3), we have in . As a consequence, The strict convexity of follows from the previous inequality and the convexity of .

Since is convex, satisfies (3.11) if and only if .

Lemma 3.7. There exist two constants such that, for all ,

Proof. Let . From (3.2) and , there exists such that As a consequence, there exist two constants such that . By using Korn's inequality (Lemma 3.3), there exists such that From Young's inequality, we have, for all , where . We set small enough such that . From inequalities (3.20), (3.21), and , we obtain which is exactly (3.18) with and .

Theorem 3.8. There exists a unique such that . Moreover, is the unique solution of (3.11).

Proof. Clearly, there exists such that . Lemma 3.7 ensures the existence of . Let be a sequence of such that . There exists an integer such that, for all , we have . Owing to Lemma 3.7, the sequence is bounded in . Since is a closed subspace of , is reflexive. Consequently, there exist and a subsequence of (still denoted ) that converges weakly to in . By Lemmas 3.5 and 3.6, is weakly lower semicontinuous; see, for instance, Corollary III.8 in [19] page 38. Then, we have and possesses at least one minimum . Since is strictly convex (Lemma 3.6), this minimum is unique. Moreover, is the unique solution of (3.11).

Spaces and are required to satisfy the inf-sup condition, see [5, 20], to ensure the existence and the uniqueness of such that satisfies the mixed formulation (3.8). The inf-sup condition is proved in [15, 21] when (or, equivalently, ). By following the proof of Proposition  5.3.2 in [22], we can easily generalise this result when ; see details in [12].

Lemma 3.9. Spaces and satisfy the inf-sup condition; that is, there exists such that

Theorem 3.10. There exists a unique couple satisfying (3.8).

Proof. Although the result is a straightforward application of Theorem  2.1 in [20] together with Theorem 3.8 and Lemma 3.9, we give all the arguments of the proof. Let and be the operators defined by where and are dual to and , respectively. From Theorem 3.8, there exists a unique such that for all , which means that . Owing to the inf-sup condition (4.1), the operator is surjective, , and is closed; see Lemma  A.40 in [15]. As a consequence, and there exists such that . Since , the pressure is necessarily unique. Eventually, there exists a unique couple satisfying or equivalently (3.8).

#### 4. Finite Element Approximation and A Priori Estimates

We assume that is a convex polygonal or polyhedral domain and is a regular mesh of parametrized by , the highest diameter of the elements of . We say that and , some finite-dimensional approximation spaces on of and , satisfy the inf-sup condition if, for all , there exists a constant such that The discrete problem is obtained by replacing the spaces and by and , respectively. It consists of finding such that The discrete similar space to is Note that is not necessarily included in . The discrete reduced problem consists of finding such that Since is a closed subspace of , Theorem 3.8 and the proof can be rewritten by replacing by .

Theorem 4.1. There exists a unique such that . Moreover, is the unique solution of (4.4).

Remark 4.2. By setting in (4.4) and by using inequality (3.2), (2.5), and Korn's inequality (Lemma 3.3), we can show that the solution of problem (4.4) satisfies where does not depend on .

From Theorem 4.1 and the inf-sup condition (4.1), we can rewrite Theorem 3.10 and its proof for the discrete mixed problem.

Theorem 4.3. If and satisfy the inf-sup condition (4.1), then there exists a unique couple satisfying (4.2).

Remark 4.4. The spaces and are two examples that satisfy the inf-sup condition (4.1) while does not satisfy (4.1); see [15].

The error analysis that follows is partly inspired from [5, 7]. We give a priori estimates for the numerical approximation of the stationary Stokes problem in Theorem 4.9. For the sake of simplicity, we suppose ; that is, the boundary Robin-Dirichlet condition is not considered; see also Remark 4.11. The nonlinearity of problem (3.8) is treated by introducing (in Lemma 4.5) a quasi-norm that depends on the solution; see [5]. The orthogonality of the error (Lemma 4.6) together with properties (3.4) and (3.5) of the function allow quasi-norm estimates to be established in Theorem 4.7. The properties of the quasi-norm given in Lemma 4.5 allow estimates with standard norms to be proved in Theorem 4.8. Eventually, these estimates together with interpolation inequalities yield to the main Theorem 4.9.

Lemma 4.5. Let be the solution of (3.8); the application is a quasi-norm of ; that is, it satisfies all properties of norms, except homogeneity. Moreover, there exists such that, for all , one has and there exists such that, for all and for all , one has

Proof. The quasi-norm properties are shown in Lemma  3.1 in [5]. Inequalities (4.7) and (4.8) result from Korn and Hölder's inequalities; see details in [12].

By setting in (3.8) it is easy to prove the next lemma.

Lemma 4.6. Let be the solution of problem (3.8) and the solution of problem (4.4), then holds for all . Moreover, if the spaces and satisfy the condition (4.1), then the solution of (4.2) satisfies for all .

Theorem 4.7. Let be the solution of (3.8) and the solution of (4.4). For all , one has where . Moreover, if the spaces and satisfy the condition (4.1), then the solution of (4.2) satisfies, for all and for all , where . The constants do not depend on and ; however, increasingly depends on .

Proof. By using, respectively, the definition (4.6) of the quasi-norm , inequality (3.4) with , and (4.9), there exists such that where . For the sake of simplicity, and are handled separately. By using inequality (3.5) with , there exists such that By using the inequality (see Lemma  2.2 in [5] or  (3.10) in [8]), with , , , and , we obtain, for all , We now use respectively Hölder's inequality, Young's inequality, and (4.7); there exist such that By setting and , we obtain By moving to the left-hand side, we obtain (4.11). From the inf-sup condition (4.1), we have, for all , From (4.10), we have, for all , From (4.19), (4.20), and (3.5) with , there exist such that where . Eventually, the previous inequality together with leads to (4.12).

Theorem 4.8. Let be the solution of (3.8), and the solution of (4.4). For all and for all , assuming , one has where . Moreover, if the spaces and satisfy the inf-sup condition (4.1), then the solution of (4.2) satisfies for all and for all , assuming , where . The constants do not depend on and ; however, and increasingly depends on .

Proof. On one hand, inequality (4.22) follows from inequalities (4.7), (4.11), and (4.8). On the other hand, (4.23) follows from (4.22) and from the following property (see (1.16), page 115 in [16]): for all and for all , there exists such that where depends on the inf-sup constant . Eventually, (4.24) follows from (4.12), (4.11), and (4.8).

Theorem 4.9. Assume that, for all , there exists a continuous operator that satisfies and a continuous operator that satisfies where is the size of the higher diameter of the elements of . Assume that and satisfy the inf-sup condition (4.1). Let be the solution of problem (3.8) and let be the solution of problem (4.2). Assume that , where , then one has where .

Proof. Apply (4.23) and (4.24) with and . By using the continuity of , (4.5), (4.26), and (4.27), there exist such that The estimate (4.28) directly follows from (4.29).

Remark 4.10. The combination for spaces and , introduced in [23], satisfies the assumptions of Theorem 4.9; see Lemma  4.20 page 190 of [15] for the inf-sup condition (4.1) and [15, 16] for the interpolation inequalities (4.26) and (4.27).

Remark 4.11. If , a similar analysis can be led by replacing the norm defined by (4.6) by

#### 5. Successive Approximations

In this section, several successive approximation algorithms are proposed for solving the nonlinearity of the discrete problem (4.4) when . For the sake of simplicity, we suppose in this section; see Remark 5.8. We present a unified scheme that contains the classical fixed point method together with Newton's method. The mesh is fixed, and we assume that the approximation spaces satisfy and . In what follows, denotes an arbitrary norm of . Since is a finite-dimensional space, all norms are equivalent. Let . We define where solves

The application is well defined. Indeed, by using, respectively, , inequalities (3.3) and (3.2), there exist such that As a consequence, the problem (5.2) is coercive. From the Lax-Milgram Theorem, see [15] page 83, there exists a unique solution of (5.2).

In what follows, denotes the solution (4.4), which is also the unique fixed point of . Assume that is given; we define iteratively a sequence , for all , by Our goal is to prove that converges to when goes to the infinity. When , we obtain the classical fixed point method, widely used to solve the nonlinearity of Glen's law; see [6, 9, 14]. When , we have an additional term in (5.2) which corresponds to Newton's method; see Remark 5.5. The case corresponds to a hybrid fixed point—Newton's method. The convergence of sequence requires several preliminary results. We compute the first derivative of in Lemma 5.1. Lemma 5.2 provides an upper bound of the first derivative. Eventually, Theorem 5.3 states the linear convergence of by using the Banach fixed point theorem. Theorem 5.7 states the second-order convergence when . By differentiating formally (5.2) at point in direction , with , we obtain the following lemma.

Lemma 5.1. Let satisfy . The application is Gâteaux differentiable at point and its derivative is given by where solves

The problems (5.2) and (5.6) have the same coercivity properties to compute (resp., ) from (resp., ). As a consequence, the problem (5.6) is well-posed by the Lax-Milgram theorem. To prove the convergence of the sequence , we look for a norm that makes a contraction at point .

Lemma 5.2. Let , and let be the fixed point of . The application satisfies where is the subordinated norm to .

Proof. Since , then (5.6), with and , is rewriten as for all . From (5.8), , and (3.3), we have By setting in (5.9), we obtain From (5.10) and Cauchy-Schwarz's inequality, we obtain Eventually, (5.7) follows from the definition of norm .

Theorem 5.3. Let , let be a given mesh of , and let be a norm of . There exist and such that if , then one has for all , and is linearly convergent to .

Proof. From Lemma 5.2, the spectral radius of is lower than constant: which is lower than 1. The theorem is then a direct application of the Banach fixed point theorem in .

Remark 5.4. In Theorem 5.3, it should be stressed that depends on . When (i.e., if we replace by ), we cannot ensure Theorem 5.3 to remain true. Nevertheless, in practise, seems to be independent of , see Section 6.

When , we have from Lemma 5.2. It suggests that the convergence of sequence is quadratic. To establish the second-order convergence, we define, for all , the application by

Let . We compute formally the first-order derivative of at point in direction :