BLAISE

The motion of a three-dimensional glacier is considered. Ice is modeled as an incompressible non-Newtonian ﬂuid. At each time step, given the shape of the glacier, a nonlinear elliptic system has to be solved in order to obtain the two components of the horizontal velocity ﬁeld. Then, the shape of the glacier is updated by solving a transport equation. Finite element techniques are used to compute the velocity ﬁeld and to solve the transport equation. Numerical results are compared to experiments on Storglaciaren (Sweden) between 1959 and 1990.


Introduction
Over the last 150 years, the Swiss Alps lost 50 % of their glaciers.Today, most of alpine glaciers are retreating, leaving visible marks on the landscape.The future evolution of glaciers is not only of concern for the tourist industry, but also for agriculture and hydro-power production.The past evolution of glaciers is an archive of climatic changes [2].
Numerical glacier modelling has become an important tool for glaciologists in order to understand the interactions between glaciers and climate.Most of the numerical methods use finite differences with fixed cartesian grids [18].
Recently, finite element techniques have been introduced in this area, see for instance [20,15], and several mathematical papers have been published in two space dimensions.The well posedness of a shallow ice model for the computation of the velocity is discussed in [10].A priori and a posteriori error estimates are derived for the corresponding finite element method in [14,8].The question of modelling errors for this model is discussed in [6].A finite element method to compute the ice velocity and the motion of the glacier is presented in [21].
The goal of this paper is to present a three-dimensional finite element model in order to compute the ice velocity and the motion of the glacier, thus extending the two-dimensional model presented in [21] to three space dimensions [23].The model is presented in the next section.Given precipitations estimates, the size and motion of a glacier is computed.Energy balance is disregarded.The velocity field is modelled as an incompressible Non-Newtonian fluid.Assuming that the horizontal extent of the ice sheet is much larger than the ice thickness, a first order analysis of the momentum equations yields to a strongly non-linear elliptic system for the two horizontal components of the velocity field.Then, mass conservation yields to a transport equation for the glacier height, with a given right-hand side corresponding to the climatic input due to snow falls or ice melting.At each time step, given the shape of the glacier, the strongly non-linear elliptic system studied in [22] has to be solved in order to obtain the two components of the horizontal velocity field.Then, the shape of the glacier is updated by solving the transport equation for the glacier height.Finite element techniques are used to compute the velocity field and to solve the transport equation.The transport equation is discretized using an Arbitrary Lagrangian Eulerian (ALE) formulation.In section 4, numerical results are compared to experiments on Storglaciaren (Sweden) between 1959 and 1990.
When comparing the present paper to the two-dimensional model presented in [21], the interested reader will note that the overall methodology is the same.However, the extension from 2D to 3D computations is not obvious, specially from the point of view of mesh generation and deformation.This issue is explained in details in Section 3.

The model
Our goal is to present a model predicting the evolution of a glacier during decades, possibly centuries given the amount of ice accumulating or disappearing on the top surface of the glacier.At time t, the unknown ice domain Ω(t) is parametrized as follows : where Λ(t) is the unknown glacier support, B is the prescribed mountain bedrock, S = B + H the unknown top surface of the glacier, see Fig. 1 for the notations.Given the snow falls or ice melting, our goal is to compute when time t is going on, the horizontal components u and v of the ice velocity and the shape of the glacier, that is the glacier support Λ and the glacier height H.
Coupling the size and motion of a glacier with the energy balance would be a significant accomplishment in the field but is beyond the scope of the present paper.We refer to [13] for a full description of models involved in glacier dynamics. z Figure 1.The notations.The ice domain at time t is Ω(t), its projection onto the xy plane is the glacier support Λ(t).The mountain bedrock is parametrized by B, the top surface by S, the glacier height is H.The function b is a given quantity corresponding to mass accumulation or ablation (snow falls or ice melting).The velocity field of ice has components u, v and w.The normal at the top surface has components n x , n y and n z .

