Numerical modelling of ice sheets and ice shelves Tony Payne Centre for Polar Observation and Modelling University of Bristol, U.K.

1
Numerical modelling of ice sheets and ice
shelves Tony PayneCentre for Polar
Observation and ModellingUniversity of Bristol,
U.K.
2
Reasons for numerical modelling

• Underlying equations cannot be solved
  analytically.
• Several different equations that are coupled.
• Irregular geometries when modelling real ice
  masses.

3
Process of numerical modelling

• Continuum mechanics description of problem
• Scaling the equations
• Discretization (quantities no longer in
  continuous form but represented at finite
  locations)
• Model coding
• Validation of the model
• Use of the model

4
Validating of numerical modelling

• Results from the numerical model are only useful
  if they can be shown to be independent of the
  discretization process
• Need to show that results do not vary with time
  step and grid sizes
• Repeat experiments using a range of grid and time
  steps
• Numerical models should be tested against any
  analytical solutions available
• EISMINT benchmarks also very useful in model
  development

5
Types of numerical modelling

• Finite element
• Great freedom in dealing with irregular
  geometries
• Fairly complicated to develop although software
  now available to do many of standard tasks
• Computationally demanding
• Finite difference (finite volume)
• Simple to use and develop
• Problems in coping with highly irregular domains

6
Basis of finite differences

• Taylor series give a relation between continuous
  and finite representations
• The series is truncated after a certain number of
  terms
• The terms omitted give rise to the truncation
  error of the approximation
• When equations implemented also incur rounding
  error (due to finite representation of number in
  memory use dp)

7
Taylor series
H(t ?t)
H(t)
t?t
t
8
Basis of finite differences
The order of the approximation is the order of
the first omitted term
9
Basis of finite differences

• In general second-order approximations are
  sufficient
• In certain circumstances more accurate higher
  orders may be needed, e.g. near abrupt
  transitions

10
Stability of finite differences

• Linear instability Courant-Friedlichs-Lewy
  constraint on grid size for explicit techniques
• Non-linear instability arises when we try to
  solve non-linear problems using linear techniques
• In general it manifests itself as high-frequency
  peaks which grow rapidly

11
Potted history of numerical modelling

• 1970s Modelling of glaciers and ice sheets
• 1- and 2-d. models
• zero-order shallow ice approximation
• Mahaffy, Budd, Oerlmans, Pollard,
• Key techniques solving non-linear parabolic
  eqns. in 1 and 2d.

12
Potted history of numerical modelling

• 1980s Ice temperature evolution
• 3d. models
• still zero-order flow models
• Jensen, Huybrechts, Ritz, Greve
• Key techniques stretched vertical coordinate

13
Potted history of numerical modelling

• 1980s Modelling of ice shelves and streams
• 2-d. models
• first-order shallow ice approximation but
  vertically integrated
• MacAyeal, Determann, Thomas
• Key techniques solving coupled non-linear
  elliptic eqns. in 2d.

14
Potted history of numerical modelling

• 1990s Coupled system evolution
• 2d. Ice shelf models
• 3d. Ice sheet, temperature models
• Marshall, Huybrechts, Hulbe
• Key techniques internal boundaries

15
Potted history of numerical modelling

• 2000s General ice flow models
• 3-d. models
• first-order shallow ice approximation with no
  vertically integration
• Pattyn, Blatter, Gudmundsson
• Key techniques solving coupled non-linear
  elliptic eqns. in 3d.

16
Stresses

• Easiest to think of as a force per unit area.
• Distinct from pressure in that they have a
  direction (are vector quantities).
• Units are Pa or bar (1 bar 100 kPa)
• Normal stresses (usually denoted by s) act
  perpendicular to a surface.
• They can be tensile (positive, extensional) or
  compressive (negative).
• Shear stresses (t) act parallel to a surface.

17
Stress tensor

• Imagine a cube.
• There are three normal stresses corresponding to
  the three axes
