# Numerical modelling of ice sheets and ice shelves

Document Sample

```					    Numerical
modelling

Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
numerics
Numerical modelling
analogy
ﬁnite differences
exact solutions
of ice sheets and ice shelves
shallow ice
sheets
solving the SIA
is it right?
Ed Bueler
application
the most basic
shallow assumption          Dept of Mathematics and Statistics
ﬂowline Stokes
and Geophysical Institute
derivation of SIA
shelves and
streams
shallow shelf
approximation (SSA)            September 2010
ﬂow-line ice shelf
numerical solution         Karthaus Summer School
scaling up

next steps
practicalities
omitted models
Numerical
modelling                                                                                                    notation
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
ﬁnite differences
exact solutions

shallow ice                                               ﬁgure modiﬁed from [Schoof, 2007]
sheets
solving the SIA
is it right?
application
•   coordinates t, x, y , z                                          (z is vertical, positive upward)
the most basic           •   subscripts for partial derivatives: ux = ∂u/∂x
shallow assumption
ﬂowline Stokes           •   h = ice surface elevation, b = bedrock elevation
derivation of SIA        •   H = ice thickness
shelves and
streams
•   T = ice temperature
shallow shelf            •   (u, v , w) = ice velocity
approximation (SSA)
ﬂow-line ice shelf
•   ρ = density of ice, ρw = density of ocean water
numerical solution       •   g = acceleration of gravity
scaling up
•   Dij = strain rate tensor
next steps
practicalities
•   stress tensors: σij = Cauchy and τij = deviatoric
omitted models

notation here is generally consistent with [Bueler and others, 2005], [Bueler and Brown, 2009],
[Fowler, 1997], [Greve and Blatter, 2009], [Schoof, 2006; 2007]
Numerical
modelling                                                                    scope
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)
slogans:
heat analogy &          • my focus is on approximating ice ﬂow
numerics
analogy
ﬁnite differences
• my example numerical codes actually work
exact solutions

shallow ice           continuum models:
sheets
solving the SIA         • shallow ice approximation (SIA) in 2D
is it right?
application
the most basic
• shallow shelf approximation (SSA) in ﬂowline (1D)
shallow assumption
ﬂowline Stokes          • mass continuity & surface/base kinematical equations
derivation of SIA

shelves and
• . . . and a little bit more
streams
shallow shelf
approximation (SSA)   numerical ideas:
ﬂow-line ice shelf
numerical solution      • ﬁnite difference schemes
scaling up

next steps              • solving algebraic systems arising in stress balances
practicalities
omitted models          • veriﬁcation
Numerical
modelling                                                          outside of scope
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
not covered here:
numerics
analogy                 • thermomechanical coupling
ﬁnite differences
exact solutions         • subglacial material/process modeling
shallow ice
sheets                  • snow/ﬁrn process modeling
solving the SIA
is it right?            • earth deformation
application
the most basic
shallow assumption
• numerical solution of Stokes equations
ﬂowline Stokes
derivation of SIA
• grounding line and calving front issues
shelves and
streams
• anisotropy in ice ﬂow
shallow shelf
approximation (SSA)     • paleoclimate modeling and “spin-up”
ﬂow-line ice shelf
numerical solution      • ﬁnite element, spectral, wavelet, multigrid, . . . methods
scaling up

next steps
• . . . etc. etc.
practicalities
omitted models
Numerical
modelling                                                        main equations
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
ﬁnite differences
• shallow ice approximation (SIA), equation (1)
exact solutions

shallow ice           • heat equation (2)   (for analogy purposes)
sheets
solving the SIA       • surface kinematical equation (3)
is it right?
application
the most basic
• basal kinematical equation (4)
shallow assumption
ﬂowline Stokes        • map-plane mass continuity equation (5)
derivation of SIA

shelves and
• Stokes equations in a ﬂow line, equations (6)
streams
shallow shelf         • shallow shelf approximation (SSA), equation (7)
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
Numerical
modelling                                                     M ATLAB/O CTAVE codes
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
• these lectures are structured around
ﬁnite differences
exact solutions
• 7 codes which solve PDEs
shallow ice
• an additional 13 codes which plot, pre-process, post-process
sheets
solving the SIA       • each PDE code is 1/2 to 1 page in length
is it right?
application           • all tested in M ATLAB   (v. 7.9) and   O CTAVE   (v. 3.2)
the most basic
shallow assumption
ﬂowline Stokes
• available as
derivation of SIA
◦ .zip or .tar.gz, from my USB memory stick, or
shelves and
streams
◦ online:
shallow shelf
approximation (SSA)
http://www.dms.uaf.edu/~bueler/karthaus/
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
Numerical
modelling                                                    Outline
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
numerics              1 ice ﬂow: an outsider’s superﬁcial view
analogy
ﬁnite differences
exact solutions

shallow ice
2 heat analogy & numerics
sheets
solving the SIA
is it right?
application
3 shallow ice sheets
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA     4 shelves and streams
shelves and
streams
shallow shelf
approximation (SSA)
5 next steps
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
Numerical
modelling                                                               my ﬁrst goal
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
ﬁnite differences     • I want to get to an equation for which I can say:
exact solutions

shallow ice
sheets
solving the SIA            numerically solve this equation, and you’ve got a usable
is it right?
application
model for . . . (a speciﬁc ice ﬂow problem)
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA
• to get to that goal: I will very quickly recall the continuum
shelves and
streams                 mechanics of ice ﬂow
shallow shelf
approximation (SSA)   • . . . from a naive outside view
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
Numerical
modelling                                            let’s treat ice in glaciers as a ﬂuid
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &        • we describe ﬂuids primarily by a velocity ﬁeld u(t, x, y , z)
numerics
analogy
ﬁnite differences
• if the ice ﬂuid were much faster-moving than it actually is1 ,
exact solutions              and if it were linearly-viscous like liquid water, then ice ﬂow
shallow ice
sheets
modeling would be more familiar to other the
solving the SIA
is it right?
climate-modeling and/or engineering computational ﬂuids
application                  people
the most basic
shallow assumption
ﬂowline Stokes
• . . . we would all use the Navier-Stokes equations, for
derivation of SIA
incompressible ﬂuids, as our ﬂow model:
shelves and
streams
shallow shelf
approximation (SSA)
·u=0                               incompressibility
ﬂow-line ice shelf                                                  2
numerical solution              ρ (ut + u ·     u) = − p + ν            u + ρg    force balance
scaling up

next steps
practicalities
omitted models

1 if   gravity were a lot stronger and/or ice were a lot weaker, for example
Numerical
modelling                     hmmm . . . does not sound like glaciology to me!
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
is numerical ice ﬂow modeling actually computational ﬂuid
approximation (SIA)
dynamics?
heat analogy &
numerics
analogy
• yes
ﬁnite differences
exact solutions         • at geophysical scale like atmosphere/ocean
shallow ice
sheets                  • . . . but it is a weird one
solving the SIA
is it right?
• consider what makes atmosphere/ocean ﬂow modeling
application
the most basic
exciting:
shallow assumption
ﬂowline Stokes
◦   turbulence
derivation of SIA            ◦   convection
shelves and                  ◦   coriolis force
streams
shallow shelf                ◦   density/salinity stratiﬁcation
approximation (SSA)
ﬂow-line ice shelf           ◦   lots of tracers & chemistry (CO2 , methane, ozone, . . . )
numerical solution
scaling up              • none of the above list is relevant to ice ﬂow!
next steps
practicalities
• so what could be interesting about the ﬂow of slow, cold,
omitted models
shallow, non-turbulent, chemically-dull, old ice?
Numerical
modelling                                           ice is a slow, shear-thinning ﬂuid
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)
• our ﬂuid is
heat analogy &
numerics
analogy                  slow:                        ρ (ut + u · u) ≈ 0
ﬁnite differences
exact solutions          non-Newtonian: viscosity ν is not constant in “τij = 2νDij ”
shallow ice
sheets
• and so the standard “full” model for ice ﬂow is shear-thinning
solving the SIA
is it right?
Stokes as follows:
application
the most basic
shallow assumption               ·u=0                                incompressibility
ﬂowline Stokes
derivation of SIA                 0=− p+                · τij + ρg   force balance
shelves and                                 n−1
streams                          Dij = Aτ         τij                ﬂow law
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
numerical solution
• equations above are true at every instant, and so
scaling up
geometry and boundary stress and ice strength
next steps
practicalities                   determine velocity instantaneously in the model
omitted models
Numerical
modelling                                                              ice is a slow ﬂuid 2
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)     • the “slow” fact, ρ (ut + u ·      u) ≈ 0, means no forces of inertia
heat analogy &
numerics
from acceleration
analogy
ﬁnite differences
• thus velocity is not a state variable of any model for evolving
exact solutions
ice sheets and glaciers
shallow ice
sheets                  • . . . it is a “diagnostic” output . . . but a very important output!
solving the SIA
is it right?
application
• an ice sheet code usually recomputes the full velocity ﬁeld at
the most basic
shallow assumption
every time, without storing it from the previous step2
ﬂowline Stokes
derivation of SIA
• to see how glacier ﬂow ﬁelds change in time we must track
shelves and
the changes to
streams
shallow shelf
◦ the location of the ice surface and base
approximation (SSA)
ﬂow-line ice shelf
◦ the force (stress) applied at the boundary
numerical solution              ◦ the ice strength, as it changes with ice temperature (for
scaling up
example)
next steps
practicalities
omitted models

2 Bob Dylan would say: you don’t need a glaciologist to remember which way
the wind blows
Numerical
modelling                                                         plane ﬂow Stokes
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
more concretely, the shear-thinning Stokes model is this:
numerics
analogy                 • in the x, z plane ﬂow case the Stokes equations say
ﬁnite differences
exact solutions

shallow ice
ux + wz = 0                      incompressibility
sheets
solving the SIA                     px = τ11,x + τ13,z         stress balance (x)
is it right?
application                         pz = τ13,x − τ11,z − ρg    stress balance (z)
the most basic
shallow assumption                            n−1
ﬂowline Stokes
ux = Aτ         τ11        ﬂow law (diagonal)
derivation of SIA
n−1
shelves and
uz + wx = 2Aτ          τ13       ﬂow law (off-diagonal)
streams
shallow shelf
approximation (SSA)     • ﬁve equations in ﬁve unknowns (u, w, p, τ11 , τ13 )
ﬂow-line ice shelf
numerical solution      • complicated enough . . . can I make it familiar by looking at a
scaling up

next steps
very-simpliﬁed situation?
practicalities
omitted models
Numerical
modelling                                                                   slab-on-a-slope
z
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &                                                                              z=
numerics                                                                                         H
• rotated coordinates:                                     ice
analogy
ﬁnite differences
ˆ           ˆ
g = g sin θ x − g cos θ z
exact solutions

shallow ice
sheets                                                                 bed                  g
solving the SIA                                                           roc
k
is it right?                                                                           z=
application
0
the most basic
• as on previous slide, but with n = 3:                                       x
shallow assumption
ﬂowline Stokes                                                                          2
derivation of SIA              ux + wz = 0                                       ux = Aτ τ11
shelves and
streams                              px = τ11,x + τ13,z + ρg sin θ      uz + wx = 2Aτ 2 τ13
shallow shelf
approximation (SSA)                  pz = τ13,x − τ11,z − ρg cos θ
ﬂow-line ice shelf
numerical solution
scaling up             • but for a slab-on-a-slope, by deﬁnition, there is no variation in x:
next steps
practicalities                          wz = 0                         0 = τ11
omitted models

τ13,z = −ρg sin θ               uz = 2Aτ 2 τ13
pz = −ρg cos θ
Numerical
modelling                                                                        slab-on-a-slope 2
Ed Bueler

ice ﬂow: a            • add some boundary conditions:
superﬁcial view
shallow ice
approximation (SIA)
w(base) = 0,            p(surface) = 0,               u(base) = u0 ,
heat analogy &
numerics              • further-simpliﬁed equations:
analogy
ﬁnite differences
exact solutions
τ13 = ρg sin θ(H − z)                  p = ρg cos θ(H − z)
shallow ice                      w =0                              τ11 = 0
sheets
solving the SIA                                                    ˆ     z
is it right?          • and
application                            u(z) = u0 + 2A(ρg sin θ)3             (H − z )3 dz
the most basic                                                       0
shallow assumption
ﬂowline Stokes                                         1
derivation of SIA                           = u0 +       A(ρg sin θ)3 H 4 − (H − z)4
2
shelves and
streams
shallow shelf                      H
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up
z

next steps
practicalities
omitted models

0                                   u0
τ13(z)                          u(z)
Numerical
modelling                                                                     slab-on-a-slope 3
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice                                                                       Horizontal Velocity (m/yr)
approximation (SIA)
0         10         20      30
0
heat analogy &
numerics
• do we believe these equations?
analogy
ﬁnite differences     • velocity on last slide (and below)
50
exact solutions
was from a formula

Depth (m)
shallow ice
sheets                • compare to observations at right                 100
solving the SIA
internal
is it right?
application
• seems there is an element of truth                                   deformation
the most basic
shallow assumption                                                       150
ﬂowline Stokes
derivation of SIA                                                              basal
motion
shelves and                                                              200
streams
bedrock
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up
Velocity proﬁle of the Athabasca
u0
next steps                               u(z)                 inclinometry. [Savage and Paterson,
practicalities
1963]
omitted models
Numerical
modelling                                                                  mass continuity
Ed Bueler

ice ﬂow: a
now we know the velocity u(z) . . . so what?
superﬁcial view
shallow ice             • suppose now that our slab has slowly varying thickness H(t, x)
approximation (SIA)

heat analogy &          • compute the vertical average of velocity, which depends on x:
numerics
analogy                                                   ˆ
ﬁnite differences                                       1 H
exact solutions
¯
u (x) =       u(x, z) dz
H 0
shallow ice
sheets
solving the SIA
is it right?           • consider change of area (really: ice
application
the most basic           volume) in the ﬁgure at right:                           M(x)
shallow assumption
ﬂowline Stokes                   ˆ x2
derivation of SIA          dA
=                  ¯      ¯
M(x) dx + u1 H1 − u2 H2
shelves and                 dt
streams
x1                                H1                         H2
shallow shelf                                                                      A
approximation (SSA)
ﬂow-line ice shelf
• in limit where width dx = x2 − x2 is
numerical solution       small we have A ≈ dx H                      ¯
u1                   ¯
u2
scaling up

next steps             • divide by dx and get
practicalities
omitted models
¯
Ht = M − (u H)x                      x1                   x2
• this is the mass continuity equation
Numerical
modelling         rough explanation of “shallow ice approximation” (SIA)
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice           • I’ll consider only u0 = 0 case in these lectures (“non-sliding SIA”)
approximation (SIA)

heat analogy &        • from slab-on-slope velocity formula,
numerics
analogy                                  ˆ H
ﬁnite differences                             1
exact solutions                    ¯
uH =         A(ρg sin θ)3 H 4 − (H − z)4             dz
0  2
shallow ice
sheets                                    1               4 5
solving the SIA                        = A(ρg sin θ)3       H
is it right?                              2               5
application
2
= A(ρg sin θ)3 H 5
the most basic
shallow assumption
ﬂowline Stokes                            5
derivation of SIA

shelves and           • note sin θ ≈ tan θ = −hx
streams
shallow shelf
approximation (SSA)
• combine with mass continuity Ht = M − (u H)x :
¯
ﬂow-line ice shelf
numerical solution
2
scaling up                                  Ht = M +         A(ρg)3 H 5 |hx |2 hx              (0)
next steps
5                        x
practicalities
omitted models
• it is pretty rough . . . but we will get back to it
Numerical
modelling                                              slow, non-Newtonian, shallow
Ed Bueler

ice ﬂow: a            • ice sheet ﬂow problems have three outstanding properties:
superﬁcial view
shallow ice               1 slow
approximation (SIA)
2 non-Newtonian
heat analogy &
numerics                  3 shallow
analogy
ﬁnite differences     • all three properties will be in all of my models
exact solutions

shallow ice           • regarding “shallow,” here in red is a to-scale cross section of
sheets
solving the SIA
Greenland at 71◦ :
is it right?                        4000
application
the most basic
shallow assumption
3000
ﬂowline Stokes
derivation of SIA

shelves and
streams                             2000
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
1000
numerical solution
scaling up

next steps
practicalities
0
omitted models

-1000
-600000   -400000   -200000       0   200000   400000
x
Numerical
modelling                    the isothermal shallow ice approximation (SIA)
Ed Bueler

ice ﬂow: a
superﬁcial view       is a model which applies reasonably well to
shallow ice
approximation (SIA)
• grounded ice sheets
heat analogy &
numerics                • on relatively uninteresting bed topography, whose ﬂow is
analogy
ﬁnite differences       • not dominated by sliding nor by liquid water at base or margin
exact solutions

shallow ice
sheets
solving the SIA
is it right?
application
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA

shelves and
streams
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
“Polaris Glacier,” northwest Greenland, photo 122, [Post and LaChapelle, 2000]
omitted models

• might apply to the above
Numerical
modelling                             isothermal shallow ice approximation 2
Ed Bueler

ice ﬂow: a
superﬁcial view
• the non-sliding isothermal SIA:
shallow ice
approximation (SIA)

heat analogy &
numerics
Ht = M +      · ΓH n+2 | h|n−1 h       (1)
analogy
ﬁnite differences
exact solutions

shallow ice
◦   generalizes slab-on-slope equation (0) earlier
sheets                    ◦   h =H +b
solving the SIA
is it right?              ◦   accumulation if M > 0, ablation if M < 0
application
the most basic
◦   recall Glen ﬂow law: Dij = A(T )τ n−1 τij
shallow assumption
ﬂowline Stokes
◦   Γ is a positive constant
derivation of SIA
• numerically solve equation (1), and you
shelves and
streams                 have a usable model for ice ﬂow in the
shallow shelf
approximation (SSA)     Barnes Ice Cap [Mahaffy, 1976]
ﬂow-line ice shelf
numerical solution
scaling up

next steps            • questions:
practicalities
omitted models            ◦ where does equation (1) come from?
◦ how to solve it numerically?
◦ how to think about it?
Numerical
modelling                                                    Outline
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
numerics              1 ice ﬂow: an outsider’s superﬁcial view
analogy
ﬁnite differences
exact solutions

shallow ice
2 heat analogy & numerics
sheets
solving the SIA
is it right?
application
3 shallow ice sheets
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA     4 shelves and streams
shelves and
streams
shallow shelf
approximation (SSA)
5 next steps
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
Numerical
modelling                                              compare to heat equation
Ed Bueler

ice ﬂow: a
superﬁcial view       • recall Newton’s law of cooling
shallow ice
approximation (SIA)
dT
heat analogy &                     = −µ(T − Tambient )
numerics                        dt
analogy
ﬁnite differences
exact solutions
where T is temperature and µ relates
shallow ice
to material and geometry of object
sheets
solving the SIA       • Newton’s law for segments of an
is it right?
application
insulated rod:
the most basic
shallow assumption         dTj           1
ﬂowline Stokes                 = −˜ Tj − (Tj−1 + Tj+1 )
µ
derivation of SIA          dt            2
shelves and                      ˜
µ
streams                        = (Tj−1 − 2Tj + Tj+1 )
shallow shelf                    2
approximation (SSA)
ﬂow-line ice shelf
numerical solution
(where µ is material constant
˜
scaling up              proportional to ∆x −2 )
next steps
practicalities        • this suggests ﬁnite difference
omitted models
approximation of a PDE:

Tt = DTxx
Numerical
modelling                                                compare to heat equation 2
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)
continuum version:
heat analogy &
numerics               • Fourier rewrote Newton’s law as heat
analogy
ﬁnite differences        ﬂux in continua: q = −k u
exact solutions

shallow ice
• by conservation of energy, allowing an
sheets                   additional source of heat f , get heat
solving the SIA
is it right?             equation:
application
the most basic
shallow assumption                ρcut = f +   · (k u)
ﬂowline Stokes
derivation of SIA
• set D = k /(ρc) and M = f /(ρc):
shelves and
streams
shallow shelf
approximation (SSA)
ut = M +      · (D   u)                   (2)
ﬂow-line ice shelf
numerical solution
scaling up              • this heat equation is a diffusive, time-dependent PDE, like the SIA
next steps                model for ice sheets
practicalities
omitted models
Numerical
modelling                                        analogy: SIA versus heat equation
Ed Bueler

ice ﬂow: a            • side-by-side comparison:
superﬁcial view
shallow ice
approximation (SIA)
SIA: H(t, x, y ) is ice thickness        heat: u(t, x, y ) is temperature
heat analogy &
numerics                                           n+2          n−1
analogy
Ht = M +         · ΓH          | h|         h   ut = M +           · (D     u)
ﬁnite differences
exact solutions

shallow ice           • we identify the diffusivity in the SIA:
sheets
solving the SIA
is it right?
application
D = ΓH n+2 | h|n−1
the most basic
shallow assumption
ﬂowline Stokes
• non-sliding shallow ice ﬂow diffuses the ice because the ﬂow
derivation of SIA
shelves and
streams
shallow shelf
approximation (SSA)   • some issues with this analogy:
ﬂow-line ice shelf
numerical solution        ◦ H = h when bed is not ﬂat . . . so what?
scaling up
◦ D depends on solution H(t, x, y ) . . . how does that complicate
next steps
practicalities
the numerical solution method?
omitted models            ◦ D → 0 at margin (H → 0) and at dome (| h| → 0) . . . so what?
• I’ll get back to these “issues”, but now let’s get our hands
dirty and numerically solve the heat equation
Numerical
modelling                                    ﬁnite differences for heat equation
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)
basic ideas of ﬁnite differences:
heat analogy &
numerics                • for differentiable f (x) and any ∆x, Taylor says
analogy
ﬁnite differences
1            1
f (x)∆x 2 + f (x)∆x 3 + . . .
exact solutions
f (x + ∆x) = f (x) + f (x)∆x +
shallow ice                                                  2            3!
sheets
solving the SIA
is it right?
• you can replace “∆x” with other expressions, e.g.:
application
the most basic                                              1            1
shallow assumption
f (x − ∆x) = f (x) − f (x)∆x +   f (x)∆x 2 − f (x)∆x 3 + . . .
ﬂowline Stokes                                              2            3!
derivation of SIA
4
shelves and                 f (x + 2∆x) = f (x) + 2f (x)∆x + 2f (x)∆x 2 + f (x)∆x 3 + . . .
streams                                                                   3
shallow shelf
approximation (SSA)
ﬂow-line ice shelf      • basic ﬁnite difference idea for differential equations: combine
numerical solution
scaling up
expressions like these to approximate derivatives
next steps              • for all of these notes, grid points have equal spacing ∆x
practicalities
omitted models
Numerical
modelling                                   ﬁnite differences for heat equation 2
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice           • partial derivative expressions, for example with u = u(t, x):
approximation (SIA)

heat analogy &                            u(t + ∆t, x) − u(t, x)
numerics                    ut (t, x) =                          + O(∆t),
analogy                                            ∆t
ﬁnite differences
u(t + ∆t, x) − u(t − ∆t, x)
exact solutions
ut (t, x) =                               + O(∆t 2 ),
shallow ice                                          2∆t
sheets
u(t, x + ∆x) − u(t, x)
solving the SIA
ux (t, x) =                          + O(∆x),
is it right?
∆x
application
u(t, x + ∆x) − u(t, x − ∆x)
+ O(∆x 2 ),
the most basic
shallow assumption          ux (t, x) =
ﬂowline Stokes                                        2∆x
derivation of SIA
u(t, x + ∆x) − 2u(t, x) + u(t, x − ∆x)
shelves and                uxx (t, x) =                                           + O(∆x 2 )
streams                                                    ∆x 2
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
• sometimes we want a derivative in-between grid points:
numerical solution
scaling up                                             u(t, x + ∆x) − u(t, x)
ux (t, x + (∆x/2)) =                           + O(∆x 2 )
next steps
practicalities
∆x
omitted models
• “+O(∆x 2 )” is better than “+O(∆x)” if ∆x is a small number
Numerical
modelling                                     explicit scheme for heat equation
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)    • recall 1D heat equation ut = Duxx
heat analogy &
numerics               • the explicit scheme using notation ujn ≈ u(tn , xj ), so
analogy
ﬁnite differences
exact solutions
ujn+1 − ujn        uj+1 − 2ujn + uj−1
n             n
shallow ice                                           =D
sheets                                      ∆t                   ∆x 2
solving the SIA
is it right?
application            • let ν = D∆t/(∆x)2 , so
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA
ujn+1 = νuj+1 + (1 − 2ν)ujn + νuj−1
n                     n

shelves and
streams
shallow shelf
approximation (SSA)   • scheme has stencil at right   →
ﬂow-line ice shelf                                                       n+1
numerical solution    • advantage over implicit (later): ujn+1 is
scaling up

next steps
determined by known quantities at time tn
n
practicalities
omitted models                                                                 j-1   j   j+1
Numerical
modelling                                      explicit scheme for heat equation 2
Ed Bueler

ice ﬂow: a
superﬁcial view                            n
• in 2D we write ujk ≈ u(tn , xj , yk )
shallow ice
approximation (SIA)

heat analogy &
• the 2D explicit scheme for the M = 0 heat equation
numerics                   ut = D(uxx + uyy ) is
analogy
ﬁnite differences
n+1   n            n         n     n        n          n     n
exact solutions              ujk − ujk          uj+1,k − 2ujk + uj−1,k   uj,k +1 − 2ujk + uj,k −1
shallow ice                              =D                            +
sheets                           ∆t                      ∆x 2                      ∆y 2
solving the SIA
is it right?
application
• or, with ν x := D∆t/(∆x)2 and ν y := D∆t/(∆y )2 ,
the most basic
shallow assumption
ﬂowline Stokes
n+1                     n         n        n            n         n
ujk = (1 − 2ν x − 2ν y )ujk + ν x uj+1,k + uj−1,k + ν y uj,k +1 + uj,k −1
derivation of SIA

shelves and
streams
shallow shelf                                                            k+1
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up                                                                k
next steps
practicalities
omitted models                           n+1
note: new value ujk is average (is it?)            k-1

of ﬁve quantities at old time tn  −→                     j-1      j           j+1
Numerical
modelling                                                                            implementation
Ed Bueler

ice ﬂow: a            function U = heat(D,J,K,dt,N)
superﬁcial view
shallow ice
approximation (SIA)   dx = 2 / J;    dy = 2 / K;
heat analogy &        [x,y] = meshgrid(-1:dx:1, -1:dy:1);
numerics              U = exp(-30*(x.*x + y.*y));
analogy
ﬁnite differences     nu_x = dt * D    / (dx*dx);        nu_y = dt * D / (dy*dy);
exact solutions
for n=1:N
shallow ice               U(2:J,2:K)   = U(2:J,2:K) + ...
sheets                        nu_x *   ( U(3:J+1,2:K) - 2 * U(2:J,2:K) + U(1:J-1,2:K) ) + ...
solving the SIA
is it right?
nu_y *   ( U(2:J,3:K+1) - 2 * U(2:J,2:K) + U(2:J,1:K-1) );
application
end
the most basic
shallow assumption    surf(x,y,U),     shading(’interp’),        xlabel x,       ylabel y,      zlabel u
ﬂowline Stokes
derivation of SIA
http://www.dms.uaf.edu/~bueler/karthaus/mfiles/heat.m
shelves and
streams
shallow shelf
approximation (SSA)     • solves ut = D(uxx + uyy ) on square −1 < x < 1, −1 < y < 1
ﬂow-line ice shelf
numerical solution      • choice: gaussian initial condition
scaling up

next steps              • “colon notation” removes loops over spatial variables
practicalities
omitted models          • to approximate u on 30 × 30 spatial grid, with D = 1 and N = 20
steps of length ∆t = 0.001,
» heat(1.0,30,30,0.001,20)
Numerical
modelling                                                                                     the look of success
Ed Bueler

ice ﬂow: a
superﬁcial view           • solving ut = D(uxx + uyy ) on 30 × 30 grid
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
approximate solution u(t, x, y ) at
ﬁnite differences
initial condition u(0, x, y )
exact solutions
t = 0.02 with ∆t = 0.001
shallow ice
sheets
solving the SIA
is it right?
application
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA          1
1
shelves and               0.8
streams                                                                               0.8
0.6
shallow shelf
approximation (SSA)                                                                   0.6
u   0.4
ﬂow-line ice shelf
u   0.4
numerical solution        0.2
scaling up                                                                            0.2
01
next steps                      0.5                                           1        01
practicalities                        0                                 0.5                 0.5                                           1
y                               0                                                                   0.5
omitted models                            -0.5                                                    0
-0.5       x
y                               0
-1 -1                                                -0.5
-0.5       x
-1 -1
Numerical
modelling                                                                                               the look of instability
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice                  • both ﬁgures are from solving ut = D(uxx + uyy ) on the same
approximation (SIA)

heat analogy &
50 × 50 grid, at same ﬁnal time and with same D, but with
numerics
analogy
slightly different time steps
ﬁnite differences
exact solutions

shallow ice
sheets
solving the SIA
is it right?
0.25                                                             0.25

application
the most basic             0.2                                                              0.2

shallow assumption
ﬂowline Stokes            0.15                                                             0.15

derivation of SIA
u    0.1                                                         u    0.1

shelves and
streams                   0.05                                                             0.05

shallow shelf
approximation (SSA)         01                                                               01

ﬂow-line ice shelf                   0.5                                           1                  0.5                                            1

0.5                                                               0.5
numerical solution                         0                                                                 0

y                                   0                            y                                    0

scaling up                                     -0.5
-0.5       x
-0.5
-0.5       x
-1 -1                                                             -1 -1

next steps
practicalities
omitted models                             D∆t                                                               D∆t
= 0.402                                                           = 0.625
∆x 2                                                              ∆x 2
Numerical
modelling                                                                          why unstable?
Ed Bueler

ice ﬂow: a
superﬁcial view
• the 1D ﬁrst-order explicit scheme in the form
shallow ice

ujn+1 = νuj+1 + (1 − 2ν)ujn + νuj−1
n                     n
approximation (SIA)

heat analogy &
numerics
analogy
ﬁnite differences          gives new value ujn+1 as an average of old values,
exact solutions

shallow ice             • but it is only an average if the middle coefﬁcient is positive:3
sheets
solving the SIA
is it right?                                                            D∆t
application                                         1 − 2ν = 1 − 2           ≥0
the most basic
shallow assumption
∆x 2
ﬂowline Stokes
derivation of SIA       • true averaging is always stable because averaged wiggles
shelves and
streams
are smaller than the previous wiggles
shallow shelf
approximation (SSA)
• “positive coefﬁcients” is a sufﬁcient stability criterion
ﬂow-line ice shelf
numerical solution      • the same thing as a time-step restriction:
scaling up

next steps
∆x 2
practicalities
∆t ≤
omitted models
2D
3 recommended     basic reference for ﬁnite differences, including stability: [Morton
and Mayers, 2005]
Numerical
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)   dx = 2 / J;    dy = 2 / K;
heat analogy &        [x,y] = ndgrid(-1:dx:1, -1:dy:1);
numerics              U = exp(-30*(x.*x + y.*y));
analogy
ﬁnite differences
t = 0.0;     count = 0;
exact solutions
while t < tf
shallow ice               dt0 = 0.25 * min(dx,dy)^2 / D;
sheets
dt = min(dt0, tf - t);
solving the SIA
is it right?
nu_x = dt * D / (dx*dx);    nu_y = dt * D / (dy*dy);
application               U(2:J,2:K) = U(2:J,2:K) + ...
the most basic                nu_x * ( U(3:J+1,2:K) - 2 * U(2:J,2:K) + U(1:J-1,2:K) ) + ...
shallow assumption
nu_y * ( U(2:J,3:K+1) - 2 * U(2:J,2:K) + U(2:J,1:K-1) );
ﬂowline Stokes
derivation of SIA
t = t + dt;    count = count + 1;
end
shelves and
streams
shallow shelf
surf(x,y,U),   shading(’interp’),          xlabel x,      ylabel y,         zlabel u
approximation (SSA)
numerical solution
scaling up

next steps              • same as heat.m except it gets the time step from:
practicalities
omitted models
D∆t          1
≤
(min{∆x, ∆y })2   4
Numerical
modelling                                           alternative instability ﬁx: implicitness
Ed Bueler

ice ﬂow: a
superﬁcial view       • explicit scheme is only “conditionally stable”,
shallow ice
approximation (SIA)
but there are methods which are stable for                    n+1

heat analogy &          any positive time step ∆t
numerics
analogy               • the simplest such is ﬁrst-order implicit →,
ﬁnite differences                                                                     n
exact solutions

shallow ice
ujn+1   −   ujn         n+1
uj+1   −   2ujn+1   +    n+1
uj−1         j-1      j       j+1
=D
sheets
solving the SIA
∆t                             ∆x 2
is it right?
application           • another is Crank-Nicolson →; instead of
O(∆t, ∆x 2 ) like ﬁrst-order explicit and implicit,
the most basic                                                                        n+1
shallow assumption
ﬂowline Stokes
derivation of SIA
Crank-Nicolson is O(∆t 2 , ∆x 2 )
shelves and
streams                                                                               n
shallow shelf
approximation (SSA)                                                                         j-1      j       j+1
ﬂow-line ice shelf
numerical solution
scaling up
• but for implicit and Crank-Nicolson methods you have to solve
next steps               systems of equations to take each step
practicalities
omitted models         •   Donald Knuth has advice for ice sheet modelers:
We should forget about small efﬁciencies, say about 97% of the time:
premature optimization is the root of all evil.
Numerical
modelling                                                  variable diffusivity and time steps
Ed Bueler

ice ﬂow: a
superﬁcial view
• recall analogy                (SIA) ↔ (heat eqn)
shallow ice
approximation (SIA)       • the SIA has a diffusivity D which varies in time and space, so
heat analogy &
numerics
by analogy:
analogy                                              ut = M +             · (D(t, x, y ) )
ﬁnite differences
exact solutions
• the explicit method is conditionally stable with the “same”
shallow ice
sheets                       stability condition if we evaluate diffusivity D(t, x, y ) at
solving the SIA              staggered grid points:
is it right?
application
the most basic                                            Dj+1/2,k (uj+1,k − uj,k ) − Dj−1/2,k (uj,k − uj−1,k )
shallow assumption                  · (D(t, x, y ) u) ≈
ﬂowline Stokes
∆x 2
derivation of SIA                                              Dj,k +1/2 (uj,k +1 − uj,k ) − Dj,k −1/2 (uj,k − uj,k −1 )
+
shelves and                                                                             ∆y 2
streams
shallow shelf
approximation (SSA)                                                               k+1
ﬂow-line ice shelf
numerical solution
scaling up
in stencil at right:
next steps                 diamonds: u                                             k
practicalities
omitted models             triangles: D
k-1

j-1           j           j+1
Numerical
modelling                                              general diffusion equation code
Ed Bueler

ice ﬂow: a            function [U,dtav] = diffusion(Lx,Ly,J,K,Dup,Ddown,Dright,Dleft,U0,tf,b)
superﬁcial view
shallow ice           dx = 2 * Lx / J;    dy = 2 * Ly / K;
approximation (SIA)
[x,y] = ndgrid(-Lx:dx:Lx, -Ly:dy:Ly);
heat analogy &        U = U0;
numerics
if nargin < 11, b = zeros(size(U0)); end
analogy
ﬁnite differences
exact solutions       t = 0.0;     count = 0;
while t < tf
shallow ice
sheets                    Dregular = max(max(Dup,Ddown),max(Dleft,Dright));
solving the SIA           maxD = max(max(Dregular));
is it right?              dt0 = 0.25 * min(dx,dy)^2 / maxD;
application               dt = min(dt0, tf - t);
the most basic
shallow assumption
mu_x = dt / (dx*dx);    mu_y = dt / (dy*dy);
ﬂowline Stokes            Ushift = U + b;
derivation of SIA         U(2:J,2:K) = U(2:J,2:K) + ...
shelves and
mu_y * Dup    .* ( Ushift(2:J,3:K+1) - Ushift(2:J,2:K)              ) - ...
streams                       mu_y * Ddown .* ( Ushift(2:J,2:K) - Ushift(2:J,1:K-1)               ) + ...
shallow shelf                 mu_x * Dright .* ( Ushift(3:J+1,2:K) - Ushift(2:J,2:K)              ) - ...
approximation (SSA)
mu_x * Dleft .* ( Ushift(2:J,2:K) - Ushift(1:J-1,2:K)               );
ﬂow-line ice shelf
numerical solution
t = t + dt;    count = count + 1;
scaling up            end
dtav = tf / count;
next steps
practicalities
http://www.dms.uaf.edu/~bueler/karthaus/mfiles/diffusion.m
omitted models

• solves abstract diffusion equation ut =             · (D(x, y )    u)
• user must supply diffusivity on staggered grid
Numerical
modelling                                                           on exact solutions
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)   • I am not quite done with the (heat) ↔ (SIA) analogy
heat analogy &
numerics              • . . . I also want to use it to address exact solutions and
analogy
ﬁnite differences       veriﬁcation
exact solutions

shallow ice
sheets                • the two senses of “solution”:
solving the SIA
is it right?              ◦ a solution u(t, x, y ) of the heat equation is a function of time
application
the most basic
and space for which ut = M + · (D u) is true
shallow assumption
ﬂowline Stokes            ◦ a solution of the heat equation is a prediction of that model
derivation of SIA
• exact solutions are exact predictions of the model, but not
shelves and
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
• if we can only approximately ﬁnd solutions of a model then
numerical solution
scaling up
knowing a few exact solutions can help test and maintain the
next steps              quality of the actual code that does the approximation . . . this
practicalities
omitted models
is veriﬁcation [Wesseling, 2001]
Numerical
modelling                                          exact solutions of heat equation
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)    • many solutions to the heat equation are known, but one is
heat analogy &
numerics                 “fundamental” to the time-dependent equation, namely the
analogy
ﬁnite differences
Green’s function
exact solutions

shallow ice
sheets
solving the SIA
is it right?
application
the most basic
• as time goes it changes shape
shallow assumption
ﬂowline Stokes
by shrinking the output
derivation of SIA       (vertical) axis and
shelves and             simultaneously lengthening the
streams
shallow shelf           input (horizontal) axis
approximation (SSA)
ﬂow-line ice shelf    • . . . but otherwise it is the same
numerical solution
scaling up              shape, so it is self-similar
next steps            • there are solutions of the SIA
increasing time →
practicalities
omitted models          just like this
Numerical
modelling                ﬁnding the Green’s function of the heat equation
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice           • the best-known way to ﬁnd the Green’s function of the heat
approximation (SIA)

heat analogy &
equation is by Fourier transform, but we use a method which
numerics
analogy
generalizes to the SIA
ﬁnite differences
exact solutions
• facts used to ﬁnd it:
shallow ice               ◦ the Green’s function starts at time t = 0 with a delta function at
sheets
solving the SIA
the origin r = 0
is it right?              ◦ the Green’s function is angularly-symmetric: u is a function of
application
the most basic              the polar coordinate r = x 2 + y 2
shallow assumption
ﬂowline Stokes
◦ calculus says: if f = f (r ) then 2 f = r −1 (rf )
derivation of SIA         ◦ thus the heat equation in 2D, without additional heat sources,
shelves and
streams
with constant diffusivity D > 0, and for a function of r , is:
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
ut = Dr −1 (rur )r
numerical solution
scaling up
• we use above facts to get Green’s function by its “similarity”
next steps
practicalities
properties
omitted models
• but we do not use the linearity of the heat equation
Numerical
modelling                      ﬁnding a “similarity” solution of heat equation
Ed Bueler

ice ﬂow: a
superﬁcial view       • function u(t, r ) is replaced by function φ of one new variable s
shallow ice
approximation (SIA)
• input scaling: s = t −β r
heat analogy &
numerics              • output scaling: u = t −α φ(s)
analogy
ﬁnite differences     • chain rule says:
exact solutions

ut = −αt −α−1 φ − βt −α−β−1 r φ ,     ur = t −α−β φ ,
shallow ice
sheets                                                                                 etc.
solving the SIA
is it right?
application           • so heat equation ut = Dr −1 (rur )r is replaced by:
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA
−αφ − βsφ = Ds−1 t −2β+1 φ + sφ
shelves and
streams
shallow shelf
• choose β = 1/2 so equation has no “t” and further simplify to
approximation (SSA)
ﬂow-line ice shelf
numerical solution                              1
scaling up                                  −     2αsφ + s2 φ    = D sφ
2
next steps
practicalities
omitted models        • choose α = 1 so that quantity in square brackets simpliﬁes to a
derivative: 2αsφ + s2 φ = s2 φ
Numerical
modelling                   ﬁnding a “similarity” solution of heat equation 2
Ed Bueler

ice ﬂow: a
• simplify and integrate, choose C = 0, integrate again:
superﬁcial view
shallow ice                                         1
approximation (SIA)                                − s2 φ = D s φ + C
heat analogy &
2
numerics                                                       s
analogy                                               φ =−       φ
ﬁnite differences
2D
2
φ(s) = Ae−s               /(4D)
exact solutions

shallow ice
sheets
solving the SIA       • return to original variables, recalling s = t −β r :
is it right?
2                            2
application
+y 2 )/(4Dt)
the most basic                       u(t, r ) = A t −1 e−r       /(4Dt)
= A t −1 e−(x
shallow assumption
ﬂowline Stokes
derivation of SIA     • it is a spreading gaussian distribution
shelves and
streams
• heat equation is linear so all time-dependent solutions of the heat
shallow shelf
approximation (SSA)
equation are convolutions with this solution . . . which is why it is
ﬂow-line ice shelf      fundamental
numerical solution
scaling up

next steps
practicalities
omitted models
Numerical
modelling                                                                             similarity solutions
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice              • conclusion from previous slides: similarity variables for heat
approximation (SIA)

heat analogy &
equation are
numerics
analogy                                input scaling                                output scaling
ﬁnite differences
−β
exact solutions
s        =          t        x,        u(t, x)        =           t −α φ(s)
shallow ice
sheets
solving the SIA
is it right?
• dimension dependence:
application                                                        1D         2D          3D
the most basic
shallow assumption
ﬂowline Stokes
input scaling (t −β )               t −1/2     t −1/2      t −1/2
derivation of SIA            output scaling (t −α )              t −1/2       t −1      t −3/2
shelves and
streams
shallow shelf         historical note: In 1905 Einstein saw
approximation (SSA)
ﬂow-line ice shelf    that the average distance traveled in
numerical solution
time t by particles in thermal motion
√

x
scaling up

next steps            goes like t. This is a microscopic
practicalities
omitted models
explanation of the similarity variable
s = t −1/2 x.
t
Numerical
modelling                                                    Outline
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
numerics              1 ice ﬂow: an outsider’s superﬁcial view
analogy
ﬁnite differences
exact solutions

shallow ice
2 heat analogy & numerics
sheets
solving the SIA
is it right?
application
3 shallow ice sheets
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA     4 shelves and streams
shelves and
streams
shallow shelf
approximation (SSA)
5 next steps
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
Numerical
modelling                                                                                     similarity solution to the SIA
Ed Bueler

ice ﬂow: a            • jump forward to 1981 . . . breaking news:                                                      P. Halfar ﬁnds the
superﬁcial view
shallow ice                  fundamental solution of the SIA!                                                 [Halfar, 1981; 1983]
approximation (SIA)

heat analogy &
• when n = 3, Halfar’s 2D solution for Glen ﬂow law has
numerics
analogy
scalings
ﬁnite differences
exact solutions
H(t, r ) = t −1/9 φ(s),                                s = t −1/18 r
shallow ice
sheets
• . . . therefore, if an ice sheet starts steep and it ﬂows without
solving the SIA              additional accumulation then its thinning really slows down as
is it right?
application                  the shape ﬂattens out
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA     7000                         7000                          7000                         7000                         7000

shelves and           6000                         6000                          6000                         6000                         6000
streams
shallow shelf         5000                         5000                          5000                         5000                         5000
approximation (SSA)
ﬂow-line ice shelf    4000                         4000                          4000                         4000                         4000
numerical solution
3000                         3000                          3000                         3000                         3000
scaling up
2000                         2000                          2000                         2000                         2000
next steps
practicalities        1000                         1000                          1000                         1000                         1000
omitted models
0                            0                             0                            0                            0
0   300   600   900          0    300   600   900          0   300   600   900          0   300   600   900          0   300   600   900

times 1, 10, 100, 1000, 10000 years
Numerical
modelling                                                    similarity solution to the SIA 2
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
ﬁnite differences
exact solutions

shallow ice
sheets
solving the SIA
is it right?
application
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA

shelves and
streams
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models

frames from t = 4 months to t = 106 years, equal spaced in exponential time
Numerical
modelling                                                   similarity solution to the SIA 3
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)
• for n = 3 the solution is:
heat analogy &
numerics
                                             3/7
analogy
1/9                             1/18           4/3
ﬁnite differences
t0            1 −              t0               r
exact solutions
H(t, r ) = H0                                                                   
shallow ice                                        t                               t               R0
sheets
solving the SIA
is it right?
application             where
3    4
the most basic
shallow assumption                                             1           7         R0
ﬂowline Stokes                                      t0 =                              7
derivation of SIA                                             18Γ          4         H0
shelves and
streams                 and H0 , R0 are central height and ice cap radius at t = t0
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
• you choose H0 and R0 , and these determine the
numerical solution      “characteristic” time t0
scaling up

next steps                             more on this in [Nye, 2000] and [Bueler and others, 2005]
practicalities
omitted models
Numerical
modelling                                                                              ﬁnding Halfar solution
Ed Bueler
the method is the same as that which found the fundamental solution of
ice ﬂow: a            the heat equation:
superﬁcial view
shallow ice             •   allow scaling of the input and output:          s = t −β r ,       H(t, r ) = t −α ϕ(s)
approximation (SIA)

heat analogy &
•   insert into Ht =          5
· ΓH | H|        2
H                  (← n = 3 case for simplicity)
numerics
analogy                 •   to eliminate t from the equation, choose 7α + 4β = 1 ; get
ﬁnite differences
exact solutions
−1         5       2
shallow ice                                          −αφ − βsϕ = s                     Γsϕ |ϕ | ϕ
sheets
solving the SIA
is it right?            •   to write as equality of derivatives, choose −α + 2β = 0 :
application
the most basic
shallow assumption                                                  2                  5       2
−βs ϕ            =    Γsϕ |ϕ | ϕ
ﬂowline Stokes
derivation of SIA

shelves and             •   integrate this and use ﬁniteness of ϕ(0), ϕ (0), and use ϕ (s) ≤ 0, to get separable
streams                     ODE
4       3
shallow shelf                                                   βs = Γϕ (−ϕ )
approximation (SSA)
ﬂow-line ice shelf      •   suppose H(t0 , 0) = H0 for some t0 > 0 to ﬁnd formula
numerical solution
scaling up                                                                                         1/3
7/3           α    7/3       7   β                 4/3
next steps                                          ϕ(s)        = (t0 H0 )         −                     s
practicalities
4   Γ
omitted models
•   write as formula for H(t, r ), on previous slide, by solving boxed equations above to get
α = 1/9, β = 1/18
•                               −β
margin location equation φ(t0 R0 ) = 0 gives formula for t0 on previous slide
Numerical
modelling                                                         evaluating Halfar solution
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)
function H = halfar(t,x,y)
heat analogy &
numerics
analogy               g = 9.81;   rho = 910.0;    secpera = 31556926;
ﬁnite differences     n = 3;
exact solutions       A = 1.0e-16/secpera;     Gamma = 2 * A * (rho * g)^3 / 5;
shallow ice
sheets                H0 = 3600;   R0 = 750e3;
solving the SIA       alpha = 1/9;   beta = 1/18;
is it right?
t0 = (beta/Gamma) * (7/4)^3 * (R0^4/H0^7);
application
the most basic
shallow assumption    r = sqrt(x.*x + y.*y);
ﬂowline Stokes        r = r / R0;   t = t / t0;
derivation of SIA     inside = max(0, 1 - (r / t^beta).^((n+1) / n));
shelves and           H = H0 * inside.^(n / (2*n+1)) / t^alpha;
streams
shallow shelf                        http://www.dms.uaf.edu/~bueler/karthaus/mfiles/halfar.m
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up              • it is easy to compute it, so . . .
next steps
practicalities          • please use it to verify SIA codes, o.k.?
omitted models
Numerical
modelling                              is the Halfar solution good for any modeling?
Ed Bueler

ice ﬂow: a
superﬁcial view       • Nye and others [2000] compared the long-time
shallow ice
approximation (SIA)       consequences of different ﬂow laws for the South Mars Polar
heat analogy &
numerics
Cap
analogy
ﬁnite differences
• they compared the long-time behavior of the corresponding
exact solutions
Halfar solutions for the different ﬂow laws, with both CO2 ice
shallow ice
sheets                  and H2 O ice parameters
solving the SIA
is it right?          • conclusions:
application
the most basic
. . . none of the three possible [CO2 ] ﬂow laws will
shallow assumption
ﬂowline Stokes                   allow a 3000-m cap, the thickness suggested by
derivation of SIA
stereogrammetry, to survive for 107 years, indicating
shelves and
streams                          that the south polar ice cap is probably not
shallow shelf
approximation (SSA)              composed of pure CO2 ice . . . the south polar cap
ﬂow-line ice shelf
numerical solution
probably consists of water ice, with an unknown
scaling up
next steps
practicalities
omitted models
• in 2008, NASA’s Phoenix lander found water ice in the north
polar cap on Mars4

4 P. Smith and 35 others, 2009. “H O at the Phoenix landing site,” Science 325(5936), 58–61.
2
Numerical
modelling                                               on degenerate diffusivity
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
• recall that the SIA is
approximation (SIA)

heat analogy &                   Ht = M +    · (D   h)   where   D = ΓH n+2 | h|n−1
numerics
analogy
ﬁnite differences     • “degeneracy” of D happens in two ways:
exact solutions

shallow ice
sheets                              why D → 0              so what?
solving the SIA
is it right?                                        H and h are continuous
domes         h→0
application
the most basic
but 2 h is singular
shallow assumption
H is continuous
ﬂowline Stokes           margins      H→0
derivation of SIA                                      but h is singular
shelves and
streams
shallow shelf
approximation (SSA)
• numerically, margins are worse than domes
ﬂow-line ice shelf
numerical solution    • degenerate diffusion equations are free boundary problems
scaling up

next steps
• transformation η = H (2n+2)/n is known to help . . . for ﬂat beds
practicalities
omitted models
• as a practical matter, numerical solvers for the SIA must
handle ut = M +      · (D   u) for very general diffusivity
Numerical
modelling                                                  computing diffusivity in SIA
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
• for numerical stability we want
ﬁnite differences
exact solutions
D = ΓH n+2 | h|n−1 on the
shallow ice
staggered grid
sheets                                                               k+1
solving the SIA       • various schemes proposed; stencil
is it right?
application
for Mahaffy [1976] scheme at right
k
the most basic
shallow assumption    • all schemes involve
ﬂowline Stokes
derivation of SIA          ◦ averaging H
k-1
shelves and                ◦ differencing h
streams                                                                    j-1   j   j+1
shallow shelf
◦ in a “balanced” way, for better
approximation (SSA)
ﬂow-line ice shelf
accuracy,
numerical solution
scaling up
to put D on the staggered grid
next steps
practicalities
omitted models
Numerical
modelling                                          SIA implementation: ﬂat bed case
Ed Bueler

function [H,dtlist] = siaflat(Lx,Ly,J,K,H0,deltat,tf)
ice ﬂow: a
superﬁcial view
shallow ice           g = 9.81;    rho = 910.0;    secpera = 31556926;
approximation (SIA)   A = 1.0e-16/secpera;    Gamma = 2 * A * (rho * g)^3 / 5;
heat analogy &        H = H0;
numerics
analogy               dx = 2 * Lx / J;   dy = 2 * Ly / K;
ﬁnite differences
N = ceil(tf / deltat);     deltat = tf / N;
exact solutions
j = 2:J;     k = 2:K;
shallow ice           nk = 3:K+1;    sk = 1:K-1;     ej = 3:J+1;               wj = 1:J-1;
sheets
solving the SIA
is it right?
t = 0;    dtlist = [];
application           for n=1:N
the most basic          Hup = 0.5 * ( H(j,nk) + H(j,k) );     Hdn = 0.5 * ( H(j,k) + H(j,sk) );
shallow assumption
Hrt = 0.5 * ( H(ej,k) + H(j,k) );     Hlt = 0.5 * ( H(j,k) + H(wj,k) );
ﬂowline Stokes
derivation of SIA
a2up = (H(ej,nk) + H(ej,k) - H(wj,nk) - H(wj,k)).^2 / (4*dx)^2 + ...
(H(j,nk) - H(j,k)).^2 / dy^2;
shelves and
streams
a2dn = (H(ej,k) + H(ej,sk) - H(wj,k) - H(wj,sk)).^2 / (4*dx)^2 + ...
shallow shelf
(H(j,k) - H(j,sk)).^2 / dy^2;
approximation (SSA)     a2rt = (H(ej,k) - H(j,k)).^2 / dx^2 + ...
ﬂow-line ice shelf              (H(ej,nk) + H(j,nk) - H(ej,sk) - H(j,sk)).^2 / (4*dy)^2;
numerical solution
a2lt = (H(j,k) - H(wj,k)).^2 / dx^2 + ...
scaling up
(H(wj,nk) + H(j,nk) - H(wj,sk) - H(j,sk)).^2 / (4*dy)^2;
next steps              Dup = Gamma * Hup.^5 .* a2up; Ddn = Gamma * Hdn.^5 .* a2dn;
practicalities
Drt = Gamma * Hrt.^5 .* a2rt; Dlt = Gamma * Hlt.^5 .* a2lt;
omitted models
t = t + deltat;      dtlist = [dtlist dtadapt];
end
http://www.dms.uaf.edu/~bueler/karthaus/mfiles/siaflat.m
Numerical
modelling                                                    the worst ice sheet
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
• adaptive time-stepping works: siaflat.m takes short time
approximation (SIA)
steps when the driving stress ρgH| h| is large, and then
heat analogy &
numerics                longer steps when it is “easier”
analogy
ﬁnite differences     • for example, a worst-case ice sheet is thick but with a very
exact solutions
bumpy surface
shallow ice
sheets
solving the SIA
• roughice.m generates initial ice sheet, which is thick and
is it right?            with a white-noise surface, then runs for 50 years
application
the most basic
shallow assumption
• initial surface, ﬁnal surface, and time-steps shown below
ﬂowline Stokes
derivation of SIA

shelves and
streams
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
Numerical
modelling                                veriﬁcation of numerical ice ﬂow codes
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)   • how do we make sure a numerical scheme and its
heat analogy &
numerics
implementation are right? some available techniques
analogy                   1   don’t make any mistakes
ﬁnite differences
exact solutions           2   “eyeball” the results when you change the code . . . if it looks
shallow ice                   right it is right
sheets
solving the SIA
3   from the beginning, build-in a comparison to an exact solution
is it right?
application
• where to get exact solutions for ice ﬂow models?
the most basic
shallow assumption
ﬂowline Stokes           similarity solutions to SIA                   [Halfar, 1983; Bueler and others, 2005]
derivation of SIA
manufactured solutions to SIA; some include   [Bueler and others 2005; 2007]
shelves and                   thermo-coupling
streams                  SSA solutions in ﬂow line and                 [Schoof, 2006; van der Veen, 1985]
shallow shelf
approximation (SSA)           cross-ﬂow cases
ﬂow-line ice shelf       Blatter solution in ﬂow line                  [Glowinski and Rappaz, 2003]
numerical solution       ﬂow line Stokes solutions in constant         [Ladyzhenskaya, 1963]
viscosity case ( 4 ψ = 0)
scaling up

next steps               manufactured solutions to Stokes              [Sargent and Fastook, 2010]
practicalities
equations with variable viscosity
omitted models
Numerical
modelling                                               verifying our SIA code siaflat.m
Ed Bueler

ice ﬂow: a
superﬁcial view                                                                       • this code calls
shallow ice
approximation (SIA)                                                                     siaflat.m
heat analogy &
numerics
• with initial condition
function [avthkerr,maxthkerr] = verifysia(J,dtyears)
analogy                                                                                 which is a t1 Halfar
ﬁnite differences     if nargin<1, J=40; end;
exact solutions       if nargin<2, dtyears=10.0; end;                                   solution
shallow ice
sheets
L = 1200e3;    dx = 2 * L / J;                                  • numerical code runs to
[x,y] = meshgrid(-L:dx:L, -L:dx:L);
solving the SIA                                                                         later time t2 and the
is it right?          t1 = 200;    t2 = 20000;    secpera = 31556926;
application           H1 = halfar(t1 * secpera,x,y);                                    compares to Halfar at t2
the most basic
shallow assumption    [H2approx,dtlist] = siaflat(L,L,J,J,H1,dtyears*secpera,...
ﬂowline Stokes                                    (t2-t1)*secpera);
derivation of SIA
map of thickness error from
H2exact = halfar(t2 * secpera,x,y);                              » verifysia(160,3.0):
shelves and           error = H2approx - H2exact;
streams
shallow shelf         avthkerr = sum(sum(abs(error)))/(J+1)^2;
approximation (SSA)   maxthkerr = max(max(abs(error)));
ﬂow-line ice shelf    erreta = H2approx.^(8/3) - H2exact.^(8/3);
numerical solution
maxeta = max(max(H2exact.^(8/3)));
scaling up

next steps             http://www.dms.uaf.edu/~bueler/karthaus/mfiles/verifysia.m
practicalities
omitted models
Numerical
modelling                                           verifying our SIA code 2
Ed Bueler

ice ﬂow: a
superﬁcial view       >> verifysia(20,24.0);
shallow ice
approximation (SIA)   errors for 20 x 20 grid:
heat analogy &
numerics
average thickness error =    22.335
analogy               maximum thickness error =    227.275   Trust but verify.
ﬁnite differences
exact solutions
>> verifysia(40,12.0);                 (Ronald Reagan)

shallow ice
errors for 40 x 40 grid:
sheets                average thickness error =    9.519
solving the SIA
is it right?          maximum thickness error =    241.040
application
the most basic
>> verifysia(80,6.0);
shallow assumption
ﬂowline Stokes
errors for 80 x 80 grid:
derivation of SIA     average thickness error =    2.937
shelves and           maximum thickness error =    159.442
streams
shallow shelf         >> verifysia(160,3.0);
approximation (SSA)
ﬂow-line ice shelf
errors for 160 x 160 grid:
numerical solution    average thickness error =    0.988     ﬁgure 2 in [Huybrechts and
scaling up
maximum thickness error =    103.456   others, 1996]
next steps
practicalities        >> verifysia(200,2.0);
omitted models
errors for 200 x 200 grid:
average thickness error =    0.827
maximum thickness error =    94.561
Numerical
modelling                                        numerical mass conservation
Ed Bueler

ice ﬂow: a
superﬁcial view       • in addition to measuring geometry errors, we might want a
shallow ice
approximation (SIA)     numerical code not to lose mass!
heat analogy &
numerics
• again, the Halfar solution is a perfect tool
analogy
ﬁnite differences
◦ because the Halfar exact solution has exact (continuum)
exact solutions             volume conservation
shallow ice
sheets
• add numerical volume calculation to verifysia.m:
solving the SIA
is it right?
>> verifysia(20)
application
vol1       = 3.96112407720492e+15
the most basic
shallow assumption      ...
ﬂowline Stokes
derivation of SIA       vol2approx = 3.96112407720492e+15
shelves and             voldiff = 1.50000000000000e+00
streams
shallow shelf         • huh? how can it be that good?
approximation (SSA)
ﬂow-line ice shelf        ◦ the Mahaffy [1976] scheme has local numerical mass
numerical solution
scaling up                  conservation
next steps                ◦ but what about the margin where H → 0?
practicalities
omitted models
◦ siaflat.m implements boundary condition “H ≥ 0”
◦ apparently this boundary condition produces global numerical
mass conservation to rounding error!
Numerical
modelling                                                                  application to the Antarctic ice sheet
Ed Bueler

ice ﬂow: a            • at Karthaus 2009, a computer project: modify siaflat.m to
superﬁcial view
shallow ice                            • allow arbitrary bed elevation b(x, y )
approximation (SIA)
• allow arbitrary mass balance M(x, y )
heat analogy &
numerics                               • apply a marine-margin condition
analogy
ﬁnite differences                and then simulate the Antarctic ice sheet
exact solutions

shallow ice
• my solution:
sheets
solving the SIA                        • add 8 lines to siaflat.m in total, creating siageneral.m
is it right?
• use the SeaRISE 5km data set for present-day Antarctica
application
the most basic                             (NetCDF)
shallow assumption
ﬂowline Stokes                         • add a page or so of code for pre- and post-processing
derivation of SIA
• do 40 ka run starting with present-day geometry
shelves and
streams                        initial (observed) surface elevation                                40 ka (modeled) surface elevation
4000                                                                4000
shallow shelf
approximation (SSA)
-2000                                               3500            -2000                                               3500
ﬂow-line ice shelf
numerical solution                                                                 3000                                                                3000
-1000                                                               -1000
scaling up                                                                         2500                                                                2500
y (km)

y (km)
next steps                        0
2000
0
2000
practicalities
1000                                                1500            1000                                                1500
omitted models
1000                                                                1000
2000                                                                2000
500                                                                 500

3000                                                                3000
0                                                                   0
-2000   -1000    0     1000   2000   3000                           -2000   -1000    0     1000   2000   3000
x (km)                                                              x (km)
Numerical
modelling                                    application to the Antarctic ice sheet 2
Ed Bueler

ice ﬂow: a
superﬁcial view                                                                            34

shallow ice
approximation (SIA)   • volume grows
32
heat analogy &
numerics
◦ levels out at about 33 × 106 km3
◦ compare present-day observed

volume (10 km )
3
analogy                                                                                    30

volume of 25 × 106 km3

6
ﬁnite differences
exact solutions

shallow ice               ◦ exponential time of about 15 ka,                               28

sheets
◦ but exponential time is much
solving the SIA                                                                            26
is it right?                longer if our model were
application
the most basic
thermocoupled, because of long                                 24
0   5000   10000   15000   20000    25000   30000   35000   40000
shallow assumption
ﬂowline Stokes

• I used A = 3 × 10−16 Pa−3 a−1
derivation of SIA

shelves and
streams                      ◦ which is an “enhancement factor” of 3, relative to the EISMINT I
shallow shelf
approximation (SSA)             [Huybrechts and others, 1996] value of A = 10−16 used earlier
ﬂow-line ice shelf
numerical solution      • an enhancement factor of 24 would make the volume level
scaling up

next steps
out at about the present-day value
practicalities
omitted models
• pop quiz: why is that not a good idea?
Numerical
modelling                                                          the most basic shallow assumption
Ed Bueler

ice ﬂow: a                                                                            air
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
• careful derivation of the SIA is next                                  ice
numerics
analogy               • makes “shallowness assumptions”
ﬁnite differences
exact solutions       • the most basic is a limitation which
shallow ice
sheets
applies to all shallow ice theoriesa                                   bedrock

solving the SIA           but not to the general Stokes
is it right?
application               theory:
the most basic
shallow assumption
ﬂowline Stokes            the surface and base of the ice are
derivation of SIA
given by differentiable functions
shelves and
streams                     z = h(t, x, y ) and z = b(t, x, y )
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
• the conﬁgurations at right, with
numerical solution
scaling up
overhang in surface, are not
next steps
allowed
practicalities
a e.g. SIA, SSA, Blatter-Pattyn, SIA+SSA hybrids, . . .
omitted models
Numerical
modelling                                  kinematic and mass continuity equations
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
• what does this “most basic shallow assumption” get you?
approximation (SIA)

heat analogy &
• answer: a map-plane mass continuity equation
numerics
analogy               • consider these three equations:
ﬁnite differences
exact solutions
◦ the surface kinematical equation
shallow ice                ◦ the base kinematical equation
sheets
solving the SIA
◦ the map-plane mass continuity equation
is it right?
application
• all three equations are important for shallow ice sheets and
the most basic
shallow assumption
shelves, but any two imply the third
ﬂowline Stokes
derivation of SIA
• the key facts we need:
shelves and               ◦ the incompressibility of ice
streams
shallow shelf
approximation (SSA)
ux + vy + wz = 0
ﬂow-line ice shelf
numerical solution
scaling up                 ◦ the Leibniz rule for differentiating integrals
next steps
ˆ   f (x)                                                          ˆ   f (x)
practicalities               d
omitted models                                h(x, y ) dy   = f (x)h(x, f (x)) − g (x)h(x, g(x)) +               hx (x, y ) dy
dx    g(x)                                                              g(x)
Numerical
modelling                       kinematic and mass continuity equations 2
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice           • let a be the surface mass balance function (a > 0 is
approximation (SIA)

heat analogy &
accumulation)
numerics
analogy               • let s be the basal melt rate function (s > 0 is basal melting)
ﬁnite differences
exact solutions       • deﬁne the map-plane ﬂux of ice,
ˆ
shallow ice
sheets                                               h
solving the SIA
is it right?                              q=                           ¯ ¯
(u, v ) dz = (u , v ) H
application                                      b
the most basic
shallow assumption
ﬂowline Stokes        • the three equations are:
derivation of SIA

shelves and
streams
shallow shelf
surface kinematical           ht = a − u h hx − v h hy + w   h
(3)
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up                    base kinematical           bt = s − u b bx − v b by + w   b
(4)
next steps
practicalities
omitted models
mass continuity     Ht = M −             · q where M = a − s     (5)
Numerical
modelling                                kinematic and mass continuity equations 3
Ed Bueler

ice ﬂow: a            for example, here is a sketch of how to show                        (3) & (4) =⇒ (5):
superﬁcial view
shallow ice
approximation (SIA)     1 recalling H = h − b and M = a − s, the difference of (3) and (4) is
heat analogy &
numerics                                   Ht = M − u h hx − v h hy + w        h
+ u b bx + v b by − w             b
(*)
analogy
ﬁnite differences
exact solutions         2 incompressibility gives difference of surface and base values of w by integration,
shallow ice                                                                                    ˆ       h
sheets
solving the SIA                          wz = −(ux + vy )     =⇒      w   h
−w      b
=−               (ux + vy ) dζ
is it right?                                                                                       b
application
the most basic             which reduces (*) to:
shallow assumption
ﬂowline Stokes                                                                                 ˆ        h
derivation of SIA
Ht = M − u h hx − v h hy + u b bx + v b by −                        (ux + vy ) dζ   (**)
shelves and                                                                                         b
streams
shallow shelf                                                          ´h
approximation (SSA)     3 on the other hand, the Leibniz rule on q =      b
(u, v ) dz gives
ﬂow-line ice shelf
numerical solution                                                                          ˆ       h
scaling up                                   · q = u h hx + v h hy − u b bx − v b by +                  (ux + vy ) dz
b
next steps
practicalities
omitted models
4 . . . which reduces (**) to
Ht = M −           ·q
Numerical
modelling                       kinematic and mass continuity equations 4
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &        • papers have lots of calculations like on the previous slide
numerics
analogy               • usually mixed in with other arguments about shallowness
ﬁnite differences
exact solutions
• most ice sheet models use the mass continuity equation to
shallow ice
sheets                  describe change in ice sheet geometry,
solving the SIA
is it right?          • . . . but they could instead use the surface kinematical
application
the most basic
equation
shallow assumption
ﬂowline Stokes        • in using a stress balance in a shallow theory, all we need to
derivation of SIA

shelves and
do is get the horizontal velocity (u, v ), which gives q and thus
streams                 the mass continuity equation
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
• for example, the derivation of the SIA from Stokes needs to
numerical solution
scaling up
only go far enough to ﬁnd an expression for the horizontal
next steps
velocity (u, v )
practicalities
omitted models
Numerical
modelling                                                                  origin of SIA
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &        where does the “shallow ice approximation” come from?
numerics
analogy                 • historically:
ﬁnite differences
exact solutions
• Fowler and Larson [1978], Morland and Johnson [1980], and
shallow ice                    Hutter [1983]
sheets
• fairly recent compared to the ﬂow law of ice (1950s) and
solving the SIA
is it right?                   Stokes linearly-viscous model (∼1860)
application
the most basic          • logically:
shallow assumption
ﬂowline Stokes               • by a “small-parameter argument”—depth-to-width ratio is
derivation of SIA

shelves and
small—to rescale the Stokes model
streams                      • analogous to an old argument giving shallow water equations
shallow shelf
approximation (SSA)            from Navier-Stokes
ﬂow-line ice shelf
numerical solution      • more precisely:
scaling up
• as sketched in the next slides, for the x, z plane ﬂow case only
next steps
practicalities
omitted models
Numerical
modelling                           Stokes equations in x, z plane ﬂow case
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &        • suppose all ﬂow is in the (screen) x, z plane
numerics
analogy               • . . . so out-of-plane velocity is zero (v = 0) and all quantities
ﬁnite differences
exact solutions         are constant in y direction (∂/∂y = 0)
shallow ice
sheets
• the symmetric strain rate tensor has many zeros
solving the SIA
is it right?
               
D11      0   D13
application
1    ∂ui   ∂uj
the most basic
shallow assumption
Dij =          +        = 0       0    0 
2    ∂xj   ∂xi
ﬂowline Stokes
derivation of SIA
D13      0   D33
shelves and
streams               • incompressibility of ice in plane ﬂow case says D has trace
shallow shelf
approximation (SSA)     zero:
ﬂow-line ice shelf
numerical solution                        ux + 0 + wz = D11 + D33 = 0
scaling up

next steps            • . . . so D33 = −D11 in this case
practicalities
omitted models
Numerical
modelling                          Stokes equations in x, z plane ﬂow case 2
Ed Bueler

ice ﬂow: a
superﬁcial view       • recall that the Cauchy stress tensor σij , with its pressure part
shallow ice
approximation (SIA)     p removed, is called “deviatoric”:
heat analogy &
numerics
analogy                    τij = σij + pδij      where         p = −(1/3)(σ11 + σ22 + σ33 )
ﬁnite differences
exact solutions

shallow ice
• the deviatoric stress tensor τij is proportional to the tensor
sheets
solving the SIA
Dij , in the isotropic (Glen law) case:
is it right?

Dij = A(T )τ n−1 τij
application
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA                                                                          1/2
shelves and
• . . . where    2τ 2 = 2(II τ )2 = τij τij   so          2     2
τ = τ11 + τ13
streams
shallow shelf
• so τij is also symmetric, and has trace zero,
approximation (SSA)
ﬂow-line ice shelf    • and τij has only two nonzero entries τ11 , τ13 :
numerical solution
scaling up                                                         
next steps                                           τ11     0 τ13
practicalities                                      0       0  0 
omitted models

τ13     0 −τ11
Numerical
modelling                         Stokes equations in x, z plane ﬂow case 3
Ed Bueler

ice ﬂow: a
superﬁcial view
• in the 3D case: n = 3, Glen, isothermal Stokes equations:
shallow ice
approximation (SIA)

heat analogy &
ux + vy + wz = 0
numerics
analogy
0=       · σ + ρg
ﬁnite differences
exact solutions                                           Dij = Aτ 2 τij
shallow ice
sheets
solving the SIA       • x, z plane ﬂow case: from simpliﬁcations on previous slides,
is it right?
application
the isothermal Stokes equations say
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA
ux + wz = 0
shelves and
streams
px = τ11,x + τ13,z
shallow shelf
approximation (SSA)                                pz = τ13,x − τ11,z − ρg            (6)
ﬂow-line ice shelf
numerical solution                                 ux = Aτ 2 τ11
scaling up

next steps                                  uz + wx = 2Aτ 2 τ13
practicalities
omitted models

• ﬁve scalar equations in ﬁve scalar unknowns (u, w, p, τ11 , τ13 )
but still fairly hard to solve . . .
Numerical
modelling                        Stokes equations in x, z plane ﬂow case 4
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)
• previous equations apply within the ice ﬂuid
heat analogy &
numerics              • . . . but boundary conditions are needed
analogy
ﬁnite differences     • at the upper surface there is no stress from the air: σij nj = 0
exact solutions

shallow ice
where n is a normal vector to the surface
sheets
solving the SIA       • in plane ﬂow case use n = (hx , 0, −1), which is gradient of
is it right?
application
F (x, z) = h(x) − z, which is constant on z = h(x)
the most basic
shallow assumption    • get stress equations at surface:
ﬂowline Stokes
derivation of SIA

shelves and                                           τ13 = (τ11 − p)hx
streams
shallow shelf
approximation (SSA)
p + τ11 + τ13 hx = 0
ﬂow-line ice shelf
numerical solution
scaling up
• we will only consider non-sliding SIA, so base boundary
next steps              conditions are
practicalities
omitted models
u=0      and    w =0
Numerical
modelling                                                                scaling the variables
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)
• following Chapter 18 of [Fowler, 1997], Mathematical Models
heat analogy &            in the Applied Sciences
numerics                    • which does a much more general job than here: 3D case with conservation of
analogy
energy, while we only do 2D isothermal case
ﬁnite differences
exact solutions             • we also miss Andrew standing on a table and singing, too
shallow ice
sheets                • consider ice sheet of typical depth d and length , and with
solving the SIA
is it right?
velocity ([u]), deviatoric stress ([τ ]), and softness ([A]) scales
application
the most basic        •    = d/ is the “aspect ratio” of the ice sheet
shallow assumption
ﬂowline Stokes        • change to new starred variables using scales:
derivation of SIA

shelves and
streams                                x = x∗                                   u = [u]u ∗
shallow shelf
approximation (SSA)                    z = d z∗                                 w = [u]w ∗
ﬂow-line ice shelf
∗                                         ∗
numerical solution                     h=dh                                    τ11 = [τ ]τ11
scaling up

next steps
H = d H∗                                            ∗
τ13 = [τ ]τ13
∗
practicalities
omitted models
A = [A]A                  p − ρg(h − z) = [τ ]p∗
Numerical
modelling                                                  scaling the variables 2
Ed Bueler

ice ﬂow: a
superﬁcial view       • we will necessarily choose relations among the scales
shallow ice
approximation (SIA)
• the choice of such relations affects the outcome and is
heat analogy &
numerics                subject to “scientiﬁc judgment” (i.e. controversy)
analogy
ﬁnite differences     • “. . . the purpose is to attain ‘properly scaled’ equations in
exact solutions

shallow ice
which the largest dimensionless parameters are numerically
sheets
solving the SIA
of order one . . . ” [Fowler, 1997]
is it right?
application
• method: we rewrite the Stokes equations and the boundary
the most basic
shallow assumption
conditions in the newly (starred) variables, then check that
ﬂowline Stokes
derivation of SIA
the scaling relations give reasonable values, then remove
shelves and
terms to get a new, more-solvable model, the SIA
streams
shallow shelf         • the scaling relations used here are:
approximation (SSA)
ﬂow-line ice shelf
numerical solution                 balance of shear stress
scaling up                                                          [τ ] = ρgd
next steps
practicalities
[u]
omitted models
balance of ﬂow law scales               = 2[A][τ ]3
d
Numerical
modelling                                                                         scaling the equations
Ed Bueler
example:
ice ﬂow: a
superﬁcial view         • recall this equation from Stokes equations (6):
shallow ice
approximation (SIA)

heat analogy &                                             px = τ11,x + τ13,z
numerics
analogy
ﬁnite differences       • replace all variables with their starred versions:
exact solutions

shallow ice                              [τ ]    ∗       dρg    ∗           [τ ]    ∗          [τ ] ∗
sheets                                          px ∗ +         hx ∗ =              τ11,x ∗ +       τ13,z ∗
solving the SIA
d
is it right?
application             • multiply by /( [τ ]):
the most basic
shallow assumption
ﬂowline Stokes                                   ∗       dρg ∗         ∗           ∗
derivation of SIA                               px ∗ +         hx ∗ = τ11,x ∗ +   τ13,z ∗
[τ ]                  d
shelves and
streams
shallow shelf           • use relation [τ ] = ρgd and = d/ :
approximation (SSA)
ﬂow-line ice shelf
∗        −2 ∗        ∗                −2 ∗
numerical solution                              px ∗ +      hx ∗   = τ11,x ∗ +           τ13,z ∗
scaling up

next steps
practicalities
• re-arrange to taste and remove stars:
omitted models
2
hx = τ13,z +            (τ11,x − px )

• . . . do this for each of the Stokes equations
Numerical
modelling                                                           scaling the equations 2
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
• Stokes system now appears this way (stars are removed):
approximation (SIA)

heat analogy &
numerics                                   ux + wz = 0
analogy
2
ﬁnite differences                                    hx = τ13,z +           (τ11,x − px )
exact solutions

shallow ice                                     pz = τ13,x − τ11,z
sheets
solving the SIA                                      1
is it right?
2      2
ux = A 2 τ11 + τ13 τ11
application                                          2
the most basic
shallow assumption
ﬂowline Stokes
2      2
uz + 2 wx = A 2 τ11 + τ13 τ13
derivation of SIA

shelves and           • this is still full Stokes!
streams
shallow shelf
approximation (SSA)
• surface conditions are
ﬂow-line ice shelf
numerical solution                                                  2
scaling up                                                τ13 =         (τ11 − p) hx
next steps
practicalities
p + τ11 + τ13 hx = 0
omitted models

• base conditions are u = 0 and w = 0 as before
Numerical
modelling                                                      checking the scales
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
• the scale for vertical velocity is already chosen as [u]
approximation (SIA)

heat analogy &
• . . . and we now identify this scale with a typical value for
numerics
analogy
accumulation: [a] = [u]
ﬁnite differences
exact solutions
• thus we have chosen these scale relations:
shallow ice
sheets                                       d
solving the SIA                          =                         [τ ] = ρgd
is it right?
application
[u]
the most basic
shallow assumption                        = 2[A][τ ]3              [a] = [u]
ﬂowline Stokes                         d
derivation of SIA

shelves and           • to check their reasonableness, we will determine other scale
streams
shallow shelf
factors in terms of , [a], and [A], which are taken from
approximation (SSA)
ﬂow-line ice shelf      observations:
numerical solution
scaling up
= 3000 km,
next steps
practicalities
[a] = 0.1 m a−1
omitted models

[A] = 2 × 10−16 Pa−3 a−1
Numerical
modelling                                                         checking the scales 2
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice           • solving for d, , [τ ], [u] gives
approximation (SIA)

heat analogy &
numerics                                               4        1/8
[a]
analogy
d=                           ≈ 3600 m
ﬁnite differences
exact solutions
2[A](ρg)3
shallow ice                                     d
sheets                                      =       ≈ 10−3
solving the SIA
is it right?
application                             [τ ] = ρgd ≈ 4 × 104 Pa = 0.4 bar
the most basic
shallow assumption
ﬂowline Stokes                          [u] = 2[A]d[τ ]3 ≈ 80 m a−1
derivation of SIA

shelves and
streams
• these are all reasonable values for large ice sheets
shallow shelf
approximation (SSA)   • it is reasonable to proceed to delete terms from the scaled
ﬂow-line ice shelf
numerical solution
equations to get the reduced (shallow) theory, the SIA
scaling up
• the actual measure of the quality of SIA is the comparison of
next steps
practicalities          its predictions to observations, when modeling various real
omitted models
ice sheets etc., not these scalings
Numerical
modelling                                                             ﬁnally get to the SIA
Ed Bueler

ice ﬂow: a
• setting = 0 in the scaled Stokes equations, get these equations
superﬁcial view         (left column) and boundary conditions (right):
shallow ice
approximation (SIA)                                                                               ∗
ux + wz = 0                                              τ13   h
=0
heat analogy &
numerics                                ∗
analogy                              hx = τ13,z                 p   h
+ τ11   h
+ τ13 h hx = 0
ﬁnite differences
exact solutions                      pz = τ13,x − τ11,z
shallow ice
sheets                                      1 2                                                   †
ux =     Aτ13 τ11                                    u   b
=0
solving the SIA
is it right?
2
†   3
application
the most basic
uz = Aτ13                                           w    b
=0
shallow assumption
ﬂowline Stokes        • combine the two marked “∗”, integrating vertically from surface
derivation of SIA

shelves and
z = h to arbitrary elevation z, to get
streams
shallow shelf                                         τ13 = −(h − z)hx
approximation (SSA)
ﬂow-line ice shelf
numerical solution
• combine the two marked “†”, integrating vertically from base z = b
scaling up              to arbitrary level z, and use the new expression for τ13 , to get
next steps                                     ˆ z                   ˆ z
3
τ13 dζ = −A(hx )3     (h − ζ)3 dζ
practicalities
omitted models
u(z) = A
b                        b
A
= − (hx )3 H 4 − (h − z)4
4
Numerical
modelling                                                          ﬁnally get to the SIA 2
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
• remember I said that it sufﬁces to ﬁnd an expression for the
approximation (SIA)
horizontal velocity, in order to state the SIA? we are there!
heat analogy &
numerics              • the previous slides have been in scaled quantities; now we recall the
analogy
ﬁnite differences       stars:
A∗ ∗
u ∗ = − (hx ∗ )3 (H ∗ )4 − (h∗ − z ∗ )4
exact solutions

shallow ice                                     4
sheets
solving the SIA       • and then unscale to ﬁnd a dimensional expression for u, using
is it right?
application
∗
hx ∗ = −1 hx , and using scaling relations:
the most basic
shallow assumption
[u]A
u = [u]u ∗ = −          −3
(hx )3 d −4 H 4 − (h − z)4
ﬂowline Stokes
derivation of SIA
4[A]
shelves and
streams                                       [u]
shallow shelf                       =−                    A(hx )3 H 4 − (h − z)4
approximation (SSA)                        4[A] 3 d 4
ﬂow-line ice shelf
numerical solution                         2[A]d(ρgd )3
scaling up                          =−                          A(hx )3 H 4 − (h − z)4
next steps
4[A] 3 d 4
practicalities                         1
omitted models                      = − A(ρg)3 (hx )3 H 4 − (h − z)4
2
Numerical
modelling                                                      ﬁnally get to the SIA 3
Ed Bueler

ice ﬂow: a
superﬁcial view                                             ´h
shallow ice
approximation (SIA)
• recall that the ﬂux of ice is q =   b
u(z) dz
heat analogy &        • so we integrate the expression for u(z):
numerics
analogy
ﬁnite differences
ˆ h
1
exact solutions
q = − A(ρg)3 (hx )3     H 4 − (h − z)4 dz
shallow ice                            2               b
sheets
solving the SIA                        1                     1
is it right?                        = − A(ρg)3 (hx )3 H 5 − (H 5 − 0)
application                            2                     5
the most basic
shallow assumption
2
ﬂowline Stokes
= − A(ρg)3 (hx )3 H 5
derivation of SIA
5
shelves and
streams
shallow shelf
• the mass continuity equation now says:
approximation (SSA)
ﬂow-line ice shelf
numerical solution                                           2
scaling up                         Ht = M − qx = M +           A(ρg)3 (hx )3 H 5       (∗)
next steps
5                     x
practicalities
omitted models
• equation (∗) matches equation (1) in the plane ﬂow case
Numerical
modelling                                              have I oversold diffusivity?
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &        • I have asserted that the default model for ice sheets, the SIA,
numerics
analogy
is a diffusion
ﬁnite differences
exact solutions
• recall the analogy:
shallow ice
sheets
solving the SIA           diffusion eqn       ↔           SIA
is it right?
application
ut = · (D u)         ↔    ht = M + · (D h)
the most basic
shallow assumption
D = D(x, y )       ↔    D = ΓH n+2 | h|n−1
ﬂowline Stokes
derivation of SIA

shelves and           • have I oversold this diffusivity analogy?
streams
shallow shelf             ◦ probably
approximation (SSA)
ﬂow-line ice shelf        ◦ . . . and I’ve acknowledged there are “issues”
numerical solution
scaling up            • but it turns out there is some robustness to my analogy
next steps                ◦ as follows on the next two slides
practicalities
omitted models
Numerical
modelling                                            have I oversold diffusivity? 2
Ed Bueler
DIFFUSIVE IDEA 1: rough beds have the effect of reducing
ice ﬂow: a
superﬁcial view       diffusivity
shallow ice
approximation (SIA)      • deﬁne the local bed topography by removing the local mean
heat analogy &
numerics
over some range λ ≈ 10 km:
analogy
ﬁnite differences
λ
exact solutions
˜
b(x, ξ) = b(x + ξ) −        b(x + ξ) dξ
shallow ice                                                    −λ
sheets
solving the SIA         • deﬁne this strange average of the local bed:
is it right?
application                                                              −1/n
the most basic                                                −(n+2)/n
shallow assumption
λ      ˜
b(x, ξ)
ﬂowline Stokes                      θ(x) =      1−                    dξ 
derivation of SIA
−λ       H(x)
shelves and
streams
shallow shelf           • using a multiple-scales analysis, Schoof [2003] says you will
approximation (SSA)
ﬂow-line ice shelf        get closer to solving the Stokes equations by making these
numerical solution
scaling up
two modiﬁcations of the SIA:
next steps                  ◦ smooth the bed,        (modelers want to do that anyway)
practicalities
omitted models
◦ but don’t lose track of the smoothed-away local bed roughness;
use it to reduce the diffusivity:

Dnew = θD      where 0 ≤ θ ≤ 1
Numerical
modelling                                           have I oversold diffusivity? 3
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
DIFFUSIVE IDEA 2: the large-scale effect of sliding in ice
approximation (SIA)
streams (addressed in next section), is diffusive
heat analogy &
numerics                • suppose that, for an ice stream modeled by the SSA equation
analogy
ﬁnite differences
exact solutions

shallow ice
2A−1/n H|ux |1/n−1 ux       − C|u|m−1 u = ρgHhx
sheets
x
solving the SIA
is it right?              we assume that the basal resistance term balances the
application
the most basic            driving stress:
shallow assumption
ﬂowline Stokes                               −C|u|m−1 u = ρgHhx
derivation of SIA

shelves and             • then the ice stream geometry evolves by a diffusion,
streams
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
ht = M +     · (D   h)     where D = C H (1/m)+1 | h|(1/m)−1
numerical solution
scaling up
• this “outer problem” is part of a matched asymptotic
next steps
practicalities            expansion that does a good job of tracking the grounding line
omitted models
in a marine ice sheet [Schoof, 2007]
Numerical
modelling                                                    Outline
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
numerics              1 ice ﬂow: an outsider’s superﬁcial view
analogy
ﬁnite differences
exact solutions

shallow ice
2 heat analogy & numerics
sheets
solving the SIA
is it right?
application
3 shallow ice sheets
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA     4 shelves and streams
shelves and
streams
shallow shelf
approximation (SSA)
5 next steps
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
Numerical
modelling            the shallow shelf approximation (SSA) stress balance
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
is a model which applies well to large ice shelves
approximation (SIA)
• . . . for parts away from grounding lines
heat analogy &
numerics
analogy
• . . . and away from calving fronts
ﬁnite differences
exact solutions

shallow ice
sheets
solving the SIA
is it right?
application
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA

shelves and
streams
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models                     edge of Ekström ice shelf, photo Hans Grobe, Polarstern expedition ANT-XX/2
Numerical
modelling                                                    SSA stress balance 2
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
ﬁnite differences
exact solutions
SSA also applies reasonably well to
shallow ice           ice streams
sheets
solving the SIA         • . . . with modest bed topography
is it right?
application             • . . . and weak bed strength
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA

shelves and
streams
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up
comment: energy conservation issues—both ice temperature and basal
next steps             melt rate—are a major part of ice stream ﬂow modeling [Raymond,
practicalities
omitted models
2000], but they are not addressed in these lectures
Numerical
modelling                                      what is, and is not, an ice stream?
Ed Bueler

ice ﬂow: a
superﬁcial view       • ice streams have fast ﬂow (50 to
shallow ice
approximation (SIA)     > 1000 m a−1 ) by sliding, with
heat analogy &            ◦ low bed resistance and
numerics
analogy                   ◦ a critical role of liquid water at
ﬁnite differences
exact solutions
bed [Clarke, 2005]
shallow ice           • “outlet glaciers” may have the
sheets
solving the SIA         same properties, but they have
is it right?
application             substantial vertical shear and
the most basic
shallow assumption      proportionally less sliding
ﬂowline Stokes
derivation of SIA     • and they have lower aspect ratio
shelves and
streams
than “true” ice streams
shallow shelf
approximation (SSA)
• and some outlet glaciers are really      Cross sections of Jakobshavns
Isbrae (a) and Whillans Ice Stream
ﬂow-line ice shelf
numerical solution
fast, e.g. Jakobshavn Isbrae             (b). Plotted without vertical
scaling up                                                       exaggeration in order to better
• thus: few simplifying assumptions        illustrate the difference between the
next steps                                                       two types. (Figure 1 in [Truffer and Echelmeyer,
practicalities          are possible for outlet glaciers         2003])
omitted models
• not clear you should use SSA for
outlet glaciers
Numerical
modelling         stream-to-shelf ﬂow line: notation again
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
ﬁnite differences
exact solutions

shallow ice
sheets
solving the SIA
is it right?
application
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA

shelves and
streams
shallow shelf
ﬁgure modiﬁed from [Schoof, 2007]
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
Numerical
modelling                                                                   the SSA equation
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
• we will only work in plane ﬂow case
numerics
analogy               • we will numerically solve only the stress balance equation
ﬁnite differences
exact solutions
which determines the velocity in an ice shelf or ice stream:
shallow ice
sheets
solving the SIA
is it right?
2A−1/n H|ux |1/n−1 ux             − C|u|m−1 u − ρgHhx = 0     (7)
x
application
the most basic
shallow assumption
ﬂowline Stokes            ◦ derived originally by Morland [1987], MacAyeal [1989]
derivation of SIA
◦ references: Weis and others [1999], Schoof [2006; 2007]
shelves and
streams
shallow shelf
• the blue term inside parentheses is the vertically-integrated
approximation (SSA)
ﬂow-line ice shelf
“longitudinal” or “membrane” stress
numerical solution
scaling up            • how to solve this equation numerically?
next steps
practicalities
• how to think about it?
omitted models
Numerical
modelling                                               the full monty, with a grounding line
Ed Bueler
here is a full ﬂow line context:
ice ﬂow: a
superﬁcial view
shallow ice
u(0) = (given)
approximation (SIA)                                                               
heat analogy &
numerics
2A−1/n H|ux |1/n−1 ux         − C|u|m−1 u − ρgHhx = 0  

x                             

analogy
ﬁnite differences                                                Ht = M − (uH)x      on 0 < x < xg
exact solutions
ρH ≥ ρw (−b) 



shallow ice                                                                       
sheets                                                                 h =H +b
solving the SIA
is it right?                                                      h, u, ux continuous at x = xg
application                                                                           
the most basic
shallow assumption
2A−1/n H|ux |1/n−1 ux           − ρg(1 − ρ/ρw )HHx = 0   

ﬂowline Stokes                                          x
derivation of SIA                                                                         on xg < x < xc
Ht = M − (uH)x 
shelves and                                                                           

streams                                                                 ρH < ρw (−b)
shallow shelf
approximation (SSA)
2A−1/n H|ux |1/n−1 ux − 2 ρ(1 − ρ/ρw )gH 2 = 0
1
ﬂow-line ice shelf
at x = xc
numerical solution                 (a condition describing movement of xc )
scaling up

next steps
practicalities          • this is the default model in the Marine Ice Sheet Model
omitted models
Intercomparison Project [MISMIP; Schoof and others, 2008]
• effectively an open problem: what is best numerical treatment of
grounding line movement in above model?
Numerical
modelling                       exact thickness and velocity for steady ice shelf
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &        • we need a more limited goal!
numerics
analogy               • ﬁrst goal: solve the steady state ice shelf in the isothermal,
ﬁnite differences
exact solutions          plane ﬂow, and Glen law case
shallow ice
sheets
• for this nice situation, there is a nice result: the thickness and
solving the SIA          velocity in the ice shelf can be completely determined5 from:
is it right?
application                    1   ice thickness Hg at the grounding line,
the most basic
shallow assumption             2   ice velocity ug at the grounding line,
ﬂowline Stokes
derivation of SIA
3   an integrable (e.g. constant) surface balance M(x)
shelves and           • we do this by-hand on the next slide
streams
shallow shelf
approximation (SSA)
• we will use the by-hand result to:
ﬂow-line ice shelf
numerical solution
◦ understand the SSA, and
scaling up                ◦ verify a numerical code
next steps
practicalities
omitted models

5 e.g.   [van der Veen, 1985]
Numerical
modelling                 exact thickness and velocity for steady ice shelf 2
Ed Bueler

ice ﬂow: a
1   using ﬂotation criterion h = (1 − ρ/ρw )H, and because of no bed
superﬁcial view           resistance, equation (7) says “derivative of something = 0”:
shallow ice
approximation (SIA)
1
heat analogy &                     2A−1/n H|ux |1/n−1 ux       −      ρg(1 − ρ/ρw )H 2       =0    (8)
numerics                                                   x        2                    x
analogy
ﬁnite differences
2   use the calving-front condition to integrate (8):
exact solutions

1
2A−1/n H|ux |1/n−1 ux −        ρg(1 − ρ/ρw )H 2 = 0
shallow ice
sheets                                                                                             (9)
solving the SIA
2
is it right?
application
3   the steady (Ht = 0) mass continuity equation can be integrated;
the most basic
shallow assumption
here M = M0 =constant but some other M(x) are allowed:
ﬂowline Stokes
derivation of SIA                               uH = M0 (x − xg ) + ug Hg ,                       (10)
shelves and
streams               4   multiply (9) by u/H, replace uH from (10), assume positive strain
shallow shelf
approximation (SSA)
rate (ux > 0), compute nth power
ﬂow-line ice shelf
numerical solution
5   get u n ux = Cs (M0 (x − xg ) + ug Hg )n , where Cs is known
scaling up
6   integrate the last result to ﬁnd u(x), then return to (10) to ﬁnd H(x):
next steps
Cs
n+1
u(x)n+1 = ug +         (M0 (x − xg ) + ug Hg )n+1 − (ug Hg )n+1 ,
practicalities
omitted models
M0
M0 (x − xg ) + ug Hg
H(x) =
u(x)
Numerical
modelling               exact thickness and velocity for steady ice shelf 3
Ed Bueler

ice ﬂow: a
superﬁcial view       • Matlab/Octave code exactshelf.m uses shelf length L = 200 km,
shallow ice
approximation (SIA)     and Hg = 500, ug = 50 m a−1
heat analogy &
numerics
• computes this shelf geometry and velocity:
analogy
ﬁnite differences
exact solutions

shallow ice
sheets
solving the SIA
is it right?
application
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA

shelves and
streams                                                      exact velocity; u = 303.8539 at x=L
shallow shelf                                      350
approximation (SSA)
300
velocity (m a-1)

ﬂow-line ice shelf
numerical solution                                 250
scaling up                                         200
next steps                                         150
practicalities                                     100
omitted models
50
0   50             100              150   200
x (km)
Numerical
modelling                                                  numerically solving the SSA
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &        • the stress balance is a two-point boundary value problem
numerics
analogy
(BVP) for the velocity u(x):
ﬁnite differences

2A−1/n H|ux |1/n−1 ux       − C|u|m−1 u − ρgHhx = 0,
exact solutions

shallow ice                                                    x
sheets
solving the SIA           left b.c.:                            u(xg ) = ug ,
is it right?
application
right b.c.:     2A−1/n H|ux |1/n−1 ux − 1 ρ(1 − ρ/ρw )gH 2 = 0 at x = xc
2
the most basic
shallow assumption
ﬂowline Stokes
• here we will prescribe the ice geometry, thus both the
derivation of SIA
thickness H(x) and the driving stress ρgHhx , and we will ﬁnd
shelves and
streams                 the velocity u(x)
shallow shelf
approximation (SSA)   • I’ll describe the numerical method for a ice shelf or ice stream
ﬂow-line ice shelf
numerical solution
scaling up
• . . . but then I’ll give a code only for an ice shelf, so C = 0 and
next steps              h = (1 − ρ/ρw )H
practicalities
omitted models
Numerical
modelling                    numerically solving the SSA stress balance 2
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
• like most such nonlinear equations, iteration is needed to
approximation (SIA)
solve this one
heat analogy &
numerics              • red term ν = A−1/n |ux |1/n−1 is the “effective viscosity”:
¯
analogy
ﬁnite differences
exact solutions

shallow ice
2A−1/n |ux |1/n−1 Hux               − C|u|m−1 u − ρgHhx = 0
sheets
x
solving the SIA
is it right?          • let
application
the most basic                         W (ux ) = 2¯H = 2A−1/n |ux |1/n−1 H
ν
shallow assumption
ﬂowline Stokes
derivation of SIA     • idea: use old effective viscosity to get new velocity solution,
shelves and             and repeat until the solution is not changing too much and/or
streams
shallow shelf           the differential equation is nearly-satisﬁed
approximation (SSA)
ﬂow-line ice shelf
• idea as algorithm: from u (k −1) , ﬁnd next values u (k ) by
numerical solution
scaling up              solving the linear problem
next steps
practicalities
(k −1)     (k )
omitted models                 W (ux          )ux          − C|u (k −1) |m−1 u (k ) − ρgHhx = 0
x
Numerical
modelling                      numerically solving the SSA stress balance 3
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)
where do you get an initial guess u (0) (x) for the velocity?
heat analogy &
numerics
analogy
• for ﬂoating ice, a velocity comes from assuming a uniform
ﬁnite differences
exact solutions           strain rate:
shallow ice
sheets
u (0) = γ(x − xg ) + ug
solving the SIA
is it right?
application
the most basic               ◦ for example: γ is the value of ux found from the calving front
shallow assumption
ﬂowline Stokes
stress imbalance and ug is the known velocity at the grounding
derivation of SIA
line
shelves and
streams                 • for grounded ice, a velocity comes from assuming the ice is
shallow shelf
approximation (SSA)       held by basal resistance only:
ﬂow-line ice shelf
numerical solution
1/m
scaling up
u (0) = −C −1 ρgHhx
next steps
practicalities
omitted models
Numerical
modelling                                   solving the “inner” linear problem
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice           • reset x interval to be 0 < x < L, instead of xg < x < xc
approximation (SIA)

heat analogy &        • abstract the problem to a linear BVP:
numerics
analogy
ﬁnite differences
exact solutions
(W (x) ux )x − α(x) u = β(x)
shallow ice
sheets                  with boundary conditions
solving the SIA
is it right?
application                                u(0) = V ,      ux (L) = γ
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA
• W (x), α(x), β(x) are all known functions in this abstract
shelves and             problem
streams
shallow shelf
approximation (SSA)
• for the nonlinear SSA equation, both W (x) and α(x) will
ﬂow-line ice shelf      come from the previous iteration
numerical solution
scaling up
• W (x) is needed on the staggered grid, for O(∆x 2 ) accuracy,
next steps
practicalities
as with time-dependent diffusion problem earlier
omitted models
• α(x), β(x) are needed on regular grid
Numerical
modelling                                    solving the “inner” linear problem 2
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)
• indices j = 1, 2, . . . , J + 1
heat analogy &
numerics
analogy
• equal-spaced grid x1 , x2 , . . . , xJ+1 , where x1 = 0 and xJ+1 = L
ﬁnite differences
exact solutions
• an approximation of the abstracted problem is:
shallow ice
sheets                          Wj+1/2 (uj+1 − uj ) − Wj−1/2 (uj − uj−1 )         ∗
solving the SIA                                                           − αj uj = βj
is it right?                                      ∆x 2
application
the most basic
shallow assumption    • u1 = V given
ﬂowline Stokes
derivation of SIA     • for right-hand boundary condition “ux (L) = γ”:
shelves and
streams                    ◦ introduce notional point xJ+2
shallow shelf              ◦
approximation (SSA)
uJ+2 − uJ
ﬂow-line ice shelf
=γ
numerical solution                                       2∆x
scaling up

next steps
◦ using equation ∗ in j = J + 1 case, eliminate uJ+2 variable
practicalities                “by-hand” [Morton and Mayers, 2005]
omitted models
Numerical
modelling                                              solving the “inner” linear problem 3
Ed Bueler

ice ﬂow: a
superﬁcial view       •   so abstract linear system has matrix-vector form “ Ax = b ”:
shallow ice
approximation (SIA)
  V 
1                                                           u1
                                                      
heat analogy &                                                                                                2
numerics                          W3/2      A22        W5/2                                 u2  β2 ∆x 
 u3  β3 ∆x 2 
analogy
          W5/2       A33                                                 
                                                                           
ﬁnite differences
 . = . 
                                                                
exact solutions
                      ..             ..               
 .   . 

                         .             .              
    .   .         
shallow ice
βJ ∆x 2 
                    WJ−1/2           AJJ   WJ+1/2         u  
J
sheets
solving the SIA
AJ+1,J    AJ+1,J+1         uJ+1      bJ+1
is it right?
application
•   with diagonal entries (j = 2, 3, . . . , J)
the most basic
shallow assumption
ﬂowline Stokes                                                                                2
derivation of SIA
Ajj = −(Wj−1/2 + Wj+1/2 + αj ∆x )
shelves and
streams               •   with special cases in last equation:
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
AJ+1,J = WJ+1/2 + WJ+3/2
numerical solution
scaling up                                                                                         2
AJ+1,J+1 = −(WJ+1/2 + WJ+3/2 + αJ+1 ∆x )
next steps
practicalities                                                                                2
omitted models
bJ+1 = −2γ∆xWJ+3/2 + βJ+1 ∆x
•   this is a tridiagonal system
Numerical
modelling                                         solving the “inner” linear problem 4
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)   function u = flowline(L,J,gamma,W,alpha,beta,V0)
heat analogy &
numerics              dx = L / J;
analogy               rhs = dx^2 * beta(:);
ﬁnite differences
rhs(1) = V0;
exact solutions
rhs(J+1) = rhs(J+1) - 2 * gamma * dx * W(J+1);
shallow ice
sheets
A = sparse(J+1,J+1);
solving the SIA
is it right?
A(1,1) = 1.0;
application           for j=2:J
the most basic          A(j,j-1:j+1) = [ W(j-1), -(W(j-1) + W(j) + alpha(j) * dx^2), W(j) ];
shallow assumption
end
ﬂowline Stokes
A(J+1,J) = W(J) + W(J+1);
derivation of SIA
A(J+1,J+1) = - (W(J) + W(J+1) + alpha(J+1) * dx^2);
shelves and
streams
shallow shelf
scale = full(max(abs(A),[],2));
approximation (SSA)   for j=1:J+1, A(j,:) = A(j,:) ./ scale(j);               end
ﬂow-line ice shelf    rhs = rhs ./ scale;
numerical solution
scaling up
u = A \ rhs;
next steps
practicalities                       http://www.dms.uaf.edu/~bueler/karthaus/mfiles/flowline.m
omitted models
Numerical
modelling                                             solving the “inner” linear problem 5
Ed Bueler

ice ﬂow: a
• before proceeding to solve nonlinear SSA problem, we can
superﬁcial view
shallow ice
test the “abstracted” code flowline.m
approximation (SIA)
• . . . by “manufacturing” solutions as in code
heat analogy &
numerics                testflowline.m (not shown)
analogy
ﬁnite differences
• results: works pretty well! converges at optimal rate
exact solutions
O(∆x 2 )
shallow ice
sheets
-3
solving the SIA                                10
is it right?
application
the most basic
shallow assumption                             10-4
ﬂowline Stokes
maximum error

derivation of SIA

shelves and                                    10-5
streams
shallow shelf
approximation (SSA)
ﬂow-line ice shelf                             10-6
numerical solution
scaling up

next steps
10-7
practicalities
convergence rate O(dx2.00065)
omitted models

10-8
10-4   10-3               10-2        10-1
dx
Numerical
modelling                                                                numerical solution to SSA
Ed Bueler

ice ﬂow: a
superﬁcial view       function [u,u0] = ssaflowline(p,J,H,b,ug,initchoice)
shallow ice
approximation (SIA)   if nargin ~= 6, error(’exactly 6 input arguments required’), end

heat analogy &        dx = p.L / J;
numerics              x = (0:dx:p.L)’;
analogy               xstag = (dx/2:dx:p.L+dx/2)’;
ﬁnite differences
exact solutions       alpha = p.C * ones(size(x));
h = H + b;
shallow ice           hx = regslope(dx,h);
sheets                beta = p.rho * p.g * H .* hx;
solving the SIA       gamma = ( 0.25 * p.A^(1/p.n) * (1 - p.rho/p.rhow) *...
is it right?                     p.rho * p.g * H(end) )^p.n;
application
u0 = ssainit(p,x,beta,gamma,initchoice);
the most basic
u = u0;
shallow assumption
ﬂowline Stokes
Hstag = stagav(H);
derivation of SIA     tol = 1.0e-14;
maxdiff = Inf;
shelves and
W = zeros(J+1,1);
streams
while maxdiff > tol
shallow shelf
uxstag = stagslope(dx,u);
approximation (SSA)
W(1:J) = 2 * p.A^(-1/p.n) * Hstag .* (abs(uxstag)).^((1/p.n)-1);
ﬂow-line ice shelf
W(J+1) = W(J);
numerical solution
scaling up              unew = flowline(p.L,J,gamma,W,alpha,beta,ug);
maxdiff = max(abs(unew-u));
next steps              u = unew;
practicalities        end
omitted models

http://www.dms.uaf.edu/~bueler/karthaus/mfiles/ssaflowline.m
Numerical
modelling           numerical thickness and velocity for steady ice shelf
Ed Bueler

ice ﬂow: a            • ﬁrst: let’s numerically solve the problem for which we know
superﬁcial view
shallow ice
approximation (SIA)
heat analogy &        • here use very coarse 10 km grid
numerics
analogy
-1
ﬁnite differences                                         case with dx = 10 km has (max error) = 3.044 m a
exact solutions                                 350
initial guess
shallow ice                                                                                 exact solution
sheets                                                                                   numerical solution
300
solving the SIA
is it right?
application
250
the most basic
shallow assumption
ﬂowline Stokes
velocity (m a )
-1

derivation of SIA
200

shelves and
streams                                         150
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
100
numerical solution
scaling up

next steps                                      50
practicalities
omitted models
0
0      50000             100000            150000        200000
x (km)
Numerical
modelling                                                                  did we make any mistakes?
Ed Bueler

ice ﬂow: a
• this does convergence analysis of ssaflowline.m:
superﬁcial view         J = [50 100 200 400 800 1600 3200];   dxkm = 200.0 ./ J;
shallow ice
approximation (SIA)
for j=1:length(J)
[av,maxerr(j)] = testshelf(J(j));   % calls ssaflowline.m
heat analogy &
numerics
end
analogy                 loglog(dxkm,maxerr,’0-’,’markersize’,16)
ﬁnite differences
exact solutions       • goal is to see velocity error shrink by rate built-into our ﬁnite
shallow ice             differences, namely at rate O(∆x 2 )
sheets
solving the SIA       • by this standard, we made no mistakes!
is it right?
application                                                     0
10
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA
10-1
shelves and
maximum error (m a-1)

streams
shallow shelf
approximation (SSA)
10-2
ﬂow-line ice shelf
numerical solution
scaling up

next steps                                                 10-3
practicalities                                                             convergence rate O(dx-2.00846)
omitted models

10-4
10-1               100                  101
dx (km)
Numerical
modelling                                                    realistic ice shelf modeling
Ed Bueler

ice ﬂow: a
• ﬂow lines (1D) are never very realistic . . .
superﬁcial view
shallow ice
• “diagnostic” (= ﬁxed geometry) ice shelf modeling in two horizontal
approximation (SIA)
variables (2D) has been quite successful
heat analogy &
numerics              • solve an elliptic PDE boundary value problem
analogy
ﬁnite differences     • velocity measurements are available for validation
exact solutions            ◦ Ross ice shelf example based on [RIGGS; Bentley, 1974] data below: model (red
shallow ice                  arrows) vs observations (black arrows)
sheets                     ◦ this is PISM but many models can do it [MacAyeal and others, 1996]
solving the SIA
is it right?               ◦ something of a mystery why regular grid methods do well!
application
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA

shelves and
streams
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
Numerical
modelling                             numerical solution of SSA: a summary
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
• SSA and Blatter [1995] stress balance equations each
approximation (SIA)
◦ determine horizontal velocity from geometry and b.c.s
heat analogy &
numerics                  ◦ are nonlinear problems: require iterative linearization
analogy
◦ each “inner” linear problem has a sparse pattern
ﬁnite differences
exact solutions                ∗ tridiagonal in SSA ﬂow-line (1D) case here
shallow ice                    ∗ much less structured in SSA 2D, Blatter 2D, and Blatter 3D cases
sheets
solving the SIA       • the Stokes stress balance is much harder because
is it right?
application             incompressibility must be handled explicitly
the most basic
shallow assumption    • a wise numerical modeler would give the sparse linear
ﬂowline Stokes
derivation of SIA       problem to a numerical linear algebra software package:
shelves and               ◦ M ATLAB/O CTAVE, LAPACK, PETSc, Trilinos, . . .
streams
shallow shelf             ◦ generally: modularize your code and test the parts separately!
approximation (SSA)
ﬂow-line ice shelf    • comment on mass continuity:
numerical solution
scaling up                                                              ¯
◦ the mass continuity PDE “Ht = M − · (uH)” is not very
next steps                                             ¯
diffusive if the source of u is SSA or Blatter,
practicalities
omitted models            ◦ and there is not much theory to guide how to solve such a
generic transport problem
Numerical
modelling           scaling-up the numerical solution of stress balances
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
• all stress balance equations used for ice ﬂow modeling are
approximation (SIA)
nonlinear equations for velocity u in the form
heat analogy &
numerics
analogy
ﬁnite differences
F(u) = 0
exact solutions

shallow ice
sheets
solving the SIA
• the “residual” F involves various partial derivatives, nonlinear
is it right?            operations, and additional ﬁelds (e.g. thickness)
application
the most basic
shallow assumption
• SSA and Blatter stress balance equations can be written in
ﬂowline Stokes
derivation of SIA
form
shelves and
A(u) u = g
streams
shallow shelf
approximation (SSA)       ◦ where A is a nonlinear function of u,
ﬂow-line ice shelf
numerical solution
◦ and A(u) acts linearly on u
scaling up
• we have already treated the SSA in the second way
next steps
practicalities        • both forms are nonlinear, elliptic PDEs
omitted models

• Newton [1669] knew a way of solving nonlinear equations!
Numerical
modelling         scaling-up the numerical solution of stress balances 2
Ed Bueler

ice ﬂow: a            • the iteration scheme we have used for SSA is “Picard”:
superﬁcial view
shallow ice
approximation (SIA)
A(u(k −1) ) u(k ) = g
heat analogy &
numerics
analogy                   ◦ this iteration converges linearly, if it converges:
u(k +1) − u(k ) ≤ γ u(k ) − u(k −1) , where 0 < γ < 1
ﬁnite differences
exact solutions

shallow ice               ◦ it can be robust, but it is not fast
sheets
solving the SIA       • Newton’s method (e.g. [Press and others, 1992]) also solves
is it right?
application             F(u) = 0 iteratively, but by linearization:
the most basic
shallow assumption
ﬂowline Stokes                              J(u(k −1) ) w = −F(u(k −1) ),
derivation of SIA

shelves and
streams
u(k ) = u(k −1) + w
shallow shelf
approximation (SSA)
ﬂow-line ice shelf      where
numerical solution
∂F(u)i
scaling up
J(u)ij =               is the Jacobian of F
next steps                                      ∂uj
practicalities
omitted models
◦ this iteration converges quadratically, if it converges:
u(k +1) − u(k ) ≤ C u(k ) − u(k −1) 2
◦ much faster!
Numerical
modelling         scaling-up the numerical solution of stress balances 3
Ed Bueler

ice ﬂow: a            • big computers still take a long time to solve SSA equations
superﬁcial view
shallow ice
when grids are ∼ 2 km for Greenland
approximation (SIA)
◦ Antarctica is 10× the area, and that much worse!
heat analogy &
numerics
analogy
• there are various numerical paradigms for discretizing (ﬁnite
ﬁnite differences
exact solutions
difference, ﬁnite element, spectral, . . . ), but choosing among
shallow ice
these is not the major scalability issue
sheets
solving the SIA
• effective scalability of 3D stress balances (e.g. Blatter &
is it right?
application
Stokes) at high spatial resolution, on large parallel machines
the most basic          requires:
shallow assumption
ﬂowline Stokes            ◦ the “inner” linear solve must be iterative–not by Gaussian
derivation of SIA
elimination—and must be preconditioned
shelves and
streams                   ◦ the “outer” nonlinear iteration must offer faster-than-linear
shallow shelf
approximation (SSA)         convergence
ﬂow-line ice shelf
◦ because the exact Jacobian is too hard to evaluate—certainly
numerical solution
scaling up                  so with coupling to other physics (e.g. thermo-)—the method
next steps                  must allow approximated Jacobians
practicalities
omitted models        • this is the promise of “preconditioned matrix-free
Newton-Krylov” (e.g. [Knoll and Keyes, 2004]) as a paradigm
for big, nasty nonlinear problems
Numerical
modelling            scaling-up the numerical solution of stress balances 4
Ed Bueler

ice ﬂow: a            • the PETSc [Balay and others,                                   Level of
superﬁcial view                                                                       Abstraction              Application Codes
shallow ice
2010] library has evolved into a
approximation (SIA)
Newton-Krylov toolkit →                                                                                               TS
heat analogy &                                                                                                                         (Time Stepping)
numerics              • my online materials turn                                                        SNES
(Nonlinear Equations Solvers)
analogy
ﬁnite differences
ssaflowline.m into a C                                                                                         KSP
PC              (Krylov Subspace Methods)
exact solutions         PETSc example                                                               (Preconditioners)

shallow ice             ssaflowline.c:
sheets                                                                                          Matrices                Vectors          Index Sets
solving the SIA
◦ parallel over an arbitrary
is it right?
BLAS               MPI
application                 number of processors
the most basic
shallow assumption        ◦ uses Newton method (or                             104
M=10
ﬂowline Stokes
Picard)                                                                                                                M=100
derivation of SIA                                                              102                                                                 M=103

shelves and
◦ try matrix-free, different                                                                                             M=104
M=105
100
streams                     preconditioners, . . .
shallow shelf
approximation (SSA)       ◦ result from                                       10-2

residual norm
ﬂow-line ice shelf          ssaflowline.c is at                               10-4
numerical solution
scaling up
lower-right:
10-6
next steps
∗ evidence of super-linear
practicalities
convergence from Newton                     10-8

omitted models                    method
∗ on coarsest grids, 10 or                    10-10

100 points, get rapid
10-12
convergence to 10−13 of                             0             1          2           3              4           5    6
Newton iteration
the initial residual norm
Numerical
modelling                                                    Outline
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
numerics              1 ice ﬂow: an outsider’s superﬁcial view
analogy
ﬁnite differences
exact solutions

shallow ice
2 heat analogy & numerics
sheets
solving the SIA
is it right?
application
3 shallow ice sheets
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA     4 shelves and streams
shelves and
streams
shallow shelf
approximation (SSA)
5 next steps
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
Numerical
modelling                                    what technical skills are needed for
Ed Bueler                                        numerical ice sheet modeling?
ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &        you need:
numerics
analogy
• comfort in a technical computing environment, usually Unix,
ﬁnite differences
exact solutions           in which you need to know:
shallow ice
sheets
◦   an editor,
solving the SIA             ◦   a compiled language (Fortran or C),
is it right?
application
◦   a scripting/prototyping language (Matlab, Python, etc.), and
the most basic
shallow assumption
◦   a version control system (Subversion, git, etc.)
ﬂowline Stokes
derivation of SIA       • willingness to read math, numerical analysis, computer
shelves and               science books
streams
shallow shelf
approximation (SSA)
• . . . and willingness to ignor some of the advice found there
ﬂow-line ice shelf
numerical solution      • know some tools for NetCDF ﬁles
scaling up

next steps
• get exposed to parallel computing
practicalities
omitted models
• but, at the end of the day: physics
Numerical
modelling                                 what technical skills are needed? 2
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)
• an important skill is to not re-invent the wheel
heat analogy &
numerics
analogy
• never re-invent the wheel for basic numerics:
ﬁnite differences
exact solutions
◦ M ATLAB, Comsol, PETSc, libmesh, Elmer, Trilinos, Triangle,
shallow ice
etc. handle fundamental numerical linear algebra, mesh
sheets                      generation, ﬁnite element assembly and solve, etc. tasks; don’t
solving the SIA
is it right?                even try to compete!
application
the most basic
◦ . . . except to write throw-away codes to help you understand
shallow assumption          numerical ideas
ﬂowline Stokes
derivation of SIA     • sometimes there is no need to re-invent ice sheet modeling:
shelves and
streams                   ◦ open source SIA-based comprehensive models: GLIMMER,
shallow shelf
approximation (SSA)
SICOPOLIS
ﬂow-line ice shelf        ◦ open source hybrid/higher-order/Stokes models: PISM,
numerical solution
scaling up                  Elmer-ice, CISM
next steps
practicalities
• glacier and ice sheet modeling is young! there is much to do!
omitted models
Numerical
modelling                    omitted models 1: thermomechanical coupling
Ed Bueler

1e-23

ice ﬂow: a
superﬁcial view
• the softness of ice “A” is not
shallow ice
approximation (SIA)
constant!                                          1e-24

heat analogy &            ◦ A = A(T ) varies by more than

A(T)
1e-25
numerics
analogy
103 in the 50 ◦ C ice temperature
Paterson-Budd form for A(T)
ﬁnite differences           range relevant to Antarctic ice                1e-26

exact solutions
◦ dissipation of ﬂow energy is
shallow ice
sheets                      signiﬁcant in conservation of                  1e-27
-50   -40   -30        -20
T (degrees C)
-10      0

solving the SIA
is it right?
energy equation
application               ◦ . . . and there is a “new” ﬂuid
the most basic
shallow assumption          instability from feedback between
ﬂowline Stokes
derivation of SIA
the above two items (lower ﬁgure;
shelves and
EISMINT II [Payne and others, 2000])
streams
shallow shelf         • ice temperature is part of the “long
approximation (SSA)
ﬂow-line ice shelf      memory” of the ice sheets for past
numerical solution
scaling up
climate
next steps            • conservation of energy is just as
practicalities
omitted models          important for fast scale dynamics
(involving sliding and hydrology) as
it is for “long memory” questions
Numerical
modelling                         omitted models 2: surface mass balance
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
ﬁnite differences
• ice sheet models need modeled surface mass balance
exact solutions

shallow ice                ◦ over paleoglacial time scales
sheets
solving the SIA
◦ and when modeling response to future changes
is it right?
application
• models include:
the most basic
shallow assumption
◦ positive degree-day (PDD) and other index models
ﬂowline Stokes            ◦ energy-balance models
derivation of SIA

shelves and           • very signiﬁcant source of uncertainty for Greenland ice sheet
streams
shallow shelf           dynamics!
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
Numerical
modelling                        omitted models 3: solid Earth deformation
Ed Bueler

ice ﬂow: a                                1
superﬁcial view       • ice has almost    3   of the density of the viscous hot rock in the
shallow ice
approximation (SIA)     Earth’s mantle,
heat analogy &
numerics              • so 1000 m of ice will depress the Earth’s crust almost 300 m
analogy
ﬁnite differences
if allowed enough time
exact solutions
• this changes bed topography and thus ice ﬂow, so Earth
shallow ice
sheets                  deformation is important to modeling ice ﬂow
solving the SIA
application
the most basic            ◦ [Peltier, 1998] = survey of observations and processes
shallow assumption
ﬂowline Stokes            ◦ [Greve, 2001] = comparison of practical models
derivation of SIA

shelves and
streams
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
Numerical
modelling                            omitted models 4: numerical Stokes equations
Ed Bueler

ice ﬂow: a
superﬁcial view             • the Stokes model itselfcan be solved numerically
shallow ice
approximation (SIA)             ◦ no shallow assumptions!
heat analogy &                  ◦ but many more equations and unknowns than SIA, SSA,
numerics
analogy                           hybrids, Blatter, . . .
ﬁnite differences
exact solutions             • requires explicit accounting for
shallow ice
sheets
◦ incompressibility = a constraint on the ﬂow
solving the SIA                 ◦ pressure = a Lagrange multiplier for that constraint
is it right?
application                 • very tough scalability issues: can you afford the loss of
the most basic
shallow assumption             resolution? loss of long-time runs?
ﬂowline Stokes
derivation of SIA           • all of the success so far is at a smaller scale (e.g. below)
shelves and
streams
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
numerical solution
Figure 7 in [Maxwell and others, 2008]:
scaling up
Athabasca Glacier: (a) modeled velocity
next steps
contour lines (m a−1 ); (b) contour lines
practicalities
derived from measurements [Raymond,
omitted models        1971]
Numerical
modelling                        omitted models 5: “higher-order” schemes
Ed Bueler

ice ﬂow: a
superﬁcial view       • both the SIA and the SSA are derived by small-parameter
shallow ice
approximation (SIA)
arguments from the Stokes equations, so . . .
heat analogy &
numerics              • is there a common shallow antecedent model?
analogy
ﬁnite differences     • Schoof and Hindmarsh [2009] answer:
exact solutions

shallow ice
◦ yes, the Blatter [1995] model is intermediate between the
sheets                      Stokes stress balance and both the SIA and SSA
solving the SIA
is it right?
application                                          Stokes
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA
other higher-order?
shelves and
streams
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
Numerical
modelling               POSTSCRIPT on free boundary value problems
Ed Bueler

ice ﬂow: a
superﬁcial view       • ice sheet/shelf modeling means free boundary problems
shallow ice
approximation (SIA)
• Hutter [1999] identiﬁes some below
heat analogy &
numerics
analogy
ﬁnite differences
exact solutions

shallow ice
sheets
solving the SIA
is it right?
application
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA

shelves and
streams
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
Numerical
modelling                                                  what is a “free boundary”?
Ed Bueler
• a free boundary for a PDE is an unknown location at which
ice ﬂow: a
superﬁcial view             there is a boundary condition
shallow ice
approximation (SIA)            ◦ the location of the free boundary must be found at the same
heat analogy &                   time as one solves the PDE problem
numerics
analogy
◦ in addition to the information which would already be present at
ﬁnite differences                a ﬁxed-location boundary, there must be enough additional
exact solutions
information at a free boundary to determine its location
shallow ice
sheets                         ◦ all free boundary problems are nonlinear, regardless of the
solving the SIA
is it right?
linearity of the PDE
application               • classic example: consider an elastic membrane attached to a
the most basic
shallow assumption
ﬂowline Stokes
wire frame and stretched over an obstacle:
derivation of SIA

shelves and
streams               constraint:
shallow shelf
approximation (SSA)                  u≥ψ
ﬂow-line ice shelf
numerical solution    in locations where u > ψ, solve:
scaling up
2
next steps                                u=0
practicalities
omitted models
where is the free boundary, and
what facts about u are true there?
Numerical
modelling                  free boundary value problem 1: polythermal ice
Ed Bueler

ice ﬂow: a            • by volume, majority of ice sheet is cold (T < 0 ◦ C)
superﬁcial view
shallow ice           • . . . but there is some ice which is temperate, where T = 0 ◦ C and
approximation (SIA)

heat analogy &
there is a positive liquid fraction within the ice matrix
numerics
analogy
• . . . so ice sheets are polythermal
ﬁnite differences
exact solutions
• boundary between cold and temperate ice is “CTS” (=
shallow ice             cold-temperate transition surface):
sheets
solving the SIA            ◦ must be found, as free boundary, when solving conservation of
is it right?
application
energy equation (“Stefan probem”)
the most basic
shallow assumption
◦ can be tracked explicitly [Greve, 1997]
ﬂowline Stokes             ◦ or treated as a level surface of the enthalpy variable
derivation of SIA
[Aschwanden and Blatter, 2009]
shelves and
streams
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models

temperate       cold
Numerical
modelling               free boundary value problem 2: for ice streaming
Ed Bueler

ice ﬂow: a
superﬁcial view
• Schoof’s [2006] insight, for diagnostic case
shallow ice
approximation (SIA)
well-posed free boundary problem
heat analogy &
SSA + (plastic till) =
numerics
analogy
for location and velocity of sliding
ﬁnite differences
exact solutions
• “plastic till” means the basal strength (resistance) is given by
shallow ice
sheets                  a yield stress τc :       τb = τc vb /|vb |
solving the SIA
is it right?          • Schoof’s scheme is a whole ice sheet form of MacAyeal’s
application
the most basic          [1989] individual ice stream models
shallow assumption
ﬂowline Stokes
derivation of SIA

shelves and
streams
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
Numerical
modelling          free boundary problem 3: for grounded ice sheet margin
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)
• side-by-side comparison, classical elastic membrane
heat analogy &              problem versus steady ice sheet problem
numerics
analogy
ﬁnite differences
constraint:                 constraint:
exact solutions

shallow ice                         u≥ψ                   h≥b       ⇐⇒        H≥0
sheets
solving the SIA
is it right?          where u > ψ, solve:         where h > b, solve steady SIA:
application
the most basic                      2
shallow assumption                      u=0            0=M+       · ΓH n+2 | h|n−1 h
ﬂowline Stokes
derivation of SIA

shelves and
streams
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
Numerical
modelling                    for SIA model, an illuminating transformation
Ed Bueler

ice ﬂow: a
superﬁcial view       • there are degenerate diffusions in other ﬁelds
shallow ice
approximation (SIA)
• (e.g. porous media where there can be wet and dry parts)
heat analogy &
numerics              • . . . whence transformations useful to ice sheet modeling
analogy
ﬁnite differences
exact solutions
• speciﬁcally one by Raviart [1970]           (n = 3 case):     η = H 8/3
shallow ice           • the ﬂat bed SIA simpliﬁes:
sheets
solving the SIA
is it right?
application
Ht = M +             · ΓH 5 | H|2     H
the most basic
shallow assumption
becomes                                 ˜
ﬂowline Stokes
derivation of SIA
→          η 3/8       =M+        · (Γ| η|2    η)
t
shelves and
streams
shallow shelf
approximation (SSA)
◦ new diffusivity (red) depends on | η| only, not η
ﬂow-line ice shelf        ◦ the margin shape is not singular in new variable!
numerical solution
scaling up
∗   η is globally continuous
next steps
∗ z = η(t, x, y ) is tangent to bed at margin
practicalities
omitted models
◦ numerical methods applied to the new equation make errors
which are as expected for “good” free boundary problems
[Bueler and others, 2005]
Numerical
modelling                                 free boundary problems: so what?
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
• why does it matter that many glaciological problems are free
ﬁnite differences
exact solutions
boundary problems in their mathematical form?
shallow ice           • the location of the free boundary may be the glaciological
sheets
solving the SIA         question
is it right?
application           • free boundaries are always locations of loss of smoothness
the most basic
shallow assumption      relative to ﬁxed boundary solutions
ﬂowline Stokes
derivation of SIA
◦ numerical error may be dominated by errors near these free
shelves and                 boundaries,
streams
shallow shelf
◦ hard to know whether model results at free boundaries are
approximation (SSA)
poor because of numerical problems or because of missing
ﬂow-line ice shelf
numerical solution          physical processes
scaling up

next steps
practicalities
omitted models
Numerical
modelling                                       The End
Ed Bueler

ice ﬂow: a
superﬁcial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
ﬁnite differences
exact solutions

shallow ice
sheets
solving the SIA
is it right?
application
the most basic
shallow assumption
ﬂowline Stokes
derivation of SIA

shelves and
streams
shallow shelf
approximation (SSA)
ﬂow-line ice shelf
numerical solution
scaling up
Thanks for bearing with me!
next steps
practicalities
omitted models

```
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
 views: 12 posted: 7/23/2011 language: Italian pages: 125