Docstoc

Numerical modelling of ice sheets and ice shelves

Document Sample
Numerical modelling of ice sheets and ice shelves Powered By Docstoc
					    Numerical
    modelling

    Ed Bueler

ice flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
numerics
                           Numerical modelling
analogy
finite 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
flowline Stokes
                               and Geophysical Institute
derivation of SIA
                             University of Alaska, Fairbanks
shelves and
streams
shallow shelf
approximation (SSA)            September 2010
flow-line ice shelf
numerical solution         Karthaus Summer School
scaling up

next steps
practicalities
omitted models
    Numerical
    modelling                                                                                                    notation
    Ed Bueler

ice flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
finite differences
exact solutions

shallow ice                                               figure modified 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
flowline 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)
flow-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
                         •   please ask about notation!       . . . there are no stupid questions

                      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 flow: a
superficial view
shallow ice
approximation (SIA)
                      slogans:
heat analogy &          • my focus is on approximating ice flow
numerics
analogy
finite 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 flowline (1D)
shallow assumption
flowline Stokes          • mass continuity & surface/base kinematical equations
derivation of SIA

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

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

ice flow: a
superficial view
shallow ice
approximation (SIA)

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

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

ice flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
finite 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
flowline Stokes        • map-plane mass continuity equation (5)
derivation of SIA

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

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

ice flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
                      • these lectures are structured around
finite 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
flowline 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/
flow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
    Numerical
    modelling                                                    Outline
    Ed Bueler

ice flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
numerics              1 ice flow: an outsider’s superficial view
analogy
finite 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
flowline Stokes
derivation of SIA     4 shelves and streams
shelves and
streams
shallow shelf
approximation (SSA)
                      5 next steps
flow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
    Numerical
    modelling                                                               my first goal
    Ed Bueler

ice flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
finite 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 specific ice flow problem)
the most basic
shallow assumption
flowline Stokes
derivation of SIA
                      • to get to that goal: I will very quickly recall the continuum
shelves and
streams                 mechanics of ice flow
shallow shelf
approximation (SSA)   • . . . from a naive outside view
flow-line ice shelf
numerical solution
scaling up

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

ice flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &        • we describe fluids primarily by a velocity field u(t, x, y , z)
numerics
analogy
finite differences
                      • if the ice fluid were much faster-moving than it actually is1 ,
exact solutions              and if it were linearly-viscous like liquid water, then ice flow
shallow ice
sheets
                             modeling would be more familiar to other the
solving the SIA
is it right?
                             climate-modeling and/or engineering computational fluids
application                  people
the most basic
shallow assumption
flowline Stokes
                      • . . . we would all use the Navier-Stokes equations, for
derivation of SIA
                             incompressible fluids, as our flow model:
shelves and
streams
shallow shelf
approximation (SSA)
                                               ·u=0                               incompressibility
flow-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 flow: a
superficial view
shallow ice
                      is numerical ice flow modeling actually computational fluid
approximation (SIA)
                      dynamics?
heat analogy &
numerics
analogy
                        • yes
finite 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 flow modeling
application
the most basic
                          exciting:
shallow assumption
flowline Stokes
                             ◦   turbulence
derivation of SIA            ◦   convection
shelves and                  ◦   coriolis force
streams
shallow shelf                ◦   density/salinity stratification
approximation (SSA)
flow-line ice shelf           ◦   lots of tracers & chemistry (CO2 , methane, ozone, . . . )
numerical solution
scaling up              • none of the above list is relevant to ice flow!
next steps
practicalities
                        • so what could be interesting about the flow of slow, cold,
omitted models
                          shallow, non-turbulent, chemically-dull, old ice?
    Numerical
    modelling                                           ice is a slow, shear-thinning fluid
    Ed Bueler

ice flow: a
superficial view
shallow ice
approximation (SIA)
                      • our fluid is
heat analogy &
numerics
analogy                  slow:                        ρ (ut + u · u) ≈ 0
finite differences
exact solutions          non-Newtonian: viscosity ν is not constant in “τij = 2νDij ”
shallow ice
sheets
                      • and so the standard “full” model for ice flow is shear-thinning
solving the SIA
is it right?
                        Stokes as follows:
application
the most basic
shallow assumption               ·u=0                                incompressibility
flowline Stokes
derivation of SIA                 0=− p+                · τij + ρg   force balance
shelves and                                 n−1
streams                          Dij = Aτ         τij                flow law
shallow shelf
approximation (SSA)
flow-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 fluid 2
    Ed Bueler