• There are six possible shear stresses

18
Stress tensor

• Subscripts denote direction for normal stresses
• Direction perpendicular to surface and then
  direction acting for shears
• However, reduces to three shears if assume body
  in equilibrium and can not rotate

19
Stress equilibrium or force balance

• In glaciology we can neglect acceleration so that
  Newtons second law reduces to a force balance
• If there were a force imbalance, acceleration
  would operate instantaneously to restore
  equilibrium
• Balance is between internal stresses and body
  forces such as gravity

20
Equations of stress equilibrium

• Assume static balance of forces by ignoring
  acceleration

21
Equations of stress equilibrium
22
Equations of stress equilibrium
?z
?y
?x
23
Equations of stress equilibrium
?z
?y
?x
24
Equations of stress equilibrium
?z
?y
?x
25
Equations of stress equilibrium

• Assume static balance of forces by ignoring
  acceleration

26
Strain rates

• Strain is a measure of the deformation of an
  object
• Typically ratio of a length before deformation to
  one after (proportional change) and has no units
  (i.e. is a ratio)
• A strain rate is simply the rate at which this
  happens.
• Denoted by epsilon dot (dot shows that it is
  rate) and has units of yr-1

27
Strain rates

• Defined as gradients in velocity
• Strain can be normal (simple extension or
  contraction)
• Or shear which changes shape
• Six strain rates that are equivalent to six
  stresses

28
Glens flow law
relates to strain rates to stress deviators

• A is related to temperature and impurities
• n is a constant equal to 3
• is effective stress (measure of overall stress
  regime)
• ? is effective strain rate

29
Definition of stress deviators

• in the case of normal stress, the deviator equals
  the stress minus the pressure (mean of three
  normal stresses)
• in the case of shear stress, the deviator simply
  equals the shear stress

30
Effective stresses

• Measure of total stress environment
• Also for strains

31
Shallow-ice approximation(s)

• Provide a rigorous simplification of the full
  stress equilibrium equations.
• Based on use of scaling and in particular on the
  aspect ratio of an ice mass (thickness to span).
• The aspect ratio is small O(0.005) hence shallow
  ice.
• Three are 3 flavours (or orders) depending on
  what effects are left in.

32
Equations of stress equilibrium

• Assume static balance of forces by ignoring
  acceleration

33
Scaling

• Replace variable with a constant scale multiplied
  by a scaled variable
• allows size of different terms to be compared
  without further analysis
• Scaled variable should always be O(1)

34
Scaling
35
Word on stress deviator scaling

• Some justification for this scale comes from the
  condition at the ice mass upper/lower boundary
• At the upper surface, stress components must
  balance

36
Glaciology is a ONE BAR subject (or 100 kPa in
S.I.)
37
Vertical equation

• Scale equations (note have dropped tildas from
  scaled variables)

38
Horizontal equations

• Scale equations (note have dropped tildas from
  scaled variables)

39
Zero-order approximation

• Retain highest order (biggest) terms in each
  equation

40
Zero-order approximation

• Retain highest order (biggest) terms in each
  equation

41
Zero-order approximation

• Easier to move back to unscaled version of
  equations now

42
Higher-order approximations

• Successively include smaller terms

43
First-order approximation

• Successively include smaller terms

44
Second-order approximation

• Successively include smaller terms

45
Quick zero-order derivation

• assume on basis of aspect ratio that
• normal stress dominates in vertical
• vertical shear dominates horizontal shear

46
Quick zero-order derivation

• Integrate vertical stress balance
• It can be shown that all normal stresses are
  equal (i.e. that the pressure is hydrostatic)

47
Quick zero-order derivation

• Substitute horizontal normal stresses with
  vertical normal stress

48
Quick zero-order derivation

• Tidy up
• Gravitational driving stress (RHS) balanced
  entirely by vertical shear (LHS)
• For whole ice column gravitational driving is
  balanced locally by basal drag alone

49
Quick zero-order approximation

