Guillaume Jouvet

A multilayer

come

Introduction

To better evaluate sea level rise [START_REF] Vaughan | Why is it hard to predict the future of ice sheets ?[END_REF] in a climate change regime, it is crucial to develop ice flows models that are mechanically complete while being computationally-tractable. Ice is known to behave like a non-Newtonian fluid, governed by Glen's flow law [START_REF] Glen | The flow law of ice[END_REF]. Thus, the velocity and the pressure of the ice satisfy non-linear Stokes equations. In practise, solving these equations requires considerable computational resources and complex meshing procedures to be implemented at the large scale, which is necessary when modelling ice sheets. As a consequence, a number of simplified models have been proposed in the last decades to lead simulations of real glaciers. Simplifications are always based on the small aspect ratio (quotient between the characteristic heights and lengths) of glaciers. More precisely, after writing the solution of the Stokes equations as an asymptotic expansion in the aspect ratio ǫ, high order terms in ǫ are neglected. Two categories of models derive from such an approach. In the first category, each order that was not neglected gives rise to a system of equations: the Shallow Ice Approximation (SIA) for ǫ 0 [START_REF] Hutter | Theoretical Glaciology[END_REF], the first order SIA for ǫ 1 , and the second order SIA for ǫ 2 [START_REF] Baral | Asymptotic theories of large-scale motion, temperature and moisture distribution in land-based polythermal ice sheets: A critical review and new develop-ments[END_REF][START_REF] Ahlkrona | A numerical study of scaling relations for non-newtonian thin-film flows with applications in ice sheet modelling[END_REF][START_REF] Egholm | Modeling the flow of glaciers in steep terrains: The integrated second-order shallow ice approximation (isosia)[END_REF], which can be solved iteratively at a much lower cost than the original Stokes problem. In contrast, the models of the second category exploit the lost of higher-order terms by eliminating unknown variables and by reducing the dimension of the mathematical model. In this last category, the first level of simplification assumes hydrostatic vertical normal stresses [START_REF] Greve | Dynamics of Ice Sheets and Glaciers[END_REF], with the simplification that the pressure variable is eliminated from the Stokes equations. However, this hydrostatic approximation is almost never used in practise. Instead, the First Order Approximation (FOA) [START_REF] Blatter | Velocity and stress fields in grounded glaciers: a simple algorithm for including deviatoric stress gradients[END_REF][START_REF] Pattyn | A new three-dimensional higher-order thermomechanical ice sheet model: Basic sensitivity, ice stream development, and ice flow across subglacial lakes[END_REF] model, which further assumes negligible horizontal derivatives of the vertical velocity compared to vertical derivatives of the horizontal velocity, is usually preferred. Compared to the hydrostatic approximation, the third component of the velocity has vanished from the FOA system. Despite the fact that the unknowns of the FOA are reduced to the horizontal components of the velocity, the FOA model is still mathematically 3D. As a consequence, solving the FOA still requires meshing complex and shallow geometries, which change in time.

To remove such complexity, the dimension of the mathematical model can be reduced by further mechanical simplifications. For instance, the Shallow Shelf Approximation (SSA) [START_REF] Morland | Unconfined ice-shelf flow[END_REF][START_REF] Macayeal | Large-scale ice flow over a viscous basal sediment: Theory and application to ice stream b, antarctica[END_REF], which accounts only for longitudinal (or membrane) stresses, is 2D after being depth-integrated. In contrast, the ǫ 0 SIA, which accounts only for vertical shear stresses, reduces to a 1D (vertical) mathematical model, independently in each column of ice. In the literature [START_REF] Schoof | Thin-film flows with wall slip: An asymptotic analysis of higher order glacier flow models[END_REF], the SSA is called a "membrane" model while the SIA is called a "lubrication" model. Those two models (SIA and SSA) are popular to describe the dynamics of ice sheets and ice shelves since the size and the complexity of the system to solve are definitely reduced compared to any 3D models.

