A moving mesh approach to an ice sheet model

Computers & Fluids 46 (2011) 381–386

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

A moving mesh approach to an ice sheet model D. Partridge ⇑, M.J. Baines School of Mathematics, Meteorology and Physics, University of Reading, Reading, United Kingdom

a r t i c l e

i n f o

Article history: Received 13 May 2010 Received in revised form 1 February 2011 Accepted 1 February 2011 Available online 26 February 2011 Keywords: Moving mesh Ice sheets Finite differences

a b s t r a c t A moving mesh approach to the numerical modelling of problems governed by nonlinear time-dependent partial differential equations (PDEs) is applied to the numerical modelling of glaciers driven by ice diffusion and accumulation/ablation. The primary focus of the paper is to demonstrate the numerics of the moving mesh approach applied to a standard parabolic PDE model in reproducing the main features of glacier flow, including tracking the moving boundary (snout). A secondary aim is to investigate waiting time conditions under which the snout moves.  2011 Elsevier Ltd. All rights reserved.

1. Introduction/Background In this paper we take a standard one-dimensional PDE model of a glacier and discretise it using a finite difference moving mesh method. The aim is to show that the numerical scheme can reproduce the features of the model and in particular handle the phenomenon of waiting time which is a known feature of glacier movement. Computational studies of glaciers are particularly challenging to the modeller. Although ice sheet models are well-established, prediction of profiles and grounding movement are infeasible analytically and difficult to achieve numerically, see Payne and Vieli [5]. Ice moves in a similar manner to a viscous fluid, though with a very high viscosity approximately 1015 times that of water [2]. However, viscous theory cannot solely be used to describe flow, since glaciers are unique in experiencing basal sliding. This can be caused in two ways, via friction where the ice makes contact with the ground as it is flowing, or geothermal heat below the surface. In order for glaciers to form they first need enough snow over the winter period to be able to survive through the summer, i.e. more accumulation of snow than is lost through melting and evaporation. This needs to be repeated over a number of successive years, and as more snow builds up, the weight increases and pressure compresses the firn (old snow) into ice. Once this ice is thick enough, gravity, amongst other forces, causes the ice to flow. This is a long, complex process which takes less time in regions where temperature changes quicker, such as the Alps and North America [2]. ⇑ Corresponding author. E-mail address: [email protected] (D. Partridge). 0045-7930/$ - see front matter  2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2011.02.001

On a global scale, ice quantities vary considerably. At present glaciers make up around 2% of the Earths water, but during an ice age this vastly increases. Either way they have a large impact on the climate system, and are becoming increasingly affected by climate change. If all this ice melted into the oceans, there would be a sea level rise of around 70 m. We are interested in glaciers for more than just the climate change reasons, as they can have a large effect on the local terrain, causing events such as landslides and flash floods. We discritise a standard PDE model of glacier movement in a moving frame of reference, on a moving mesh, using a local mass balance principle to define a velocity in order to move the mesh. We note that, as in other nonlinear diffusion problems, glaciers experience a waiting time before they begin to move. We suggest a mechanism whereby waiting ends and the snout moves. Finally, consideration is given to ways the model may be extended, and the impacts that these extensions may have on the results we have obtained, leading the way to potential further work to be undertaken on the problem. One of the main concepts to take into consideration when modelling glaciers is the idea of mass balance, and where on the glacier mass is gained or lost. Generally, near the source of the glacier, the accumulation of snow is greater than the ablation (melting/evaporation), so the mass increases. Further away the ablation becomes greater than the accumulation, and the mass decreases. However ice can build up in the lower zone due to ice flow coming from the glacier’s upper zone. The front-most end of the glacier is known as the snout, which rarely moves straight away; it waits until the velocity behind it is great enough to push it down the mountain. It is this feature which is of special interest.

382

D. Partridge, M.J. Baines / Computers & Fluids 46 (2011) 381–386

2. Model description A standard PDE model for glaciers was proposed by Oerlemans [3] in 1984, with a flat bed occupying the region x 2 [0,b(t)] as shown in Fig. 1. Let H(x, t) represent the thickness of the ice. At the ends of this domain we have the boundary conditions, @H ¼0 @x at the fixed point x = 0, and H = 0 at the moving boundary x = b(t). In one dimension the continuity equation for ice can be written as

@H @ðHuÞ ¼ þ sðxÞ; @t @x