ice flow: a
superficial view
shallow ice
approximation (SIA)     • the “slow” fact, ρ (ut + u ·      u) ≈ 0, means no forces of inertia
heat analogy &
numerics
                           from acceleration
analogy
finite 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 field at
the most basic
shallow assumption
                          every time, without storing it from the previous step2
flowline Stokes
derivation of SIA
                        • to see how glacier flow fields change in time we must track
shelves and
                          the changes to
streams
shallow shelf
                                ◦ the location of the ice surface and base
approximation (SSA)
flow-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 flow Stokes
    Ed Bueler

ice flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
                      more concretely, the shear-thinning Stokes model is this:
numerics
analogy                 • in the x, z plane flow case the Stokes equations say
finite 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
flowline Stokes
                                    ux = Aτ         τ11        flow law (diagonal)
derivation of SIA
                                               n−1
shelves and
                              uz + wx = 2Aτ          τ13       flow law (off-diagonal)
streams
shallow shelf
approximation (SSA)     • five equations in five unknowns (u, w, p, τ11 , τ13 )
flow-line ice shelf
numerical solution      • complicated enough . . . can I make it familiar by looking at a
scaling up

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

ice flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &                                                                              z=
numerics                                                                                         H
                      • rotated coordinates:                                     ice
analogy
finite 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
flowline 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 θ
flow-line ice shelf
numerical solution
scaling up             • but for a slab-on-a-slope, by definition, 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 flow: a            • add some boundary conditions:
superficial view
shallow ice
approximation (SIA)
                               w(base) = 0,            p(surface) = 0,               u(base) = u0 ,
heat analogy &
numerics              • further-simplified equations:
analogy
finite 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
flowline Stokes                                         1
derivation of SIA                           = u0 +       A(ρg sin θ)3 H 4 − (H − z)4
                                                       2
shelves and
streams
shallow shelf                      H
approximation (SSA)
flow-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 flow: a
superficial view
shallow ice                                                                       Horizontal Velocity (m/yr)
approximation (SIA)
                                                                              0         10         20      30
                                                                          0
heat analogy &
numerics
                      • do we believe these equations?
analogy
finite 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
flowline Stokes
derivation of SIA                                                              basal
                                                                               motion
shelves and                                                              200
streams
                                                                                        bedrock
shallow shelf
approximation (SSA)
flow-line ice shelf
numerical solution
scaling up
                                                              Velocity profile of the Athabasca
                                    u0
                                                              Glacier, Canada, derived from
next steps                               u(z)                 inclinometry. [Savage and Paterson,
practicalities
                                                              1963]
omitted models
    Numerical
    modelling                                                                  mass continuity
    Ed Bueler

ice flow: a
                      now we know the velocity u(z) . . . so what?
superficial 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                                                   ˆ
finite 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 figure at right:                           M(x)
shallow assumption
flowline 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)
flow-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 flow: a
superficial 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
finite 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
flowline 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 :
                                                               ¯
flow-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 flow: a            • ice sheet flow problems have three outstanding properties:
superficial view
shallow ice               1 slow
approximation (SIA)
                          2 non-Newtonian
heat analogy &
numerics                  3 shallow
analogy
finite 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
flowline Stokes
derivation of SIA

shelves and
streams                             2000
shallow shelf
approximation (SSA)
flow-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 flow: a
superficial 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 flow is
analogy
finite 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
flowline Stokes
derivation of SIA

shelves and
streams
shallow shelf
approximation (SSA)
flow-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 flow: a
superficial view
                      • the non-sliding isothermal SIA:
shallow ice
approximation (SIA)

heat analogy &
numerics
                                        Ht = M +      · ΓH n+2 | h|n−1 h       (1)
analogy
finite 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 flow law: Dij = A(T )τ n−1 τij
shallow assumption
flowline Stokes
                          ◦   Γ is a positive constant
derivation of SIA
                      • numerically solve equation (1), and you
shelves and
streams                 have a usable model for ice flow in the
shallow shelf
approximation (SSA)     Barnes Ice Cap [Mahaffy, 1976]
flow-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 flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
numerics              1 ice flow: an outsider’s superficial view
analogy
finite 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
flowline Stokes
derivation of SIA     4 shelves and streams
shelves and
streams
shallow shelf
approximation (SSA)
                      5 next steps
flow-line ice shelf
numerical solution
scaling up

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

ice flow: a
superficial view       • recall Newton’s law of cooling
shallow ice
approximation (SIA)
                                dT