The velocity field
For the time scales involved (from several years to a century), glacier ice is behaving like a fluid, not like a solid.Ice can be modeled as an incompressible viscous fluid and time-dependent effects can be neglected.Thus, given the ice domain Ω(t), a stationary problem has to be solved in order to obtain the ice velocity.Several models are available and the case of a shallow glacier is considered here.More precisely, it is assumed that the horizontal extent of the ice sheet is much larger than the ice thickness.In [3], a scale analysis is performed for such a glacier, the scaling parameter being the ratio between the ice thickness and the length (or width) of the glacier.For instance in small valley glaciers this ratio ranges between 0.01 and 0.1 (eventually more in steep regions).Below, we briefly present the corresponding first order model in three space dimensions, see [3] for details.A mathematical and numerical analysis of this model is provided in [22].
It should be noted that this first order model requires only the horizontal components u and v of the velocity field to be computed.Indeed, the horizontal component w is not contained in the set of equations anymore, although it can be recovered from the mass conservation equation.The fact that w is absent from the set of equations is due to the fact that w occurs only in terms which are of second order magnitude with respect to the small parameter (the ratio between the ice thickness and the length or width of the glacier).Thus, the feedback of the horizontal component w to the horizontal components u and v is small, or negligible in first order accuracy.
At time t, the equations governing the horizontal components u(t), v(t) : Ω(t) → R of the velocity field are as follows: Here µ is the ice viscosity, which is implicitly defined by the relation where A, T 0 and n are positive coefficients and s is given by On the mountain bedrock z = B(x, y), (x, y) ∈ Λ(t), zero Dirichlet conditions are imposed: It should be noted that the above no-slip boundary condition is too simple to reproduce complex phenomena such as glacier surges [12,13].Indeed, in temperate glaciers, due to water flowing between the bedrock and the ice sheet, ice may slip along the bedrock.In the 2D simulation of Gries glacier [21], the velocity on the bedrock was tuned so that it was fitting the measured velocity on the top surface.Here, the measured velocity on the top surface is available only in a small portion of the glacier (Storglaciaren), thus a different strategy is adopted.No-slip boundary conditions are imposed on the mountain bedrock and the parameter A involved in (2.3) is tuned so that the computed velocity best fits the measured one.In this subsection, the two components u and v of the horizontal velocity are assumed to be known in the ice domain Ω(t) and an additional equation governing the evolution of the glacier height H is derived.Consider, as in Fig. 2 Let ū, v be the vertical average of the horizontal velocities u, v defined by

A model for the glacier height
(2.10) Then, since the domain D is an arbitrary subset of the glacier support Λ(t), (2.8) leads to The evolution of the ice domain Ω(t) is known whenever the glacier support Λ(t) and the glacier thickness H(t) are known.The glacier thickness is governed by the transport equation (2.11).In order to determine Λ(t) we need to provide the additional equation (2.12)

Model recapitulation
Our mathematical model is then the following.Given the glacier bedrock B, the initial glacier support Λ(0), the initial glacier height H(0), given positive parameters A, T 0 , n, given the function b, find the horizontal velocities u, v satisfying (2.1) (2.2), the glacier height H(t) satisfying (2.11) and the glacier support Λ(t) satisfying (2.12).

Algorithm
The numerical procedure used to solve the above mathematical problem is based on our physical understanding of the problem and is the following.
At each time step, given an approximation of the ice domain Ω(t) we proceed as follows.
• Compute the new glacier height H by solving (2.11) with an Arbitrary Lagrangian Eulerian (ALE) finite element method.
The three steps of this decoupling algorithm are presented in details hereafter.At initial time Λ(0) and are known so that the initial ice domain Ω(0) = {(x, y, z) ∈ R 3 ; (x, y) ∈ Λ(0); B(x, y) ≤ z ≤ B(x, y) + H(x, y, 0)} is known.Let ∆t be the time step, t m = m∆t, m = 0, 1, 2, . . . .Consider the time iteration number m and let Λ m , H m : Λ m → R and Ω m be approximations of Λ(t m ), H(t m ) : Λ(t m ) → R and Ω(t m ), respectively.The goal is : • to compute a new approximation • to compute an approximation Ω m+1 of the ice domain Ω(t m+1 ).In this subsection, the domain Λ m and the function H m : Λ m → R are known and we want to compute an approximation u m , v m of the horizontal velocity components u(t m ), v(t m ) : Ω(t m ) → R. Let T h (Λ m ) be a regular triangulation of Λ m into triangles with size less than h.We assume that H m : Λ m → R is a continuous piecewise linear function on each triangle of T h (Λ m ).Let further B m be a piecewise linear interpolation of B on this same triangulation and let S m = B m + H m .The ice domain Ω m is then defined by