with sdx the stress term, and parameters A and n taken from Glen’s flow law, an established general law for steady state ice deformation [6]. From Van Der Veen [9] the driving stress is given by

sdx ¼ qgH

u¼

 n 2AH n n n @H qgH : nþ2 @x

u ¼ cHnþ1

Z

bðtÞ

Hðx; tÞdx ¼ hðtÞ; say:

ð2Þ

0

From (1), using Leibniz’s integral rule, and applying the boundary conditions

d dt

Z

bðtÞ

Z

bðtÞ

@H dbðtÞ dx þ HðbðtÞ; tÞ @t dt 0 Z bðtÞ Z bðtÞ @ ¼ sðxÞdx ½Hudx þ @x 0 0 Z bðtÞ Z bðtÞ bðtÞ sðxÞdx ¼ sðxÞdx; ¼ ½Hu0 þ

Hðx; tÞdx ¼

0

0

ð3Þ

0

the physical equivalent of which states that any change in the integral of ice thickness over the whole glacier, or equivalently any change in the ice volume, is due only to the snow term, which represents the net accumulation/ablation of snow over the whole glacier.

ð6Þ

The parameters A, n, q and g are set as constant to simplify the model, giving

2.1. Mass balance An important property concerns the integral of the ice thickness over the whole domain (the volume), i.e.

ð5Þ

with q the ice density, g representing gravity, and h equal to ice thickness plus the surface elevation. On a flat bed there is no surface elevation so we may put h = H. From (4) and (5) we get an equation for the depth integrated horizontal velocity

ð1Þ

where H is the ice thickness, s(x) = sa(x)  sb(x), with sa the accumulation rate of snow and sb the basal melting rate.

@h ; @x



@H @x

n ð7Þ

;

where c is a single positive constant parameter. Expressing the velocity in this form may present problems when dealing with the boundary condition H = 0 at x = b(t), apparently giving a zero velocity at the right-hand boundary and resulting in a glacier that will never move, which we know physically is not the case. However it is perfectly possible for u to be non-zero as long as Hnþ1 Hnx is finite, which requires Hx to be infinite. For the most part though we are not concerned with physical values for the variables, but more with the numerical behaviour of the moving mesh approach, hence x and H are non-dimensionae ¼ 103 H. We drop the tilde notation for lised, with ~ x ¼ 106 x and H convenience. From Roberts [7] we set c = 0.000022765 in standard SI units, and n = 3. Substituting the velocity into Eq. (1) we get the model equation

@H @ h 5 3i ¼c H Hx þ sðxÞ; @t @x

ð8Þ

which incorporates nonlinear diffusion and a source term. In this paper we set sb(x)  0, making s(x) = sa(x), although the non-zero basal melting case is considered in [4].

2.2. Model velocity 3. Snout behaviour The model velocity u is defined as the mean depth integrated horizontal velocity, and assuming lamellar flow it is given by [9]

2AH n u¼ s ; n þ 2 dx

ð4Þ

From (7), with n = 3 we derive the useful form

u ¼ cðH4=3 Hx Þ3 ¼ 

27 c½ðH7=3 Þx 3 : 343

ð9Þ

When expressing the velocity in this manner it is interesting to substitute an expression for H that has the right general shape and satisfies the boundary conditions, i.e.

H ¼ ð1  x2 Þa

s(x)

ð10Þ

where a > 0, for which 7

7a

H3 ¼ ð1  x2 Þ 3 H(x,t)

ðH7=3 Þx ¼ 2x:

7a 7a ð1  x2 Þ 3 1 : 3

ð11Þ

The velocity (9) then has some interesting properties as x ? 1, depending on the value of a.

0

x(t) Fig. 1. One-dimensional domain.

b(t)

Case 1 :

7a > 1; ) ðH7=3 Þx is zero 3

ð12Þ

Case 2 :

7a < 1; ) ðH7=3 Þx is infinite 3

ð13Þ

Case 3 :

7a ¼ 1; ) ðH7=3 Þx is finite 3

ð14Þ

383

D. Partridge, M.J. Baines / Computers & Fluids 46 (2011) 381–386

In case 1, from Eq. (9), the initial velocity of the snout of the glacier is zero, and it remains stationary. In case 2 the velocity is infinite, which is not physical. In case 3 we get a finite velocity value at the snout when a ¼ 37, and the right-hand boundary moves. From this analysis we can expect H to be asymptotically of the form H = (1  x2)3/7 for (1  x) small at the moment of initial movement, i.e. if H approaches case 3 asymptotically at the snout then the boundary will move.

