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

Authors: Guillaume Jouvet; Jacques Rappaz
Journal: Advances in Numerical Analysis
 (2011)

---

## 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, anda 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.


---

## Body

## 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, calledfirst-order approximation [2]. This model is obtained by rewriting the Stokes equations into a dimensionless form and by dropping all terms of order 𝒪(ϵ2), 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 [6–8].Unlike thefirst-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, 6–8], 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 Section2. 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Ω⊂ℝd, with d=2 or 3. Ice can be considered as an incompressible non-Newtonian fluid with negligible inertial effects [11]. It follows that the velocity u and the pressure p of ice solve the stationary nonlinear Stokes problem in Ω:
(2.1)-2div⁡(με(u))+∇p=f,div⁡(u)=0,
where ε(u)=(1/2)(∇u+∇uT) denotes the rate of strain tensor, μ the viscosity of ice, and f the gravity force. Here above, the viscosity μ depends on |ε(u)|:=ε(u):ε(u) and is defined by the regularised Glen's flow law [11]. More precisely, for a given velocity field u, the viscosity μ satisfies the following nonlinear equation:
(2.2)12μ=A(τ0n-1+(2μs)n-1),
where s=|ε(u)|, A is a positive parameter, n≥1 is Glen's exponent, and τ0>0 is a small regularization parameter which prevents infinite viscosity for zero strain (τ0=0 in the original Glen's law [1]). When n=1, 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, n 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: ΓN, ΓR, and ΓD, referring to circumstances (i), (ii), and (iii), respectively. We assume throughout that Ω is bounded, its boundaries ΓN and ΓR, are 𝒞1 and ΓD≠∅. We consider the free surface condition:
(2.3)2με(u)⋅n-pn=0,on⁡ΓN,
where n is the unit outward normal vector along the boundary of the domain Ω. We apply the nonlinear sliding condition [10, 13, 14]:
(2.4)u⋅n=0,(2με(u)⋅n)⋅ti=-αu⋅ti,i=1,d-1on⁡ΓR,
where {ti}i=1,d-1 are the orthogonal vectors tangent to the boundary ΓR, that is, t1 when d=2 and t1,t2 when d=3. Here above, α=α(|u|) is the sliding coefficient that is given by
(2.5)α(t)=c(t+t0)1/n-1,
where t=|u| is the Euclidean norm of u, n is Glen's exponent, c is a positive parameter, and t0>0 is a small parameter which prevents infinite α for zero velocity. The no-sliding condition writes
(2.6)u=0,on⁡ΓD.Note that the conditions applied on boundariesΓN, ΓR, and ΓD are Neumann, Robin-Dirichlet, and Dirichlet conditions, respectively. When n=1 (Newtonian flow) and ΓR=∅, the problem (2.1) with boundary conditions (2.3), (2.6) has already been widely studied; see, for instance, [15–17].