heat analogy &                     = −µ(T − Tambient )
numerics                        dt
analogy
finite 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
flowline Stokes                 = −˜ Tj − (Tj−1 + Tj+1 )
                                   µ
derivation of SIA          dt            2
shelves and                      ˜
                                 µ
streams                        = (Tj−1 − 2Tj + Tj+1 )
shallow shelf                    2
approximation (SSA)
flow-line ice shelf
numerical solution
                        (where µ is material constant
                                ˜
scaling up              proportional to ∆x −2 )
next steps
practicalities        • this suggests finite difference
omitted models
                        approximation of a PDE:

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

ice flow: a
superficial view
shallow ice
approximation (SIA)
                      continuum version:
heat analogy &
numerics               • Fourier rewrote Newton’s law as heat
analogy
finite differences        flux 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)
flowline Stokes
derivation of SIA
                        • set D = k /(ρc) and M = f /(ρc):
shelves and
streams
shallow shelf
approximation (SSA)
                                                  ut = M +      · (D   u)                   (2)
flow-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 flow: a            • side-by-side comparison:
superficial 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)
finite 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
flowline Stokes
                      • non-sliding shallow ice flow diffuses the ice because the flow
derivation of SIA
                        is down the surface gradient
shelves and
streams
shallow shelf
approximation (SSA)   • some issues with this analogy:
flow-line ice shelf
numerical solution        ◦ H = h when bed is not flat . . . 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                                    finite differences for heat equation
    Ed Bueler

ice flow: a
superficial view
shallow ice
approximation (SIA)
                      basic ideas of finite differences:
heat analogy &
numerics                • for differentiable f (x) and any ∆x, Taylor says
analogy
finite 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 + . . .
flowline 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)
flow-line ice shelf      • basic finite 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                                   finite differences for heat equation 2
    Ed Bueler

ice flow: a
superficial 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
finite 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) =
flowline 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)
flow-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 flow: a
superficial 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
finite 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
flowline 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   →
flow-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 flow: a
superficial 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
finite 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
flowline 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)
flow-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 five quantities at old time tn  −→                     j-1      j           j+1
    Numerical
    modelling                                                                            implementation
    Ed Bueler

ice flow: a            function U = heat(D,J,K,dt,N)
superficial 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
finite 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
flowline 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
flow-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 flow: a
superficial 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
finite 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
flowline Stokes
derivation of SIA          1
                                                                                       1
shelves and               0.8
streams                                                                               0.8
                          0.6
shallow shelf
approximation (SSA)                                                                   0.6
                      u   0.4
flow-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 flow: a
superficial view
shallow ice                  • both figures are from solving ut = D(uxx + uyy ) on the same
approximation (SIA)

heat analogy &
                                 50 × 50 grid, at same final time and with same D, but with
numerics
analogy
                                 slightly different time steps
finite 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
flowline 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


flow-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 flow: a
superficial view
                        • the 1D first-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
finite 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 coefficient 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
flowline 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 coefficients” is a sufficient stability criterion
flow-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 finite differences, including stability: [Morton
                      and Mayers, 2005]
    Numerical
    modelling                        adaptive implementation: guaranteed stability
    Ed Bueler

ice flow: a
superficial view
                      function U = heatadapt(D,J,K,tf)
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
finite 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) );
flowline 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)
flow-line ice shelf                   http://www.dms.uaf.edu/~bueler/karthaus/mfiles/heatadapt.m
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 fix: implicitness
    Ed Bueler

ice flow: a
superficial 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 first-order implicit →,
finite 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 first-order explicit and implicit,
the most basic                                                                        n+1
shallow assumption
flowline 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
flow-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 efficiencies, say about 97% of the time:
                                      premature optimization is the root of all evil.
    Numerical
    modelling                                                  variable diffusivity and time steps
    Ed Bueler

ice flow: a
superficial 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 ) )
finite 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) ≈
flowline 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
flow-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 flow: a            function [U,dtav] = diffusion(Lx,Ly,J,K,Dup,Ddown,Dright,Dleft,U0,tf,b)
superficial 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
finite 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);
flowline 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)               );
flow-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 flow: a
superficial 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
finite differences       verification
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
flowline 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
streams                 exact predictions about nature
shallow shelf
approximation (SSA)
flow-line ice shelf
                      • if we can only approximately find 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 verification [Wesseling, 2001]
    Numerical
    modelling                                          exact solutions of heat equation
    Ed Bueler

ice flow: a
superficial 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
finite 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
flowline Stokes
                        by shrinking the output