In practise, the vertical shear components of the stress tensor are significant where ice is grounded while the longitudinal components are dominant in floating parts. As a consequence, both components must be combined if one wants to model the entire domain (grounded and floating). This has motivated the construction of "hybrid" models, which account for both kind of stresses, while being mathematically 2D. The simplest hybrid model consists of the linear combination SIA+SSA, which consists of summing the velocities of each model [START_REF] Bueler | Shallow shelf approximation as a "sliding law" in a thermomechanically coupled ice sheet model[END_REF]. Unfortunately, this model does not include the simultaneous coupling between the vertical shear and the longitudinal stresses. As a result, this model can not capture the 3D ice flows that occur in deep and narrow valleys or in the vicinity of grounding lines (Pattyn and others, 2013;[START_REF] Jouvet | An adaptive newton multigrid method for a model of marine ice sheets[END_REF]. In contrast, the L1L2 [START_REF] Hindmarsh | A numerical comparison of approximations to the stokes equations used in ice sheet and glacier modeling[END_REF] or some variants like the ones proposed in [START_REF] Pollard | A Coupled Ice-Sheet/Ice-Shelf/Sediment Model Applied to[END_REF], [START_REF] Schoof | Thin-film flows with wall slip: An asymptotic analysis of higher order glacier flow models[END_REF] or in [START_REF] Goldberg | A variationally derived, depthintegrated approximation to a higher-order glaciological flow model[END_REF] include the vertical shear stress in the computation of the effective viscosity of the SSA. All these hybrid models have in common that they solve a single non-linear elliptic 2D problem, and that the velocity profile is reconstructed a posteriori via an implicit relation [START_REF] Schoof | Thin-film flows with wall slip: An asymptotic analysis of higher order glacier flow models[END_REF][START_REF] Winkelmann | The potsdam parallel ice sheet model (PISM-PIK) part 1: Model description[END_REF][START_REF] Cornford | Adaptive mesh, finite volume modeling of marine ice sheets[END_REF]. The hierarchy of the aforementioned models is drawn in Figure 1.

Stokes 3D FOA 3D SIA 1D SIA 1D +SSA 2D Multilayer 2D "L1L2-like" 2D SSA 2D
Hybrid Models In this paper, a new hybrid model generalizing the SSA is introduced. The SSA model assumes a verticallyconstant velocity profile, such that it only accounts for the longitudinal components of the stress while neglect-ing the vertical components. To recover these components, the velocity profile of the new model partitioned and then assumed to be vertically piecewise-constant.

Inspired by an ocean model [START_REF] Audusse | A multilayer saint-venant system with mass exchanges for shallow water flows. derivation and numerical validation[END_REF], this approach consists of seeing the thickness of ice as a pile of thin layers, which can spread out, tighten and can slide over each others. Similarly to the SSA, the model is obtained by integrating vertically the FOA model, but locally on each layer. The boundary terms appearing when integrating give rise to interlayer tangential stresses. To account for sliding between layers, only the vertical shear components in the local frame are kept. These components redefine the interlayer tangential stresses. The final model consists of a system of 2D non-linear elliptic equations of size the number of layers. This multilayer1 model naturally generalizes the SSA, which corresponds to the 1-layer case of the model. By construction, the multilayer model is expected to be as accurate as any higher-order model when increasing the number of vertical layers. However, in contrast to the FOA or the Stokes model, the multilayer model is mathematically 2D, and then much easier to solve. Moreover, any solver that has been developed for the SSA can be extended to solve the multilayer model. To test the multilayer model against other higher-order models, the 2D velocity fields of two benchmark experiments are computed for the experiment B of the ISMIP-HOM project (Pattyn and others, 2008) and for one real mountain glacier.

This paper is organized as follows. The model is derived in Section 2. Then, two numerical methods to solve the resulting system are described in Section 3. Finally, the numerical results are reported in Section 4.

Model derivation

In this section, a generic 3D system of ice sheet and ice shelf is considered. The most complex ice flow model, non-linear Stokes, and the First Order Approximation (FOA), which is a simplification, are described in Section 2.1 and 2.2, respectively. Then an integration procedure derives the multilayer model in Section 2.3. The redefinition of the interlayer stress terms is reported in Section 2.4 while the boundary conditions are rewritten in the multilayer setting in Section 2.5. The multilayer model is reformulated in a vectorial way in Section 2.6. Finally, the vectorial problem is rewritten as a variational problem and a minimisation problem in Section 2.7.

Let V be a three-dimensional domain of ice, which is defined by

V = {(x, y, z), s(x, y) ≤ z ≤ s(x, y)}, (1)
where (x, y) denote the horizontal coordinates, z denotes the vertical coordinate, s(x, y) and s(x, y) are the elevations of the lower and upper ice surfaces. Call b(x, y) the elevation of the bedrock. Note that s = b holds where ice is grounded and s > b where ice is floating. The flotation of ice is driven by Archimedes principle,

s = max b, - ρ ρ w h , (2)
where h := s -s is the ice thickness and ρ and ρ w are the densities of ice and water, respectively ; see Figure 2. Relation (2) says that if the buoyancy -ρ w gb is less than the force exerted by ice ρgh, then ice is grounded, otherwise ice is floating and ρ/ρ w of the ice thickness is below sea level.

Bedrock Water

Water Level The boundary of V is divided into the upper interface

Γ s = {(x, y, z), z = s(x, y)}, (3)
the lower interface

Γ f ∪ Γ n ∪ Γ l = {(x, y, z), z = s(x, y)} (4)
and a possible vertical ice-cliff at the calving front Γ c , see Figure 2. At the lower interface, ice might be frozen to the ground, sliding on the ground or floating on the water. For this reason, the cases are distinguished as follows: Γ f denotes the non-sliding part, Γ n the sliding grounded part, and Γ l the floating part, respectively.

Additionally Ω and Ω k denote the projections of V and its boundaries Γ k on the horizontal plane (Oxy) for k ∈ {s, f, n, l, c}.

In what follows, the velocity field of the ice fluid in V is denoted by u = (u x , u y , u z ), the pressure field is denoted by p, the derivative for variable i is denoted by ∂ i (i ∈ {x, y, z}), and Einstein summation is adopted.

Stokes

The Stokes model consists of the momentum conservation equation and the incompressibility condition:

-∂ j σ ij = ρg i , in V, (5) ∂ i u i = 0, in V, (6)
where σ ij and g i are the components of the Cauchy stress tensor and (g x , g y , g z ) = (0, 0, -g), where g is the gravitational constant. Call τ ij the components of the deviatoric stress tensor defined by

σ ij = τ ij -pδ ij , (7)
where δ ij is the Kronecker symbol. Glen's flow law [START_REF] Glen | The flow law of ice[END_REF], which describes the mechanical behaviour of ice, consists of the following non-linear viscosity relation:

τ ij = 2µ εij , (8)
where εij denotes the components of the strain rate tensor defined by εij = 1 2 (∂ j u i + ∂ i u j ), µ is the viscosity defined by

µ = 1 2 A -1 n 1 2 εij εji 1 2 ( 1 n -1) , (9)
A > 0 and n ≥ 1 are two constants called the rate factor and Glen's exponent, respectively. In reality, A is not constant since it depends on ice temperature [START_REF] Greve | Dynamics of Ice Sheets and Glaciers[END_REF]. However, for simplicity, it is assumed in this paper that the ice is isothermal.

The boundary conditions that supplement (5), (6) (7), (8), and (9) are the following. No force applies on the ice-air interface,

σ ij n j = 0, on Γ s , (10)
where

n = (n x , n y , n z ) T = (-∂ x s, -∂ y s, 1) T (11)
is an outer normal vector along Γ s . Along the lower surface interface, the no-slip condition is

u i = 0, on Γ f , (12)
the non-linear friction condition reads [START_REF] Hutter | Theoretical Glaciology[END_REF]:

u i n i = 0, on Γ n , (13)
τ ij n j t k i = -C| u| 1 m -1 u i t k i , on Γ n , (14)
for k ∈ {x, y}, where m > 0, C = C(x, y) > 0, and the condition on the floating interface reads [START_REF] Greve | Dynamics of Ice Sheets and Glaciers[END_REF]:

σ ij n i = ρ w g z n i , on Γ l , (15)
for k ∈ {x, y}, where

n = (n x , n y , n z ) T = (∂ x s, ∂ y s, -1) T 1 + (∂ x s) 2 + (∂ y s) 2 , ( 16
)
t x = (t x x , t x y , t x z ) T = (1, 0, ∂ x s) T , (17)
t y = (t y x , t y y , t y z ) T = (0, 1, ∂ y s) T , (18)
are the outward normal unit vector and two orthogonal vectors tangent to the boundaries Γ n and Γ l , respectively. Finally, the same condition applies at the calving front below sea level:

σ ij n i = ρ w g min(z, 0) n i , on Γ c . ( 19
)
where (n x , n y , n z ) = (n x , n y , 0) is an outer normal vector to Γ c .

First Order Approximation (FOA)

Call ǫ = [h]/[x] the aspect ratio of V , where [h] and [x] denote its typical height and length. A dimensionless scaling [START_REF] Blatter | Velocity and stress fields in grounded glaciers: a simple algorithm for including deviatoric stress gradients[END_REF][START_REF] Schoof | Thin-film flows with wall slip: An asymptotic analysis of higher order glacier flow models[END_REF] shows that

∂ j σ zj = ∂ z σ zz + O(ǫ 2 ) in V, (20)
σ zj n j = σ zz n z + O(ǫ 2 ), on Γ s ∪ Γ f ∪ Γ n ∪ Γ l , (21)
τ zj n j = τ zz n z + O(ǫ 2 ), on Γ s ∪ Γ f ∪ Γ n ∪ Γ l . (22) and εiz = 1 2 ∂ z u i + O(ǫ 2 ), i ∈ {x, y} in V. ( 23
)
From now on, the remainders O(ǫ 2 ) in ( 20), ( 21), ( 22), and ( 23) are neglected such that the Stokes problem and its boundary conditions simplify. Indeed, using (20), the third equation of (5) simplifies into

∂ z σ zz = ρg, in V, (24)
while, using (21), the third equation of (10) becomes

σ zz = 0, on Γ s . ( 25
)
Integrating vertically ( 24) with (25) yields

σ zz = τ zz -p = -ρg(s -z), in V. ( 26
)
By ( 6) and (8), it follows that

p = ρg(s -z) -τ xx -τ yy , in V. ( 27
)
Thus p can be eliminated from the two first equations of (5):

∂ x (2τ xx + τ yy ) + ∂ y σ xy + ∂ z σ xz = ρg∂ x s, (28) ∂ x σ xy + ∂ y (2τ yy + τ xx ) + ∂ z σ yz = ρg∂ y s, (29)
and from the stress-free boundary condition (10):

(2τ xx + τ yy )n x + (σ xy )n y + σ xz n z = 0, on Γ s , ( 30
) (σ xy )n x + (2τ yy + τ xx )n y + σ yz n z = 0, on Γ s . (31)
Using ( 16), ( 17) and ( 18), the simplification due to ( 22), the friction condition ( 14) becomes:

(2τ xx + τ yy )n x + (σ xy )n y + σ xz n z = -C| u| 1 m -1 u i t x i , on Γ n , (32) (σ yx )n x + (2τ yy + τ xx )n y + σ yz n z = -C| u| 1 m -1 u i t y i , on Γ n . (33)
Using ( 21), ( 26) and s = -ρ ρw h, which derives from the floating condition (2), the third equation of ( 15) becomes

σ zz = τ zz -p = ρ w gs = -ρgh, on Γ l . ( 34
)
Again using ( 27) to eliminate the pressure from the two first equations of (15) and using ( 34) imply

(2τ xx + τ yy )n x + (σ xy )n y + τ xz n z = 0, on Γ l , ( 35
) (σ yx )n x + (2τ yy + τ xx )n y + τ yz n z = 0, on Γ l . (36)
Similarly, the condition at the calving front (19) becomes:

(2τ xx + τ yy )n x + (σ xy )n y = (ρ w g min(z, 0) + ρg (s -z)) n x , on Γ c , (37)
(σ yx )n x + (2τ yy + τ xx )n y = (ρ w g min(z, 0) + ρg (s -z)) n y , on Γ c . (38)
In addition, (23) says that the horizontal derivatives of the vertical velocities are small compared to the vertical derivatives of the horizontal velocities. Consequently, using and the incompressibility ( 6), (9) becomes:

µ = 1 2 A -1 n 1 2 (∂ x u x ) 2 + 1 2 (∂ y u y ) 2 + 1 2 (∂ x u x + ∂ y u y ) 2 + 1 4 (∂ y u x + ∂ x u y ) 2 + 1 4 (∂ z u x ) 2 + 1 4 (∂ z u y ) 2 1 2 ( 1 n -1)
.

(39)

Depth integration over layers

The domain of ice is now divided in the vertical direction into L layers of thickness h 1 , ..., h L such that l=1,...,L

h l = h, (40)
see Figure 3. Call s l = s + h l the elevation of the upper surface of the layer l for l = 0, ..., L, with the convention h 0 = 0. The choice of this vertical division is discussed later. The derivation of the SSA model is based on the assumption of a constant velocity profile [START_REF] Macayeal | Large-scale ice flow over a viscous basal sediment: Theory and application to ice stream b, antarctica[END_REF]. Instead, here u is assumed to be vertically piecewise constant, equal to u l (x, y) on layer l:

u(x, y, z) = l=1,...,L u l (x, y)1 (s l-1 ,s l ] (z), (41)
where 1 I (z) equals 1 if z ∈ I and 0 otherwise, see Figure 3. The discontinuities of the velocity lead to undefined stresses between the layers. The redefinition of such stresses is addressed in Section 2.4. Consider an arbitrary layer indexed by l ∈ {1, ..., L}. Using Leibnitz's rule, the integration of ( 28) vertically over the layer l yields

u l u l+1 h L h l+1 h l h 1 s l-1 s l s l+1 s L = s s 0 = s
2∂ x s l s l-1 τ xx dz -2[τ xx ] z=s l ∂ x s l + 2[τ xx ] z=s l-1 ∂ x s l-1 +∂ x s l s l-1 τ yy dz -[τ yy ] z=s l ∂ x s l + [τ yy ] z=s l-1 ∂ x s l-1 +∂ y s l s l-1 σ xy dz -[σ xy ] z=s l ∂ y s l + [σ xy ] z=s l-1 ∂ y s l-1 +[σ xz ] z=s l -[σ xz ] z=s l-1 = ρgh l ∂ x s,
where [•] z=s l (resp. [•] z=s l-1 ) stands for the limit z → s l (resp. z → s l-1 ) with z < s l (resp. z > s l-1 ). From ( 7), ( 8), ( 9) and using the fact that u = u l is constant in the layer l, the latter becomes:

2∂ x µ l h l 2∂ x u l x + ∂ y u l y + ∂ y µ l h l ∂ y u l x + ∂ x u l y + Σ l,0 x,y + Σ l,-1 x,y = ρgh l ∂ x s, (42)
where

µ l = 1 2 A -1 n 1 2 ∂ x u l x 2 + 1 2 ∂ y u l y 2 (43) + 1 2 ∂ x u l x + ∂ y u l y 2 + 1 4 ∂ y u l x + ∂ x u l y 2 1 2 ( 1 n -1)
, and, for k = 0, -1,

Σ l,k x,y = (-1) -k α l (44) (2τ xx + τ yy )n l+k x + σ xy n l+k y + σ xz n l+k z z=s l+k ,
where

n l = (n l x , n l y , n l z ) T = (∂ x s l , ∂ y s l , -1) T α l , (45)
is the outer normal unit vector to the upper boundary of layer l, and

α l = 1 + (∂ x s l ) 2 + (∂ y s l ) 2 . ( 46
)
Similarly, integrating (29) vertically over the layer l ∈ {1, ..., L} leads to

2∂ y µ l h l 2∂ y u l y + ∂ x u l x + ∂ x µ l h l ∂ y u l x + ∂ x u l y + Σ l,0 y,x + Σ l,-1 y,x = ρgh l ∂ y s. ( 47
)
Define the vector

Σ l,k := Σ l,k x,y Σ l,k y,x
, fork = 0, -1.

Using the simplification due to ( 22), one can verify that Σ l,0 and Σ l,-1 correspond to the tangential stresses at the limit of the top and the bottom of the layer l:

Σ l,k = (-1) -k α l τ ij n l+k i t x,l+k j τ ij n l+k i t y,l+k j z=s l+k (48)
where n l is defined by (45) and

t x,l = (t x,l x , t x,l y , t x,l z ) T = (1, 0, ∂ x s l ) T (49)
t y,l = (t y,l x , t y,l y , t y,l z ) T = (0, 1, ∂ y s l ) T (50)

are two orthogonal vectors tangent to the upper boundary of layer l.

On the one hand, the continuity of the stress across the layers implies:

Σ l,0 = -Σ l+1,-1 , ∀l = 1, ..., L -1, (51)
thus that the layers are coupled. On the other hand, assuming that the layers are oriented in the direction of the flow, the layers can slide on each other such that only the shear components are significant at the interface between layers. Based on this statement, Σ l,0 and -Σ l+1,-1 are redefined by the quantity in equation ( 68) in the next section.

Interlayer stresses

From ( 48), ( 8) and ( 9), the k-th component of Σ l,0 and -Σ l+1,-1 should be equal to

A -1 n α l 1 2 εij εji 1 2 ( 1 n -1) εij n l j t k,l i . (52)
However, ε is not defined between the layers because of the discontinuity of the velocity field (41). In order to make (52) meaningful and to redefine Σ l,0 and -Σ l+1,-1 , several hypothesis are made. The key assumption is that the multilayer vertical splitting ( 40) is chosen such that the layers are aligned with the direction of the flow:

(∂ j u i )n l i = 0. ( 53
)
Since the layers can slide on each other, the vertical shear components dominate in the stress expressed in the local frame tangential to the layer boundary. As a consequence, the longitudinal stresses in the local frame are set zero:

εij t k,l j t k,l i = 0, ∀(k, k) ∈ {x, y}. (54)
On the one hand, (53) implies

εij n l j = 1 2 (∂ i u j )n l i , (55)
while, on the other hand, (54) yields

1 2 εij εji = ( εij n l j ) 2 . ( 56
)
Additionally, (∂ i u j )n l i is the derivative of u in the direction n l :

(∂ i u j )n l i = ∂ z u j , ( 57
)
where z is the local variable defined by z = n l x x + n l y y + n l z z. In (57), the derivative with respect to z (according to the direction orthogonal to the interface) is approximated by the finite difference:

∂ z u k 2 = u l+1 k -u l k hl+1 + hl , k ∈ {x, y, z}. (58)
It should be stressed that (58) depends on the local coordinate z and not on z. However, (53) implies that the Taylor expansion of u l k at z has no first-order term,

u l k (x, y, z) = u l k (z) + O([u l ]δ 2 ǫ 2 /L 2 ), (59)
for k ∈ {x, y, z}. By neglecting the remainder, u l k is considered as locally constant in the plan-orthogonal to n l , so (58) holds in the primary variables (x, y, z) too. In addition, calling hl the thickness function in the local frame, the Taylor expansion of hl and hl+1 at z yields

hl = 1 + ∇s l • ∇h l 1 + |∇s l | 2 h l + O([h l ]δ 3 ǫ/L), (60) hl+1
= 1 -∇s l • ∇h l+1 1 + |∇s l | 2 h l+1 + O([h l+1 ]δ 3 ǫ/L), (61)
where δ = [s]/[x] and [s] denotes the scales for variations in surface elevation, see Figure 4. In the following, the O(δ 2 ) terms are kept in ( 60) and ( 61), but the higherorder remainders are neglected. Finally, using ( 52), ( 55), ( 56), ( 57), ( 58), ( 59), ( 60) and ( 61), the tangential stress is redefined by

S l k = A -1 n α l α l β l 1 n u l+1 -u l 1 n -1 2 u l+1 -u l • t l,k (62) Zoom + rotation Local frame Global frame t x n i t x s l hl+1 n h l+1 h l j O([h l+1 ]δ 3 ǫ/L) 1 √ 1+|∇s l | 2 h l+1 -∇s l •∇h l √ 1+|∇s l | 2 h l+1
β l = (h l + h l+1 ) + ∇s l • (h l ∇h l -h l+1 ∇h l+1 ), ( 63
)
which simplifies into

β l = h l + h l+1 , (64)
if the multilayer splitting (40) is chosen uniform, i.e. if h l = h/L. Finally, (53) implies

(u l+1 i -u l i )n l i = 0,
so that the third component can be eliminated:

u l+1 -u l 2 2 = u l+1 x -u l x u l+1 y -u l y T M l u l+1 x -u l x u l+1 y -u l y , (65)
u l+1 -u l • t l,x t l,y = M l u l+1 x -u l x u l+1 y -u l y , (66)
where

M l = 1 + (∂ x s l ) 2 (∂ x s l )(∂ y s l ) (∂ x s l )(∂ y s l ) 1 + (∂ y s l ) 2 . ( 67
)
By combining (62), ( 65), (66), S l becomes:

S l = A -1 n α l α l β l 1 n × u l+1 x -u l x u l+1 y -u l y ( 1 n -1) M l M l u l+1 x -u l x u l+1 y -u l y , (68)
where the norm | • | M l derives from the scalar product

|(u x , u y )| 2 M l = (u x , u y )M l (u x , u y ) T . (69)
However, S l can be more simply redefined by neglecting O(δ 2 ) terms in (68). Doing so yields

α l = 1, M l = I, (70)
and

S l = A -1 n × u l+1 x -u l x h l+1 + h l 2 + u l+1 y -u l y h l+1 + h l 2 1 2 ( 1 n -1) × u l+1 x -u l x h l+1 + h l , u l+1 y -u l y h l+1 + h l . ( 71
)
The multilayer models, which are based on (71), instead of ( 68), are later labelled with an asterisk in exponent, e.g. 16-layer * . The approximation of the stress tensor in ( 71) is similar to the one involved when deriving the SIA [START_REF] Greve | Dynamics of Ice Sheets and Glaciers[END_REF]; see also Appendix B. In fact, S l defined in (71) corresponds to the vertical stresses in the global frame, assuming the other components to be zero. If the surface gradients are small (i.e. if δ 2 is small), then the global frame is close enough to the local frame so that S l can be redefined by (71) instead of (68). In contrast, if these gradients are not negligible, the slope at the layer interface must be accounted, for and so S l is better redefined by (68); see Section 4.

Boundary conditions

The boundary condition at the top of the highest layer, i.e. Σ L,0 , is first considered. Returning to (44) and the free-stress condition (30) (31), it follows that

Σ L,0 = 0. ( 72
)
The boundary conditions at the bottom of the lowest layer, i.e. Σ 1,-1 , are now considered. On Ω f , no sliding occurs below the lowest layer and the approach of Section 2.4 can be applied, i.e. Σ 1,-1 is redefined by

S 0 = A -1 n α 0 α 0 β 0 1 n u 1 x u 1 y 1 n -1 M 0 M 0 u 1 x u 1 y . (73)
On the sliding part Ω n , the conditions ( 32) and ( 33) with (44) become

Σ 1,-1 = C α 0 | u 1 | 1 m -1 u 1 i t 0,x i u 1 i t 0,y i on Ω n . ( 74
)
It remains to rewrite (74) in the horizontal velocity components. Since the velocity field is tangential to ( t 0,x , t 0,y ) (condition (13)), it follows, like in Section 2.4,

Σ 1,-1 = C α 0 u 1 x u 1 y 1 m -1 M 0 M 0 u 1 x u 1 y on Ω n . (75)
Finally, the conditions ( 35) and ( 36) on the floating part Ω l become

Σ 1,-1 = 0, on Ω l . ( 76
)
At the calving front Ω c , (37), ( 38) are integrated from s l-1 to s l , and the following boundary condition is obtained:

2µ l h l 2∂ x u l x + ∂ y u l y n x + µ l h l ∂ y u l x + ∂ x u l y n y = F l n x , (77) µ l h l ∂ y u l x + ∂ x u l y n x + 2µ l h l 2∂ y u l y + ∂ x u l x n y = F l n y , ( 78
)
where

F l = 1 2 ρg (s -s l ) 2 -(s -s l+1 ) 2 + 1 2 ρ w g (min(s l+1 , 0)) 2 -(min(s l , 0)) 2 . ( 79
)

Vectorial reformulation

From now on, Σ l,0 and -Σ l+1,-1 defined by ( 48) are replaced by S l defined by (68) (or (71) for the multilayer * model). Additionally, the notation u = (u x , u y ) is used in order to rewrite ( 42) and ( 47) in the compact form

-A -1 n ∇ • h l D(u l ) 1 n -1 * [D(u l ) + tr(D(u l ))I] -S l + S l-1 = -ρgh l ∇s, (80)
where D(u) := 1 2 ∇u + ∇u T denotes the strain-rate of u, ∇, ∇• denote gradient and divergence, respectively, with respect to the horizontal variables (x, y), tr is the trace operator, I the identity second order tensor and |•| * denotes the norm |X| * := (X, X) * associated with the scalar product defined by (X, Y ) * := 1 2 (tr(XY ) + tr(X)tr(Y )) .

Now, ( 73), ( 75) and ( 76) can be summarized by

S 0 = A -1 n α 0 α 0 β 0 1 n u 1 1 n -1 M 0 M 0 u 1 × 1 Ω f + Cα 0 u 1 1 m -1 M 0 (M 0 u 1 ) × 1 Ωn , (81)
where 1 R (x, y) equals to 1 if (x, y) ∈ R and 0 otherwise. Finally, the multilayer solution (u 1 , ..., u L ) solves the following 2 × 2-block tridiagonal system of equations:

-A -1 n ∇ • h L D(u L ) 1 n -1 * [D(u L ) + tr(D(u L ))I] + A -1 n α L-1 α L-1 β L-1 1 n u L -u L-1 1 n -1 M L-1 M L-1 (u L -u L-1 ) = -ρgh L ∇s, (82)
for all l ∈ {2, ..., L -1}:

-A -1 n ∇ • h l D(u l ) 1 n -1 * [D(u l ) + tr(D(u l ))I] + A -1 n α l-1 α l-1 β l-1 1 n u l -u l-1 1 n -1 M l-1 × M l-1 (u l -u l-1 ) + A -1 n α l α l β l 1 n u l -u l+1 1 n -1 M l M l (u l -u l+1 ) = -ρgh l ∇s, (83)
and

-A -1 n ∇ • h 1 D(u 1 ) 1 n -1 * [D(u 1 ) + tr(D(u 1 ))I] + A -1 n α 0 α 0 β 0 1 n u 1 1 n -1 M 0 M 0 u 1 × 1 Ω f + Cα 0 u 1 1 m -1 M 0 (M 0 u 1 ) × 1 Ωn (84) + A -1 n α 1 α 1 β 1 1 n u 1 -u 2 1 n -1 M 1 M 1 (u 1 -u 2 ) = -ρgh 1 ∇s.
At the calving front, (77) (78) can be rewritten as:

A -1 n h l |D(u l )| 1 n -1 * [D(u l ) + tr(D(u l ))I] • n := F l n, (85)
where n denotes an horizontal outward normal vector to Ω c . It is interesting to notice the similarity of the system (82) (83) (84) with the vectorial equation of the SSA [START_REF] Macayeal | Large-scale ice flow over a viscous basal sediment: Theory and application to ice stream b, antarctica[END_REF][START_REF] Schoof | A variational approach to ice stream flow[END_REF][START_REF] Jouvet | An adaptive newton multigrid method for a model of marine ice sheets[END_REF], which corresponds to the 1-layer model (i.e. when L = 1). When neglecting the O(δ 2 ) components in the friction term of (84), it reduces to the common expression |u 1 | 1 m -1 u 1 × 1 Ωn , see e.g. [START_REF] Cornford | Adaptive mesh, finite volume modeling of marine ice sheets[END_REF]. In order to get a better overview of the multilayer model, the 2-layer * model is rewritten in Appendix A in a much simpler setting, which assumes a 2D flow in the vertical x -z plane, such that no dependency on the traverse coordinate y occurs. Additionally, the multilayer model is rewritten in Appendix B in the "infinite parallel-sided slab" setting. In particular, it is shown that the multilayer solution (101) converges to the exact solution of the Stokes equations when refining the vertical multilayer splitting (40). Interestingly, this convergence does not hold with the FOA. This shows that the reconstruction of the interlayer stress of Section 2.4 recovers some mechanical components that are neglected in the FOA.

Variational and minimisation problems

To analyse the multilayer system (82) (83) (84) and to implement a finite element method, it must be rewritten as a variational problem. For that, the equations of the system (82) (83) (84) are multiplied by a test function v l , summed, integrated-by-parts using (85), to give

A -1 n l=1,...,L Ω h l |D(u l )| 1 n -1 * (D(u l ), D(v l )) * dΩ + A -1 n Ω f α 0 α 0 β 0 1 n u 1 1 2 ( 1 n -1) M 0 M 0 u 1 • v 1 dΩ + C Ωn α 0 u 1 1 m -1 M 0 M 0 u 1 • v 1 dΩ + A -1 n l=2,...,L Ω α l-1 α l-1 β l-1 1 n u l -u l-1 1 n -1 M l-1 M l-1 (u l -u l-1 ) • (v l -v l-1 )dΩ = -ρg Ω l=1,...,L h l ∇s • v l dΩ + Ωc l=1,...,L F l n • v l dS. ( 86
)
One can verify that ( 86) is as the Euler-Lagrange equation DJ (u), v = 0 where u = (u 1 , ..., u L ) and v = (v 1 , ..., v L ), for the functional

J (u) = A -1 n 1 n + 1 l=1,...,L Ω h l |D(u l )| 1 n +1 * dΩ (87) + A -1 n 1 n + 1 Ω f α 0 α 0 β 0 1 n u 1 1 2 ( 1 n +1) M 0 dΩ (88) + C 1 m + 1 Ωn α 0 |u 1 | 1 2 ( 1 m +1) M 0 dΩ (89) + A -1 n 1 n + 1 l=2,...,L Ω α l-1 α l-1 β l-1 1 n (90) × u l -u l-1 1 n +1 M l-1 dΩ + ρg Ω l=1,...,L h l ∇s • u l dΩ (91) - Ωc l=1,...,L F l n • u l dS. ( 92
)
More precisely, one can show that solving (86) is equivalent to solving the minimisation problem:

Find u s.t. J (u) ≤ J (v), ∀v. (93)
Like the SSA, minimization problem (93) consists of a vector p-Laplace problem with p = 1 + 1 n < 2. However, unlike the SSA, this p-Laplace problem is additionally vectorial over the interlayer terms (90), which couples the layers. Terms (90) can be seen as a penalisation of jumps for the piecewise-constant velocity profile, as in formulations of the Discontinuous Galerkin method [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications[END_REF]. While the p-Laplace term (87) corresponds to the longitudinal stresses, the penalisation (90) corresponds to the vertical shear stress components. In addition, the terms ( 88) and ( 89) apply only on the lowest layer and correspond to the no-slip and the sliding conditions. Finally, (91) represents the gravitational forces while (92) corresponds to the balance between the ice pressure and the hydrostatic sea water pressure.

One can show that the functional J is convex, strongly continuous in W 1,1+ 1 n (Ω) and therefore weakly lower semi-continuous [START_REF] Schoof | A variational approach to ice stream flow[END_REF][START_REF] Schoof | Coulomb friction and other sliding laws in a higher order glacier flow model[END_REF]. Additionally, J is coercive if the h l are uniformly lower bounded by a positive constant [START_REF] Schoof | A variational approach to ice stream flow[END_REF][START_REF] Schoof | Coulomb friction and other sliding laws in a higher order glacier flow model[END_REF]. It follows the existence of a minimizer of J , provided this last assumption is satisfied. In addition, J is strictly convex and the minimizer is unique if Ω n ∪ Ω f has a positive measure.

For convenience, the functional J corresponding to the 2-layers * model for a simplified 2D flow in the vertical x -z plane is rewritten in Appendix A.

Numerical method

In this section, two methods to solve the multilayer system (82) (83) (84) with the boundary condition (85) are described. Like for the SSA [START_REF] Jouvet | An adaptive newton multigrid method for a model of marine ice sheets[END_REF], the variational problem (86), or equivalently the minimisation problem (93), are more appropriate than the strong form (82) (83) (84) to implement the finite element method. Consider a regular triangulation of Ω ⊂ R 2 and the finite element space, which is spanned by the continuous, linear functions on each triangle. Since this finite element space can be identified to R I×2 , where I is the number of nodes of the mesh, the Ritz-Galerkin approximation of the minimisation problem (93) reads:

Find U ∈ R L×I×2 s.t. J h (U ) ≤ J h (V ), ∀V ∈ R L×I×2 ,
(94) where the vector U contains the nodal values of the approximation of u and J h : R L×I×2 -→ R ∪ {+∞} is strictly convex, coercive and lower semi-continuous. Here J h is either J h = J defined by ( 87) or an approximation of J by numerical quadratures. The simplest strategy for solving (94) is to use a Gauß-Seidel method. This method consists of minimizing successively J h in each coordinate directions (l, i, k) ∈ L×I ×2. Since these minimisations are scalar, there exist many methods (bisection, Newton, ...) to perform these tasks. There are two main ways to loop over the indices (l, i, k) ∈ L×I ×2: looping first over the layer indices l and then over the horizontal node indices i or the opposite. Opting for one or for the other strategy leads to two different methods, which are successively described in the two next sections.

Method 1

In the first method, the approximation sequence

U ν = (U 1 ν , ..., U L ν ), (where U l ν ∈ R I×2
), which is initialized by U 0 , is defined recursively by taking the solutions of the successive minimisation problems:

J U 1 ν+1 , U 2 ν , ..., U L ν ≤ J V, U 2 ν , ..., U L ν , ∀V, ... J ., U l-1 ν+1 , U l ν+1 , U l+1 ν , . ≤ J ., U l-1 ν+1 , V, U l+1 ν , . , ∀V, ... J U 1 ν+1 , U 2 ν+1 , ..., U L ν+1 ≤ J U 1 ν+1 , U 2 ν+1 , ..., V , ∀V,
where the V are taken in R I×2 . Thus, this first method consists of solving the multilayer system ( 82), ( 83) or (84) layer by layer, e.g. from the lowest one l = 1 to the highest one l = L, using the old solution U ν in the interface terms, and to iterate. It remains to define a method to solve each individual minimisation problem. Since each of these problems are similar to the one resulting from the SSA, one can use the solvers that was developed for the SSA. The non-linear Gauß-Seidel method is well-known to suffer from poor convergence rates. To improve the convergence, one can use the method described in [START_REF] Jouvet | An adaptive newton multigrid method for a model of marine ice sheets[END_REF], which combines a Newton-type acceleration with a linear geometric multigrid method for solving the correction step. Additionally, it might take a lot of iterations to reach the final solution when initializing U 0 by zero. To get a better initialisation, one can apply a "coarse-to-fine" strategy, which consists of first using the 1-layer model, copying the solution on two layers, using the 2-layer model, ect, until to reach the prescribed number of layers. As said above, one can naturally upgrade any SSA solver to a multilayer solver with this first method. However, computing the ice flows of the "infinite parallel-sided slab" for which the exact solution is known, see Appendix B, shows that the convergence severely deteriorates when the number of layers grows. Indeed, the number of iterations needed to reach a given accuracy increases exponentially with respect to the number of layers. This can be justified as follows: further efforts are needed to transfer some informations between the lowest and the highest layers if they are many layers. This phenomena is even more pronounced if the vertical coupling dominates. In contrast, if the vertical coupling is slight, then the convergence can be reached after a single loop. As a consequence, this first method is wellsuited for computing the ice flows of a system of ice sheet / ice shelf since the longitudinal stresses dominate the vertical ones. However, this method is not efficient for computing the ice flows of smaller mountain glaciers, in which the shear stresses significantly increase the vertical coupling between layers. For these types of ice flows, it is recommended to use the second method, which is described in the next section.

Method 2

In contrast with the first method, the second one first loop over the node indices i, i.e. the sequence U ν = (U 1 ν , ..., U I ν ), (where U i ν ∈ R L×2 ), which is initialized by U 0 , is defined recursively by taking the solutions of the successive minimisation problems:

J U 1 ν+1 , U 2 ν , ..., U I ν ≤ J V 1 , U 2 ν , ..., U I ν , ∀V, ... J ., U i-1 ν+1 , U i ν+1 , U i+1 ν , . ≤ J ., U i-1 ν+1 , V i , U i+1 ν , . , ∀V, ... J U 1 ν+1 , U 2 ν+1 , ..., U I ν+1 ≤ J U 1 ν+1 , U 2 ν+1 , ..., V I , ∀V.
where the V are taken in R L×2 . For simplicity, the Gauß-Seidel method is written above, but this one can be better replaced by the Newton multigrid method described in [START_REF] Jouvet | An adaptive newton multigrid method for a model of marine ice sheets[END_REF]. In any case, it remains to define a method to solve each columnwise minimisation problem (of size R L×2 ). For that, a non-linear Gauß-Seidel method, which minimizes successively from the lowest to the highest layer, can be used. However, since one never needs more than around ten layers (L ∼ 10) in practise, those minimisation problems are of small size and then can be more efficiently solved by Newton's method. The convergence of this second method is expected to be much better compared to the first method when solving strongly vertically coupled ice flows.

Validation

The numerical methods presented in this section were validated against exact solutions. On the one hand, the 1-layer model was tested with manufactured solutions similar to the ones used in [START_REF] Schoof | A variational approach to ice stream flow[END_REF]. However, these solutions could not be used to validate the multilayer model, which involves additional interlayer terms. For this reason, these terms were tested separately with the exact solution of the "infinite parallel-sided slab" built in Appendix B.

Results

In this section, the ice flows of two different applications are computed with several multilayer and higher-order models to perform a comparison study. In contrast with Section 2, the solutions of Section 4 are restricted to a simple setting: the flow is 2D in the vertical x -z plane (no y dependency), no shelf part and no sliding condition on the bedrock are considered. The first application concerns the experiment B of the ISMIP-HOM project (Pattyn and others, 2008) while the second one concerns a vertical section of Gries glacier, Switzerland. In both applications, the multilayer solutions are computed using the method 2 of Section 3 with a uniform division of the ice thickness, i.e. h l = h/L = H. On the one hand, the horizontal segment Ω, which supports the glacier, is uniformly divided into 64 segments to generate a 1D mesh. This 1D mesh is used to compute the velocity field with the 1-layer, 2-layer, 4-layer, 8-layer, 16-layer * and the 16-layer models (the L-layer * corresponds to the L-layer model, but with the simplified redefinition of the interlayer tangential stress (71) instead of ( 68), see Section 2.4). On the other hand, a triangular 2D mesh is built by extruding 16 vertical layers of the 1D mesh between the lower and the upper surfaces. This 2D mesh is used to compute the velocity fields of the FOA and the Stokes models. The numerical convergence of the 16layer * , 16-layer, the FOA and the Stokes solutions was assessed by checking the discrepancy between the solution and the one obtained by doubling the horizontal resolution or the number of layers. The following physical parameters were used: ρ = 910 kg m -3 , ρ w = 1000 kg m -3 , n = 3, and g = 9.81 m s -2 .

ISMIP-HOM B

In this section, the results of experiment B of the ISMIP-HOM project (Pattyn and others, 2008) are reported. In this experiment, the geometry was defined by

s(x) = -x tan(0.5 • ), s(x) = s(x) -1000 + 500 sin (2πx/L x ) .
for x ∈ Ω = [0, L x ], the no-slip condition (12) was prescribed on the bedrock, periodic boundary conditions were used at the boundaries of Ω and the rate factor was taken equal to A = 3.17 × 10 -24 Pa -3 s -1 , see (Pattyn and others, 2008) for further details. Figure 5 displays the solutions of all models with L x = 10, 40 and 160 km. Figure 5 (top) shows that the multilayer models converge when increasing the number of layers, such that the 8 and the 16-layer solutions are nearly indistinguishable. Additionally, this convergence does not depend on the wavelengths L x in the boundary conditions. As a matter of fact, the multilayer model solutions get closer to the FOA and the Stokes solutions when increasing the number of layers, however, do not converge to any of those solutions. For instance, the 16-layer solution is above the FOA solution by about 10% when L x = 10 km, however, this disagreement decreases for the higher wavelengths (L x = 40 km and L x = 160 km). Interestingly, when L x = 10 km, the 16-layer matches well the Stokes solution in the first half of the domain, but moves slightly away in the second part, where a depression occurs in the bedrock. However, the FOA and Stokes solutions are identical for the higher wavelengths while the 16-layer one shows a discrepancy in the depression area. Here, the deficiency of the multilayer model is attributed to the choice of the vertical division (40). Indeed, the uniform splitting leads to nearly flat layers close to the top. However, the streamlines of the Stokes velocity field (not shown) are significantly more curved at the top of the depression. This violates the assumption of alignment between the ice flows and the layers, on which the multilayer model is based. Finally, the simplified 16-layer * model leads to higher disagreements with the reference solutions (FOA and Stokes) than the 16-layer model at small wavelengths (i.e. when the surface gradient are not negligible) while no differences between the 16-layer and the 16-layer * solutions are observed for larger wavelengths. This is due to the fact that, in contrast with the 16-layer model, the 16layer * neglects O(δ 2 ) terms in the redefinition of the interlayer stress (71), where δ is the slope of the layers. This proves that these terms are relevant and must be kept where the glacier is steep.

Gries glacier

Gries glacier is a 5-km-long glacier situated in Switzerland with an ice thickness ranging from 0 to about 200 meters. The geometry of the glacier along a central flow-line is available in 1961 from [START_REF] Kirner | Modélisation mathématique et simulation numérique des phénomènes dynamiques et thermiques apparaissant dans un glacier[END_REF] such that the section of the glacier can be represented in the x -z plan. The velocity field of Gries glacier was computed with the no-slip condition (12) on the bedrock and A = 2.5 × 10 -24 Pa -3 s -1 as a rate factor. For comparison purposes, Figure 6 displays the horizontal velocity fields obtained with the 16-layer * , the 16-layer, the FOA, and the Stokes models, respectively.

As a first result, the 16-layer * and the 16-layer models show significantly different results in the steepest part, while the discrepancy is invisible elsewhere. Taking the Stokes solution as a reference (since it is free of any simplification), it appears that the 16-layer * model overes- timates the solution in this area. As in Section 4.1, this confirms the relevance of the O(δ 2 ) terms in the redefinition of the interlayer stress (71). On the other hand, the FOA and the Stokes solutions also show some discrepancies in the steepest part, which are due to the simplifications O(ǫ 2 ) of Section 2.2. Interestingly, the 16-layer solution matches the FOA solution very well. Like in Section 4.1, the discrepancy between the multilayer and the Stokes solutions is attributed to a mismatch between the alignment of the layers and the real direction taken by the ice flows. Fcombi

Discussion

The numerical results of Section 4 have shown that the multilayer model, which is 2D-vectorial, can be used as an alternative of the 3D higher-order models. Three aspects of this work will be further investigated in the future: the mathematical model, the numerical model and the applications, each one being described in the three next paragraphs.

In most of existing ice flow models, the mechanical modelling strictly precedes the numerical modelling. It is sometimes necessary to go back to the physics after the discretization, e.g. when preconditioning the linear systems [START_REF] Brown | Composable linear solvers for multiphysics[END_REF]. As in [START_REF] Audusse | A multilayer saint-venant system with mass exchanges for shallow water flows. derivation and numerical validation[END_REF], the common order of modellings was here partly inverted since the vertical discretization comes first followed by the mechanical modelling and the horizontal numerical discretization. A question naturally arises from this approach: Can the multilayer model be derived from a vertical semi-discretization of the Stokes equations by a Discontinuous Galerkin method [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications[END_REF] or a finite volume method with an extruded mesh? Indeed, it is intriguing to notice that the interlayer terms (Eq. ( 90)) look similar to penalisation terms for the jumps of the velocity between layers. This statement calls for several model extensions. One of them consists of increasing the order of the vertical discretization, by allowing the velocity field to be described by vertically piecewise polynomials. However, it is not clear if this approach keeps intact the mathematical structure of the multilayer model, which takes strongly advantage of the reformulation as a minimisation problem. All other existing hybrid models have this structure while the velocity profile is given by an implicit relation involving an integral for "L1L2-like" model [START_REF] Hindmarsh | A numerical comparison of approximations to the stokes equations used in ice sheet and glacier modeling[END_REF][START_REF] Schoof | Thin-film flows with wall slip: An asymptotic analysis of higher order glacier flow models[END_REF][START_REF] Goldberg | A variationally derived, depthintegrated approximation to a higher-order glaciological flow model[END_REF] or by a polynomial for the isothermal SIA+SSA [START_REF] Winkelmann | The potsdam parallel ice sheet model (PISM-PIK) part 1: Model description[END_REF]. As a consequence, a natural way to extend the multilayer model is to include any of the "L1L2-like" models individually on each layer. In this paper the ice flows were assumed to be isothermal for simplicity. However, coupling the multilayer model with a thermal model should be considered in practical applications.

In Section 3, a first numerical method was presented to solve the multilayer system using a solver for the SSA and an iterative method that loops over layers. This method is suited to compute the ice flows of a shallow system of ice sheet / ice shelf, in which the horizontal stresses dominate. However, its convergence rate decreases dramatically when computing strongly vertically coupled ice flows, like these of mountain glaciers. To overcome this problem, a second method, which reverses the order of loops (horizontal over nodes or vertical over layers), is expected to converge must faster for such an application. To verify this, a comparison study of the numerical performance will be addressed in a future work. The second method presented in this paper is still subject to improvements to make its convergence mesh-independent and inherent to the type of flows. For instance, one could apply the non-linear multigrid solver [START_REF] Jouvet | An adaptive newton multigrid method for a model of marine ice sheets[END_REF] directly in the 3D structure of the minimisation problem (93). Another crucial as-pect of the multilayer model concerns the choice of the vertical subdivision and its impact on the results. This splitting must be done in such way that the ice flows and the layers are approximatively aligned. The choice of a uniform splitting likely contributed to the discrepancies observed between the multilayer and the Stokes solutions of Section 4. The optimality of the multilayer splitting will investigated in a future work.

In this paper, the multilayer model was derived for a general system of ice flows, however, the solutions computed in Section 4 were restricted to a simple setting: 2D geometry, no shelf part and no sliding condition on the bedrock. To complete this work, further comparisons with other higher-order models will be performed. In particular, the multilayer model will be used to simulate the idealized marine ice sheets of the project MISMIP [START_REF] Schoof | MISMIP: Marine ice sheet model intercomparison project[END_REF] and MISMIP 3D (Pattyn and others, 2013). A special attention will be given on the ability of the model to reproduce the motion of the grounding line indicated by other higher-order models. In such an application, the simultaneous use of several multilayer models will be useful: the multilayer model can be used in the vicinity of the grounding line while the 1-layer model (the SSA) model can be used elsewhere. In addition, the multilayer model will be also tested to simulate 3D mountain glaciers like in [START_REF] Jouvet | Numerical simulation of Rhonegletscher from 1874 to 2100[END_REF][START_REF] Jouvet | Modelling the retreat of grosser Aletschgletscher in a changing climate[END_REF] in diagnostic and prognostic ways, such that the computational and the mechanical performances of the multilayer model can be compared with those of the Stokes model.

A 2-layer * flow-line model for a marine ice sheet

In this appendix, the multilayer model of Section 2.6 is rewritten in a simplified setting, which depicts a 2D marine ice sheet like the one of Figure 2. More precisely, the following assumptions are made:

i) The flow is 2D in the vertical x -z plane (no y dependency), Ω = [0, x c ] is an interval, and x c is the abscissa of the calving front.

ii) The bedrock b decreases with respect to x such that there are a grounded ice sheet on the left-hand-side Ω = [0, x g ] and a floating ice shelf on the righthand-side Ω = [x g , x c ] of the domain Ω.

iii) The non-linear sliding law (13), ( 14) applies on the grounded area [0, x g ] while the floating condition (15) applies elsewhere.