• Normally used in combination with Glens flow law

50
Quick zero-order approximation

• Shown that vertical shear stresses are only
  important ones

51
Quick zero-order approximation

• Integrate from bed to some height within ice to
  get velocity

52
Quick zero-order approximation

• For A assumed constant with depth

53
Simple glacier flow model 1-d. equations
54
Simple glacier flow model 1-d. equations

• Ice thickness (H) evolution
• Zero order model of ice flow (u)
• Mass balance model (b)
• Here simplified by assuming flat bed

55
Simple glacier flow model 1-d. equations

• Mass balance model (b)
• Nonlinear (D) parabolic equation
• Here simplified by assuming flat bed

56
Glacier model spatial discretization

• Regular 1-d. grid
• Staggered so that fluxes evaluated halfway
  between cells

GRID i-1 i-½ i i½ i1
H,s,b ? ? ?
?s/?x ? ?
D(1) ? ? ? ? ?
D(2) ? ?
57
Glacier model spatial discretization

• Substitute FD fluxes

58
Glacier model spatial discretization

• Forms a tridiagonal matrix

59
Glacier model explicit time discretization

• First-order, forward approximation to time
  derivative
• Disadvantages first-order error (Taylors) and
  not very stable (CFD limit)

i-1 i-½ i i½ i1
time ? ? ? ? ?
time dt ?
60
Glacier model fully implicit discretization

• First-order, backward approximation to time
  derivative
• still first-order error
• but now far more stable
• however need to evaluate flows a new time step

i-1 i-½ i i½ i1
time ?
time dt ? ? ? ? ?
61
Glacier model implicit solution

• Use direct Gaussian solver on tridiagonal (NR)
• No additional computing cost compared to explicit
• Stable as long as diagonally dominant (i.e.
  always in this case)

62
Glacier model semi-implicit discretization

• Second-order, centred approximation to time
  derivative (Crank-Nicolson)
• need to evaluate flows at both time steps

i-1 i-½ i i½ i1
time ? ? ? ? ?
time dt ? ? ? ? ?
63
Glacier model non-linear instability

• Diffusivity (D) is a highly non-linear function
  of thickness (power 5) and surface slope (power
  2)
• Good type problem for methods later on
• Options
• Ignore the problem (just use D from previous time
  step)
• Picard iteration
• Under-relaxation, Hindmarsh schemes
• Newton-Raphson iteration

64
Glacier model Picard iteration
Solve for new S given D
Calculate D from H, s
Check for convergence
DONE
65
Glacier model relaxation and Hindmarsh

• Non-linear instability manifests itself as
  high-frequency oscillation as numerical solution
  over corrects past the true solution but never
  finds it
• Under-relaxation reduces the correction applied
  in each iteration and makes convergence more
  likely
• The scheme published by Hindmarsh and Payne
  (1996) recognises the (chaotic) oscillatory
  behaviour and aims for the mid-point

H
iterations
66
Glacier model Newton-Raphson solver

• Recognises the non-linear nature of the problem
  fully
• Very fast convergence
• More complicated to programme but MAPLE etc makes
  evaluation of Jacobian easy

67
Glacier model final comments

• Can insert expressions for basal slip but these
  typically make assumption about all gravitational
  driving balanced locally
• Strictly only appropriate at grid scales 10-20
  times ice thickness
• Easily coupled to simple mass balance climate
  models (also isostasy)
• Very successful in explaining features of Ice
  Ages
• VERIFY numerics by comparing the Nye-Vialov
  steady-state profiles

68
Examples

• This type of model was very useful in
  understanding causes of Ice Ages and relation to
  Milankovitch forcing

snowline
69
Examples

• Coupled via snowline to simple climate models.
• Used to explain disparity between forcing and
  response

70
Examples

• Range of experiments showed that 100 kyr cycles
  could be reproduced if extra physics added to
  cause fast deglaciation

Pollard, 1982
71
Simple ice sheet model 2-d. equations