Computation of the velocity field
and is meshed into tetrahedrons.For stability purposes, all the triangulations T h (Λ m ) and all the tetrahedral meshes T h (Ω m ), m = 0, 1, 2, ... have the same topology.Moreover, the vertices of T h (Ω m ) are aligned along vertical segments with endpoints (x m j , y m j , B m j ) and (x m j , y m j , B m j +H m j ), where (x m j , y m j ) are the vertices of triangulation T h (Λ m ) and B m j = B m (x m j , y m j ), H m j = H m (x m j , y m j ), see Fig. 3 and 4.
for all test functions (ϕ, ψ), continuous, piecewise linear on the tetrahedrons of T h (Ω m ), vanishing on the bedrock of the glacier.Here µ m is the ice viscosity which is implicitly defined from the derivatives of u m , v m as in (2.3) (2.4).It should be noted that existence of a solution to (3.2) and convergence when the mesh size h goes to zero has been addressed in [22].
The difficulty is due to the fact that (3.2) is a nonlinear problem since µ m depends implicitly on the derivatives of the unknown velocities u m , v m .Picard's iterative method is used to solve the nonlinear system, no under-relaxation being necessary to reach convergence, see the proof in [23].
Since the vertices of T h (Ω m ) are vertically aligned, their projection onto the xy plane coincide with the vertices of T h (Λ m ), which facilitates the computation of the vertically averaged horizontal velocity components ū, v defined by (2.9) (2.10).Using a trapezoidal quadrature formula we can compute approximations ūm j and vm j of ū(x m j , y m j , t m ) and v(x m j , y m j , t m ) as follows: Here N j is the number of vertices in the mesh T h (Ω m ) aligned on top of a vertex (x m j , y m j ) in the triangulation T h (Λ m ) (see Fig. 3) and ∆z m j = H m (x m j , y m j )/N j is the vertical mesh spacing.The functions ūm , vm : Λ m → R are then defined as where ϕ m j (x, y) are the continuous, piecewise linear shape functions associated with the triangulation T h (Λ m ).

Computation of the new glacier support Λ m+1
In order to determine the new glacier support Λ m+1 we proceed in two steps: first, we move the vertices lying on the boundary of Λ m , then we move the internal vertices of Λ m .

Moving the vertices lying on the boundary of Λ m .
At time t, let (x(t), y(t)) be a point lying on the boundary of glacier support Λ(t).We compute an approximation of x(t + ∆t), y(t + ∆t) using the first order expansion In order to compute the ẋ(t), ẏ(t), we take the derivative of (2.12) with respect to time and obtain x + d 2 y = 1, be the direction in which the vertices located on the boundary of Λ(t) will be moved.We have: ẋ(t) = δd x and ẏ(t) = δd y , where δ is the velocity at (x(t), y(t)) in direction d.From (3.5) we obtain Now, let (x m j , y m j ) be a vertex of the triangulation T h (Λ m ) lying on the boundary of Λ m .We compute the position of this vertex at time t m+1 as follows: where δ m j is an approximation of δ(x m j , y m j , t m ) defined by (3.6).Since ū, v and H are continuous, piecewise linear on each triangle K of T h (Λ m ), their space derivatives are piecewise constant, so that a reconstructed gradient has to be computed at the vertices (x m j , y m j ) in order to compute δ m j .This reconstructed gradient is obtained by averaging the gradient over all triangles containing (x m j , y m j ).It now remains to choose the direction d m j = ((d x ) m j , (d y ) m j ) T along which the boundary vertex (x m j , y m j ) is moved.The most obvious choice is to move the points towards the direction normal to the boundary, that is d m j proportional to −∇H m (x m j , y m j ).In practice, this choice is not the best one.Indeed, when the boundary of Λ m has a strong curvature, then two neighboring vertices may intersect, which leads to a non conforming mesh, see Fig. 5.In order to remedy this problem, we modify the direction in which we move the boundary vertices.If the curvature at a given vertex (x m j , y m j ) is larger than the curvature at its neighbors, then the direction d m j is expected to be aligned with the normal at (x m j , y m j ).On the other side, if the curvature at a given vertex (x m j , y m j ) is smaller than the curvature at its neighbors, then the direction d m j is expected to be aligned with the normal at its neighbors.The following algorithm is used to obtain the direction d m j at each vertex (x m j , y m j ) lying on the boundary of Λ m .We start with a direction normal to the boundary.
where ∇H m j is the reconstructed gradient of H m at vertex (x m j , y m j ).Then, the new direction is computed as a linear combination of its own direction vector d m j and the directions of the two neighbors, d m j−1 , d m j+1 , weighted by a coefficient c m j depending on the curvature.This first step propagates the direction of a node with large curvature only to its two immediate neighbors, which might not be sufficient to avoid crossover.We apply the above step several times to achieve the desired result : and then normalize the vector to one.The coefficient c m j depends on the angle θ m j reported in Fig. 6 and we set c j = π/θ m j .
Figure 6.Definition of the angle θ m j .