iv) The multilayer splitting is uniform and made of two layers, i.e. if h l = h/2 for l = 1, 2.

v) In the 2-layer * model, the interlayer tangential stress is redefined by (71) instead of (68), see Section 2.4.

In this simplified setting, the system (82), ( 83), (84) becomes: Find (u 1 , u 2 ) (the index x is removed for clarity) such that

-2A -1 n ∂ x h 2 ∂ x u 2 1 n -1 ∂ x u 2 (95) + A -1 n u 2 -u 1 h 1 n -1 u 2 -u 1 h = -ρg h 2 ∂ x s,
and

-2A -1 n ∂ x h 2 ∂ x u 1 1 n -1 ∂ x u 1 (96) + C|u 1 | 1 m -1 u 1 1 [0,xg] + A -1 n u 1 -u 2 h 1 n -1 u 1 -u 2 h = -ρg h 2 ∂ x s,
while at the calving front, the condition (85) becomes:

2A -1 n h 2 |∂ x u l | 1 n -1 ∂ x u l := F l . (97)
Finally, the variational formulation (86) becomes B Multilayer exact solution for the "infinite parallel-sided slab"

2A -1 n l=1,2 xc 0 h 2 |∂ x u l | 1 n -1 (∂ x u l )(∂ x v l )dx (98) + C xg 0 |u 1 | 1 m -1 (u 1 )(v 1 )dx + A -1 n xc 0 u 2 -u 1 h 1 n -1 u 2 -u 1 h (v 2 -v 1 )dx = -ρg
The "infinite parallel-sided slab" simplified setting [START_REF] Greve | Dynamics of Ice Sheets and Glaciers[END_REF] relies on the following assumptions:

i) The flow is 2D in the vertical x -z plane (no y dependency) and Ω = [0, x l ] is an interval.

ii) The bedrock slope ∂ x b and the ice thickness h = H are constants, such that ∂ x b = ∂ x s.

iii) The no-slip condition (12) applies everywhere on the bedrock.

iv) The bedrock is everywhere above sea level such that there is no floating ice shelf.

v) Periodic boundary conditions are prescribed at x = 0 and x = x l in order to simulate an infinite horizontal domain.

Since the geometry shows no x variation, the solution (u 1 , ..., u L ) is independent of x such that all derivatives can be removed in the multilayer system (82) (83) (84), which rewrites: Find (u 1 , ..., u L ) (the index x is removed for clarity) s.t.

α 2 α 2 u L -u L-1 2AH 1 n = -ρgH(∂ x b) ... ... α 2 α 2 u l -u l-1 2AH 1 n -α 2 α 2 u l+1 -u l 2AH 1 n = -ρgH(∂ x b) ... ... α 2 α 2 u 1 AH 1 n -α 2 α 2 u 2 -u 1 2AH 1 n = -ρgH(∂ x b),
where α = 1 + (∂ x b) 2 and α = 1 for the multilayer and the simplified multilayer * models, respectively. It is easy to verify that the solution of the system above is