• 2-d. version of previous non-linear parabolic
  equation
• Very similar derivation
• 2-d. nature of problem increases complexity of
  numerics

72
Simple ice sheet model spatial discretization

• Number of choices for how quantities calculated
  on staggered grid
• Problem originally studied in GCM by Arakawa
• C-grid traditionally used by glaciologists
  although B may have advantages

73
Simple ice sheet model spatial discretization

• Problems with evaluating all derivatives on
  C-grid
• B-grid more symmetrical
• ? flux
• ? thickness

i-1 i-½ i i½ i1
j1 ? ? ?
j½ ?
j ? ? ? ? ?
j-½ ?
j-1 ? ? ?
i-1 i-½ i i½ i1
j1 ? ? ?
j½ ? ?
j ? ? ?
j-½ ? ?
j-1 ? ? ?
74
Simple ice sheet model implicit solution

• No longer have simple tridiagonal matrix
• This means direct (Gaussian) methods of solution
  will be very, very expensive
• Need to look towards iterative methods

75
Simple ice sheet model matrix solvers

• The traditional method (i.e., explicit) uses a
  point iteration
• The convergence of the technique can be improved
  by over-relaxation (SOR)
• Could use line inversion (alternating direct
  implicit, ADI)

76
Simple ice sheet model matrix solvers

• Conjugate gradient methods are ideally suited
• Commonly available as library routines that can
  be easily incorporated
• Numerical recipes or SLAP library (which offers
  many preconditioners)
• All use sparse matrix storage
• Multigrid techniques now popular in many other
  fields but untried in glaciology

77
Examples

• Type of model was work horse of modelling
  community and has produced many of the
  simulations on which previous IPCC reports are
  based (eg Huybrechts and de Wolde 1999)

78
Reasons for interest in temperature

• Determines whether water present at the bed
• ice frozen to bed exerts greater traction
• meltwater allows lubrication and reduced traction
• all known ice streams have water at the bed
• Determines the softness of ice
• warm is deforms more rapidly
• very large effect over natural temp. range

79
Ice temperature evolution equations

• Vertical diffusion
• Vertical advection
• Horizontal advection
• Dissipation
• Local rate of change
• Upper boundary condition air temperature
• Lower boundary condition geothermal heat flux

80
Potential coupling to ice flow?

• Effect of thicker ice?
• Effect of increased driving stresses?
• Effect of faster velocities?
• Effect of warmer temperatures?

81
Ice temperature evolution numerics

• Treat as column model with horizontal terms as
  corrections
• 1-d. diffusion equation can be solved using
  previous methods
• Additional techniques needed
• Boundary conditions
• Advection terms
• Stretched vertical coordinate

82
Temperature evolution upper boundary

• Dirichlet boundary condition
• Simply use air temperature at k 1

83
Temperature evolution lower boundary

• von Neumann or flux boundary condition
• Use boundary condition to supply extra equation
  for n1
• Substitute this into original approximation at n

84
Ice temperature evolution advection

• Number of options for first derivatives
  associated with advection
• Centred, second order
• Non-centred, first order
• Non-centred, second order
• Centred derivative is unconditional unstable
  (solution splits)
• Non-centred, first order introduces excessive
  artificial diffusion

85
Ice temperature evolution advection

• Non-centred, second order normally produces
  satisfactory results if used as part of an
  upwinding scheme

86
Ice temperature evolution stretched grid

• Ice masses have highly irregular vertical
  profiles
• This causes problems if a regular grid is used in
  the vertical
• Points will jump in and out of the domain as the
  ice mass thins and thickens
• Have implicitly been ignoring this problem in
  spatial evolution of ice mass

87
Ice temperature evolution stretched grid

• Grid stretching is a technique that ensures the
  grid always fits to the local ice thickness
• Also ensures that boundaries always lie on grid
  points (hence boundary conditions easier to
  employ)
• However does add the complexity of equations