derivation of SIA       (vertical) axis and
shelves and             simultaneously lengthening the
streams
shallow shelf           input (horizontal) axis
approximation (SSA)
flow-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                finding the Green’s function of the heat equation
    Ed Bueler

ice flow: a
superficial view
shallow ice           • the best-known way to find 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
finite differences
exact solutions
                      • facts used to find 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
flowline 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)
flow-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                      finding a “similarity” solution of heat equation
    Ed Bueler

ice flow: a
superficial 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
finite 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
flowline 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)
flow-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 simplifies to a
                        derivative: 2αsφ + s2 φ = s2 φ
    Numerical
    modelling                   finding a “similarity” solution of heat equation 2
    Ed Bueler

ice flow: a
                      • simplify and integrate, choose C = 0, integrate again:
superficial view
shallow ice                                         1
approximation (SIA)                                − s2 φ = D s φ + C
heat analogy &
                                                    2
numerics                                                       s
analogy                                               φ =−       φ
finite 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
flowline 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
flow-line ice shelf      fundamental
numerical solution
scaling up

next steps
practicalities
omitted models
    Numerical
    modelling                                                                             similarity solutions
    Ed Bueler

ice flow: a
superficial view
shallow ice              • conclusion from previous slides: similarity variables for heat
approximation (SIA)

heat analogy &
                           equation are
numerics
analogy                                input scaling                                output scaling
finite 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
flowline 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)
flow-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 flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
numerics              1 ice flow: an outsider’s superficial view
analogy
finite 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
flowline Stokes
derivation of SIA     4 shelves and streams
shelves and
streams
shallow shelf
approximation (SSA)
                      5 next steps
flow-line ice shelf
numerical solution
scaling up

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

ice flow: a            • jump forward to 1981 . . . breaking news:                                                      P. Halfar finds the
superficial view
shallow ice                  fundamental solution of the SIA!                                                 [Halfar, 1981; 1983]
approximation (SIA)

heat analogy &
                      • when n = 3, Halfar’s 2D solution for Glen flow law has
numerics
analogy
                             scalings
finite 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 flows without
solving the SIA              additional accumulation then its thinning really slows down as
is it right?
application                  the shape flattens out
the most basic
shallow assumption
flowline 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)
flow-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 flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
finite differences
exact solutions

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

shelves and
streams
shallow shelf
approximation (SSA)
flow-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 flow: a
superficial view
shallow ice
approximation (SIA)
                      • for n = 3 the solution is:
heat analogy &
numerics
                                                                                                             3/7
analogy
                                                         1/9                             1/18           4/3
finite 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
flowline 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)
flow-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                                                                              finding Halfar solution
    Ed Bueler
                      the method is the same as that which found the fundamental solution of
ice flow: a            the heat equation:
superficial 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
finite 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ϕ |ϕ | ϕ
flowline Stokes
derivation of SIA

shelves and             •   integrate this and use finiteness of ϕ(0), ϕ (0), and use ϕ (s) ≤ 0, to get separable
streams                     ODE
                                                                        4       3
shallow shelf                                                   βs = Γϕ (−ϕ )
approximation (SSA)
flow-line ice shelf      •   suppose H(t0 , 0) = H0 for some t0 > 0 to find 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 flow: a
superficial view
shallow ice
approximation (SIA)
                      function H = halfar(t,x,y)
heat analogy &
numerics
analogy               g = 9.81;   rho = 910.0;    secpera = 31556926;
finite 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);
flowline 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)
flow-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 flow: a
superficial view       • Nye and others [2000] compared the long-time
shallow ice
approximation (SIA)       consequences of different flow laws for the South Mars Polar
heat analogy &
numerics
                          Cap
analogy
finite differences
                      • they compared the long-time behavior of the corresponding
exact solutions
                        Halfar solutions for the different flow 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 ] flow laws will
shallow assumption
flowline 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
flow-line ice shelf
numerical solution
                                 probably consists of water ice, with an unknown
scaling up
                                 admixture of dust.
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 flow: a
superficial 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
finite 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
flowline Stokes           margins      H→0
derivation of SIA                                      but h is singular
shelves and
streams
shallow shelf
approximation (SSA)
                      • numerically, margins are worse than domes
flow-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 flat 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 flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
                      • for numerical stability we want
finite 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
flowline 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)
flow-line ice shelf
                             accuracy,
numerical solution
scaling up
                        to put D on the staggered grid
next steps
practicalities
omitted models
    Numerical
    modelling                                          SIA implementation: flat bed case
    Ed Bueler

                      function [H,dtlist] = siaflat(Lx,Ly,J,K,H0,deltat,tf)