u l = -2 A α 2(n+1) (ρg(∂ x b)) n H (100)
× 1 2 (HL) n + (H(L -1)) n + ... + (H(L -l + 1)) n .

Therefore, this exact solution can be used to check the convergence of the methods of Section 3. Note that (100) consists of the integration of

∂ z u = -2 A α 2(n+1) [ρg(∂ x b)(s -z)] n , (101)
with the rectangle formula as a quadrature rule on any layer. Interestingly, the multilayer solution (101) equals the exact solution of the Stokes system [START_REF] Greve | Dynamics of Ice Sheets and Glaciers[END_REF], while the multilayer * (when α = 1) solution (101) equals the SIA solution [START_REF] Greve | Dynamics of Ice Sheets and Glaciers[END_REF])

∂ z u = -2A[ρg(∂ x b)(s -z)] n . (102)

Figure 1 :

 1 Figure 1: Overview of the hierarchy of ice flows models: Stokes, First Order Approximation (FOA), Shallow Ice Approximation (SIA), Multilayer model, "L1L2-like" and Shallow Shelf Approximation (SSA). The dimension of the mathematical model is indicated in exponent.

FFigure 2 :

 2 Figure 2: Section of an ice sheet and an ice shelf, with notation.

Figure 3 :

 3 Figure 3: Multilayer splitting of the ice thickness.