## 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 alls∈ℝ+, there exists a unique μ=μ(s)∈ℝ+ satisfying (2.2). The function s→μ(s) is 𝒞∞(0,+∞) and decreasing. There exist D1,D2,D3,D4>0 such that:
(3.1)D1(1+s+t)1-1/n(s-t)≤sμ(s)-tμ(t)≤D2(1+s+t)1-(1/n)(s-t),∀s≥t≥0,(3.2)D1(1+s)1-1/n≤μ(s)≤D2(1+s)1-1/n,∀s≥0,(3.3)1nμ(s)≤μ(s)+sμ′(s),∀s>0,(3.4)(μ(|ξ|)ξ-μ(|η|)η):(ξ-η)≥D3(1+|ξ|+|ξ-η|)(1/n-1)|ξ-η|2,∀ξ,η∈Rd×d,(3.5)|μ(|ξ|)ξ-μ(|η|)η|≤D4(1+|ξ|+|ξ-η|)(1/n-1)|ξ-η|,∀ξ,η∈Rd×d.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 s. Inequalities (3.4) and (3.5) result from inequality (3.1), Lemma  2.1 in [5] and inequality (1/2)(|ξ|+|η|)≤|ξ|+|ξ-η|≤2(|ξ|+|η|). 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:
(3.6)V:={v∈[W1,r(Ω)]d,v=0on⁡ΓD,v⋅n=0on⁡ΓR},Q:=Lr′(Ω),
where
(3.7)r:=1+1n,r′:=n+1
are conjugate exponents and n is Glen's exponent. By using (3.2), we have μ(s)s≤Csr-1 for all s>0. Then, if u∈V, we have μ(|ε(u)|)ε(u)∈[Lr′(Ω)]d×d. By using the trace inequality ∥v∥Lr(ΓR)≤∥v∥W1-1/r,r(ΓR)≤C∥v∥W1,r(Ω) for all v∈[W1,r(Ω)]d, see [19 page 197], we obtain (u·ti)∈Lr(ΓR),i=1,d-1. Similarly, we can show α(|u|)(u·ti)∈Lr′(ΓR),i=1,d-1. 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 (u,p)∈V×Q such that
(3.8)2∫Ωμ(|ε(u)|)ε(u):ε(v)dV+∑i=1,d-1∫ΓRα(|u|)(u⋅ti)(v⋅ti)dS,-∫Ωpdiv⁡(v)dV+∫Ωqdiv⁡(u)dV=∫Ωg⋅vdV,∀(v,q)∈V×Q
is meaningful.Remark 3.2.
IfΓN=∅, pressure p in (3.8) is defined up to a constant. In that case, Q=Lr′ is replaced by Q=L0r′={q∈Lr′,∫ΩqdV=0}. Moreover, if n=1, then r=2 and μ is constant and if ΓR=∅, then the (linear) problem (3.8) is well posed; see, for instance, [15–17].The next lemma states the equivalence of norms∥|ε(·)|∥Lr and ∥·∥W1,r on space V.Lemma 3.3 (Korn's inequality).
IfΓD≠∅ and if 1<γ<∞, then there exists a constant C>0 such that
(3.9)‖v‖W1,γ≤C‖|ε(v)|‖Lγ,
for all v∈W1,γ(Ω) such that v=0 on ΓD.Proof.
We apply Corollary  4.1 in [18] (F being the identity matrix) and Lemma  3.1 page 40 in [16].We consider the divergence-free velocity space:(3.10)Vdiv⁡:={v∈[W1,r(Ω)]d,div⁡(v)=0,v=0on⁡ΓD,v⋅n=0on⁡ΓR}.InVdiv, the pressure field p vanishes of the variational formulation (3.8). The reduced formulation consists then of finding u∈Vdiv such that
(3.11)2∫Ωμ(|ε(u)|)ε(u):ε(v)dV+∑i=1,d-1∫ΓRα(|u|)(u⋅ti)(v⋅ti)dS=∫Ωg⋅vdV,∀v∈Vdiv⁡.
To transform problem (3.11) into a minimisation problem, we introduce the functional
(3.12)J(u):=∫ΩM(|ε(u)|)dV+12∫ΓRN(|u|)dS-∫Ωu⋅fdV,
where
(3.13)M(x):=∫0xsμ(s)ds,N(x):=∫0xtα(t)dt.
The functional J is Gâteaux differentiable, and its first derivative DJ, at point u∈Vdiv, in direction v∈Vdiv, is given by
(3.14)〈DJ(u),v〉=2∫Ωμ(|ε(u)|)ε(u):ε(v)dV+∑i=1,d-1∫ΓRα(|u|)(u⋅ti)(v⋅ti)dS-∫Ωg⋅vdV.
Clearly, any minimiser of J in Vdiv 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 J in Lemma 3.5, the strict convexity of J in Lemma 3.6, and the coercivity (in the sense of (3.18)) of J in Lemma 3.7. The continuity of J requires the following result (Lemma  4 in [3])Lemma 3.4.
LetO be a measurable set of ℝd and  f,g∈Lr(O), then one has the following inequality:
(3.15)∫O||f|r-|g|r|dV≤r‖|f|+|g|‖Lr(O)r-1‖f-g‖Lr(O).Lemma 3.5.
The functionalJ is  ∥·∥W1,r-continuous.Proof.
By using (3.2), (2.5), and 1-1/n=2-r, we have, for all u,v∈Vdiv(3.16)|M(|ε(u)|)-M(|ε(v)|)|=|∫|ε(v)||ε(u)|sμ(s)ds|≤|D1∫|ε(v)||ε(u)|sr-1ds|=D1||ε(u)|r-|ε(v)|rr|,|N(|u|)-N(|v|)|=|∫|v||u|tα(t)dt|≤|c∫|v||u|tr-1dt|=c||u|r-|v|rr|.
These two inequalities together with Lemma 3.4 imply the ∥·∥W1,r-continuity of J.Lemma 3.6.
The functionalJ is strictly convex on V.Proof.
Clearly,M′(s)=sμ(s) and M′′(s)=sμ′(s)+μ(s). From (3.3), we have M′′(s)=sμ′(s)+μ(s)≥(1/n)μ(s)>0 if s>0, and then M is strictly convex. Since M is an increasing function, M(|·|) is strictly convex. In the same way, we can show that N(|·|) is strictly convex by using (2.5). Let u,v∈Vdiv satisfying u≠v and θ∈(0,1). From Korn's inequality (Lemma 3.3), we have ε(u)≠ε(v) in Lr. As a consequence,
(3.17)∫ΩM(|θε(u)+(1-θ)ε(v)|)dV<θ∫ΩM(|ε(u)|)dV+(1-θ)∫ΩM(|ε(v)|)dV.
The strict convexity of J follows from the previous inequality and the convexity of N(|·|).SinceJ is convex, u∈Vdiv satisfies (3.11) if and only if J(u)≤J(v),∀v∈Vdiv.Lemma 3.7.
There exist two constantsD1,D2>0 such that, for all  u∈V,
(3.18)J(u)≥D1‖u‖W1,rr-D2.Proof.
Letu∈V. From (3.2) and 1-1/n=2-r, there exists C0>0 such that
(3.19)M(|ε(u)|)≥∫0|ε(u)|C0s(1+s)2-rds=∫0|ε(u)|2/2C0(1+2t)2-rdt≥|ε(u)|22C0(1+|ε(u)|)r-2.
As a consequence, there exist two constants C1,C2>0 such that M(|ε(u)|)≥C1(1+|ε(u)|)r-C2. By using Korn's inequality (Lemma 3.3), there exists C3>0 such that
(3.20)∫ΩM(|ε(u)|)dV≥C3‖u‖W1,rr-C2∫ΩdV.
From Young's inequality, we have, for all δ>0,
(3.21)∫Ω|u⋅f|dV≤∫Ω(1r′δr′|f|r′+δrr|u|r)dV=C4δr′∫ΩdV+C5δr‖u‖Lrr,
where C4,C5>0. We set δ small enough such that C3-δrC5>0. From inequalities (3.20), (3.21), and N≥0, we obtain
(3.22)J(u)=∫Ω(M(|ε(u)|)-u⋅f)dV+∫ΓR12N(|u|)dS≥C3‖u‖W1,rr-C2∫ΩdV-C4δr′∫ΩdV-C5δr‖u‖W1,rr,
which is exactly (3.18) with D1:=C3-C5δr and D2:=(C2+C4/δr′)∫ΩdV.Theorem 3.8.
There exists a uniqueu∈Vdiv such that J(u)=inf{J(v);v∈Vdiv}. Moreover, u is the unique solution of (3.11).Proof.
Clearly, there existsu̅∈Vdiv such that J(u̅)<+∞. Lemma 3.7 ensures the existence of m̅=inf{J(v);v∈Vdiv}. Let {uν} be a sequence of Vdiv such that limν→∞J(uν)=m̅. There exists an integer K such that, for all ν>K, we have m̅+1>J(uν). Owing to Lemma 3.7, the sequence {uν} is bounded in Vdiv. Since Vdiv is a closed subspace of V, Vdiv is reflexive. Consequently, there exist u∈Vdiv and a subsequence of {uν} (still denoted {uν}) that converges weakly to u in Vdiv. By Lemmas 3.5 and 3.6, J is weakly lower semicontinuous; see, for instance, Corollary III.8 in [19] page 38. Then, we have
(3.23)m̅=lim⁡inf⁡ν→+∞J(uν)≥J(u)≥m̅,
and J possesses at least one minimum u∈Vdiv. Since J is strictly convex (Lemma 3.6), this minimum is unique. Moreover, u is the unique solution of (3.11).SpacesV and Q are required to satisfy the inf-sup condition, see [5, 20], to ensure the existence and the uniqueness of p∈Q such that (u,p) satisfies the mixed formulation (3.8). The inf-sup condition is proved in [15, 21] when ΓD=∂Ω (or, equivalently, ΓN=ΓR=∅). By following the proof of Proposition  5.3.2 in [22], we can easily generalise this result when ΓR∪ΓN≠∅; see details in [12].Lemma 3.9.
SpacesV and Q satisfy the inf-sup condition; that is, there exists C>0 such that
(3.24)C<inf⁡q∈Qsup⁡v∈V∫Ωqdiv⁡(v)dV‖q‖Lr′‖v‖W1,r.Theorem 3.10.
There exists a unique couple(u,p)∈(V,Q) 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 A:V→V′ and B:V→Q′ be the operators defined by
(3.25)〈Au,v〉=2∫Ωμ(|ε(u)|)ε(u):ε(v)dV+∑i=1,d-1∫ΓRα(|u|)(u⋅ti)(v⋅ti)dS,∀v∈V,〈Bu,q〉=∫Ωqdiv⁡udV,∀q∈Q,
where V′ and Q′ are dual to V and Q, respectively. From Theorem 3.8, there exists a unique u∈kerB such that 〈Au-f,v〉=0 for all v∈kerB, which means that Au-f∈(kerB)⊥. Owing to the inf-sup condition (4.1), the operator B:V→Q′ is surjective, kerBT=∅, and ℛ(BT) is closed; see Lemma  A.40 in [15]. As a consequence, Au-f∈(kerB)⊥=ℛ(BT)¯=ℛ(BT) and there exists p∈Q such that Au-f=BTp. Since kerBT=∅, the pressure p is necessarily unique. Eventually, there exists a unique couple (u,p)∈V×Q satisfying
(3.26)Au-BTp=f,Bu=0,
or equivalently (3.8).


## 4. Finite Element Approximation andA Priori Estimates

We assume thatΩ is a convex polygonal or polyhedral domain and 𝒯h is a regular mesh of Ω¯ parametrized by h, the highest diameter of the elements of 𝒯h. We say that Vh⊂V and Qh⊂Q, some finite-dimensional approximation spaces on 𝒯h of V and Q, satisfy the inf-sup condition if, for all κ∈(1,∞), there exists a constant Ch>0 such that
(4.1)Ch<inf⁡qh∈Qhsup⁡vh∈Vh∫Ωqhdiv⁡(vh)dV‖qh‖Lκ′‖vh‖W1,κ.
The discrete problem is obtained by replacing the spaces V and Q by Vh and Qh, respectively. It consists of finding (uh,ph)∈(Vh,Qh) such that
(4.2)2∫Ωμ(|ε(uh)|)ε(uh):ε(vh)dV+∑i=1,d-1∫ΓRα(|uh|)(uh⋅ti)(vh⋅ti)dS,-∫Ωphdiv⁡vhdV+∫Ωqhdiv⁡uhdV=∫Ωf⋅vhdV,∀(vh,qh)∈(Vh,Qh).
The discrete similar space to Vdiv is (4.3)Vdiv⁡,h={vh∈Vh;∫Ωqhdiv⁡(vh)dV=0;∀qh∈Qh}.
Note that Vdiv,h is not necessarily included in Vdiv. The discrete reduced problem consists of finding uh∈Vdiv,h such that
(4.4)2∫Ωμ(|ε(uh)|)ε(uh):ε(vh)dV+∑i=1,d-1∫ΓRα(|uh|)(uh⋅ti)(vh⋅ti)dS=∫Ωg⋅vhdV,∀vh∈Vdiv⁡,h.
Since Vdiv,h is a closed subspace of V, Theorem 3.8 and the proof can be rewritten by replacing Vdiv by Vdiv,h.Theorem 4.1.
There exists a uniqueuh∈Vdiv,h such that J(uh)=inf{J(vh);vh∈Vdiv,h}. Moreover, uh is the unique solution of (4.4).Remark 4.2.
By settingvh=uh in (4.4) and by using inequality (3.2), (2.5), and Korn's inequality (Lemma 3.3), we can show that the solution uh of problem (4.4) satisfies
(4.5)‖uh‖W1,r≤C‖f‖Lr′,
where C>0 does not depend on uh.From Theorem4.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.
IfVh and  Qh satisfy the inf-sup condition (4.1), then there exists a unique couple  (uh,ph)∈(Vh,Qh) satisfying (4.2).Remark 4.4.
The spaces[ℙ1/Bulle]d-ℙ1 and [ℙ2]d-ℙ1 are two examples that satisfy the inf-sup condition (4.1) while ℙ1-ℙ1 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 ΓR=∅; 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(u,p) be the solution of (3.8); the application
(4.6)v⟶‖|v|‖:=∫Ω|ε(v)|2(1+|ε(u)|+|ε(v)|)2-rdV
is a quasi-norm of V; that is, it satisfies all properties of norms, except homogeneity. Moreover, there exists D1>0 such that, for all v∈W1,r(Ω), one has
(4.7)‖v‖W1,r2≤D1[1+‖u‖W1,r+‖v‖W1,r]2-r‖|v|‖2,
and there exists D2>0 such that, for all κ∈[r,2] and for all v∈W1,κ(Ω), one has
(4.8)‖|v|‖2≤D2‖v‖W1,κκ.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 settingv=vh∈Vdiv,h in (3.8) it is easy to prove the next lemma.Lemma 4.6.
Let(u,p)∈(V,Q) be the solution of problem (3.8) and uh∈Vdiv,h the solution of problem (4.4), then
(4.9)∫Ω2(μ(|ε(u)|)ε(u)-μ(|ε(uh)|)ε(uh)):ε(vh)dV-∫Ωdiv⁡(vh)(p-qh)dV=0
holds for all (vh,qh)∈(Vdiv,h,Qh). Moreover, if the spaces Vh and Qh satisfy the inf-sup condition (4.1), then the solution (uh,ph) of (4.2) satisfies
(4.10)∫Ω2(μ(|ε(u)|)ε(u)-μ(|ε(uh)|)ε(uh)):ε(vh)dV-∫Ωdiv⁡(vh)(p-ph)dV=0,
for all vh∈Vh.Theorem 4.7.
Let(u,p)∈(V,Q) be the solution of (3.8) and uh∈Vdiv,h the solution of (4.4). For all (vh,qh)∈Vdiv,h×Qh, one has
(4.11)‖|u-uh|‖≤D1[1+‖uh‖W1,r+‖vh‖W1,r](2-r)/2(‖|u-vh|‖+‖p-qh‖Lr′),
where D1>0. Moreover, if the spaces Vh and Qh satisfy the inf-sup condition (4.1), then the solution (uh,ph) of (4.2) satisfies, for all κ∈[r,2] and for all qh∈Qh,
(4.12)‖p-ph‖Lκ′≤D2(‖|u-uh|‖2/κ′+‖p-qh‖Lκ′),
where D2>0. The constants D1,D2 do not depend on uh and vh; however, D2 increasingly depends on (Ch)-1.Proof.
By using, respectively, the definition (4.6) of the quasi-norm |||·|||, inequality (3.4) with 1-1/n=2-r, and (4.9), there exists C1>0 such that
(4.13)‖|u-uh|‖2≤C1{∫Ω2(μ(|ε(u)|)ε(u)-μ(|ε(uh)|)ε(uh)):ε(u-uh)dV}=C1{∫Ω2(μ(|ε(u)|)ε(u)-μ(|ε(uh)|)ε(uh)):ε(u-vh)dV︸:=A1+∫Ωdiv⁡(vh-uh)(p-qh)dV︸:=A2},
where (vh,qh)∈Vdiv,h×Qh. For the sake of simplicity, A1 and A2 are handled separately. By using inequality (3.5) with 1-1/n=2-r, there exists C2>0 such that
(4.14)|A1|≤∫Ω2|μ(|ε(u)|)ε(u)-μ(|ε(uh)|)ε(uh)||ε(u-vh)|dV≤C2∫Ω|ε(u-uh)||ε(u-vh)|(1+|ε(u)|+|ε(u-uh)|)2-rdV.
By using the inequality (see Lemma  2.2 in [5] or  (3.10) in [8]),
(4.15)(1+a+t)-ets≤α(1+a+t)-et2+α-1(1+a+s)-es2,∀a,t,s≥0,∀α∈(0,1],e∈(0,1),
with a=|ε(u)|, t=|ε(u-uh)|, s=|ε(u-vh)|, and e=2-r, we obtain, for all α∈[0,1],
(4.16)|A1|≤C2{α‖|u-uh|‖2+α-1‖|u-vh|‖2}.
We now use respectively Hölder's inequality, Young's inequality, and (4.7); there exist C3,C4,C5>0 such that
(4.17)|A2|≤‖div⁡(vh-uh)‖Lr‖p-qh‖Lr′≤12C3β‖vh-uh‖W1,r2+12β-1‖p-qh‖Lr′2≤C4β[1+‖u‖W1,r+‖vh-uh‖W1,r]2-r‖|vh-uh|‖2+12β-1‖p-qh‖Lr′2≤C5β[1+‖uh‖W1,r+‖vh‖W1,r]2-r(‖|u-uh|‖2+‖|u-vh|‖2)+12β-1‖p-qh‖Lr′2.
By setting α=1/(4C2C1) and β=1/(4C5C1[1+∥uh∥W1,r+∥vh∥W1,r]2-r), we obtain
(4.18)‖|u-uh|‖2≤C1{|A1|+|A2|}≤12‖|u-uh|‖2+4C22C12‖|u-vh|‖2+14‖|u-vh|‖2+2C12C5[1+‖uh‖W1,r+‖vh‖W1,r]2-r‖p-qh‖Lr′(Ω)2.
By moving (1/2)|||u-uh|||2 to the left-hand side, we obtain (4.11). From the inf-sup condition (4.1), we have, for all qh∈Qh,
(4.19)Ch‖qh-ph‖Lκ′<sup⁡vh∈Vh⁡∫Ω(qh-ph)div⁡(vh)dV‖vh‖W1,κ.
From (4.10), we have, for all vh∈Vh,
(4.20)∫Ωdiv⁡(vh)(qh-ph)dV=∫Ω2(μ(|ε(u)|)ε(u)-μ(|ε(uh)|)ε(uh)):ε(vh)dV+∫Ωdiv⁡(vh)(qh-p)dV.
From (4.19), (4.20), and (3.5) with 1-1/n=2-r, there exist C6,C7>0 such that
(4.21)Ch‖qh-ph‖Lκ′<sup⁡vh∈Vh|∫Ω2(μ(|ε(u)|)ε(u)-μ(|ε(uh)|)ε(uh)):ε(vh)dV|‖vh‖W1,κ+sup⁡vh∈Vh|∫Ωdiv⁡(vh)(qh-p)dV|‖vh‖W1,κ≤C6{‖2(μ(|ε(u)|)ε(u)-μ(|ε(uh)|)ε(uh))‖Lκ′+‖qh-p‖Lκ′}≤C7{∫Ω(1+|ε(u)|+|ε(u-uh)|)(r-2)κ′|ε(u-uh)|κ′dV}1/κ′+C6‖qh-p‖Lκ′≤C7C8{∫Ω(1+|ε(u)|+|ε(u-uh)|)(r-2)|ε(u-uh)|2dV}1/κ′+C6‖qh-p‖Lκ′≤C7‖|u-uh|‖2/κ′+C6‖qh-p‖Lκ′,
where C8:=∥(1+|ε(u)|+|ε(u-uh)|)(r-2)|ε(u-uh)|2-κ∥L∞1/κ<1. Eventually, the previous inequality together with ∥p-ph∥Lκ′≤∥qh-ph∥Lκ′+∥p-qh∥Lκ′ leads to (4.12).Theorem 4.8.
Let(u,p)∈(V,Q) be the solution of (3.8), and uh∈Vdiv,h the solution of (4.4). For all (vh,qh)∈Vdiv,h×Qh and for all κ∈[r,2], assuming u∈W1,κ(Ω), one has
(4.22)‖u-uh‖W1,r≤D1[1+‖uh‖W1,r+‖vh‖W1,r]2-r(‖u-vh‖W1,κκ/2+‖p-qh‖Lr′),
where D1>0. Moreover, if the spaces Vh and Qh satisfy the inf-sup condition (4.1), then the solution (uh,ph) of (4.2) satisfies for all (vh,qh)∈Vh×Qh and for all κ∈[r,2], assuming u∈W1,κ(Ω),
(4.23)‖u-uh‖W1,r≤D2[1+‖uh‖W1,r+‖vh‖W1,r]2-r(‖u-vh‖W1,κκ/2+‖p-qh‖Lr′),(4.24)‖p-ph‖Lκ′≤D3[1+‖uh‖W1,r+‖vh‖W1,r](2-r)/κ′×{(‖u-vh‖W1,κκ/2+‖p-qh‖Lr′)2/κ′+‖p-qh‖Lκ′},
where D2,D3>0. The constants D1,D2,D3 do not depend on uh and vh; however, D2 and D3 increasingly depends on (Ch)-1.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 v∈Vdiv and for all wh∈Vh, there exists vh∈Vdiv,h such that
(4.25)‖v-vh‖W1,κ≤C‖v-wh‖W1,κ,
where C depends on the inf-sup constant Ch. Eventually, (4.24) follows from (4.12), (4.11), and (4.8).Theorem 4.9.
Assume that, for allκ∈[r,2], there exists a continuous operator πh:[W2,κ]d→Vh that satisfies
(4.26)‖u-πh(u)‖W1,κ≤Ch‖u‖W2,κ,∀u∈[W2,κ]d,
and a continuous operator ρh:W1,κ′→Qh that satisfies
(4.27)‖p-ρh(p)‖Lκ′≤Ch‖p‖W1,κ′,∀p∈W1,κ′,
where h is the size of the higher diameter of the elements of 𝒯h. Assume that Vh and Qh satisfy the inf-sup condition (4.1). Let (u,p) be the solution of problem (3.8) and let (uh,ph) be the solution of problem (4.2). Assume that (u,p)∈([W2,κ]d,W1,κ′), where κ∈[r,2], then one has
(4.28)‖u-uh‖W1,r+(‖p-ph‖Lκ′)κ′/2≤Dhκ/2,
where D=D(∥u∥W2,κ,∥p∥W1,r′,(Ch)-1)>0.Proof.
Apply (4.23) and (4.24) with vh=πh(u) and qh=ρh(p). By using the continuity of πh, (4.5), (4.26), and (4.27), there exist C1,C2,C3,C4>0 such that
(4.29)‖u-uh‖W1,r≤C1[1+‖f‖Lr′+‖u‖W1,r]2-r(h(κ/2)+h)≤C2h(κ/2),‖p-ph‖Lκ′≤C3[1+‖f‖Lr′+‖u‖W1,r](2-r)/κ′{(h(κ/2)+h)(2/κ′)+h}≤C4h(κ/κ′).
The estimate (4.28) directly follows from (4.29).Remark 4.10.
The combination[ℙ1/Bulle]d-ℙ1 for spaces Vh and Qh, 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ΓR≠∅, a similar analysis can be led by replacing the norm defined by (4.6) by(4.30)v⟶‖|v|‖:=∫Ω|ε(v)|2(1+|ε(u)|+|ε(v)|)2-rdV+∫ΓR|v|2(1+|u|+|v|)2-rdS.


## 5. Successive Approximations

In this section, several successive approximation algorithms are proposed for solving the nonlinearity of the discrete problem (4.4) when n>1. For the sake of simplicity, we suppose ΓR=∅ 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 𝒯h is fixed, and we assume that the approximation spaces satisfy Vh⊂V∩[W1,∞(Ω)]d and Qh⊂Q∩L∞(Ω). In what follows, ∥·∥ denotes an arbitrary norm of Vdiv,h. Since Vdiv,h is a finite-dimensional space, all norms are equivalent. Let γ∈[0,1]. We define
(5.1)E:{Vdiv⁡,h⟶Vdiv⁡,h,ũh⟶w̃h,
where w̃h∈Vdiv,h solves
(5.2)2∫Ωμ(|ε(ũh)|)ε(w̃h):ε(vh)dV+2γ∫Ωμ′(|ε(ũh)|)|ε(ũh)|(ε(ũh):ε(w̃h-ũh))(ε(ũh):ε(vh))dV=∫Ωg⋅vhdV,∀vh∈Vdiv⁡,h.The applicationE is well defined. Indeed, by using, respectively, μ′<0, inequalities (3.3) and (3.2), there exist C1,C2>0 such that
(5.3)2∫Ωμ(|ε(ũh)|)|ε(w̃h)|2dV+2γ∫Ωμ′(|ε(ũh)|)|ε(ũh)|(ε(ũh):ε(w̃h))2dV≥2∫Ωμ(|ε(ũh)|)|ε(w̃h)|2dV-2γ(1-1/n)∫Ωμ(|ε(ũh)|)|ε(w̃h)|2dV≥2C1(1-γ(1-1/n))(1+‖ε(ũh)‖L∞)1-1/n∫Ω|ε(w̃h)|2dV≥C2‖w̃h‖2.
As a consequence, the problem (5.2) is coercive. From the Lax-Milgram Theorem, see [15] page 83, there exists a unique solution w̃h∈Vdiv,h of (5.2).In what follows,uh denotes the solution (4.4), which is also the unique fixed point of E. Assume that uh,0 is given; we define iteratively a sequence uh,k, for all k≥1, by
(5.4)uh,k+1=E(uh,k).
Our goal is to prove that uh,k converges to uh when k goes to the infinity. When γ=0, we obtain the classical fixed point method, widely used to solve the nonlinearity of Glen's law; see [6, 9, 14]. When γ=1, we have an additional term in (5.2) which corresponds to Newton's method; see Remark 5.5. The case γ∈(0,1) corresponds to a hybrid fixed point—Newton's method. The convergence of sequence uh,k requires several preliminary results. We compute the first derivative of E in Lemma 5.1. Lemma 5.2 provides an upper bound of the first derivative. Eventually, Theorem 5.3 states the linear convergence of uh,k by using the Banach fixed point theorem. Theorem 5.7 states the second-order convergence when γ=1. By differentiating formally (5.2) at point ũh in direction u̅h, with w̃h=E(ũh), we obtain the following lemma.Lemma 5.1.
Let(ũh,w̃h) satisfy E(ũh)=w̃h. The application E is Gâteaux differentiable at point ũh and its derivative is given by
(5.5)DE(ũh):{Vdiv⁡,h⟶Vdiv⁡,h,u̅h⟶w̅h,
where w̅h solves
(5.6)2∫Ωμ(|ε(ũh)|)ε(w̅h):ε(vh)dV+2∫Ωμ′(|ε(ũh)|)|ε(ũh)|(ε(ũh):ε(u̅h))(ε(w̃h):ε(vh))dV+γ[2∫Ωμ′′(|ε(ũh)|)|ε(ũh)|-μ′(|ε(ũh)|)|ε(ũh)|3(ε(ũh):ε(u̅h))×(ε(ũh):ε(w̃h-ũh))(ε(ũh):ε(vh))dVμ′(|ε(ũh)|)|ε(ũh)|]+γ[2∫Ωμ′(|ε(ũh)|)|ε(ũh)|(ε(u̅h):ε(w̃h-ũh))(ε(ũh):ε(vh))dV]+γ[2∫Ωμ′(|ε(ũh)|)|ε(ũh)|(ε(ũh):ε(w̅h-u̅h))(ε(ũh):ε(vh))dV]+γ[2∫Ωμ′(|ε(ũh)|)|ε(ũh)|(ε(ũh):ε(w̃h-ũh))(ε(u̅h):ε(vh))dV]=0,∀vh∈Vdiv⁡,h.The problems (5.2) and (5.6) have the same coercivity properties to compute w̅h (resp., w̃h) from u̅h (resp., ũh). As a consequence, the problem (5.6) is well-posed by the Lax-Milgram theorem. To prove the convergence of the sequence uh,k, we look for a norm that makes E a contraction at point uh.Lemma 5.2.
Letγ∈[0,1], and let  uh be the fixed point of  E. The application  DE(uh) satisfies
(5.7)‖|DE(uh)|‖μ≤(1-γ)(1-1/n)1-(1-1/n)γ<1,
where |||·|||μ is the subordinated norm to ∥·∥μ:=∫Ωμ(|ε(uh)|)|ε(·)|2dV.Proof.
SinceE(uh)=uh, then (5.6), with ũh=uh and w̃h=uh, is rewriten as
(5.8)∫Ωμ(|ε(uh)|)ε(w̅h):ε(vh)dV+(1-γ)[∫Ωμ′(|ε(uh)|)|ε(uh)|(ε(uh):ε(u̅h))(ε(uh):ε(vh))dV]+γ[∫Ωμ′(|ε(uh)|)|ε(uh)|(ε(uh):ε(w̅h))(ε(uh):ε(vh))dV]=0,
for all vh∈Vdiv,h. From (5.8), μ′<0, and (3.3), we have
(5.9)∫Ωμ(|ε(uh)|)ε(w̅h):ε(vh)dV≤-(1-γ)[∫Ωμ′(|ε(uh)|)|ε(uh)||ε(u̅h)||ε(vh)|dV]-γ[∫Ωμ′(|ε(uh)|)|ε(uh)||ε(w̅h)||ε(vh)|dV]≤(1-γ)(1-1n)[∫Ωμ(|ε(uh)|)|ε(u̅h)||ε(vh)|dV]+γ(1-1n)[∫Ωμ(|ε(uh)|)|ε(w̅h)||ε(vh)|dV].
By setting vh=w̅h in (5.9), we obtain
(5.10)‖w̅h‖μ2≤(1-γ)(1-1n)∫Ωμ(|ε(uh)|)|ε(u̅h)||ε(w̅h)|dV+γ(1-1n)‖w̅h‖μ2.
From (5.10) and Cauchy-Schwarz's inequality, we obtain
(5.11)‖w̅h‖μ2≤[(1-γ)(1-1/n)1-(1-1/n)γ]‖u̅h‖μ‖w̅h‖μ.
Eventually, (5.7) follows from the definition of norm |||·|||μ.Theorem 5.3.
Letγ∈[0,1], let 𝒯h be a given mesh of Ω, and let ∥·∥ be a norm of Vdiv,h. There exist δ>0 and C>0 such that if ∥uh,0-uh∥<δ, then one has
(5.12)‖uh,k-uh‖≤C[(1-γ)(1-1/n)1-(1-1/n)γ]k‖uh,0-uh‖,
for all k≥0, and uh,k is linearly convergent to uh.Proof.
From Lemma5.2, the spectral radius of DE(uh) is lower than constant:
(5.13)[(1-γ)(1-1/n)1-(1-1/n)γ],
which is lower than 1. The theorem is then a direct application of the Banach fixed point theorem in Vdiv,h.Remark 5.4.
In Theorem5.3, it should be stressed that δ depends on h. When h→0 (i.e., if we replace uh by u), we cannot ensure Theorem 5.3 to remain true. Nevertheless, in practise, δ seems to be independent of h, see Section 6.Whenγ=1, we have DE(uh)=0 from Lemma 5.2. It suggests that the convergence of sequence uh,k is quadratic. To establish the second-order convergence, we define, for all vh∈Vdiv,h, the application ℱ(·,vh):Vdiv,h→ℝ by
(5.14)F(ũh;vh):=2∫Ωμ(|ε(ũh)|)ε(ũh):ε(vh)dV-∫Ωf⋅vhdV.Letũh,u̅h,u̿h∈Vdiv,h. We compute formally the first-order derivative of ℱ at point ũh in direction u̅h: (5.15)〈DF(ũh;vh),u̅h〉=2∫Ωμ(|ε(ũh)|)ε(u̅h):ε(vh)dV+2∫Ωμ′(|ε(ũh)|)|ε(ũh)|(ε(ũh):ε(u̅h))(ε(ũh):ε(vh))dV,
and the second-order derivative of ℱ at point ũh in direction (u̅h,u̿h): (5.16)〈D2F(ũh;vh),u̅h,u̿h〉=2∫Ωμ′(|ε(ũh)|)|ε(ũh)|(ε(ũh):ε(u̿h))(ε(u̅h):ε(vh))dV+2∫Ωμ′′(|ε(ũh)|)|ε(ũh)|-μ′(|ε(ũh)|)|ε(ũh)|3×(ε(ũh):ε(u̿h))(ε(ũh):ε(u̅h))(ε(ũh):ε(vh))dV+2∫Ωμ′(|ε(ũh)|)|ε(ũh)|(ε(u̿h):ε(u̅h))(ε(ũh):ε(vh))dV+2∫Ωμ′(|ε(ũh)|)|ε(ũh)|(ε(ũh):ε(u̅h))(ε(u̿h):ε(vh))dV.Remark 5.5.
Ifγ=1,we have, from the definition of E (5.2),ofℱ (5.14), and of Dℱ (5.15),
(5.17)uh,k+1=E(uh,k)⟺F(uh,k;vh)+〈DF(uh,k;vh),uh,k+1-uh,k〉=0,∀vh∈Vdiv⁡,h,
which highlights Newton's method.Lemma 5.6.
The following inequalities hold, for allũh,u̅h,u̿h∈Vdiv,h:
(5.18)2n∫Ωμ(|ε(ũh)|)|ε(u̅h)|2dV≤〈DF(ũh;u̅h),u̅h〉,(5.19)-〈D2F(ũh;u̅h),u̿h,u̿h〉≤∫Ω(8|μ′(|ε(ũh)|)|+2|μ′′(|ε(ũh)|)||ε(ũh)|)|ε(u̿h)|2|ε(u̅h)|dV.Proof.
Inequality (5.18) follows from (5.15) and (3.3), while inequality (5.19) directly follows from (5.16). Computational details are given in [12].Theorem 5.7.
Supposeγ=1, let 𝒯h be a given mesh of Ω, and let ∥·∥ be a norm of Vdiv,h. There exist δ>0 and C>0 such that if ∥uh,0-uh∥<δ, then one has
(5.20)‖uh-uh,k+1‖≤C‖uh-uh,k‖2,
for all k≥0, and uh,k is quadratically convergent to uh.Proof.
Owing to Theorem5.3, there exists δ>0 such that if ∥uh,0-uh∥<δ, then ∥|ε(uh,k-uh)|∥L∞→0 when k→∞. As a consequence, there exists K1 such that ∥|ε(uh,k-uh)|∥L∞<δ for all k>K1. By writing the Taylor expansion of ℱ at point uh,k, there exists ûh,k∈Vdiv,h such that
(5.21)‖|ε(ûh,k-uh,k)|‖L∞<‖|ε(uh-uh,k)|‖L∞<δ,(5.22)F(uh;vh)=F(uh,k;vh)+〈DF(uh,k;vh),uh-uh,k〉+12〈D2F(ûh,k;vh),uh-uh,k,uh-uh,k〉.
Since uh solves (4.4), then ℱ(uh;vh)=0 in (5.22). By setting vh=uh-uh,k+1, we obtain, from (5.22) and (5.17),
(5.23)〈DF(uh,k;uh-uh,k+1),uh-uh,k+1〉=-12〈D2F(ûh,k;uh-uh,k+1),uh-uh,k,uh-uh,k〉.
Thanks to (3.2), there exists C0=D1(1+∥ε(uh)∥L∞)1/n-1 such that 0<C0≤μ(|ε(uh)|). As a consequence, since ∥|ε(uh,k-uh)|∥L∞→0 and ∥|ε(ûh,k-uh)|∥L∞→0 when k→∞ and μ∈𝒞∞(ℝ+*), there exist C1,C2>0 and K2≥K1 such that, for all k>K2,
(5.24)C1<μ(|ε(uh,k)|),(5.25)(8|μ′(|ε(ûh,k)|)|+2|μ′′(|ε(ûh,k)|)||ε(ûh,k)|)<C2.
By applying (5.18) with ũh=uh,k and u̅h=uh-uh,k+1 and (5.24), we obtain
(5.26)C1∫Ω|ε(uh-uh,k+1)|2dV≤〈DF(uh,k;uh-uh,k+1),uh-uh,k+1〉.
By applying (5.19) with ũh=ûh,k, u̅h=uh-uh,k+1, and u̿h=uh-uh,k and (5.25), we obtain:
(5.27)-〈D2F(ûh,k;uh-uh,k+1),uh-uh,k,uh-uh,k〉≤C2∫Ω|ε(uh-uh,k)|2|ε(uh-uh,k+1)|dV.
By combining (5.23), (5.26), and (5.27), we obtain
(5.28)C1∫Ω|ε(uh-uh,k+1)|2dV≤C2∫Ω|ε(uh-uh,k)|2|ε(uh-uh,k+1)|dV.
By using Cauchy-Schwarz's inequality and the equivalence of norms, there exists C3>0 such that
(5.29)‖|ε(uh-uh,k+1)|‖L22≤C3‖|ε(uh-uh,k)|‖L22‖|ε(uh-uh,k+1)|‖L2.
Clearly, (5.20) follows from (5.29).Remark 5.8.
IfΓR≠∅, the nonlinear Robin-Dirichlet condition can be handled in the same way as for the viscosity function. In that case, we modify the application E by adding
(5.30)∑i=1,d-1(∫ΓRα(|ũh|)(w̃h⋅ti)(vh⋅ti)dS+γ∫ΓRα′(|ũh|)|ũh|(ũh⋅(w̃h-ũh))(ũh⋅ti)(vh⋅ti)dS)
to the left-hand side of (5.2). Theorems 5.3 and 5.7 can be easily extended to this case.


## 6. Numerical Results

In this section, numerical experiences are performed in two dimensions (d=2) to validate the results of Theorems 4.9, 5.3, and 5.7. An exact solution of the Stokes problem (2.1) in the square Ω=[0,1]2 is considered in the pure Dirichlet case, that is, ∂Ω=ΓD. Let
(6.1)u(x,y)=((x(1-x))θ+1(y(1-y))θ(1-2y)-(x(1-x))θ(y(1-y))θ+1(1-2x)),θ∈[1,2],
be a divergence-free velocity field that vanishes on the boundaries of Ω. Let p(x,y)=xy-1/4 be the pressure field such that ∫ΩpdV=0. The right-hand side term f is deduced by computing (2.1). The parameters involving in Glen's law are n=2, τ0=0.1 bar and A=0.1bar-2  a-1. Since n=2, the function μ defined by (2.2) is given explicitly, when s>0, by
(6.2)μ(s)=((2Aτ0)2+82As-2Aτ0)(42As).
Numerical solutions are obtained after several successive approximations uh,k, as described in Section 5. Each uh,k corresponds to a unique ph,k. The algorithm is initialised by (uh,0,ph,0):=(0,0). Each linearised problem is solved by using the finite element open source code Freefem++; see [24]. As spaces Vh and Qh, we opt for the combination [ℙ1/Bulle]d-ℙ1 that satisfy the inf-sup condition (4.1) and the interpolation properties (4.26) and (4.27); see Remark 4.10. Six Delaunay unstructured regular meshes 𝒯h of the square Ω=[0,1]2 are generated with various resolutions h. Since the Dirichlet condition u=0 is applied on the whole boundary ∂Ω=ΓD, a penalisation term is added in the variational formulation to constrain the pressure average to be close to zero. For all norms ∥·∥, the error between u and uh,k has two components:
(6.3)‖u-uh,k‖≤‖u-uh‖+‖uh-uh,k‖,
where uh is the exact solution of the nonlinear discrete problem. The convergence of the first component ∥u-uh∥ with respect to h is the concern of Theorem 4.9, while the convergence of the second component ∥uh-uh,k∥ with respect to k is the concern of Theorems 5.3 and 5.7. Let k̅ be an integer large enough such that uh and uh,k̅ can be confused, that is, such that ∥uh-uh,k̅∥≪∥u-uh∥. To check the convergence of the second component, we compute the following error:
(6.4)Eku=‖∇uh,k̅-∇uh,k‖Lr‖∇u‖Lr,
where the norm ∥·∥Lr is evaluated by using the trapezoidal rule. For a fixed h, Theorem 5.3 states the linear convergence of Eku that depends on constant (5.13) when γ∈[0,1) and the quadratic convergence when γ=1. Three values of γ are considered: γ=0 to test the fixed point algorithm, γ=0.5 to test the hybrid method, and γ=1 to test Newton's method. Figure 1 displays Eku according to k for each method: γ∈{0,0.5,1}, and for two different meshes. The recorded orders of convergence are consistent with Theorem 5.3: Newton's method (γ=1) converges quadratically, the fixed point method and the hybrid method (γ<1) converge linearly, and convergence is faster for bigger γ and then smaller constant (5.13). Newton's method is especially very efficient: in our example, only 3 iterations are needed against 8 for the fixed point algorithm to obtain the same accuracy of the numerical solution. Figure 1 also shows that the convergence of Eku with respect to k is not affected by any mesh refinement, as noticed in [3]. Moreover, the addition of supplementary terms in the Stokes system does not increase significantly the computational time for solving the linear system with a direct method.The errorEku with respect to k for each method (fixed point: γ=0, hybrid: γ=0.5, Newton: γ=1). The error Eku obtained with the coarsest mesh is displayed on (a) and the error Eku obtained with the finest mesh is displayed on (b).
(a)(b)The estimate of Theorem4.9 is now tested by computing the following errors:
(6.5)Ehu=‖∇u-∇uh,k̅‖Lr‖∇u‖Lr,Ehp=‖p-ph,k̅‖Lκ′‖p‖Lκ′.
We can change the regularity of u by changing the parameter θ in (6.1) from 2 to 1.34. Indeed, if θ=2, then u∈𝒞∞(Ω), while if θ=1.34, then u∉[W2,2(Ω)]2, but u∈[W2,r(Ω)]2, where r=3/2. In any case p∈𝒞∞(Ω). Figure 2 displays Ehu and Ehp with respect to h in both cases θ=2 and θ=1.34. The estimate (4.28) anticipates Ehu=𝒪(h), Ehp=𝒪(h) if θ=2, and Ehu=𝒪(h3/4), Ehp=𝒪(h1/2) if θ=1.34. In both cases, the observed order of convergence for Ehu and Ehp is close to one, which is greater or equal to the estimate. It suggests the nonoptimality of estimate (4.28) in the nonregular case, as noticed in [7] for a comparable problem.The errorsEhu and Ehp with respect to h with θ=2 (a) and with θ=1.34 (b). Regression straight lines have been drawn from the last three points of each set of six recorded errors.
(a)(b)


## 7. Conclusions and Perspectives

We have proved the existence and the uniqueness of a weak solution of a nonlinear Stokes problem that describes the motion of glaciers. We have also proved the convergence of the finite element approximation and givena priori error estimates. New successive approximation algorithms have been proposed to solve the system nonlinearity and have been proved to be convergent. When implementing Newton's method, both theoretical and numerical studies have shown the efficiency of this method in comparison with the classical fixed point method.Two extensions of our work should be investigated in future research. First,a posteriori estimates could be an aspect to be developed in order to implement an adaptive mesh procedure. Second, the presented Stokes model could benefit from recent improvements of the basal sliding description with Coulomb-type laws [10].



---

*Source: 164581-2012-01-12.xml*