ice flow: a
superficial 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;
finite 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) );
flowline 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 + ...
flow-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
                        [H,dtadapt] = diffusion(Lx,Ly,J,K,Dup,Ddn,Drt,Dlt,H,deltat);
                        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 flow: a
superficial 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
finite 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, final surface, and time-steps shown below
flowline Stokes
derivation of SIA

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

next steps
practicalities
omitted models
    Numerical
    modelling                                verification of numerical ice flow codes
    Ed Bueler

ice flow: a
superficial 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
finite 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 flow models?
the most basic
shallow assumption
flowline 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 flow line and                 [Schoof, 2006; van der Veen, 1985]
shallow shelf
approximation (SSA)           cross-flow cases
flow-line ice shelf       Blatter solution in flow line                  [Glowinski and Rappaz, 2003]
numerical solution       flow 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 flow: a
superficial 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
finite 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,...
flowline 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)));
flow-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 flow: a
superficial 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.
finite 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
flowline 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)
flow-line ice shelf
                      errors for 160 x 160 grid:
numerical solution    average thickness error =    0.988     figure 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 flow: a
superficial 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
finite 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      ...
flowline 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)
flow-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 flow: a            • at Karthaus 2009, a computer project: modify siaflat.m to
superficial 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
finite 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
flowline 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
flow-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 flow: a
superficial 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
finite 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
flowline Stokes
                            advection time-scale                                                                            t (a)



                        • 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
flow-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 flow: a                                                                            air
superficial view
shallow ice
approximation (SIA)

heat analogy &
                      • careful derivation of the SIA is next                                  ice
numerics
analogy               • makes “shallowness assumptions”
finite 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
flowline 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)
flow-line ice shelf
                      • the configurations 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 flow: a
superficial 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:
finite 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
flowline Stokes
derivation of SIA
                      • the key facts we need:
shelves and               ◦ the incompressibility of ice
streams
shallow shelf
approximation (SSA)
                                                                ux + vy + wz = 0
flow-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 flow: a
superficial 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)
finite differences
exact solutions       • define the map-plane flux 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
flowline 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)
flow-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 flow: a            for example, here is a sketch of how to show                        (3) & (4) =⇒ (5):
superficial 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
finite 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
flowline 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
flow-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 flow: a
superficial 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
finite 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
flowline 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)
flow-line ice shelf
                      • for example, the derivation of the SIA from Stokes needs to
numerical solution
scaling up
                        only go far enough to find an expression for the horizontal
next steps
                        velocity (u, v )
practicalities
omitted models
    Numerical
    modelling                                                                  origin of SIA
    Ed Bueler

ice flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &        where does the “shallow ice approximation” come from?
numerics
analogy                 • historically:
finite differences
exact solutions
                            • Fowler and Larson [1978], Morland and Johnson [1980], and
shallow ice                    Hutter [1983]
sheets
                            • fairly recent compared to the flow law of ice (1950s) and
solving the SIA
is it right?                   Stokes linearly-viscous model (∼1860)
application
the most basic          • logically:
shallow assumption
flowline 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
flow-line ice shelf
numerical solution      • more precisely:
scaling up
                            • as sketched in the next slides, for the x, z plane flow case only
next steps
practicalities
omitted models
    Numerical
    modelling                           Stokes equations in x, z plane flow case
    Ed Bueler

ice flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &        • suppose all flow is in the (screen) x, z plane
numerics
analogy               • . . . so out-of-plane velocity is zero (v = 0) and all quantities
finite 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
flowline Stokes
derivation of SIA
                                                             D13      0   D33
shelves and
streams               • incompressibility of ice in plane flow case says D has trace
shallow shelf
approximation (SSA)     zero:
flow-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 flow case 2
    Ed Bueler

ice flow: a
superficial 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 )
finite 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
flowline 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)
flow-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 flow case 3
    Ed Bueler

ice flow: a
superficial 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
finite differences
exact solutions                                           Dij = Aτ 2 τij
shallow ice
sheets
solving the SIA       • x, z plane flow case: from simplifications on previous slides,
is it right?
application
                        the isothermal Stokes equations say
the most basic
shallow assumption
flowline 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)
flow-line ice shelf
numerical solution                                 ux = Aτ 2 τ11
scaling up

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

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

ice flow: a
superficial view
shallow ice
approximation (SIA)
                      • previous equations apply within the ice fluid
heat analogy &
numerics              • . . . but boundary conditions are needed
analogy
finite 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 flow 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:
flowline Stokes
derivation of SIA