Figure 4 :

 4 Figure 4: Representation of the global frame (left) and the local frame (right).

  Figure 5: ISMIP-HOM Experiment B results for L x = 10, 40 and 160 km, from left to right, respectively. Top: surface horizontal velocities of the 1-layer, 2-layer, 4-layer, 8-layer, and the 16-layer models. Middle: surface horizontal velocities of the 16-layer * , the 16-layer, the FOA and the Stokes models. Bottom: horizontal velocity field of the 16-layer model. For convenience, the horizontal domain is stretched.

Figure 6 :

 6 Figure 6: The horizontal velocity field across a section of Gries glacier in 1961 obtained with the 16-layer * , the 16-layer, the FOA, and the Stokes model, respectively. The unit is meters per years and the continuous lines correspond to level lines every 10 meters per year. For convenience, the horizontal domain is stretched by a factor 5.

F

  s)v l dx + l=1,2 F l (x c )v l (x c ),and the functional J rewritesJ (u 1 , u 2 ) = A -l (x c )u l (x c ).

In the literature(Hindmarsh,

[START_REF] Hindmarsh | A numerical comparison of approximations to the stokes equations used in ice sheet and glacier modeling[END_REF][START_REF] Egholm | Modeling the flow of glaciers in steep terrains: The integrated second-order shallow ice approximation (isosia)[END_REF], the terminology "multilayer" is sometimes used as a synonym of "hybrid", as defined in the paper

Acknowledgements

The author wishes to acknowledge Ralf Kornhuber for support and Ed Bueler for his detailed and helpful comments to the manuscript.

* Supported by the Deutsche Forschungsgemeinschaft (project KL 1806 5-1)