88
Ice temperature evolution stretched grid

• Need to add terms to deal with deformed grid
• This is easy for terms in z

89
Ice temperature evolution stretched grid

• But trickier for other terms which in t, x, and y
• If this is done with each term then an apparent
  velocity for the grid can be found

90
Ice temperature evolution stretched grid

• An irregularly spaced grid is often used in the
  vertical, stretched coordinate
• Points are concentrated near the bed (most
  deformation)
• This does increase truncation error (non centred)
• Better alternative is to use non-linear coordinate

91
Ice temperature evolution tidying up

• Horizontal velocities still come from zero order
  shallow ice model
• Vertical velocities found using incompressibility
  condition
• Numerical integration from bed to surface
  (Trapezoidal Rule, NR)
• Test of integration accuracy from surface
  kinematic boundary condition

92
Ice temperature evolution uses

• Work coupling melt rates to simple local water
  storage model

93
Ice temperature evolution uses

• Interaction between flow, temperature and basal
  water can lead to interesting results

94
Ice streams and thermomechanics

• Coupled evolution of ice sheet flow, form and
  temperature
• Activation front of steep surface slopes triggers
  warming and sliding
• Eventual cooling through cold ice advection and
  loss of dissipation

95
Heinrich events

• Marshall and Clark (1997) simulation of Heinrich
  events with 3-d. thermomechanical model

96
EISMINT thermomechanical results

• Comparison of 10 ice sheet models with
  thermomechanical coupling (Payne and others 2000)
• Results show sensitivity details of numerics
• Hindmarsh (2004) uses normal modes to show that
  while feedback exists it is relatively weak
  (hence numerical problems)

97
Ice temperature evolution uses

• Also for Antarctica

log10(U)
98
Ice shelf models derivation

• Still assume that vertical balance is dominated
  by normal stress gradient
• Show remainder of derivation for 1-d. confined
  shelf

99
Ice shelf models derivation

• Repeat zero-order proceedure
• Left with additional stress deviator term

100
Ice shelf models derivation

• Integrate through ice thickness
• Apply Leibniz rule

101
Ice shelf models derivation

• Subsititute boundary conditions for stresses

102
Ice shelf models derivation

• Perform remaining integrations
• Terms cancel etc.

103
Ice shelf models derivation

• Two 2-d. vertically integrated equations
• also used to study ice streams with inclusion of
  a basal drag term
• Linear slip law normally used
• Up to this point the derivation is general

104
Ice shelf models derivation

• Normally solved by replacing stress deviators
  with strain rates
• Inverse form of Glens flow law
• Define an effective viscosity (f)

105
Ice shelf models derivation

• Problem is that need vertically averaged
  quantities for stresses and velocities etc.
• No simple relation between these averages if
  vertical shear is present
• Must assume (e.g., MacAyeal) that no variation in
  the vertical
• This limits application to true ice shelves or
  shelfy ice streams

106
Ice shelf models derivation

• Effective viscosity can also be found in terms of
  strain rates

107
Ice shelf models derivation

• Final equations are coupled non-linear elliptical
  equations in u and v
• Non-linearity enters via effective viscosity

108
Ice shelf models solution

• Solve for u given v, and vice versa
• Use previous techniques for 2-d. parabolic
  equations (conjugate gradients)
• Deal with non-linearity using previous techniques
  also

109
Ice shelf models boundary conditions

• Kinematic specify a velocity either from
• observations (if modelling shelf in isolation) or
• zero-order model (if coupling to an ice sheet
  model)
• Dynamic specify a stress
• Appropriate to front of ice shelf
• Balance of forces with displaced water

H
h
110
Ice shelf models boundary conditions

• Dynamic boundary condition messy
• Greatly simplified if implemented so that shelf
  front aligned along x or y
• Problem when modelling irregular shaped shelves
• Use an artificial shelf with arbitrarily low
  thickness to extend shelf to domain edge

111
Ice shelf models thickness evolution