Computation is performed with initial conditions (10) with a set to 1. The snow term is approximated for all time by the linear function (as in Van Der Veen [9])

sðxÞ ¼ eð1  dxÞ;

ð18Þ

where d and e are the snow parameters, set to be 0.5 and 0.05, respectively. The model is run for a sufficient length of time for the boundary to wait, then move. The initial mesh is chosen to be evenly spaced.

4. Numerical method 4.3. Explicit time-stepping 4.1. Moving framework-subdomain mass balance Eq. (8) is generally impossible to solve analytically, so we seek a numerical approximation via a mesh. Since the problem involves a moving boundary, a natural description is to use a moving framework, for which we define a two-stage method, firstly by using a velocity v(x, t) to move an arbitrary point ^ x of the glacier [1], then updating the ice thickness using a mass balance principle. To define this velocity we assume that the mass balance principle (3) of the model holds in any moving subdomain ½0; ^ xðtÞ of [0, b(t)]. In physical terms this velocity v is such that the partial ice volume changes only due to the local accumulation/ablation of snow. In equation form we assume that

d dt

Z

^xðtÞ

Hðx; tÞdx ¼

0

Z

^xðtÞ

sðxÞdx

ð15Þ

0

for each subdomain ð0; ^xðtÞÞ, consistent with the model property (3). We now show that a consequence of (15) is that v(x, t) = u(x, t). By Leibniz’s integral rule, making use of (8) and the boundary conditions given in Section 2,

^x Z ^xðtÞ Z ^xðtÞ d @H d^xðtÞ dx þ Hð^xðtÞ; tÞ Hðx; tÞdx ¼  dt 0 @t dt  0 0 Z ^xðtÞ d^xðtÞ 5 3 ¼ cH Hx j^x þ sðxÞdx þ Hð^xðtÞ; tÞ : dt 0

To advance the node positions ^ xi ðtÞ from the velocity (17) we use an explicit Euler scheme. Letting k denote the time discretisation level,

h  i3 k ^xkþ1  ^xki i H7=3 ; ¼ c2 i x Dt

ð19Þ

Therefore, dropping the hat notation for convenience, at each time step we update each mesh point by:

xkþ1 ¼ xki  c2 Dt i

h  i3 k Hi7=3 ;

ð20Þ

x

estimating ðH7=3 i Þx using an upwind difference scheme. Since we do not seek high levels of accuracy Euler time-stepping is sufficient, provided the time step is suitably small to ensure stability. To determine the updated ice thickness we use the same timestepping scheme on Eq. (15). Note that the limits have been chosen to give an incremental form.

hR

xjþ1 xj1

ikþ1 hR ik xjþ1 Hdx  xj1 Hdx

Dt

¼

"Z

xjþ1

#k sdx :

xj1

Using the midpoint rule we obtain the approximation

ð16Þ

      kþ1 kþ1 xjþ1  xkþ1  xkjþ1  xkj1 Hkj ¼ Dt xkjþ1  xkj1 skj ; j1 Hj

Therefore the assumption (15) is equivalent to

giving

 Z ^xðtÞ  d^xðtÞ dx ¼ 0 cH5 H3x þ Hð^xðtÞ; tÞ dt 0

Hkþ1 j

  xkjþ1  xkj1  ðHkj þ Dtskj Þ: ¼ kþ1  x xkþ1 jþ1 j1

ð21Þ

which, since ^xðtÞ is arbitrary, gives

cH5 H3x þ Hð^xðtÞ; tÞ

d^xðtÞ ¼ 0: dt

Hence the velocity and we obtain

v¼

v ¼ d^x=dt

4.4. Convergence considerations is driven only by the diffusion term

d^xðtÞ cH5 H3x ¼ ¼ cH4 H3x ¼ c2 ½ðH7=3 Þx 3 dt H

ð17Þ

where c2 = (3/7)3c. Note that this velocity is the same as taking v to be the model velocity (7) at each of the nodes. Reversing the argument implies that the assumption (15) and the velocity (17) are equivalent. Following [8] it is reasonable to make the assumption that v = u at the boundary irrespective of the moving framework since the model velocity coincides with the ice velocity at all times. 4.2. Numerical method-finite differences To form our finite difference solution in the first part of our twostage method we discretise (17), updating the mesh positions accordingly. The second part uses Eq. (15) to update H. This cycle is performed for every time step.