shelves and                                           τ13 = (τ11 − p)hx
streams
shallow shelf
approximation (SSA)
                                         p + τ11 + τ13 hx = 0
flow-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 flow: a
superficial 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
finite 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
flowline 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 ∗
flow-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 flow: a
superficial 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 “scientific judgment” (i.e. controversy)
analogy
finite 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
flowline 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)
flow-line ice shelf
numerical solution                 balance of shear stress
scaling up                                                          [τ ] = ρgd
next steps
                                    and pressure gradient
practicalities
                                                                    [u]
omitted models
                                balance of flow law scales               = 2[A][τ ]3
                                                                     d
    Numerical
    modelling                                                                         scaling the equations
    Ed Bueler
                      example:
ice flow: a
superficial view         • recall this equation from Stokes equations (6):
shallow ice
approximation (SIA)

heat analogy &                                             px = τ11,x + τ13,z
numerics
analogy
finite 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
flowline 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)
flow-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 flow: a
superficial view
shallow ice
                      • Stokes system now appears this way (stars are removed):
approximation (SIA)

heat analogy &
numerics                                   ux + wz = 0
analogy
                                                                        2
finite 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
flowline 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
flow-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 flow: a
superficial 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]
finite 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]
flowline 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)
flow-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 flow: a
superficial view
shallow ice           • solving for d, , [τ ], [u] gives
approximation (SIA)

heat analogy &
numerics                                               4        1/8
                                                        [a]
analogy
                                         d=                           ≈ 3600 m
finite 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
flowline 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
flow-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                                                             finally get to the SIA
    Ed Bueler

ice flow: a
                      • setting = 0 in the scaled Stokes equations, get these equations
superficial 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
finite 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
flowline 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)
flow-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                                                          finally get to the SIA 2
    Ed Bueler

ice flow: a
superficial view
shallow ice
                      • remember I said that it suffices to find 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
finite differences       stars:
                                                A∗ ∗
                                        u ∗ = − (hx ∗ )3 (H ∗ )4 − (h∗ − z ∗ )4
exact solutions

shallow ice                                     4
sheets
solving the SIA       • and then unscale to find 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
flowline 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
flow-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                                                      finally get to the SIA 3
    Ed Bueler

ice flow: a
superficial view                                             ´h
shallow ice
approximation (SIA)
                      • recall that the flux of ice is q =   b
                                                                 u(z) dz
heat analogy &        • so we integrate the expression for u(z):
numerics
analogy
finite 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
flowline 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)
flow-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 flow case
    Numerical
    modelling                                              have I oversold diffusivity?
    Ed Bueler

ice flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &        • I have asserted that the default model for ice sheets, the SIA,
numerics
analogy
                        is a diffusion
finite 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
flowline Stokes
derivation of SIA

shelves and           • have I oversold this diffusivity analogy?
streams
shallow shelf             ◦ probably
approximation (SSA)
flow-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 flow: a
superficial view       diffusivity
shallow ice
approximation (SIA)      • define the local bed topography by removing the local mean
heat analogy &
numerics
                           over some range λ ≈ 10 km:
analogy
finite differences
                                                               λ
exact solutions
                                        ˜
                                        b(x, ξ) = b(x + ξ) −        b(x + ξ) dξ
shallow ice                                                    −λ
sheets
solving the SIA         • define this strange average of the local bed:
is it right?
application                                                              −1/n
the most basic                                                −(n+2)/n
shallow assumption
                                              λ      ˜
                                                     b(x, ξ)
flowline 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)
flow-line ice shelf        get closer to solving the Stokes equations by making these
numerical solution
scaling up
                          two modifications 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 flow: a
superficial 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
finite 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
flowline 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)
flow-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 flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
numerics              1 ice flow: an outsider’s superficial view
analogy
finite 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
flowline Stokes
derivation of SIA     4 shelves and streams
shelves and
streams
shallow shelf
approximation (SSA)
                      5 next steps
flow-line ice shelf
numerical solution
scaling up

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

ice flow: a
superficial 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
finite differences
exact solutions

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

shelves and
streams
shallow shelf
approximation (SSA)
flow-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 flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
finite 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
flowline Stokes
derivation of SIA

shelves and
streams
shallow shelf
approximation (SSA)
flow-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 flow 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 flow: a
superficial view       • ice streams have fast flow (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
finite 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
flowline 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
flow-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 flow line: notation again
    Ed Bueler

ice flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
finite differences
exact solutions

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

shelves and
streams
shallow shelf
                         figure modified from [Schoof, 2007]
approximation (SSA)
flow-line ice shelf
numerical solution
scaling up

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

ice flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
                      • we will only work in plane flow case
numerics
analogy               • we will numerically solve only the stress balance equation
finite 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
flowline 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)
flow-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 flow line context:
ice flow: a
superficial 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
finite 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   
                                                                                      