• Surface elevation can be found from buoyancy
• Ice thickness evolution can no longer be coupled
  with velocity calculation
• Must be solved separately

112
Ice shelf models thickness evolution

• Number of alternatives
• Convert velocities to diffusivities and use old
  method
• Solve using simple schemes for hyperbolic
  equations (e.g. staggered leapfrog etc)
• Use more complex transport scheme (e.g.
  semi-lagrangian methods)
• A satisfactory general solver has yet to be found

113
Ice shelf model algorithm
Calculate new gravitational driving stress and
boundary conditions
Calculate for f given u, v
Solve for u given f, v
Solve for v given f, u (old)
Calculate new thickness distribution
Check for convergence
114
Antarctic mass balance
Davis and others (2005) use ERS alimetry to
determine change in surface elevation over last
decade
115
Antarctic mass balance
Comparison with ERA-40 climate reanalysis for
1980 to 2001 for precipitation
116
Interest in Pine Island Glacier

• PIG has the largest discharge (66 Gt yr-1) of all
  WAIS ice streams
• with Thwaites Glacier, it drains 40 of the WAIS
• little studied in comparison to Siple Coast ice
  streams

117
Interest in Pine Island Glacier

• grounding line retreated 8 km between 1992 and
  1994
• implies ice thinning at the grounding line of the
  order of 3.5 m yr-1
• radar altimetry shows widespread thinning
• thinning pattern extends 150 km from grounding
  line
• thinning maps on to template of fast flowing
  section of ie stream

118
Interest in Pine Island Glacier

• thinning unlikely to be related to snowfall
  variation
• hypothesized causes of thinning
• internal flow mechanics of ice stream (surging?)
• long-term response to climate change (LGM?)
• recent collapse of ice shelf and/or change in
  grounding

119
Outline of model

• Vertically-integrated MacAyeal model includes
• momentum balances in x and y
• assumes vertical shear minimal
• viscous flow law
• dynamic b.c. at shelf front
• Prognostic thickness evolution but based on
  perturbations to ice flow

• Same domain and grid
• Ice surface now two dimensional
• Time steps 0.01 yr for H and 0.1 yr for u/v
  (fast wave speeds)

120
Used to study Pine Island Glacier
Shallow ice model
Kinematic b.c. (velocity)
Dynamic b.c. (stress)
Thickness evolution throughout
MacAyeal stream
MacAyeal shelf
121
Results surface lowering in m after 0 yr
122
Results surface lowering in m after 2 yr
123
Results surface lowering in m after 4 yr
124
Results surface lowering in m after 7 yr
125
Results surface lowering in m after 10 yr
126
Results surface lowering in m after 15 yr
127
Results surface lowering in m after 20 yr
128
Results surface lowering in m after 35 yr
129
Results surface lowering in m after 50 yr
130
Results accumulated thinning after 150 yr.
Total thinning produced by a range of ?2 changes
near the grounding line
131
General ice-flow models derivation

• Both the ice sheet and shelf models are limited
  in their applicability by the assumptions that
  are made in their derivation
• Ice sheet models assume local stress balance, and
  are only strictly applicable at 10-20 time ice
  thickness
• Ice shelf models assume no vertical shear
• Many problems in glaciology lie between these two
  extremes, e.g.
• Many (most) ice streams
• Onset areas and shear margins
• Valley glaciers
• Ice divides (coring locations)

132
General ice-flow models derivation

• This motivates development of general models
• Approach is similar to ice shelf model but
  vertical shear terms are not discarded
• Stretched coordinates are again used
• But second derivatives in x introduce much
  complexity

133
General model geometry and boundary conditions
zero traction
zero velocity
force balance
viscous slip law ? ?2u
134
Dangers of numerical modelling

• Tendency to treat model as a black box or a
  surrogate for reality.
• Cannot simply use the model says without
  proving that the model solves the underlying
  equations satisfactorily and that processes being
  described are realistic consequences of the
  equations.