Moving the internal vertices of Λ m .
Once we have moved the boundary nodes of T h (Λ m ), we move the interior nodes without changing mesh topology and obtain the new mesh T h (Λ m+1 ).To this end, the Laplacian method (see for instance [19,24]) is used.The method consists in solving the following system : where X m+1 and Y m+1 are the new mesh coordinates and ∆ m is the Laplacian operator defined on the mesh prior to deformation.We use Dirichlet boundary conditions in order to impose the coordinates of the nodes on the new border: for all vertices (x m j , y m j ) lying on the boundary of T h (Λ m ).An example of mesh deformation is reported in Fig. 7 for Storglaciären, Sweden, the mesh deformation being exaggerated.It should be noted that the Laplacian method for moving the grid points exhibits poor behavior when non-convex geometries are encountered.More elaborated strategies are available for avoiding this problem, see [7].The use of such techniques has not been necessary in the framework of our model.

Computation of the new glacier height H m+1
Our goal is now to compute an approximation H m+1 : Λ m+1 → R of H(t m+1 ), the solution of (2.11), given an approximation H m : Λ m → R. Since both H m and H m+1 are computed on different domains, an Arbitrary Lagrangian Eulerian (ALE) formulation (see for instance [4]) will be used.
Our finite element discretization of (3.11) then consists in finding Ĥm+1 : Λ = Λ m+1 → R, continuous, piecewise linear, that is such that Ĥm+1 = 0 on the boundary of Λ = Λ m+1 and such that ) for all continuous, piecewise linear functions φ vanishing on the boundary of Λ = Λ m+1 .The additional diffusion term has been added here for stability purposes, the classical SUPG or GLS method being not diffusive enough in our context.Following [11], the parameter is equal to the local mesh size times the mean local velocity.Also, the vertically averaged velocities ûm , vm : Λ = Λ m+1 → R are defined by where ūm j and vm j are defined in (3.3).

Computation of the new ice domain Ω m+1
Given the new glacier support Λ m+1 , given the new glacier height H m+1 : Λ m+1 → R, the new ice domain Ω m+1 is defined by where B m+1 is a piecewise linear interpolation of the glacier bedrock B on the triangulation T h (Λ m+1 ).For stability purposes, the tetrahedral mesh T h (Ω m+1 ) has the same topology than the previous one, T h (Ω m ).Moreover, the vertices of the tetrahedral mesh T h (Ω m+1 ) are aligned on top of those of the triangulation T h (Λ m+1 ), as in Fig. 3 and 4.