As there is no exact solution available a reference solution is taken to check convergence. We choose a final time (=200), which comfortably allows the boundary to wait, then move. Here we use the solution with 321 mesh points as our reference. Table 1 shows that both the position of the moving boundary point and the L2 norm of the ice thickness do indeed converge towards the reference solution. We also estimate the rate of convergence using the reference solution. For a single boundary position we consider the relative

Table 1 Table of values for varying number of mesh points, with solutions converging towards the reference value. N

xb

kHk/N

21 41 81 161 Ref.

2.7173 2.7110 2.7083 2.7073 2.7068

.837005 .599566 .424638 .302860 .214531

384

D. Partridge, M.J. Baines / Computers & Fluids 46 (2011) 381–386

away. It is here that the diffusion term begins to dominate the flow. Similarly, the velocity (Fig. 3b) builds up more on the boundary initially, then as the outward movement speeds up, the peak velocity starts to return to the centre of the glacier. This is due to the glacier entering a region with less accumulation of snow, slowing down the end points.

Table 2 Error analysis showing decreasing errors for increasing mesh points.       ref   ref kHref  HNk/kHrefk N xb  xN b =xb  21 41 81 161

.003879 .001552 .000554 .000185

.00088069 .00053414 .00030948 .00014125

5. Waiting time

error, which in Table 2 is shown to be decreasing. Assuming e / Np we find that p  3/2. In the case of the ice thickness we also consider the relative error with respect to the reference. From Table 2 we again see that the error decreases for increasing mesh points, and in this case we find p  2.

We now discuss in more detail the mechanism whereby the snout waits and then moves. We define b0 as the initial location of the boundary, i.e. b(0). Since initially H > 0 for x < b0 and H = 0 at x = b0 " t, the function H(x) in the vicinity of x = b0 is of the form

HðxÞ ¼ ðb0  xÞa gðxÞ 4.5. Results The model is run with 51 mesh points (Dx = 0.02), with a time step Dt = 0.005. The time step is chosen because it satisfies the stability criteria [9]

"

Dt 6 min

ðDxÞ

2

ð23Þ

to leading order in (b0  x), where a > 0 is the highest power of (b0  x) in H(x) and g(x) > 0 therefore has a finite derivative at x = b0. Note that for the initial conditions in the computations g(x) = (1 + x)a. Hence to leading order in (b0  x),

HðxÞ7=3 ¼ ðb0  xÞ7a=3 GðxÞ

# ð22Þ

4cH5 H2x

Varying the parameter a in the initial conditions shows the snout profile behaviour of Section 3. In Fig. 2a we see that when a = 3/7 the gradient is effectively infinite at the boundary, while a comparison case of a > 3/7 shows a finite gradient. Similarly, looking at the initial velocities for each of the two cases we see in Fig. 2b that the boundary does not move when a > 3/7, while in the comparison case it does. Also of note is that the peak velocity is at the snout of the glacier under the required condition, compared with near the centre for the stationary boundary case. This suggests that over time the peak velocity moves towards the boundary before the glacier begins to move. We may also look at the results over a period of time. For the experiments shown, the value of a has been set to 3/7 and 60,000 steps were taken, or 300 time units. The results are plotted every 5000 steps. From Fig. 3a, the change in ice thickness is initially dominated by the build-up of snow. Even though the boundary is now moving, it is not until a larger amount of snow has fallen that the speed notably increases and the boundary moves further

(a)

ð24Þ

where G(x) = (g(x))7/3 > 0 also with a finite derivative at x = b0. The velocity (17) then takes the form

v ðxÞ ¼ ca3 ðb0  xÞ7a3 fGðxÞg3

ð25Þ

to leading order in (b0  x). When a = 3/7 the velocity reduces to

v ðxÞ ¼ cð3=7Þ3 fGðxÞg3 :

ð26Þ

We see from (25) that v ðb0 Þ ¼ limx!b0 v ðxÞ ¼ 0 when a > 3/7, while from (26) v(b0) is non-zero when a = 3/7. The discrepancy indicates a discontinuity in v(x) at x = b0 as a ? 3/7. Fig. 4 shows a schematic plot of v(x), normalised by c{G(x)}3, against x for small (b0  x) and various a, showing how the discontinuity forms. Since the velocity must remain a continuous function of x, the limit a = 3/7 cannot be attained and the snout must move. We now consider the transition in time of v ð^ x; tÞ to a = 3/7 from above for small ðb0  ^ xÞ under the subdomain mass balance assumption (15). To leading order in ðb0  ^ xÞ the time-varying form of (24) is

