[next] [tail] [up]

3.1 Ice Thickness Evolution

The evolution of the ice thickness, H, stems from the continuity equation and can be expressed as

∂H- -- ∂t = - ∇ ⋅(uH )+ B,
(3.1)

where u is the vertically averaged ice velocity, B is the surface mass balance and ∇ is the horizontal gradient operator (Payne and Dongelmans, 1997).

For large–scale ice sheet models, the shallow ice approximation is generally used. This approximation states that bedrock and ice surface slopes are assumed sufficiently small so that the normal stress components can be neglected (Hutter, 1983). The horizontal shear stresses (τxz and τyz) can thus be approximated by

∂s- τxz(z) = - ρg(s- z)∂x , ∂s- τyz(z) = - ρg(s- z)∂y ,
(3.2)

where ρ is the density of ice, g the acceleration due to gravity and s = H + h the ice surface.

Strain rates ˙εij of polycrystalline ice are related to the stress tensor by the non–linear flow law:

( ) ˙ε = 1 ∂ui + ∂uz- = A (T *)τ(n-1)τ i = x,y, iz 2 ∂z ∂i * iz
(3.3)

where τ* is the effective shear stress defined by the second invariant of the stress tensor, n the flow law exponent and A the temperature–dependent flow law coefficient. T* is the absolute temperature corrected for the dependence of the melting point on pressure (T* = T + 8.7 ⋅ 10-4(H + h - z), T in Kelvin, Huybrechts, 1986). The parameters A and n have to be found by experiment. n is usually taken to be 3. A depends on factors such as temperature, crystal size and orientation, and ice impurities. Experiments suggest that A follows the Arrhenius relationship:

A (T*) = fae-Q ∕RT *,
(3.4)

where a is a temperature–independent material constant, Q is the activation energy for creep and R is the universal gas constant (Paterson, 1994). f is a tuning parameter used to ‘speed–up’ ice flow and accounts for ice impurities and the development of anisotropic ice fabrics (Payne, 1999; Tarasov and Peltier, 1999, 2000; Peltier et al., 2000).

Integrating (3.4) with respect to z gives the horizontal velocity profile:

∫z u(z)- u (h) = - 2(ρg)n|∇s |n-1∇s A(s- z)ndz, h
(3.5)

where u(h) is the basal velocity (sliding velocity). Integrating (3.5) again with respect to z gives an expression for the vertically averaged ice velocity:

∫s∫z uH = - 2(ρg)n|∇s |n- 1∇s A(s- z)ndzdz′. h h
(3.6)

The vertical ice velocity stems from the conservation of mass for an incompressible material:

∂ux-+ ∂uy+ ∂uz= 0. ∂x ∂y ∂z
(3.7)

Integrating (3.7) with respect to z gives the vertical velocity distribution of each ice column:

∫z w (z) = - ∇ ⋅u(z)dz + w(h), h
(3.8)

with lower, kinematic boundary condition

w (h ) = ∂h-+ u(h)⋅∇h + S, ∂t
(3.9)

where S is the melt rate at the ice base given by Equation (3.55). The upper kinematic boundary is given by the surface mass balance and must satisfy:

∂s- w(s) = ∂t + u(s)⋅∇s + B.
(3.10)

3.1.1 Numerical Grid

The continuous equations descrining ice physics have to be discretised in order to be solved by a computer (which is inherently finite). This section describes the finite–difference grids employed by the model.

Horizontal Grid

The modelled region (x ∈ [0,Lx], y ∈ [0,Ly]) is discretised using a regular grid so that xi = (i- 1)Δx for i ∈ [1,N] (and similarly for yj). The model uses two staggered horizontal grids in order to improve stability. Both grids use the same grid spacing, Δx and Δy, but are off-set by half a grid (see Fig. 3.1).

pict

Figure 3.1: Horizontal Grid.

Quantities calculated on the (r,s)–grid are denoted with a tilde, i.e. F˜. Quantities are transformed between grids by averaging over the surrounding nodes, i.e. a quantity in the (i,j)–grid becomes in the (r,s) grid:

˜ Fr,s = ˜ Fi+12,j+12 = 1 4(Fi,j + Fi+1,j + Fi+1,j+1 + Fi,j+1) (3.11a)
and similarly for the reverse transformation:
Fi,j = Fr-12,s-12 = 1 4(F˜r-1,s-1 + ˜Fr,s-1 + ˜Fr,s + F˜r-1,s) (3.11b)

In general, horizontal velocities and associated quantities like the diffusivity are calculated on the (r,s) grid, ice thickness, temperatures and vertical velocities are calculated on the (i,j)–grid.

Horizontal gradients are calculated on the (r,s)–grid, i.e. surface gradients are:

( ) ∂s- ∂xr,s = ˜s r,sx = si+1,j --si,j +-si+1,j+1 --si,j+1 2Δx (3.12a)
( ) ∂s- ∂yr,s = ˜s r,sy = si,j+1 --si,j +-si+1,j+1 --si+1,j 2Δy (3.12b)
Ice thickness gradients, H˜r,sx and H˜r,sy, are formed similarly. Gradients in the (r,s)–grid are formed in a similar way, e.g.
( ) ∂u- x ˜ur,s-1---˜ur--1,s-1 +-˜ur,s --˜ur-1,s ∂x i,j = ui,j = 2Δx
(3.13)

Periodic Boundary Conditions

The model can be run with horizontal periodic boundary conditions, i.e. the western edge of the modelled region is joined with the eastern edge. Figure 3.2 illustrates the numeric grid when the model is run in torus mode.

pict

Figure 3.2: A row of the numeric grid when the model is used in torus mode. Circles indicate points in (i,j)–grid and squares indicate points in the (r,s)–grid. Points with the same colour are logically the same.

These boundary conditions are enforced by exchanging points for the temperature and vertical velocity calculations. The ice thicknesses are calculated explicitly at the ghostpoints.

σ–Coordinate System

The vertical coordinate, z, is scaled by the ice thickness analogous to the s–coordinate in numerical weather simulations (e.g. Holton, 1992). A new vertical coordinate, σ, is introduced so that the ice surface is at σ = 0 and the ice base at σ = 1 (see Fig. 3.3), i.e.

s--z- σ = H .
(3.14)

pict

Figure 3.3: Vertical scaling of the ice sheet model. The vertical axis is scaled to unity. The horizontal coordinates are not changed.

The derivatives of a function f in (x,y,z,t) become in the new (˜x,ỹ,σ,˜t ) system:

∂f ∂x- = ∂f ∂x˜ + 1 H-Δ˜x∂f ∂σ-, (3.15a)
∂f --- ∂y = ∂f --- ∂y˜ + 1 -- HΔỹ∂f --- ∂σ, (3.15b)
∂f- ∂t = ∂f- ∂˜t + -1 HΔ˜t ∂f- ∂σ, (3.15c)
∂f- ∂z = -1- H∂f- ∂σ, (3.15d)
where the geometric factors, Δ˜x, Δỹ and Δ˜t , are defined by
Δ˜x = ( ∂s ∂H ) ∂˜x-- σ∂-˜x, (3.16a)
Δỹ = ( ) ∂s-- σ∂H- ∂˜y ∂y˜, (3.16b)
Δ˜t = ( ) ∂s- ∂H- ∂˜t - σ∂ ˜t. (3.16c)
The integral of z becomes in the σ–coordinate system:
z σ ∫ ∫ fdz = - H f dσ h 1
(3.17)

The vertical coordinate is discretised using an irregular grid spacing to reflect the fact that ice flow is more variable at the bottom of the ice column. In the vertical the index k is used.

3.1.2 Ice Sheet Equations in σ–Coordinates

The horizontal velocity, Equation (3.5), becomes in the σ–coordinate system

σ n n+1 n-1 ∫ n u (σ) = - 2(ρg) H |∇s | ∇s A σ dσ +u (1) 1
(3.18)

and the vertically averaged velocity

∫1 uH = H udσ+ u (1)H 0
(3.19)

The vertical velocity, Equation (3.8), becomes

σ∫ ( ) w (σ) = - ∂u-⋅(∇s - σ∇H )+ H ∇ ⋅u dσ +w (1) 1 ∂ σ
(3.20)

and lower boundary condition

∂h w (1) = ∂t-+ u(1)⋅∇h + S.
(3.21)

3.1.3 Calculating the Horizontal Velocity and the Diffusivity

Horizontal velocity and diffusivity calculations are split up into two parts:

u(σ) = c∇s + u(1) (3.22a)
D = H ∫ 01cdσ (3.22b)
q = D∇s + Hu(1) (3.22c)
with
c(σ) = -2(ρg)nHn+1|∇s|n-1 ∫ 1σAσndσ (3.22d)

Quantities u and D are found on the velocity grid. Integrating from the ice base (k = N - 1), the discretised quantities become

˜cr,s,N = 0
(3.23a)

˜cr,s,k = - 2(ρg)nHn+1 ((˜sx )2 + (˜sy )2)n-21 r,s r,s r,s k ( ) ∑ Ar,s,κ-+-Ar,s,κ+1- σ-κ+1-+-σκ n 2 2 (σκ+1 - σκ) κ=N-1
(3.23b)
N∑- 1 D˜r,s = Hr,s ˜cr,s,k +-˜cr,s,k+1(σk+1 - σk ) k=0 2
(3.23c)