Validation of implementation
Figure 8. Validation of implementation for a simple test case.L 2 error for the glacier height H with respect to the mesh size h when setting ∆t = h.Diamonds : the glacier height is computed using eq.( 25); crosses : the glacier height is computed using using SUPG; dots : slope proportional to h.
The analysis and simulation of the horizontal velocities u, v has already been validated in [22].Our goal is to validate the coupling between the velocity and the shape of the glacier.The numerical procedure is tested with an exact solution.We choose the horizontal velocity components u(x, y, z) = yz and v(x, y, z) = xz and the glacier height The glacier support Λ(t) is then the elliptic domain centered at the origin and with large and small axis equal to 1 + t and 1, respectively.Then we set the corresponding right-hand sides in (2.1) (2.2) and (2.11).We run the simulation from time t = 0 to t = t M = 0.5.Let e be the L 2 error for the glacier height H at final time : In Figure 8 we have reported the error e with respect to the mesh size h when setting the time step ∆t = h.Numerical results show that the L 2 error for the glacier height H at final time is of order h.Also, we have compared the formulation (3.13) to the classical SUPG formulation of [5].The SUPG scheme is more accurate but both methods seem to display the same convergence order.However, the SUPG method has not been used for the numerical simulation of Storglaciaren hereafter since it produces oscillations of the free surface which prevents robust simulations to be performed.balance for each year, the summer value (negative) is simply removed to the winter value (positive).The glacier geometry was measured for the years 1959, 1969, 1980 and 1990 [16, 17].We run our simulation starting with the glacier geometry of 1959, see Fig. 9, the mesh having 1200 vertices.For the missing mass balance b of 1959 and 1960 we use the ones of 1961.The time step is half a year, we stop the simulation in 1990 (the CPU time is less than one hour on a Pentium M processor 2.13GHz) and compare the obtained glacier shape to the measured one.The shape of the glacier during the simulation is shown in Fig. 10, the projection of the meshes onto the xy plane are reported in Fig. 11.The difference between computed and measured glacier thicknesses of years 1969, 1980 and 1990 are reported in Fig. 12 and is less than 10 meters (less than 5%), except for a small region located close to the highest point of the glacier.These results are very similar those obtained in [1], where a finite difference method based on the algorithm due to [9] was used.

Figure 2 .
Figure 2. Mass conservation into a vertical cylinder with base D in the xy plane.The normal to ∂D has component ν x and ν y in the xy plane.The amount of mass due to snow accumulation or melting on the top surface is b.
, a vertical cylinder with base D. Mass conservation into this cylinder writes y)∈∂D z=S(x,y,t) z=B(x,y)(uν x + vν y )dzds(x, y) = (x,y)∈D b(x, y, t)dxdy,where b is a given quantity corresponding to mass accumulation or ablation (snow falls or ice melting) that can be measured by glaciologists.Since S = B + H, the divergence theorem applied to domain D yields

Figure 3 .
Figure 3. Notations for the finite element tetrahedral mesh of Ω m .The glacier support Λ m is meshed into triangles with vertices (x m j , y m j ).Then, the vertices of the 3D mesh are aligned on top of those of the glacier support Λ m .

Figure 4 .
Figure 4.An example of finite element meshes.Top : mesh of the glacier support Λ m into triangles.Bottom : mesh of the ice glacier domain Ω m into tetrahedrons.

Figure 5 .
Figure 5. Moving the vertices lying on the boundary of Λ m .Left : moving vertices in the direction normal to the boundary may result in crossover in the case of large curvature.Right : changing the direction prevents crossover.

Figure 7 .
Figure 7. Mesh deformation of the support of Storglaciären, Sweden.Top : The mesh T h (Λ m ) of the glacier support Λ m .Middle : the vertices of the mesh T h (Λ m ) lying on the boundary of Λ m are moved (the deformation is exaggerated).Bottom : the internal vertices of T h (Λ m ) are moved according to the Laplacian method, which yields the new mesh T h (Λ m+1 ).

Figure 11 .
Figure 11.Meshes of Storglaciären, Sweden, view from above.The triangulation T h (Λ m ) of the glacier support Λ m is reported.Top left : year 1960, top right : 1970, bottom left : 1980, bottom right 1990.

Figure 12 .
Figure 12.Differences between measured and computed ice thickness.for years 1969 (top left), 1980 (top right) and 1990 (bottom left).Positive values indicate that the measured ice thickness is larger that the computed one.The caption is in meters and %.