flowline 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
flow-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 flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &        • we need a more limited goal!
numerics
analogy               • first goal: solve the steady state ice shelf in the isothermal,
finite differences
exact solutions          plane flow, 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,
flowline 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:
flow-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 flow: a
                      1   using flotation criterion h = (1 − ρ/ρw )H, and because of no bed
superficial 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
finite 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:
flowline 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
flow-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 find u(x), then return to (10) to find 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 flow: a
superficial 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
finite differences
exact solutions

shallow ice
sheets
solving the SIA
is it right?
application
the most basic
shallow assumption
flowline 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)




flow-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 flow: a
superficial 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):
finite 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
flowline 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 find
shelves and
streams                 the velocity u(x)
shallow shelf
approximation (SSA)   • I’ll describe the numerical method for a ice shelf or ice stream
flow-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 flow: a
superficial 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
finite 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
flowline 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-satisfied
approximation (SSA)
flow-line ice shelf
                      • idea as algorithm: from u (k −1) , find 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 flow: a
superficial view
shallow ice
approximation (SIA)
                      where do you get an initial guess u (0) (x) for the velocity?
heat analogy &
numerics
analogy
                        • for floating ice, a velocity comes from assuming a uniform
finite 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
flowline 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:
flow-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 flow: a
superficial 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
finite 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
flowline 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
flow-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 flow: a
superficial 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
finite 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
flowline 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
flow-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 flow: a
superficial 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                                                 
                                                                                                             
finite 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
flowline 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)
flow-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 flow: a
superficial 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(:);
finite 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
flowline 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
flow-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 flow: a
                      • before proceeding to solve nonlinear SSA problem, we can
superficial 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
finite 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
flowline Stokes
                               maximum error




derivation of SIA

shelves and                                    10-5
streams
shallow shelf
approximation (SSA)
flow-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 flow: a
superficial 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)’;
finite 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
flowline 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);
flow-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 flow: a            • first: let’s numerically solve the problem for which we know
superficial view
shallow ice
approximation (SIA)
                        exact answer
heat analogy &        • here use very coarse 10 km grid
numerics
analogy
                                                                                                          -1
finite 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
flowline Stokes
                              velocity (m a )
                             -1




derivation of SIA
                                                200

shelves and
streams                                         150
shallow shelf
approximation (SSA)
flow-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 flow: a
                      • this does convergence analysis of ssaflowline.m:
superficial 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)
finite differences
exact solutions       • goal is to see velocity error shrink by rate built-into our finite
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
flowline Stokes
derivation of SIA
                                                           10-1
shelves and
                                   maximum error (m a-1)




streams
shallow shelf
approximation (SSA)
                                                           10-2
flow-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 flow: a
                      • flow lines (1D) are never very realistic . . .
superficial view
shallow ice
                      • “diagnostic” (= fixed 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
finite 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
flowline Stokes
derivation of SIA

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

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

ice flow: a
superficial 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
finite differences
exact solutions                ∗ tridiagonal in SSA flow-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
flowline 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)
flow-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 flow: a
superficial view
shallow ice
                      • all stress balance equations used for ice flow modeling are
approximation (SIA)
                        nonlinear equations for velocity u in the form
heat analogy &
numerics
analogy
finite 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 fields (e.g. thickness)
application
the most basic
shallow assumption
                      • SSA and Blatter stress balance equations can be written in
flowline Stokes
derivation of SIA
                        form
shelves and
                                                   A(u) u = g
streams
shallow shelf
approximation (SSA)       ◦ where A is a nonlinear function of u,
flow-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 flow: a            • the iteration scheme we have used for SSA is “Picard”:
superficial 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
finite 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
flowline 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)
flow-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 flow: a            • big computers still take a long time to solve SSA equations
superficial 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 (finite
finite differences
exact solutions
                        difference, finite 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
flowline 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
flow-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 flow: a            • the PETSc [Balay and others,                                   Level of
superficial 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
finite 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
flowline 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
flow-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 flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
numerics              1 ice flow: an outsider’s superficial view
analogy
finite 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
flowline Stokes
derivation of SIA     4 shelves and streams
shelves and
streams
shallow shelf
approximation (SSA)
                      5 next steps
flow-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 flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &        you need:
numerics
analogy
                        • comfort in a technical computing environment, usually Unix,
finite 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.)
flowline 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
flow-line ice shelf
numerical solution      • know some tools for NetCDF files
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 flow: a
superficial 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:
finite differences
exact solutions
                          ◦ M ATLAB, Comsol, PETSc, libmesh, Elmer, Trilinos, Triangle,
shallow ice
                            etc. handle fundamental numerical linear algebra, mesh
sheets                      generation, finite 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
flowline 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
flow-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 flow: a
superficial 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)
finite differences           range relevant to Antarctic ice                1e-26