Hðx; tÞ ¼ ðb0  ^xÞa gð^x; tÞ

(b) 1.5

ð27Þ

x 10−5

1

1

0.6

V

H

0.8

0.4

0.5

0.2

0

0

0.1

0.2

0.3

0.4

0.5

x

0.6

0.7

0.8

0.9

1

0

0

0.2

0.4

X

0.6

0.8

1

Fig. 2. Initial values for a = 1 (solid) and a = 3/7 (dashed). (a) The ice thickness with evidence of the infinite gradient in the 3/7ths case while (b) shows initial velocity, with an obvious difference at the right-hand boundary.

385

D. Partridge, M.J. Baines / Computers & Fluids 46 (2011) 381–386

(a)

(b)

0.012

5

0.01

4

0.008

3

0.006

v

H

6

2

0.004

1

0.002

0

0

0.5

1

1.5

x

2

2.5

3

0

3.5

0

0.5

1

1.5

x

2

2.5

3

3.5

Fig. 3. Numerical results using explicit time-stepping. (a) The evolution of the ice thickness every 12000 steps and (b) the corresponding evolution of velocity every 1000 steps.

0.09

and thus on the trajectory given by (29)

0.08

a3 ¼

0.07

1 1 ðb0  ^xÞ10 v ðtÞ c ð2wðtÞÞ7

Since d^ x=dt ¼ v ðtÞ, differentiating a3 with respect to time along this trajectory,

0.06 0.05

v

3a2

0.04 α=3/7 α=49/112 α=25/56 α=13/28 α=1/2 α=15/28

0.03 0.02 0.01 0 0.9

0.92

0.94

x

0.96

0.98

1

Fig. 4. Plots of v against x for a = 15/28, 1/2, 13/28, 25/56, 49/112, 3/7 showing formation of the discontinuity.

where ^x ¼ xðtÞ, while that of (25) is

v ð^x; tÞ ¼ ca

3

ðb0  ^xÞ

7a3

3

fGð^x; tÞg :

ð28Þ

We consider the evolution of the velocity under the subdomain mass balance assumption (15) applied to the interval ð^x; b0 Þ, where ðb0  ^ xÞ is small. Let h(t) be the positive mass in the triangle consisting of points ð^x; Hðx; tÞÞ; ð^x; 0Þ and the fixed point (b0, 0). To leading order in ðb0  ^xÞ, therefore,

1 wðtÞ ¼ ðb0  ^xÞHðx; tÞ: 2

ð29Þ

da 1 10 ðb0  ^xÞ9 v ðtÞ2 ¼ c ð2hðtÞÞ7 dt

to leading order in ðb0  ^ xÞ, assuming that dw/dt is of the same spatial order as w(t) and dv/dt is of the same spatial order as v(t). Thus da/dt is strictly negative. Hence, for sufficiently small ðb0  ^xÞ; a is strictly decreasing and when a reaches 3/7 the boundary moves. It is possible to interpret v(x, t) in terms of (curved) characteristics. Although v is given explicitly by (17) it is useful to consider a PDE satisfied by v. Differentiating v we get

v x ¼ 3cH4 H2x Hxx  4cH3 H4x v t ¼ 3cH4 H2x Hxt  4cH3 H3x Ht

ð31Þ ð32Þ

To keep Eq. (32) in terms of space derivatives, we can substitute Eq. (17) into the original Eq. (1), with u replaced by v, giving

Ht ¼ s  ðHv Þx ¼ s  Hx v  Hv x ;

ð33Þ

which we can differentiate to get

Htx ¼ sx  Hxx v  2Hx v x  Hv xx :

ð34Þ

The expressions for Ht and Htx can be substituted into (32) to give

v t ¼ 3cH5 H2x v xx  10vv x  v x v  3cH4 H2x sx  4cH3 H3x s;

Hence, from (27), to leading order in ðb0  ^xÞ,

and finally rearranged in the form of a Burgers-like equation with additional source terms,

1 ðb0  ^xÞaþ1 gð^x; tÞ ¼ wðtÞ 2

v t þ 11vv x ¼ 3cH5 H2x v xx  3cH4 H2x sx  4cH3 H3x s:

so that on the trajectory given by (29) from (28) the velocity at time t can be written

v ðtÞ ¼ ca3 ðb0  ^xÞ7a3

2wðtÞ

!7

Hence

dv ¼ 3cH5 H2x v xx  3cH4 H2x sx  4cH3 H3x s dt

ð36Þ

along the (curved) characteristics given by

ðb0  xÞaþ1

¼ ca3 ð2wðtÞÞ7 ðb0  ^xÞ10

ð35Þ

ð30Þ

dx ¼ 11v ðx; tÞ dt

ð37Þ

386

D. Partridge, M.J. Baines / Computers & Fluids 46 (2011) 381–386

where H is obtained from dH/dt = Ht + (dx/dt)Hx using (37). The presence of the dissipative terms in dv/dt prevents a shock forming in the interior but from (37) the characteristics become infinitely steep as v(x, t) ? 0 at the boundary x = b0, consistent with the formation of the discontinuity in v(x, t) at the snout. 6. Discussion We have seen that the numerics of a moving mesh method based on local mass balance qualitatively imitates a standard one-dimensional glacier model. We also analysed the conditions required for the snout to move, which required an infinite slope condition to be met asymptotically at the boundary in order for the velocity at this point to be non-zero. For an initial condition in the form of the quadratic function (10) with a > 3/7, we were able to simulate the qualitative waiting time behaviour of the glacier as the power of a evolved to 3/7. However the model (8) to which we have applied our numerics has a number of limitations, most notably that since the velocity is proportional to the negative slope at the snout the glacier cannot retreat. In real glacier systems, it is observed that glaciers also break-up, where large sections of ice completely break off from the main part and drift away or melt. This also presents problems for a depthaveraged model such as the one considered here with the ice thickness vanishing at points other than at the boundary. We may however see a positive Hx gradient in the interior, which would generate a negative velocity. A mechanism to use a negative interior velocity to allow retreat and break-up to occur needs investigating. A number of other avenues also await exploration. The model we have been analysing throughout is only a 1D depth-averaged model, so a logical next step is to take the model into more dimensions. This could mean one of two things. The first is to consider a fully 3D model of the glacier instead of depth averaging. When the velocity varies with height we should see quite different movement which would require a far more complex model with mesh

points moving in the vertical in addition to the x direction. We would expect the top and bottom of the glacier to be moving quicker than the middle section, and the impact of basal sliding could then be analysed more effectively. Alternatively we could keep the depth averaged vertically and consider the domain in the (x, y) plane, as opposed to a cross-section in the x-plane as is modelled currently. This will require additional boundary conditions at the sides of the glacier, the type depending on whether the boundary is a solid wall (no flux condition) or whether the glacier just curves to the ground (H = 0). The model itself then takes the form

Ht ðx; y; tÞ ¼ r:½Hðx; y; tÞ5 rHðx; y; tÞ3  þ sðx; yÞ:

ð38Þ

A viable 2D numerical model using the same moving mesh approach is possible using finite element approximation. A further objective in this ongoing work is to introduce the concepts of Data Assimilation. Here the aim is to improve the model accuracy using observations (mostly taken remotely), and to track the moving boundary location. References [1] Budd C, Huang W, Russell R. Adaptivity with moving grids. Acta Numer 2009;18:111–241. [2] Hambrey M, Alean J. Glaciers. Cambridge University Press; 1992. [3] Oerlemans J. [last checked 23.11.09]. [4] Partridge D. Analysis and computation of a simple glacier model using moving grids. Master’s thesis, Department of Mathematics, University of Reading; 2009. [5] Payne AJ, Vieli A. Assessing the ability of numerical ice sheet models to simulate grounding line migration. J Geophys Res 2005;110. [6] Rae B, Irving D, Hubbard B, McKinley J. Preliminary investigations of centrifuge modelling of polycrystalline ice deformation. Ann Glaciol 2000;31. [7] Roberts R. Modelling glacier flow. Master’s thesis, Department of Mathematics, University of Reading; 2007. [8] Tezduyar T, Beh M. A new strategy for finite element computations involving moving boundaries and interfaces. The deforming-spatial-domain/space-time procedure: I. The concept and the preliminary numerical tests. Comput Methods Appl Mech Eng 1992;94. [9] Van Der Veen CJ. Fundamentals of glacier dynamics. Taylor and Francis; 1999.