Expressions for ui,j,k and qi,j are straight forward.

3.1.4 Solving the Ice Thickness Evolution Equation

Equation (3.1) can be rewritten as a diffusion equation, with non–linear diffusion coefficient D:

∂H- = - ∇ ⋅D ∇s + B = - ∇ ⋅q+ B ∂t
(3.24)

This non–linear partial differential equation can be linearised by using the diffusion coefficient from the previous time step. The diffusion coefficient is calculated on the (r,s)–grid, i.e. staggered in both x and y direction. Figure 3.4 illustrates the staggered grid. Using finite differences, the fluxes in x direction, qx become

qi+12,jx = -1 - 2(D˜r,s + D˜r,s-1)s - s -i+1,j---i,j Δx (3.25a)
qi-12,jx = -1 - 2(D˜r-1,s + D˜r-1,s-1)s - s -i,j---i--1,j Δx (3.25b)
and the fluxes in y direction
qi,j+1 2y = -1 2(D˜r,s + D˜r-1,s)si,j+1 --si,j Δy (3.25c)
qi,j-1 2y = -1 2(D˜r,s-1 + D˜r-1,s-1)si,j --si,j--1 Δy. (3.25d)

pict

Figure 3.4: Illustration of the staggered grid used to calculate ice thicknesses, diffusivities and mass fluxes.

ADI Scheme

The alternating–direction implicit method (ADI) uses the concept of operator splitting where Equation (3.24) is first solved in the x–direction and then in the y–direction, (Press et al., 1992). The time step Δt is devided into two time steps Δt∕2. The descretised version of Equation (3.24) becomes (Huybrechts, 1986):

2Ht+12 - Hti,j -i,j-------- Δt = - 1 1 qx,i+t+1,2j - qxi,t-+ 12,j --2-------2-- Δx -qyi,t,j+1- qyi,,jt- 1 ----2------2- Δy + Bi,j (3.26a)
2 1 Hti+,j1- Ht+i,j2 ----Δt------ = - 1 1 qx,i+t+1,2j - qxi,t-+ 12,j --2--Δx---2-- -qyi,t,j++11- qyi,,jt+-1 1 ----2Δy----2- + Bi,j (3.26b)
Gathering all t + 1 2 terms on the left side, Equation (3.26a) can be expressed as a tri–diagonal set of equations for each row j:
t+12 t+ 12 t+12 - αi,jHi-1,j + (1- βi,j)H i,j - γi,jHi+1,j = δi,j
(3.27)

with

αi,j = D˜r--1,s +-D˜r--1,s-1 4Δx2Δt (3.28a)
βi,j = -˜Dr,s + 2˜Dr-1,s +D˜r -1,s-1 ------------2---------- 4ΔxΔt = -(αi,j + γi,j) (3.28b)
γi,j = ˜ ˜ Dr,s +-Dr2,s-1 4ΔxΔt (3.28c)
and the RHS,
( ) δ = Ht - -Δt- qy,t 1 - qy,t1 + Δt-B + α h - β h +γ h . i,j i,j 2Δy i,j+2 i,j-2 2 i,j i,j i-1,j i,j i,j i,j i+1,j
(3.28d)

A similar tri–diagonal system is found for each column, i of Equation (3.26b).

Linearised Semi–Implicit Scheme

Using the Crank–Nicolson scheme, the semi–implicit temporal discretisation of (3.24) is then:

Ht+1 - Ht qx,t+11- qx,t+11 qy,t+11- qy,t+11 -i,j-----i,j-= -i+-2,j---i-2,j-+ -i,j+-2---i,j-2- Δt 2Δx 2Δy qxi+,t1,j - qxi-,t1,j qy,i,tj+ 1- qyi,,tj- 1 + ---22Δx---2--+ ----22Δy----2-+Bi,j
(3.29)
The superscripts t and t+1 indicate at what time the ice thickness H is evaluated. Collecting all Ht+1 terms of (3.29) on the LHS and moving all other terms to the RHS we can rewrite (3.29) as
t+1 t+1 t+1 t+1 t+1 - αi,jHi-1,j - βi,jH i+1,j - γi,jHi,j-1 - δi,jHi,j+1 +(1 - εi,j)Hi,j = ζi,j
(3.30)

with the RHS,