exact solutions
                          ◦ dissipation of flow energy is
shallow ice
sheets                      significant 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” fluid
the most basic
shallow assumption          instability from feedback between
flowline Stokes
derivation of SIA
                            the above two items (lower figure;
shelves and
                             EISMINT II [Payne and others, 2000])
streams
shallow shelf         • ice temperature is part of the “long
approximation (SSA)
flow-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 flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
finite 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
flowline Stokes            ◦ energy-balance models
derivation of SIA

shelves and           • very significant source of uncertainty for Greenland ice sheet
streams
shallow shelf           dynamics!
approximation (SSA)
flow-line ice shelf
numerical solution
scaling up

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

ice flow: a                                1
superficial 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
finite differences
                        if allowed enough time
exact solutions
                      • this changes bed topography and thus ice flow, so Earth
shallow ice
sheets                  deformation is important to modeling ice flow
solving the SIA
is it right?          • read?:
application
the most basic            ◦ [Peltier, 1998] = survey of observations and processes
shallow assumption
flowline Stokes            ◦ [Greve, 2001] = comparison of practical models
derivation of SIA

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

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

ice flow: a
superficial 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, . . .
finite differences
exact solutions             • requires explicit accounting for
shallow ice
sheets
                                ◦ incompressibility = a constraint on the flow
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?
flowline 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)
flow-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 flow: a
superficial 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
finite 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
flowline Stokes
derivation of SIA
                                                                other higher-order?
shelves and
streams
shallow shelf
approximation (SSA)
flow-line ice shelf
numerical solution
scaling up

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

ice flow: a
superficial view       • ice sheet/shelf modeling means free boundary problems
shallow ice
approximation (SIA)
                      • Hutter [1999] identifies some below
heat analogy &
numerics
analogy
finite differences
exact solutions

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

shelves and
streams
shallow shelf
approximation (SSA)
flow-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 flow: a
superficial 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
finite differences                a fixed-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
flowline Stokes
                            wire frame and stretched over an obstacle:
derivation of SIA

shelves and
streams               constraint:
shallow shelf
approximation (SSA)                  u≥ψ
flow-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 flow: a            • by volume, majority of ice sheet is cold (T < 0 ◦ C)
superficial 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
finite 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]
flowline 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)
flow-line ice shelf
numerical solution
scaling up

next steps
practicalities
omitted models
                                 Canadian-type                  Scandinavian-type

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

ice flow: a
superficial 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
finite 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
flowline Stokes
derivation of SIA

shelves and
streams
shallow shelf
approximation (SSA)
flow-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 flow: a
superficial view
shallow ice
approximation (SIA)
                          • side-by-side comparison, classical elastic membrane
heat analogy &              problem versus steady ice sheet problem
numerics
analogy
finite 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
flowline Stokes
derivation of SIA

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

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

ice flow: a
superficial view       • there are degenerate diffusions in other fields
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
finite differences
exact solutions
                      • specifically one by Raviart [1970]           (n = 3 case):     η = H 8/3
shallow ice           • the flat bed SIA simplifies:
sheets
solving the SIA
is it right?
application
                                          Ht = M +             · ΓH 5 | H|2     H
the most basic
shallow assumption
                                     becomes                                 ˜
flowline Stokes
derivation of SIA
                                        →          η 3/8       =M+        · (Γ| η|2    η)
                                                           t
shelves and
streams
shallow shelf
approximation (SSA)
                          ◦ new diffusivity (red) depends on | η| only, not η
flow-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 flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
                      • why does it matter that many glaciological problems are free
finite 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 fixed boundary solutions
flowline 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
flow-line ice shelf
numerical solution          physical processes
scaling up

next steps
practicalities
omitted models
    Numerical
    modelling                                       The End
    Ed Bueler

ice flow: a
superficial view
shallow ice
approximation (SIA)

heat analogy &
numerics
analogy
finite differences
exact solutions

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

shelves and
streams
shallow shelf
approximation (SSA)
flow-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