t t t t t ζi,j = αi,jHi-1,j + βi,jH i+1,j + γi,jHi,j-1 +δi,jH i,j+1 + (1 + εi,j)H i,j + 2(αi,jhi- 1,j + βi,jhi+1,j + γi,jhi,j-1 + δi,jhi,j+1 + εi,jhi,j)+ Bi,jΔt
(3.31)
with the elements of the sparse matrix
αi,j = ˜ ˜ Dr--1,s +-Dr-1,s--1 4Δx2Δt (3.32a)
βi,j = D˜r,s +-˜Dr,s--1 4Δx2Δt (3.32b)
γi,j = D˜r,s-1 + ˜Dr-1,s- 1 ------4Δy2-------Δt (3.32c)
δi,j = ˜ ˜ Dr,s +-Dr2-1,s 4ΔyΔt (3.32d)
εi,j = -(αi,j + βi,j + γi,j + δi,j) (3.32e)

This matrix equation is solved using an iterative matrix solver for non-symmetric sparse matrices. The solver used here is the bi–conjugate gradient method with incomplete LU decomposition preconditioning provided by the SLAP package.

Non–Linear Scheme

The non–linearity of Equation (3.24) arises from the dependance of D on s. A non–linear scheme for (3.24) can be formulated using Picard iteration, which consists of two iterations: an outer, non–linear and an inner, linear equation. The scheme is started off with the diffusivity from the previous time step, i.e.

(0),t+1 t D = D
(3.33a)

and Equation (3.30) becomes

- α(ξ),t+1Ht+1 - β(ξ),t+1H (ξ+1),t+1- γ(ξ),t+1H (ξ+1),t+1 i,j i- 1,j i,j i+1,j i,j i,j- 1 - δ(i,ξ)j,t+1H (ξi,+j+11),t+1+ (1- ε(ξi,)j,t+1)H (iξ,j+1),t+1 = ζ(0i,)j,t
(3.33b)
Equation (3.33b) is iterated over ξ until the maximum ice thickness residual is smaller than some threshold:
( || (ξ+1),t+1 (ξ),t+1||) max |H - H | < Hres
(3.34)

pict

Figure 3.5: Flow diagram showing how the linearised solver (on the left) and the non–linear solver work. The inner, linear iteration is contained within the box labeled “calculate new ice distribution”.

3.1.5 Calculating Vertical Velocities

Grid Velocity

The vertical grid moves as a consequence of using a σ–coordinate system. The grid velocity is

( ) wgrid(σ) = ∂s+ u ⋅∇s - σ ∂H- + u⋅∇H ∂t ∂t
(3.35)

The numerical implementation of Equation (3.35) is straight–forward.

Vertical Velocity

The discretised version of the vertical velocity equation (3.20) is slightly more compilicated because the horizontal velocities are calculated on the (r,s) grid. The vertical velocity at the ice base is wi,j,N = wi,j,Ngrid - bi,j, where bi,j is the basal melt rate. Integrating from the bottom, the vertical velocity is then

∑1 { ( ux + ux vy + vy ) wi,j,k = - Hi,j -i,j,k----i,j,k+1+ -i,j,k----i,j,k+1 (σk+1 - σk) ˜k=N -1 2 2 ( x 1 x) + (˜ui,j,k+1 - ˜ui,j,k) ˜si,j - 2(σk+1 + σk)˜Hi,j ( )} + (˜vi,j,k+1 - ˜vi,j,k) ˜syi,j - 1(σk+1 +σk)H˜yi,j + wi,j,N 2
(3.36)

with the weighted ice thickness

H = 4Hi,j +-2(Hi-1,j-+-Hi+1,j +-Hi,j--1 +-Hi,j+1) i,j 16 + Hi-1,j-1-+Hi+1,j-1 +-Hi+1,j+1-+-Hi-1,j+1 16

This scheme produces vertical velocities at the ice divide which are too small. The vertical velocities on the ice surface are given by the upper kinematic boundary condition, Equation (3.10). Equation (3.36) can be corrected with:

w*i,j,k = wi,j,k - (1- σk)(wi,j,k - wsi,j),
(3.37)

where wsi,j is the vertical velocity at the ice surface given by (3.10). Figure 3.6 shows the different vertical velocities at the ice surface.

pict

Figure 3.6: Vertical ice surface velocities of the EISMINT-1 moving margin experiment.

The difference between the vertical velocities calculated by the model and the vertical velocities given by (3.10) at the ice margin are due to the fact that temperatures and velocities are only calculated when the ice is thicker than a certain threshold value which is not met at the ice margin.

Figure 3.7 shows vertical profiles of the vertical velocity at the ice divide and a point half–way between the divide and the domain margin. A corresponding temperature profile is also shown since the vertical velocity determines the vertical temperature advection (see Section 3.2.4).

pict

Figure 3.7: Vertical velocity and temperature distribution for columns at the ice divide and a point half–way between the divide and the domain margin.

[next] [front] [up]