VIEWS: 12 PAGES: 125 POSTED ON: 7/23/2011
Numerical modelling Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & numerics Numerical modelling analogy ﬁnite differences exact solutions of ice sheets and ice shelves shallow ice sheets solving the SIA is it right? Ed Bueler application the most basic shallow assumption Dept of Mathematics and Statistics ﬂowline Stokes and Geophysical Institute derivation of SIA University of Alaska, Fairbanks shelves and streams shallow shelf approximation (SSA) September 2010 ﬂow-line ice shelf numerical solution Karthaus Summer School scaling up next steps practicalities omitted models Numerical modelling notation Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & numerics analogy ﬁnite differences exact solutions shallow ice ﬁgure modiﬁed from [Schoof, 2007] sheets solving the SIA is it right? application • coordinates t, x, y , z (z is vertical, positive upward) the most basic • subscripts for partial derivatives: ux = ∂u/∂x shallow assumption ﬂowline Stokes • h = ice surface elevation, b = bedrock elevation derivation of SIA • H = ice thickness shelves and streams • T = ice temperature shallow shelf • (u, v , w) = ice velocity approximation (SSA) ﬂow-line ice shelf • ρ = density of ice, ρw = density of ocean water numerical solution • g = acceleration of gravity scaling up • Dij = strain rate tensor next steps practicalities • stress tensors: σij = Cauchy and τij = deviatoric omitted models • 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 ﬂow: a superﬁcial view shallow ice approximation (SIA) slogans: heat analogy & • my focus is on approximating ice ﬂow numerics analogy ﬁnite differences • my example numerical codes actually work exact solutions shallow ice continuum models: sheets solving the SIA • shallow ice approximation (SIA) in 2D is it right? application the most basic • shallow shelf approximation (SSA) in ﬂowline (1D) shallow assumption ﬂowline Stokes • mass continuity & surface/base kinematical equations derivation of SIA shelves and • . . . and a little bit more streams shallow shelf approximation (SSA) numerical ideas: ﬂow-line ice shelf numerical solution • ﬁnite difference schemes scaling up next steps • solving algebraic systems arising in stress balances practicalities omitted models • veriﬁcation Numerical modelling outside of scope Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & not covered here: numerics analogy • thermomechanical coupling ﬁnite differences exact solutions • subglacial material/process modeling shallow ice sheets • snow/ﬁrn process modeling solving the SIA is it right? • earth deformation application the most basic shallow assumption • numerical solution of Stokes equations ﬂowline Stokes derivation of SIA • grounding line and calving front issues shelves and streams • anisotropy in ice ﬂow shallow shelf approximation (SSA) • paleoclimate modeling and “spin-up” ﬂow-line ice shelf numerical solution • ﬁnite element, spectral, wavelet, multigrid, . . . methods scaling up next steps • . . . etc. etc. practicalities omitted models Numerical modelling main equations Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & numerics analogy ﬁnite differences • shallow ice approximation (SIA), equation (1) exact solutions shallow ice • heat equation (2) (for analogy purposes) sheets solving the SIA • surface kinematical equation (3) is it right? application the most basic • basal kinematical equation (4) shallow assumption ﬂowline Stokes • map-plane mass continuity equation (5) derivation of SIA shelves and • Stokes equations in a ﬂow line, equations (6) streams shallow shelf • shallow shelf approximation (SSA), equation (7) approximation (SSA) ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models Numerical modelling M ATLAB/O CTAVE codes Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & numerics analogy • these lectures are structured around ﬁnite differences exact solutions • 7 codes which solve PDEs shallow ice • an additional 13 codes which plot, pre-process, post-process sheets solving the SIA • each PDE code is 1/2 to 1 page in length is it right? application • all tested in M ATLAB (v. 7.9) and O CTAVE (v. 3.2) the most basic shallow assumption ﬂowline Stokes • available as derivation of SIA ◦ .zip or .tar.gz, from my USB memory stick, or shelves and streams ◦ online: shallow shelf approximation (SSA) http://www.dms.uaf.edu/~bueler/karthaus/ ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models Numerical modelling Outline Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & numerics 1 ice ﬂow: an outsider’s superﬁcial view analogy ﬁnite differences exact solutions shallow ice 2 heat analogy & numerics sheets solving the SIA is it right? application 3 shallow ice sheets the most basic shallow assumption ﬂowline Stokes derivation of SIA 4 shelves and streams shelves and streams shallow shelf approximation (SSA) 5 next steps ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models Numerical modelling my ﬁrst goal Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & numerics analogy ﬁnite differences • I want to get to an equation for which I can say: exact solutions shallow ice sheets solving the SIA numerically solve this equation, and you’ve got a usable is it right? application model for . . . (a speciﬁc ice ﬂow problem) the most basic shallow assumption ﬂowline Stokes derivation of SIA • to get to that goal: I will very quickly recall the continuum shelves and streams mechanics of ice ﬂow shallow shelf approximation (SSA) • . . . from a naive outside view ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models Numerical modelling let’s treat ice in glaciers as a ﬂuid Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & • we describe ﬂuids primarily by a velocity ﬁeld u(t, x, y , z) numerics analogy ﬁnite differences • if the ice ﬂuid were much faster-moving than it actually is1 , exact solutions and if it were linearly-viscous like liquid water, then ice ﬂow shallow ice sheets modeling would be more familiar to other the solving the SIA is it right? climate-modeling and/or engineering computational ﬂuids application people the most basic shallow assumption ﬂowline Stokes • . . . we would all use the Navier-Stokes equations, for derivation of SIA incompressible ﬂuids, as our ﬂow model: shelves and streams shallow shelf approximation (SSA) ·u=0 incompressibility ﬂow-line ice shelf 2 numerical solution ρ (ut + u · u) = − p + ν u + ρg force balance scaling up next steps practicalities omitted models 1 if gravity were a lot stronger and/or ice were a lot weaker, for example Numerical modelling hmmm . . . does not sound like glaciology to me! Ed Bueler ice ﬂow: a superﬁcial view shallow ice is numerical ice ﬂow modeling actually computational ﬂuid approximation (SIA) dynamics? heat analogy & numerics analogy • yes ﬁnite differences exact solutions • at geophysical scale like atmosphere/ocean shallow ice sheets • . . . but it is a weird one solving the SIA is it right? • consider what makes atmosphere/ocean ﬂow modeling application the most basic exciting: shallow assumption ﬂowline Stokes ◦ turbulence derivation of SIA ◦ convection shelves and ◦ coriolis force streams shallow shelf ◦ density/salinity stratiﬁcation approximation (SSA) ﬂow-line ice shelf ◦ lots of tracers & chemistry (CO2 , methane, ozone, . . . ) numerical solution scaling up • none of the above list is relevant to ice ﬂow! next steps practicalities • so what could be interesting about the ﬂow of slow, cold, omitted models shallow, non-turbulent, chemically-dull, old ice? Numerical modelling ice is a slow, shear-thinning ﬂuid Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) • our ﬂuid is heat analogy & numerics analogy slow: ρ (ut + u · u) ≈ 0 ﬁnite differences exact solutions non-Newtonian: viscosity ν is not constant in “τij = 2νDij ” shallow ice sheets • and so the standard “full” model for ice ﬂow is shear-thinning solving the SIA is it right? Stokes as follows: application the most basic shallow assumption ·u=0 incompressibility ﬂowline Stokes derivation of SIA 0=− p+ · τij + ρg force balance shelves and n−1 streams Dij = Aτ τij ﬂow law shallow shelf approximation (SSA) ﬂow-line ice shelf numerical solution • equations above are true at every instant, and so scaling up geometry and boundary stress and ice strength next steps practicalities determine velocity instantaneously in the model omitted models Numerical modelling ice is a slow ﬂuid 2 Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) • the “slow” fact, ρ (ut + u · u) ≈ 0, means no forces of inertia heat analogy & numerics from acceleration analogy ﬁnite differences • thus velocity is not a state variable of any model for evolving exact solutions ice sheets and glaciers shallow ice sheets • . . . it is a “diagnostic” output . . . but a very important output! solving the SIA is it right? application • an ice sheet code usually recomputes the full velocity ﬁeld at the most basic shallow assumption every time, without storing it from the previous step2 ﬂowline Stokes derivation of SIA • to see how glacier ﬂow ﬁelds change in time we must track shelves and the changes to streams shallow shelf ◦ the location of the ice surface and base approximation (SSA) ﬂow-line ice shelf ◦ the force (stress) applied at the boundary numerical solution ◦ the ice strength, as it changes with ice temperature (for scaling up example) next steps practicalities omitted models 2 Bob Dylan would say: you don’t need a glaciologist to remember which way the wind blows Numerical modelling plane ﬂow Stokes Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & more concretely, the shear-thinning Stokes model is this: numerics analogy • in the x, z plane ﬂow case the Stokes equations say ﬁnite differences exact solutions shallow ice ux + wz = 0 incompressibility sheets solving the SIA px = τ11,x + τ13,z stress balance (x) is it right? application pz = τ13,x − τ11,z − ρg stress balance (z) the most basic shallow assumption n−1 ﬂowline Stokes ux = Aτ τ11 ﬂow law (diagonal) derivation of SIA n−1 shelves and uz + wx = 2Aτ τ13 ﬂow law (off-diagonal) streams shallow shelf approximation (SSA) • ﬁve equations in ﬁve unknowns (u, w, p, τ11 , τ13 ) ﬂow-line ice shelf numerical solution • complicated enough . . . can I make it familiar by looking at a scaling up next steps very-simpliﬁed situation? practicalities omitted models Numerical modelling slab-on-a-slope z Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & z= numerics H • rotated coordinates: ice analogy ﬁnite differences ˆ ˆ g = g sin θ x − g cos θ z exact solutions shallow ice sheets bed g solving the SIA roc k is it right? z= application 0 the most basic • as on previous slide, but with n = 3: x shallow assumption ﬂowline Stokes 2 derivation of SIA ux + wz = 0 ux = Aτ τ11 shelves and streams px = τ11,x + τ13,z + ρg sin θ uz + wx = 2Aτ 2 τ13 shallow shelf approximation (SSA) pz = τ13,x − τ11,z − ρg cos θ ﬂow-line ice shelf numerical solution scaling up • but for a slab-on-a-slope, by deﬁnition, there is no variation in x: next steps practicalities wz = 0 0 = τ11 omitted models τ13,z = −ρg sin θ uz = 2Aτ 2 τ13 pz = −ρg cos θ Numerical modelling slab-on-a-slope 2 Ed Bueler ice ﬂow: a • add some boundary conditions: superﬁcial view shallow ice approximation (SIA) w(base) = 0, p(surface) = 0, u(base) = u0 , heat analogy & numerics • further-simpliﬁed equations: analogy ﬁnite differences exact solutions τ13 = ρg sin θ(H − z) p = ρg cos θ(H − z) shallow ice w =0 τ11 = 0 sheets solving the SIA ˆ z is it right? • and application u(z) = u0 + 2A(ρg sin θ)3 (H − z )3 dz the most basic 0 shallow assumption ﬂowline Stokes 1 derivation of SIA = u0 + A(ρg sin θ)3 H 4 − (H − z)4 2 shelves and streams shallow shelf H approximation (SSA) ﬂow-line ice shelf numerical solution scaling up z next steps practicalities omitted models 0 u0 τ13(z) u(z) Numerical modelling slab-on-a-slope 3 Ed Bueler ice ﬂow: a superﬁcial view shallow ice Horizontal Velocity (m/yr) approximation (SIA) 0 10 20 30 0 heat analogy & numerics • do we believe these equations? analogy ﬁnite differences • velocity on last slide (and below) 50 exact solutions was from a formula Depth (m) shallow ice sheets • compare to observations at right 100 solving the SIA internal is it right? application • seems there is an element of truth deformation the most basic shallow assumption 150 ﬂowline Stokes derivation of SIA basal motion shelves and 200 streams bedrock shallow shelf approximation (SSA) ﬂow-line ice shelf numerical solution scaling up Velocity proﬁle of the Athabasca u0 Glacier, Canada, derived from next steps u(z) inclinometry. [Savage and Paterson, practicalities 1963] omitted models Numerical modelling mass continuity Ed Bueler ice ﬂow: a now we know the velocity u(z) . . . so what? superﬁcial view shallow ice • suppose now that our slab has slowly varying thickness H(t, x) approximation (SIA) heat analogy & • compute the vertical average of velocity, which depends on x: numerics analogy ˆ ﬁnite differences 1 H exact solutions ¯ u (x) = u(x, z) dz H 0 shallow ice sheets solving the SIA is it right? • consider change of area (really: ice application the most basic volume) in the ﬁgure at right: M(x) shallow assumption ﬂowline Stokes ˆ x2 derivation of SIA dA = ¯ ¯ M(x) dx + u1 H1 − u2 H2 shelves and dt streams x1 H1 H2 shallow shelf A approximation (SSA) ﬂow-line ice shelf • in limit where width dx = x2 − x2 is numerical solution small we have A ≈ dx H ¯ u1 ¯ u2 scaling up next steps • divide by dx and get practicalities omitted models ¯ Ht = M − (u H)x x1 x2 • this is the mass continuity equation Numerical modelling rough explanation of “shallow ice approximation” (SIA) Ed Bueler ice ﬂow: a superﬁcial view shallow ice • I’ll consider only u0 = 0 case in these lectures (“non-sliding SIA”) approximation (SIA) heat analogy & • from slab-on-slope velocity formula, numerics analogy ˆ H ﬁnite differences 1 exact solutions ¯ uH = A(ρg sin θ)3 H 4 − (H − z)4 dz 0 2 shallow ice sheets 1 4 5 solving the SIA = A(ρg sin θ)3 H is it right? 2 5 application 2 = A(ρg sin θ)3 H 5 the most basic shallow assumption ﬂowline Stokes 5 derivation of SIA shelves and • note sin θ ≈ tan θ = −hx streams shallow shelf approximation (SSA) • combine with mass continuity Ht = M − (u H)x : ¯ ﬂow-line ice shelf numerical solution 2 scaling up Ht = M + A(ρg)3 H 5 |hx |2 hx (0) next steps 5 x practicalities omitted models • it is pretty rough . . . but we will get back to it Numerical modelling slow, non-Newtonian, shallow Ed Bueler ice ﬂow: a • ice sheet ﬂow problems have three outstanding properties: superﬁcial view shallow ice 1 slow approximation (SIA) 2 non-Newtonian heat analogy & numerics 3 shallow analogy ﬁnite differences • all three properties will be in all of my models exact solutions shallow ice • regarding “shallow,” here in red is a to-scale cross section of sheets solving the SIA Greenland at 71◦ : is it right? 4000 application the most basic shallow assumption 3000 ﬂowline Stokes derivation of SIA shelves and streams 2000 shallow shelf approximation (SSA) ﬂow-line ice shelf 1000 numerical solution scaling up next steps practicalities 0 omitted models -1000 -600000 -400000 -200000 0 200000 400000 x Numerical modelling the isothermal shallow ice approximation (SIA) Ed Bueler ice ﬂow: a superﬁcial view is a model which applies reasonably well to shallow ice approximation (SIA) • grounded ice sheets heat analogy & numerics • on relatively uninteresting bed topography, whose ﬂow is analogy ﬁnite differences • not dominated by sliding nor by liquid water at base or margin exact solutions shallow ice sheets solving the SIA is it right? application the most basic shallow assumption ﬂowline Stokes derivation of SIA shelves and streams shallow shelf approximation (SSA) ﬂow-line ice shelf numerical solution scaling up next steps practicalities “Polaris Glacier,” northwest Greenland, photo 122, [Post and LaChapelle, 2000] omitted models • might apply to the above Numerical modelling isothermal shallow ice approximation 2 Ed Bueler ice ﬂow: a superﬁcial view • the non-sliding isothermal SIA: shallow ice approximation (SIA) heat analogy & numerics Ht = M + · ΓH n+2 | h|n−1 h (1) analogy ﬁnite differences exact solutions shallow ice ◦ generalizes slab-on-slope equation (0) earlier sheets ◦ h =H +b solving the SIA is it right? ◦ accumulation if M > 0, ablation if M < 0 application the most basic ◦ recall Glen ﬂow law: Dij = A(T )τ n−1 τij shallow assumption ﬂowline Stokes ◦ Γ is a positive constant derivation of SIA • numerically solve equation (1), and you shelves and streams have a usable model for ice ﬂow in the shallow shelf approximation (SSA) Barnes Ice Cap [Mahaffy, 1976] ﬂow-line ice shelf numerical solution scaling up next steps • questions: practicalities omitted models ◦ where does equation (1) come from? ◦ how to solve it numerically? ◦ how to think about it? Numerical modelling Outline Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & numerics 1 ice ﬂow: an outsider’s superﬁcial view analogy ﬁnite differences exact solutions shallow ice 2 heat analogy & numerics sheets solving the SIA is it right? application 3 shallow ice sheets the most basic shallow assumption ﬂowline Stokes derivation of SIA 4 shelves and streams shelves and streams shallow shelf approximation (SSA) 5 next steps ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models Numerical modelling compare to heat equation Ed Bueler ice ﬂow: a superﬁcial view • recall Newton’s law of cooling shallow ice approximation (SIA) dT heat analogy & = −µ(T − Tambient ) numerics dt analogy ﬁnite differences exact solutions where T is temperature and µ relates shallow ice to material and geometry of object sheets solving the SIA • Newton’s law for segments of an is it right? application insulated rod: the most basic shallow assumption dTj 1 ﬂowline Stokes = −˜ Tj − (Tj−1 + Tj+1 ) µ derivation of SIA dt 2 shelves and ˜ µ streams = (Tj−1 − 2Tj + Tj+1 ) shallow shelf 2 approximation (SSA) ﬂow-line ice shelf numerical solution (where µ is material constant ˜ scaling up proportional to ∆x −2 ) next steps practicalities • this suggests ﬁnite difference omitted models approximation of a PDE: Tt = DTxx Numerical modelling compare to heat equation 2 Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) continuum version: heat analogy & numerics • Fourier rewrote Newton’s law as heat analogy ﬁnite differences ﬂux in continua: q = −k u exact solutions shallow ice • by conservation of energy, allowing an sheets additional source of heat f , get heat solving the SIA is it right? equation: application the most basic shallow assumption ρcut = f + · (k u) ﬂowline Stokes derivation of SIA • set D = k /(ρc) and M = f /(ρc): shelves and streams shallow shelf approximation (SSA) ut = M + · (D u) (2) ﬂow-line ice shelf numerical solution scaling up • this heat equation is a diffusive, time-dependent PDE, like the SIA next steps model for ice sheets practicalities omitted models Numerical modelling analogy: SIA versus heat equation Ed Bueler ice ﬂow: a • side-by-side comparison: superﬁcial view shallow ice approximation (SIA) SIA: H(t, x, y ) is ice thickness heat: u(t, x, y ) is temperature heat analogy & numerics n+2 n−1 analogy Ht = M + · ΓH | h| h ut = M + · (D u) ﬁnite differences exact solutions shallow ice • we identify the diffusivity in the SIA: sheets solving the SIA is it right? application D = ΓH n+2 | h|n−1 the most basic shallow assumption ﬂowline Stokes • non-sliding shallow ice ﬂow diffuses the ice because the ﬂow derivation of SIA is down the surface gradient shelves and streams shallow shelf approximation (SSA) • some issues with this analogy: ﬂow-line ice shelf numerical solution ◦ H = h when bed is not ﬂat . . . so what? scaling up ◦ D depends on solution H(t, x, y ) . . . how does that complicate next steps practicalities the numerical solution method? omitted models ◦ D → 0 at margin (H → 0) and at dome (| h| → 0) . . . so what? • I’ll get back to these “issues”, but now let’s get our hands dirty and numerically solve the heat equation Numerical modelling ﬁnite differences for heat equation Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) basic ideas of ﬁnite differences: heat analogy & numerics • for differentiable f (x) and any ∆x, Taylor says analogy ﬁnite differences 1 1 f (x)∆x 2 + f (x)∆x 3 + . . . exact solutions f (x + ∆x) = f (x) + f (x)∆x + shallow ice 2 3! sheets solving the SIA is it right? • you can replace “∆x” with other expressions, e.g.: application the most basic 1 1 shallow assumption f (x − ∆x) = f (x) − f (x)∆x + f (x)∆x 2 − f (x)∆x 3 + . . . ﬂowline Stokes 2 3! derivation of SIA 4 shelves and f (x + 2∆x) = f (x) + 2f (x)∆x + 2f (x)∆x 2 + f (x)∆x 3 + . . . streams 3 shallow shelf approximation (SSA) ﬂow-line ice shelf • basic ﬁnite difference idea for differential equations: combine numerical solution scaling up expressions like these to approximate derivatives next steps • for all of these notes, grid points have equal spacing ∆x practicalities omitted models Numerical modelling ﬁnite differences for heat equation 2 Ed Bueler ice ﬂow: a superﬁcial view shallow ice • partial derivative expressions, for example with u = u(t, x): approximation (SIA) heat analogy & u(t + ∆t, x) − u(t, x) numerics ut (t, x) = + O(∆t), analogy ∆t ﬁnite differences u(t + ∆t, x) − u(t − ∆t, x) exact solutions ut (t, x) = + O(∆t 2 ), shallow ice 2∆t sheets u(t, x + ∆x) − u(t, x) solving the SIA ux (t, x) = + O(∆x), is it right? ∆x application u(t, x + ∆x) − u(t, x − ∆x) + O(∆x 2 ), the most basic shallow assumption ux (t, x) = ﬂowline Stokes 2∆x derivation of SIA u(t, x + ∆x) − 2u(t, x) + u(t, x − ∆x) shelves and uxx (t, x) = + O(∆x 2 ) streams ∆x 2 shallow shelf approximation (SSA) ﬂow-line ice shelf • sometimes we want a derivative in-between grid points: numerical solution scaling up u(t, x + ∆x) − u(t, x) ux (t, x + (∆x/2)) = + O(∆x 2 ) next steps practicalities ∆x omitted models • “+O(∆x 2 )” is better than “+O(∆x)” if ∆x is a small number Numerical modelling explicit scheme for heat equation Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) • recall 1D heat equation ut = Duxx heat analogy & numerics • the explicit scheme using notation ujn ≈ u(tn , xj ), so analogy ﬁnite differences exact solutions ujn+1 − ujn uj+1 − 2ujn + uj−1 n n shallow ice =D sheets ∆t ∆x 2 solving the SIA is it right? application • let ν = D∆t/(∆x)2 , so the most basic shallow assumption ﬂowline Stokes derivation of SIA ujn+1 = νuj+1 + (1 − 2ν)ujn + νuj−1 n n shelves and streams shallow shelf approximation (SSA) • scheme has stencil at right → ﬂow-line ice shelf n+1 numerical solution • advantage over implicit (later): ujn+1 is scaling up next steps determined by known quantities at time tn n practicalities omitted models j-1 j j+1 Numerical modelling explicit scheme for heat equation 2 Ed Bueler ice ﬂow: a superﬁcial view n • in 2D we write ujk ≈ u(tn , xj , yk ) shallow ice approximation (SIA) heat analogy & • the 2D explicit scheme for the M = 0 heat equation numerics ut = D(uxx + uyy ) is analogy ﬁnite differences n+1 n n n n n n n exact solutions ujk − ujk uj+1,k − 2ujk + uj−1,k uj,k +1 − 2ujk + uj,k −1 shallow ice =D + sheets ∆t ∆x 2 ∆y 2 solving the SIA is it right? application • or, with ν x := D∆t/(∆x)2 and ν y := D∆t/(∆y )2 , the most basic shallow assumption ﬂowline Stokes n+1 n n n n n ujk = (1 − 2ν x − 2ν y )ujk + ν x uj+1,k + uj−1,k + ν y uj,k +1 + uj,k −1 derivation of SIA shelves and streams shallow shelf k+1 approximation (SSA) ﬂow-line ice shelf numerical solution scaling up k next steps practicalities omitted models n+1 note: new value ujk is average (is it?) k-1 of ﬁve quantities at old time tn −→ j-1 j j+1 Numerical modelling implementation Ed Bueler ice ﬂow: a function U = heat(D,J,K,dt,N) superﬁcial view shallow ice approximation (SIA) dx = 2 / J; dy = 2 / K; heat analogy & [x,y] = meshgrid(-1:dx:1, -1:dy:1); numerics U = exp(-30*(x.*x + y.*y)); analogy ﬁnite differences nu_x = dt * D / (dx*dx); nu_y = dt * D / (dy*dy); exact solutions for n=1:N shallow ice U(2:J,2:K) = U(2:J,2:K) + ... sheets nu_x * ( U(3:J+1,2:K) - 2 * U(2:J,2:K) + U(1:J-1,2:K) ) + ... solving the SIA is it right? nu_y * ( U(2:J,3:K+1) - 2 * U(2:J,2:K) + U(2:J,1:K-1) ); application end the most basic shallow assumption surf(x,y,U), shading(’interp’), xlabel x, ylabel y, zlabel u ﬂowline Stokes derivation of SIA http://www.dms.uaf.edu/~bueler/karthaus/mfiles/heat.m shelves and streams shallow shelf approximation (SSA) • solves ut = D(uxx + uyy ) on square −1 < x < 1, −1 < y < 1 ﬂow-line ice shelf numerical solution • choice: gaussian initial condition scaling up next steps • “colon notation” removes loops over spatial variables practicalities omitted models • to approximate u on 30 × 30 spatial grid, with D = 1 and N = 20 steps of length ∆t = 0.001, » heat(1.0,30,30,0.001,20) Numerical modelling the look of success Ed Bueler ice ﬂow: a superﬁcial view • solving ut = D(uxx + uyy ) on 30 × 30 grid shallow ice approximation (SIA) heat analogy & numerics analogy approximate solution u(t, x, y ) at ﬁnite differences initial condition u(0, x, y ) exact solutions t = 0.02 with ∆t = 0.001 shallow ice sheets solving the SIA is it right? application the most basic shallow assumption ﬂowline Stokes derivation of SIA 1 1 shelves and 0.8 streams 0.8 0.6 shallow shelf approximation (SSA) 0.6 u 0.4 ﬂow-line ice shelf u 0.4 numerical solution 0.2 scaling up 0.2 01 next steps 0.5 1 01 practicalities 0 0.5 0.5 1 y 0 0.5 omitted models -0.5 0 -0.5 x y 0 -1 -1 -0.5 -0.5 x -1 -1 Numerical modelling the look of instability Ed Bueler ice ﬂow: a superﬁcial view shallow ice • both ﬁgures are from solving ut = D(uxx + uyy ) on the same approximation (SIA) heat analogy & 50 × 50 grid, at same ﬁnal time and with same D, but with numerics analogy slightly different time steps ﬁnite differences exact solutions shallow ice sheets solving the SIA is it right? 0.25 0.25 application the most basic 0.2 0.2 shallow assumption ﬂowline Stokes 0.15 0.15 derivation of SIA u 0.1 u 0.1 shelves and streams 0.05 0.05 shallow shelf approximation (SSA) 01 01 ﬂow-line ice shelf 0.5 1 0.5 1 0.5 0.5 numerical solution 0 0 y 0 y 0 scaling up -0.5 -0.5 x -0.5 -0.5 x -1 -1 -1 -1 next steps practicalities omitted models D∆t D∆t = 0.402 = 0.625 ∆x 2 ∆x 2 Numerical modelling why unstable? Ed Bueler ice ﬂow: a superﬁcial view • the 1D ﬁrst-order explicit scheme in the form shallow ice ujn+1 = νuj+1 + (1 − 2ν)ujn + νuj−1 n n approximation (SIA) heat analogy & numerics analogy ﬁnite differences gives new value ujn+1 as an average of old values, exact solutions shallow ice • but it is only an average if the middle coefﬁcient is positive:3 sheets solving the SIA is it right? D∆t application 1 − 2ν = 1 − 2 ≥0 the most basic shallow assumption ∆x 2 ﬂowline Stokes derivation of SIA • true averaging is always stable because averaged wiggles shelves and streams are smaller than the previous wiggles shallow shelf approximation (SSA) • “positive coefﬁcients” is a sufﬁcient stability criterion ﬂow-line ice shelf numerical solution • the same thing as a time-step restriction: scaling up next steps ∆x 2 practicalities ∆t ≤ omitted models 2D 3 recommended basic reference for ﬁnite differences, including stability: [Morton and Mayers, 2005] Numerical modelling adaptive implementation: guaranteed stability Ed Bueler ice ﬂow: a superﬁcial 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 ﬁnite differences t = 0.0; count = 0; exact solutions while t < tf shallow ice dt0 = 0.25 * min(dx,dy)^2 / D; sheets dt = min(dt0, tf - t); solving the SIA is it right? nu_x = dt * D / (dx*dx); nu_y = dt * D / (dy*dy); application U(2:J,2:K) = U(2:J,2:K) + ... the most basic nu_x * ( U(3:J+1,2:K) - 2 * U(2:J,2:K) + U(1:J-1,2:K) ) + ... shallow assumption nu_y * ( U(2:J,3:K+1) - 2 * U(2:J,2:K) + U(2:J,1:K-1) ); ﬂowline Stokes derivation of SIA t = t + dt; count = count + 1; end shelves and streams shallow shelf surf(x,y,U), shading(’interp’), xlabel x, ylabel y, zlabel u approximation (SSA) ﬂow-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 ﬁx: implicitness Ed Bueler ice ﬂow: a superﬁcial view • explicit scheme is only “conditionally stable”, shallow ice approximation (SIA) but there are methods which are stable for n+1 heat analogy & any positive time step ∆t numerics analogy • the simplest such is ﬁrst-order implicit →, ﬁnite differences n exact solutions shallow ice ujn+1 − ujn n+1 uj+1 − 2ujn+1 + n+1 uj−1 j-1 j j+1 =D sheets solving the SIA ∆t ∆x 2 is it right? application • another is Crank-Nicolson →; instead of O(∆t, ∆x 2 ) like ﬁrst-order explicit and implicit, the most basic n+1 shallow assumption ﬂowline Stokes derivation of SIA Crank-Nicolson is O(∆t 2 , ∆x 2 ) shelves and streams n shallow shelf approximation (SSA) j-1 j j+1 ﬂow-line ice shelf numerical solution scaling up • but for implicit and Crank-Nicolson methods you have to solve next steps systems of equations to take each step practicalities omitted models • Donald Knuth has advice for ice sheet modelers: We should forget about small efﬁciencies, say about 97% of the time: premature optimization is the root of all evil. Numerical modelling variable diffusivity and time steps Ed Bueler ice ﬂow: a superﬁcial view • recall analogy (SIA) ↔ (heat eqn) shallow ice approximation (SIA) • the SIA has a diffusivity D which varies in time and space, so heat analogy & numerics by analogy: analogy ut = M + · (D(t, x, y ) ) ﬁnite differences exact solutions • the explicit method is conditionally stable with the “same” shallow ice sheets stability condition if we evaluate diffusivity D(t, x, y ) at solving the SIA staggered grid points: is it right? application the most basic Dj+1/2,k (uj+1,k − uj,k ) − Dj−1/2,k (uj,k − uj−1,k ) shallow assumption · (D(t, x, y ) u) ≈ ﬂowline Stokes ∆x 2 derivation of SIA Dj,k +1/2 (uj,k +1 − uj,k ) − Dj,k −1/2 (uj,k − uj,k −1 ) + shelves and ∆y 2 streams shallow shelf approximation (SSA) k+1 ﬂow-line ice shelf numerical solution scaling up in stencil at right: next steps diamonds: u k practicalities omitted models triangles: D k-1 j-1 j j+1 Numerical modelling general diffusion equation code Ed Bueler ice ﬂow: a function [U,dtav] = diffusion(Lx,Ly,J,K,Dup,Ddown,Dright,Dleft,U0,tf,b) superﬁcial view shallow ice dx = 2 * Lx / J; dy = 2 * Ly / K; approximation (SIA) [x,y] = ndgrid(-Lx:dx:Lx, -Ly:dy:Ly); heat analogy & U = U0; numerics if nargin < 11, b = zeros(size(U0)); end analogy ﬁnite differences exact solutions t = 0.0; count = 0; while t < tf shallow ice sheets Dregular = max(max(Dup,Ddown),max(Dleft,Dright)); solving the SIA maxD = max(max(Dregular)); is it right? dt0 = 0.25 * min(dx,dy)^2 / maxD; application dt = min(dt0, tf - t); the most basic shallow assumption mu_x = dt / (dx*dx); mu_y = dt / (dy*dy); ﬂowline Stokes Ushift = U + b; derivation of SIA U(2:J,2:K) = U(2:J,2:K) + ... shelves and mu_y * Dup .* ( Ushift(2:J,3:K+1) - Ushift(2:J,2:K) ) - ... streams mu_y * Ddown .* ( Ushift(2:J,2:K) - Ushift(2:J,1:K-1) ) + ... shallow shelf mu_x * Dright .* ( Ushift(3:J+1,2:K) - Ushift(2:J,2:K) ) - ... approximation (SSA) mu_x * Dleft .* ( Ushift(2:J,2:K) - Ushift(1:J-1,2:K) ); ﬂow-line ice shelf numerical solution t = t + dt; count = count + 1; scaling up end dtav = tf / count; next steps practicalities http://www.dms.uaf.edu/~bueler/karthaus/mfiles/diffusion.m omitted models • solves abstract diffusion equation ut = · (D(x, y ) u) • user must supply diffusivity on staggered grid Numerical modelling on exact solutions Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) • I am not quite done with the (heat) ↔ (SIA) analogy heat analogy & numerics • . . . I also want to use it to address exact solutions and analogy ﬁnite differences veriﬁcation exact solutions shallow ice sheets • the two senses of “solution”: solving the SIA is it right? ◦ a solution u(t, x, y ) of the heat equation is a function of time application the most basic and space for which ut = M + · (D u) is true shallow assumption ﬂowline Stokes ◦ a solution of the heat equation is a prediction of that model derivation of SIA • exact solutions are exact predictions of the model, but not shelves and streams exact predictions about nature shallow shelf approximation (SSA) ﬂow-line ice shelf • if we can only approximately ﬁnd solutions of a model then numerical solution scaling up knowing a few exact solutions can help test and maintain the next steps quality of the actual code that does the approximation . . . this practicalities omitted models is veriﬁcation [Wesseling, 2001] Numerical modelling exact solutions of heat equation Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) • many solutions to the heat equation are known, but one is heat analogy & numerics “fundamental” to the time-dependent equation, namely the analogy ﬁnite differences Green’s function exact solutions shallow ice sheets solving the SIA is it right? application the most basic • as time goes it changes shape shallow assumption ﬂowline Stokes by shrinking the output derivation of SIA (vertical) axis and shelves and simultaneously lengthening the streams shallow shelf input (horizontal) axis approximation (SSA) ﬂow-line ice shelf • . . . but otherwise it is the same numerical solution scaling up shape, so it is self-similar next steps • there are solutions of the SIA increasing time → practicalities omitted models just like this Numerical modelling ﬁnding the Green’s function of the heat equation Ed Bueler ice ﬂow: a superﬁcial view shallow ice • the best-known way to ﬁnd the Green’s function of the heat approximation (SIA) heat analogy & equation is by Fourier transform, but we use a method which numerics analogy generalizes to the SIA ﬁnite differences exact solutions • facts used to ﬁnd it: shallow ice ◦ the Green’s function starts at time t = 0 with a delta function at sheets solving the SIA the origin r = 0 is it right? ◦ the Green’s function is angularly-symmetric: u is a function of application the most basic the polar coordinate r = x 2 + y 2 shallow assumption ﬂowline Stokes ◦ calculus says: if f = f (r ) then 2 f = r −1 (rf ) derivation of SIA ◦ thus the heat equation in 2D, without additional heat sources, shelves and streams with constant diffusivity D > 0, and for a function of r , is: shallow shelf approximation (SSA) ﬂow-line ice shelf ut = Dr −1 (rur )r numerical solution scaling up • we use above facts to get Green’s function by its “similarity” next steps practicalities properties omitted models • but we do not use the linearity of the heat equation Numerical modelling ﬁnding a “similarity” solution of heat equation Ed Bueler ice ﬂow: a superﬁcial view • function u(t, r ) is replaced by function φ of one new variable s shallow ice approximation (SIA) • input scaling: s = t −β r heat analogy & numerics • output scaling: u = t −α φ(s) analogy ﬁnite differences • chain rule says: exact solutions ut = −αt −α−1 φ − βt −α−β−1 r φ , ur = t −α−β φ , shallow ice sheets etc. solving the SIA is it right? application • so heat equation ut = Dr −1 (rur )r is replaced by: the most basic shallow assumption ﬂowline Stokes derivation of SIA −αφ − βsφ = Ds−1 t −2β+1 φ + sφ shelves and streams shallow shelf • choose β = 1/2 so equation has no “t” and further simplify to approximation (SSA) ﬂow-line ice shelf numerical solution 1 scaling up − 2αsφ + s2 φ = D sφ 2 next steps practicalities omitted models • choose α = 1 so that quantity in square brackets simpliﬁes to a derivative: 2αsφ + s2 φ = s2 φ Numerical modelling ﬁnding a “similarity” solution of heat equation 2 Ed Bueler ice ﬂow: a • simplify and integrate, choose C = 0, integrate again: superﬁcial view shallow ice 1 approximation (SIA) − s2 φ = D s φ + C heat analogy & 2 numerics s analogy φ =− φ ﬁnite differences 2D 2 φ(s) = Ae−s /(4D) exact solutions shallow ice sheets solving the SIA • return to original variables, recalling s = t −β r : is it right? 2 2 application +y 2 )/(4Dt) the most basic u(t, r ) = A t −1 e−r /(4Dt) = A t −1 e−(x shallow assumption ﬂowline Stokes derivation of SIA • it is a spreading gaussian distribution shelves and streams • heat equation is linear so all time-dependent solutions of the heat shallow shelf approximation (SSA) equation are convolutions with this solution . . . which is why it is ﬂow-line ice shelf fundamental numerical solution scaling up next steps practicalities omitted models Numerical modelling similarity solutions Ed Bueler ice ﬂow: a superﬁcial view shallow ice • conclusion from previous slides: similarity variables for heat approximation (SIA) heat analogy & equation are numerics analogy input scaling output scaling ﬁnite differences −β exact solutions s = t x, u(t, x) = t −α φ(s) shallow ice sheets solving the SIA is it right? • dimension dependence: application 1D 2D 3D the most basic shallow assumption ﬂowline Stokes input scaling (t −β ) t −1/2 t −1/2 t −1/2 derivation of SIA output scaling (t −α ) t −1/2 t −1 t −3/2 shelves and streams shallow shelf historical note: In 1905 Einstein saw approximation (SSA) ﬂow-line ice shelf that the average distance traveled in numerical solution time t by particles in thermal motion √ x scaling up next steps goes like t. This is a microscopic practicalities omitted models explanation of the similarity variable s = t −1/2 x. t Numerical modelling Outline Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & numerics 1 ice ﬂow: an outsider’s superﬁcial view analogy ﬁnite differences exact solutions shallow ice 2 heat analogy & numerics sheets solving the SIA is it right? application 3 shallow ice sheets the most basic shallow assumption ﬂowline Stokes derivation of SIA 4 shelves and streams shelves and streams shallow shelf approximation (SSA) 5 next steps ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models Numerical modelling similarity solution to the SIA Ed Bueler ice ﬂow: a • jump forward to 1981 . . . breaking news: P. Halfar ﬁnds the superﬁcial view shallow ice fundamental solution of the SIA! [Halfar, 1981; 1983] approximation (SIA) heat analogy & • when n = 3, Halfar’s 2D solution for Glen ﬂow law has numerics analogy scalings ﬁnite differences exact solutions H(t, r ) = t −1/9 φ(s), s = t −1/18 r shallow ice sheets • . . . therefore, if an ice sheet starts steep and it ﬂows without solving the SIA additional accumulation then its thinning really slows down as is it right? application the shape ﬂattens out the most basic shallow assumption ﬂowline Stokes derivation of SIA 7000 7000 7000 7000 7000 shelves and 6000 6000 6000 6000 6000 streams shallow shelf 5000 5000 5000 5000 5000 approximation (SSA) ﬂow-line ice shelf 4000 4000 4000 4000 4000 numerical solution 3000 3000 3000 3000 3000 scaling up 2000 2000 2000 2000 2000 next steps practicalities 1000 1000 1000 1000 1000 omitted models 0 0 0 0 0 0 300 600 900 0 300 600 900 0 300 600 900 0 300 600 900 0 300 600 900 times 1, 10, 100, 1000, 10000 years Numerical modelling similarity solution to the SIA 2 Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & numerics analogy ﬁnite differences exact solutions shallow ice sheets solving the SIA is it right? application the most basic shallow assumption ﬂowline Stokes derivation of SIA shelves and streams shallow shelf approximation (SSA) ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models frames from t = 4 months to t = 106 years, equal spaced in exponential time Numerical modelling similarity solution to the SIA 3 Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) • for n = 3 the solution is: heat analogy & numerics 3/7 analogy 1/9 1/18 4/3 ﬁnite differences t0 1 − t0 r exact solutions H(t, r ) = H0 shallow ice t t R0 sheets solving the SIA is it right? application where 3 4 the most basic shallow assumption 1 7 R0 ﬂowline Stokes t0 = 7 derivation of SIA 18Γ 4 H0 shelves and streams and H0 , R0 are central height and ice cap radius at t = t0 shallow shelf approximation (SSA) ﬂow-line ice shelf • you choose H0 and R0 , and these determine the numerical solution “characteristic” time t0 scaling up next steps more on this in [Nye, 2000] and [Bueler and others, 2005] practicalities omitted models Numerical modelling ﬁnding Halfar solution Ed Bueler the method is the same as that which found the fundamental solution of ice ﬂow: a the heat equation: superﬁcial view shallow ice • allow scaling of the input and output: s = t −β r , H(t, r ) = t −α ϕ(s) approximation (SIA) heat analogy & • insert into Ht = 5 · ΓH | H| 2 H (← n = 3 case for simplicity) numerics analogy • to eliminate t from the equation, choose 7α + 4β = 1 ; get ﬁnite differences exact solutions −1 5 2 shallow ice −αφ − βsϕ = s Γsϕ |ϕ | ϕ sheets solving the SIA is it right? • to write as equality of derivatives, choose −α + 2β = 0 : application the most basic shallow assumption 2 5 2 −βs ϕ = Γsϕ |ϕ | ϕ ﬂowline Stokes derivation of SIA shelves and • integrate this and use ﬁniteness of ϕ(0), ϕ (0), and use ϕ (s) ≤ 0, to get separable streams ODE 4 3 shallow shelf βs = Γϕ (−ϕ ) approximation (SSA) ﬂow-line ice shelf • suppose H(t0 , 0) = H0 for some t0 > 0 to ﬁnd formula numerical solution scaling up 1/3 7/3 α 7/3 7 β 4/3 next steps ϕ(s) = (t0 H0 ) − s practicalities 4 Γ omitted models • write as formula for H(t, r ), on previous slide, by solving boxed equations above to get α = 1/9, β = 1/18 • −β margin location equation φ(t0 R0 ) = 0 gives formula for t0 on previous slide Numerical modelling evaluating Halfar solution Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) function H = halfar(t,x,y) heat analogy & numerics analogy g = 9.81; rho = 910.0; secpera = 31556926; ﬁnite differences n = 3; exact solutions A = 1.0e-16/secpera; Gamma = 2 * A * (rho * g)^3 / 5; shallow ice sheets H0 = 3600; R0 = 750e3; solving the SIA alpha = 1/9; beta = 1/18; is it right? t0 = (beta/Gamma) * (7/4)^3 * (R0^4/H0^7); application the most basic shallow assumption r = sqrt(x.*x + y.*y); ﬂowline Stokes r = r / R0; t = t / t0; derivation of SIA inside = max(0, 1 - (r / t^beta).^((n+1) / n)); shelves and H = H0 * inside.^(n / (2*n+1)) / t^alpha; streams shallow shelf http://www.dms.uaf.edu/~bueler/karthaus/mfiles/halfar.m approximation (SSA) ﬂow-line ice shelf numerical solution scaling up • it is easy to compute it, so . . . next steps practicalities • please use it to verify SIA codes, o.k.? omitted models Numerical modelling is the Halfar solution good for any modeling? Ed Bueler ice ﬂow: a superﬁcial view • Nye and others [2000] compared the long-time shallow ice approximation (SIA) consequences of different ﬂow laws for the South Mars Polar heat analogy & numerics Cap analogy ﬁnite differences • they compared the long-time behavior of the corresponding exact solutions Halfar solutions for the different ﬂow laws, with both CO2 ice shallow ice sheets and H2 O ice parameters solving the SIA is it right? • conclusions: application the most basic . . . none of the three possible [CO2 ] ﬂow laws will shallow assumption ﬂowline Stokes allow a 3000-m cap, the thickness suggested by derivation of SIA stereogrammetry, to survive for 107 years, indicating shelves and streams that the south polar ice cap is probably not shallow shelf approximation (SSA) composed of pure CO2 ice . . . the south polar cap ﬂow-line ice shelf numerical solution probably consists of water ice, with an unknown scaling up 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 ﬂow: a superﬁcial view shallow ice • recall that the SIA is approximation (SIA) heat analogy & Ht = M + · (D h) where D = ΓH n+2 | h|n−1 numerics analogy ﬁnite differences • “degeneracy” of D happens in two ways: exact solutions shallow ice sheets why D → 0 so what? solving the SIA is it right? H and h are continuous domes h→0 application the most basic but 2 h is singular shallow assumption H is continuous ﬂowline Stokes margins H→0 derivation of SIA but h is singular shelves and streams shallow shelf approximation (SSA) • numerically, margins are worse than domes ﬂow-line ice shelf numerical solution • degenerate diffusion equations are free boundary problems scaling up next steps • transformation η = H (2n+2)/n is known to help . . . for ﬂat beds practicalities omitted models • as a practical matter, numerical solvers for the SIA must handle ut = M + · (D u) for very general diffusivity Numerical modelling computing diffusivity in SIA Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & numerics analogy • for numerical stability we want ﬁnite differences exact solutions D = ΓH n+2 | h|n−1 on the shallow ice staggered grid sheets k+1 solving the SIA • various schemes proposed; stencil is it right? application for Mahaffy [1976] scheme at right k the most basic shallow assumption • all schemes involve ﬂowline Stokes derivation of SIA ◦ averaging H k-1 shelves and ◦ differencing h streams j-1 j j+1 shallow shelf ◦ in a “balanced” way, for better approximation (SSA) ﬂow-line ice shelf accuracy, numerical solution scaling up to put D on the staggered grid next steps practicalities omitted models Numerical modelling SIA implementation: ﬂat bed case Ed Bueler function [H,dtlist] = siaflat(Lx,Ly,J,K,H0,deltat,tf) ice ﬂow: a superﬁcial view shallow ice g = 9.81; rho = 910.0; secpera = 31556926; approximation (SIA) A = 1.0e-16/secpera; Gamma = 2 * A * (rho * g)^3 / 5; heat analogy & H = H0; numerics analogy dx = 2 * Lx / J; dy = 2 * Ly / K; ﬁnite differences N = ceil(tf / deltat); deltat = tf / N; exact solutions j = 2:J; k = 2:K; shallow ice nk = 3:K+1; sk = 1:K-1; ej = 3:J+1; wj = 1:J-1; sheets solving the SIA is it right? t = 0; dtlist = []; application for n=1:N the most basic Hup = 0.5 * ( H(j,nk) + H(j,k) ); Hdn = 0.5 * ( H(j,k) + H(j,sk) ); shallow assumption Hrt = 0.5 * ( H(ej,k) + H(j,k) ); Hlt = 0.5 * ( H(j,k) + H(wj,k) ); ﬂowline Stokes derivation of SIA a2up = (H(ej,nk) + H(ej,k) - H(wj,nk) - H(wj,k)).^2 / (4*dx)^2 + ... (H(j,nk) - H(j,k)).^2 / dy^2; shelves and streams a2dn = (H(ej,k) + H(ej,sk) - H(wj,k) - H(wj,sk)).^2 / (4*dx)^2 + ... shallow shelf (H(j,k) - H(j,sk)).^2 / dy^2; approximation (SSA) a2rt = (H(ej,k) - H(j,k)).^2 / dx^2 + ... ﬂow-line ice shelf (H(ej,nk) + H(j,nk) - H(ej,sk) - H(j,sk)).^2 / (4*dy)^2; numerical solution a2lt = (H(j,k) - H(wj,k)).^2 / dx^2 + ... scaling up (H(wj,nk) + H(j,nk) - H(wj,sk) - H(j,sk)).^2 / (4*dy)^2; next steps Dup = Gamma * Hup.^5 .* a2up; Ddn = Gamma * Hdn.^5 .* a2dn; practicalities Drt = Gamma * Hrt.^5 .* a2rt; Dlt = Gamma * Hlt.^5 .* a2lt; omitted models [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 ﬂow: a superﬁcial view shallow ice • adaptive time-stepping works: siaflat.m takes short time approximation (SIA) steps when the driving stress ρgH| h| is large, and then heat analogy & numerics longer steps when it is “easier” analogy ﬁnite differences • for example, a worst-case ice sheet is thick but with a very exact solutions bumpy surface shallow ice sheets solving the SIA • roughice.m generates initial ice sheet, which is thick and is it right? with a white-noise surface, then runs for 50 years application the most basic shallow assumption • initial surface, ﬁnal surface, and time-steps shown below ﬂowline Stokes derivation of SIA shelves and streams shallow shelf approximation (SSA) ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models Numerical modelling veriﬁcation of numerical ice ﬂow codes Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) • how do we make sure a numerical scheme and its heat analogy & numerics implementation are right? some available techniques analogy 1 don’t make any mistakes ﬁnite differences exact solutions 2 “eyeball” the results when you change the code . . . if it looks shallow ice right it is right sheets solving the SIA 3 from the beginning, build-in a comparison to an exact solution is it right? application • where to get exact solutions for ice ﬂow models? the most basic shallow assumption ﬂowline Stokes similarity solutions to SIA [Halfar, 1983; Bueler and others, 2005] derivation of SIA manufactured solutions to SIA; some include [Bueler and others 2005; 2007] shelves and thermo-coupling streams SSA solutions in ﬂow line and [Schoof, 2006; van der Veen, 1985] shallow shelf approximation (SSA) cross-ﬂow cases ﬂow-line ice shelf Blatter solution in ﬂow line [Glowinski and Rappaz, 2003] numerical solution ﬂow line Stokes solutions in constant [Ladyzhenskaya, 1963] viscosity case ( 4 ψ = 0) scaling up next steps manufactured solutions to Stokes [Sargent and Fastook, 2010] practicalities equations with variable viscosity omitted models Numerical modelling verifying our SIA code siaflat.m Ed Bueler ice ﬂow: a superﬁcial view • this code calls shallow ice approximation (SIA) siaflat.m heat analogy & numerics • with initial condition function [avthkerr,maxthkerr] = verifysia(J,dtyears) analogy which is a t1 Halfar ﬁnite differences if nargin<1, J=40; end; exact solutions if nargin<2, dtyears=10.0; end; solution shallow ice sheets L = 1200e3; dx = 2 * L / J; • numerical code runs to [x,y] = meshgrid(-L:dx:L, -L:dx:L); solving the SIA later time t2 and the is it right? t1 = 200; t2 = 20000; secpera = 31556926; application H1 = halfar(t1 * secpera,x,y); compares to Halfar at t2 the most basic shallow assumption [H2approx,dtlist] = siaflat(L,L,J,J,H1,dtyears*secpera,... ﬂowline Stokes (t2-t1)*secpera); derivation of SIA map of thickness error from H2exact = halfar(t2 * secpera,x,y); » verifysia(160,3.0): shelves and error = H2approx - H2exact; streams shallow shelf avthkerr = sum(sum(abs(error)))/(J+1)^2; approximation (SSA) maxthkerr = max(max(abs(error))); ﬂow-line ice shelf erreta = H2approx.^(8/3) - H2exact.^(8/3); numerical solution maxeta = max(max(H2exact.^(8/3))); scaling up next steps http://www.dms.uaf.edu/~bueler/karthaus/mfiles/verifysia.m practicalities omitted models Numerical modelling verifying our SIA code 2 Ed Bueler ice ﬂow: a superﬁcial view >> verifysia(20,24.0); shallow ice approximation (SIA) errors for 20 x 20 grid: heat analogy & numerics average thickness error = 22.335 analogy maximum thickness error = 227.275 Trust but verify. ﬁnite differences exact solutions >> verifysia(40,12.0); (Ronald Reagan) shallow ice errors for 40 x 40 grid: sheets average thickness error = 9.519 solving the SIA is it right? maximum thickness error = 241.040 application the most basic >> verifysia(80,6.0); shallow assumption ﬂowline Stokes errors for 80 x 80 grid: derivation of SIA average thickness error = 2.937 shelves and maximum thickness error = 159.442 streams shallow shelf >> verifysia(160,3.0); approximation (SSA) ﬂow-line ice shelf errors for 160 x 160 grid: numerical solution average thickness error = 0.988 ﬁgure 2 in [Huybrechts and scaling up maximum thickness error = 103.456 others, 1996] next steps practicalities >> verifysia(200,2.0); omitted models errors for 200 x 200 grid: average thickness error = 0.827 maximum thickness error = 94.561 Numerical modelling numerical mass conservation Ed Bueler ice ﬂow: a superﬁcial view • in addition to measuring geometry errors, we might want a shallow ice approximation (SIA) numerical code not to lose mass! heat analogy & numerics • again, the Halfar solution is a perfect tool analogy ﬁnite differences ◦ because the Halfar exact solution has exact (continuum) exact solutions volume conservation shallow ice sheets • add numerical volume calculation to verifysia.m: solving the SIA is it right? >> verifysia(20) application vol1 = 3.96112407720492e+15 the most basic shallow assumption ... ﬂowline Stokes derivation of SIA vol2approx = 3.96112407720492e+15 shelves and voldiff = 1.50000000000000e+00 streams shallow shelf • huh? how can it be that good? approximation (SSA) ﬂow-line ice shelf ◦ the Mahaffy [1976] scheme has local numerical mass numerical solution scaling up conservation next steps ◦ but what about the margin where H → 0? practicalities omitted models ◦ siaflat.m implements boundary condition “H ≥ 0” ◦ apparently this boundary condition produces global numerical mass conservation to rounding error! Numerical modelling application to the Antarctic ice sheet Ed Bueler ice ﬂow: a • at Karthaus 2009, a computer project: modify siaflat.m to superﬁcial view shallow ice • allow arbitrary bed elevation b(x, y ) approximation (SIA) • allow arbitrary mass balance M(x, y ) heat analogy & numerics • apply a marine-margin condition analogy ﬁnite differences and then simulate the Antarctic ice sheet exact solutions shallow ice • my solution: sheets solving the SIA • add 8 lines to siaflat.m in total, creating siageneral.m is it right? • use the SeaRISE 5km data set for present-day Antarctica application the most basic (NetCDF) shallow assumption ﬂowline Stokes • add a page or so of code for pre- and post-processing derivation of SIA • do 40 ka run starting with present-day geometry shelves and streams initial (observed) surface elevation 40 ka (modeled) surface elevation 4000 4000 shallow shelf approximation (SSA) -2000 3500 -2000 3500 ﬂow-line ice shelf numerical solution 3000 3000 -1000 -1000 scaling up 2500 2500 y (km) y (km) next steps 0 2000 0 2000 practicalities 1000 1500 1000 1500 omitted models 1000 1000 2000 2000 500 500 3000 3000 0 0 -2000 -1000 0 1000 2000 3000 -2000 -1000 0 1000 2000 3000 x (km) x (km) Numerical modelling application to the Antarctic ice sheet 2 Ed Bueler ice ﬂow: a superﬁcial view 34 shallow ice approximation (SIA) • volume grows 32 heat analogy & numerics ◦ levels out at about 33 × 106 km3 ◦ compare present-day observed volume (10 km ) 3 analogy 30 volume of 25 × 106 km3 6 ﬁnite differences exact solutions shallow ice ◦ exponential time of about 15 ka, 28 sheets ◦ but exponential time is much solving the SIA 26 is it right? longer if our model were application the most basic thermocoupled, because of long 24 0 5000 10000 15000 20000 25000 30000 35000 40000 shallow assumption ﬂowline Stokes 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 ﬂow-line ice shelf numerical solution • an enhancement factor of 24 would make the volume level scaling up next steps out at about the present-day value practicalities omitted models • pop quiz: why is that not a good idea? Numerical modelling the most basic shallow assumption Ed Bueler ice ﬂow: a air superﬁcial view shallow ice approximation (SIA) heat analogy & • careful derivation of the SIA is next ice numerics analogy • makes “shallowness assumptions” ﬁnite differences exact solutions • the most basic is a limitation which shallow ice sheets applies to all shallow ice theoriesa bedrock solving the SIA but not to the general Stokes is it right? application theory: the most basic shallow assumption ﬂowline Stokes the surface and base of the ice are derivation of SIA given by differentiable functions shelves and streams z = h(t, x, y ) and z = b(t, x, y ) shallow shelf approximation (SSA) ﬂow-line ice shelf • the conﬁgurations at right, with numerical solution scaling up overhang in surface, are not next steps allowed practicalities a e.g. SIA, SSA, Blatter-Pattyn, SIA+SSA hybrids, . . . omitted models Numerical modelling kinematic and mass continuity equations Ed Bueler ice ﬂow: a superﬁcial view shallow ice • what does this “most basic shallow assumption” get you? approximation (SIA) heat analogy & • answer: a map-plane mass continuity equation numerics analogy • consider these three equations: ﬁnite differences exact solutions ◦ the surface kinematical equation shallow ice ◦ the base kinematical equation sheets solving the SIA ◦ the map-plane mass continuity equation is it right? application • all three equations are important for shallow ice sheets and the most basic shallow assumption shelves, but any two imply the third ﬂowline Stokes derivation of SIA • the key facts we need: shelves and ◦ the incompressibility of ice streams shallow shelf approximation (SSA) ux + vy + wz = 0 ﬂow-line ice shelf numerical solution scaling up ◦ the Leibniz rule for differentiating integrals next steps ˆ f (x) ˆ f (x) practicalities d omitted models h(x, y ) dy = f (x)h(x, f (x)) − g (x)h(x, g(x)) + hx (x, y ) dy dx g(x) g(x) Numerical modelling kinematic and mass continuity equations 2 Ed Bueler ice ﬂow: a superﬁcial view shallow ice • let a be the surface mass balance function (a > 0 is approximation (SIA) heat analogy & accumulation) numerics analogy • let s be the basal melt rate function (s > 0 is basal melting) ﬁnite differences exact solutions • deﬁne the map-plane ﬂux of ice, ˆ shallow ice sheets h solving the SIA is it right? q= ¯ ¯ (u, v ) dz = (u , v ) H application b the most basic shallow assumption ﬂowline Stokes • the three equations are: derivation of SIA shelves and streams shallow shelf surface kinematical ht = a − u h hx − v h hy + w h (3) approximation (SSA) ﬂow-line ice shelf numerical solution scaling up base kinematical bt = s − u b bx − v b by + w b (4) next steps practicalities omitted models mass continuity Ht = M − · q where M = a − s (5) Numerical modelling kinematic and mass continuity equations 3 Ed Bueler ice ﬂow: a for example, here is a sketch of how to show (3) & (4) =⇒ (5): superﬁcial view shallow ice approximation (SIA) 1 recalling H = h − b and M = a − s, the difference of (3) and (4) is heat analogy & numerics Ht = M − u h hx − v h hy + w h + u b bx + v b by − w b (*) analogy ﬁnite differences exact solutions 2 incompressibility gives difference of surface and base values of w by integration, shallow ice ˆ h sheets solving the SIA wz = −(ux + vy ) =⇒ w h −w b =− (ux + vy ) dζ is it right? b application the most basic which reduces (*) to: shallow assumption ﬂowline Stokes ˆ h derivation of SIA Ht = M − u h hx − v h hy + u b bx + v b by − (ux + vy ) dζ (**) shelves and b streams shallow shelf ´h approximation (SSA) 3 on the other hand, the Leibniz rule on q = b (u, v ) dz gives ﬂow-line ice shelf numerical solution ˆ h scaling up · q = u h hx + v h hy − u b bx − v b by + (ux + vy ) dz b next steps practicalities omitted models 4 . . . which reduces (**) to Ht = M − ·q Numerical modelling kinematic and mass continuity equations 4 Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & • papers have lots of calculations like on the previous slide numerics analogy • usually mixed in with other arguments about shallowness ﬁnite differences exact solutions • most ice sheet models use the mass continuity equation to shallow ice sheets describe change in ice sheet geometry, solving the SIA is it right? • . . . but they could instead use the surface kinematical application the most basic equation shallow assumption ﬂowline Stokes • in using a stress balance in a shallow theory, all we need to derivation of SIA shelves and do is get the horizontal velocity (u, v ), which gives q and thus streams the mass continuity equation shallow shelf approximation (SSA) ﬂow-line ice shelf • for example, the derivation of the SIA from Stokes needs to numerical solution scaling up only go far enough to ﬁnd an expression for the horizontal next steps velocity (u, v ) practicalities omitted models Numerical modelling origin of SIA Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & where does the “shallow ice approximation” come from? numerics analogy • historically: ﬁnite differences exact solutions • Fowler and Larson [1978], Morland and Johnson [1980], and shallow ice Hutter [1983] sheets • fairly recent compared to the ﬂow law of ice (1950s) and solving the SIA is it right? Stokes linearly-viscous model (∼1860) application the most basic • logically: shallow assumption ﬂowline Stokes • by a “small-parameter argument”—depth-to-width ratio is derivation of SIA shelves and small—to rescale the Stokes model streams • analogous to an old argument giving shallow water equations shallow shelf approximation (SSA) from Navier-Stokes ﬂow-line ice shelf numerical solution • more precisely: scaling up • as sketched in the next slides, for the x, z plane ﬂow case only next steps practicalities omitted models Numerical modelling Stokes equations in x, z plane ﬂow case Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & • suppose all ﬂow is in the (screen) x, z plane numerics analogy • . . . so out-of-plane velocity is zero (v = 0) and all quantities ﬁnite differences exact solutions are constant in y direction (∂/∂y = 0) shallow ice sheets • the symmetric strain rate tensor has many zeros solving the SIA is it right? D11 0 D13 application 1 ∂ui ∂uj the most basic shallow assumption Dij = + = 0 0 0 2 ∂xj ∂xi ﬂowline Stokes derivation of SIA D13 0 D33 shelves and streams • incompressibility of ice in plane ﬂow case says D has trace shallow shelf approximation (SSA) zero: ﬂow-line ice shelf numerical solution ux + 0 + wz = D11 + D33 = 0 scaling up next steps • . . . so D33 = −D11 in this case practicalities omitted models Numerical modelling Stokes equations in x, z plane ﬂow case 2 Ed Bueler ice ﬂow: a superﬁcial view • recall that the Cauchy stress tensor σij , with its pressure part shallow ice approximation (SIA) p removed, is called “deviatoric”: heat analogy & numerics analogy τij = σij + pδij where p = −(1/3)(σ11 + σ22 + σ33 ) ﬁnite differences exact solutions shallow ice • the deviatoric stress tensor τij is proportional to the tensor sheets solving the SIA Dij , in the isotropic (Glen law) case: is it right? Dij = A(T )τ n−1 τij application the most basic shallow assumption ﬂowline Stokes derivation of SIA 1/2 shelves and • . . . where 2τ 2 = 2(II τ )2 = τij τij so 2 2 τ = τ11 + τ13 streams shallow shelf • so τij is also symmetric, and has trace zero, approximation (SSA) ﬂow-line ice shelf • and τij has only two nonzero entries τ11 , τ13 : numerical solution scaling up next steps τ11 0 τ13 practicalities 0 0 0 omitted models τ13 0 −τ11 Numerical modelling Stokes equations in x, z plane ﬂow case 3 Ed Bueler ice ﬂow: a superﬁcial view • in the 3D case: n = 3, Glen, isothermal Stokes equations: shallow ice approximation (SIA) heat analogy & ux + vy + wz = 0 numerics analogy 0= · σ + ρg ﬁnite differences exact solutions Dij = Aτ 2 τij shallow ice sheets solving the SIA • x, z plane ﬂow case: from simpliﬁcations on previous slides, is it right? application the isothermal Stokes equations say the most basic shallow assumption ﬂowline Stokes derivation of SIA ux + wz = 0 shelves and streams px = τ11,x + τ13,z shallow shelf approximation (SSA) pz = τ13,x − τ11,z − ρg (6) ﬂow-line ice shelf numerical solution ux = Aτ 2 τ11 scaling up next steps uz + wx = 2Aτ 2 τ13 practicalities omitted models • ﬁve scalar equations in ﬁve scalar unknowns (u, w, p, τ11 , τ13 ) but still fairly hard to solve . . . Numerical modelling Stokes equations in x, z plane ﬂow case 4 Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) • previous equations apply within the ice ﬂuid heat analogy & numerics • . . . but boundary conditions are needed analogy ﬁnite differences • at the upper surface there is no stress from the air: σij nj = 0 exact solutions shallow ice where n is a normal vector to the surface sheets solving the SIA • in plane ﬂow case use n = (hx , 0, −1), which is gradient of is it right? application F (x, z) = h(x) − z, which is constant on z = h(x) the most basic shallow assumption • get stress equations at surface: ﬂowline Stokes derivation of SIA shelves and τ13 = (τ11 − p)hx streams shallow shelf approximation (SSA) p + τ11 + τ13 hx = 0 ﬂow-line ice shelf numerical solution scaling up • we will only consider non-sliding SIA, so base boundary next steps conditions are practicalities omitted models u=0 and w =0 Numerical modelling scaling the variables Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) • following Chapter 18 of [Fowler, 1997], Mathematical Models heat analogy & in the Applied Sciences numerics • which does a much more general job than here: 3D case with conservation of analogy energy, while we only do 2D isothermal case ﬁnite differences exact solutions • we also miss Andrew standing on a table and singing, too shallow ice sheets • consider ice sheet of typical depth d and length , and with solving the SIA is it right? velocity ([u]), deviatoric stress ([τ ]), and softness ([A]) scales application the most basic • = d/ is the “aspect ratio” of the ice sheet shallow assumption ﬂowline Stokes • change to new starred variables using scales: derivation of SIA shelves and streams x = x∗ u = [u]u ∗ shallow shelf approximation (SSA) z = d z∗ w = [u]w ∗ ﬂow-line ice shelf ∗ ∗ numerical solution h=dh τ11 = [τ ]τ11 scaling up next steps H = d H∗ ∗ τ13 = [τ ]τ13 ∗ practicalities omitted models A = [A]A p − ρg(h − z) = [τ ]p∗ Numerical modelling scaling the variables 2 Ed Bueler ice ﬂow: a superﬁcial view • we will necessarily choose relations among the scales shallow ice approximation (SIA) • the choice of such relations affects the outcome and is heat analogy & numerics subject to “scientiﬁc judgment” (i.e. controversy) analogy ﬁnite differences • “. . . the purpose is to attain ‘properly scaled’ equations in exact solutions shallow ice which the largest dimensionless parameters are numerically sheets solving the SIA of order one . . . ” [Fowler, 1997] is it right? application • method: we rewrite the Stokes equations and the boundary the most basic shallow assumption conditions in the newly (starred) variables, then check that ﬂowline Stokes derivation of SIA the scaling relations give reasonable values, then remove shelves and terms to get a new, more-solvable model, the SIA streams shallow shelf • the scaling relations used here are: approximation (SSA) ﬂow-line ice shelf numerical solution balance of shear stress scaling up [τ ] = ρgd next steps and pressure gradient practicalities [u] omitted models balance of ﬂow law scales = 2[A][τ ]3 d Numerical modelling scaling the equations Ed Bueler example: ice ﬂow: a superﬁcial view • recall this equation from Stokes equations (6): shallow ice approximation (SIA) heat analogy & px = τ11,x + τ13,z numerics analogy ﬁnite differences • replace all variables with their starred versions: exact solutions shallow ice [τ ] ∗ dρg ∗ [τ ] ∗ [τ ] ∗ sheets px ∗ + hx ∗ = τ11,x ∗ + τ13,z ∗ solving the SIA d is it right? application • multiply by /( [τ ]): the most basic shallow assumption ﬂowline Stokes ∗ dρg ∗ ∗ ∗ derivation of SIA px ∗ + hx ∗ = τ11,x ∗ + τ13,z ∗ [τ ] d shelves and streams shallow shelf • use relation [τ ] = ρgd and = d/ : approximation (SSA) ﬂow-line ice shelf ∗ −2 ∗ ∗ −2 ∗ numerical solution px ∗ + hx ∗ = τ11,x ∗ + τ13,z ∗ scaling up next steps practicalities • re-arrange to taste and remove stars: omitted models 2 hx = τ13,z + (τ11,x − px ) • . . . do this for each of the Stokes equations Numerical modelling scaling the equations 2 Ed Bueler ice ﬂow: a superﬁcial view shallow ice • Stokes system now appears this way (stars are removed): approximation (SIA) heat analogy & numerics ux + wz = 0 analogy 2 ﬁnite differences hx = τ13,z + (τ11,x − px ) exact solutions shallow ice pz = τ13,x − τ11,z sheets solving the SIA 1 is it right? 2 2 ux = A 2 τ11 + τ13 τ11 application 2 the most basic shallow assumption ﬂowline Stokes 2 2 uz + 2 wx = A 2 τ11 + τ13 τ13 derivation of SIA shelves and • this is still full Stokes! streams shallow shelf approximation (SSA) • surface conditions are ﬂow-line ice shelf numerical solution 2 scaling up τ13 = (τ11 − p) hx next steps practicalities p + τ11 + τ13 hx = 0 omitted models • base conditions are u = 0 and w = 0 as before Numerical modelling checking the scales Ed Bueler ice ﬂow: a superﬁcial view shallow ice • the scale for vertical velocity is already chosen as [u] approximation (SIA) heat analogy & • . . . and we now identify this scale with a typical value for numerics analogy accumulation: [a] = [u] ﬁnite differences exact solutions • thus we have chosen these scale relations: shallow ice sheets d solving the SIA = [τ ] = ρgd is it right? application [u] the most basic shallow assumption = 2[A][τ ]3 [a] = [u] ﬂowline Stokes d derivation of SIA shelves and • to check their reasonableness, we will determine other scale streams shallow shelf factors in terms of , [a], and [A], which are taken from approximation (SSA) ﬂow-line ice shelf observations: numerical solution scaling up = 3000 km, next steps practicalities [a] = 0.1 m a−1 omitted models [A] = 2 × 10−16 Pa−3 a−1 Numerical modelling checking the scales 2 Ed Bueler ice ﬂow: a superﬁcial view shallow ice • solving for d, , [τ ], [u] gives approximation (SIA) heat analogy & numerics 4 1/8 [a] analogy d= ≈ 3600 m ﬁnite differences exact solutions 2[A](ρg)3 shallow ice d sheets = ≈ 10−3 solving the SIA is it right? application [τ ] = ρgd ≈ 4 × 104 Pa = 0.4 bar the most basic shallow assumption ﬂowline Stokes [u] = 2[A]d[τ ]3 ≈ 80 m a−1 derivation of SIA shelves and streams • these are all reasonable values for large ice sheets shallow shelf approximation (SSA) • it is reasonable to proceed to delete terms from the scaled ﬂow-line ice shelf numerical solution equations to get the reduced (shallow) theory, the SIA scaling up • the actual measure of the quality of SIA is the comparison of next steps practicalities its predictions to observations, when modeling various real omitted models ice sheets etc., not these scalings Numerical modelling ﬁnally get to the SIA Ed Bueler ice ﬂow: a • setting = 0 in the scaled Stokes equations, get these equations superﬁcial view (left column) and boundary conditions (right): shallow ice approximation (SIA) ∗ ux + wz = 0 τ13 h =0 heat analogy & numerics ∗ analogy hx = τ13,z p h + τ11 h + τ13 h hx = 0 ﬁnite differences exact solutions pz = τ13,x − τ11,z shallow ice sheets 1 2 † ux = Aτ13 τ11 u b =0 solving the SIA is it right? 2 † 3 application the most basic uz = Aτ13 w b =0 shallow assumption ﬂowline Stokes • combine the two marked “∗”, integrating vertically from surface derivation of SIA shelves and z = h to arbitrary elevation z, to get streams shallow shelf τ13 = −(h − z)hx approximation (SSA) ﬂow-line ice shelf numerical solution • combine the two marked “†”, integrating vertically from base z = b scaling up to arbitrary level z, and use the new expression for τ13 , to get next steps ˆ z ˆ z 3 τ13 dζ = −A(hx )3 (h − ζ)3 dζ practicalities omitted models u(z) = A b b A = − (hx )3 H 4 − (h − z)4 4 Numerical modelling ﬁnally get to the SIA 2 Ed Bueler ice ﬂow: a superﬁcial view shallow ice • remember I said that it sufﬁces to ﬁnd an expression for the approximation (SIA) horizontal velocity, in order to state the SIA? we are there! heat analogy & numerics • the previous slides have been in scaled quantities; now we recall the analogy ﬁnite differences stars: A∗ ∗ u ∗ = − (hx ∗ )3 (H ∗ )4 − (h∗ − z ∗ )4 exact solutions shallow ice 4 sheets solving the SIA • and then unscale to ﬁnd a dimensional expression for u, using is it right? application ∗ hx ∗ = −1 hx , and using scaling relations: the most basic shallow assumption [u]A u = [u]u ∗ = − −3 (hx )3 d −4 H 4 − (h − z)4 ﬂowline Stokes derivation of SIA 4[A] shelves and streams [u] shallow shelf =− A(hx )3 H 4 − (h − z)4 approximation (SSA) 4[A] 3 d 4 ﬂow-line ice shelf numerical solution 2[A]d(ρgd )3 scaling up =− A(hx )3 H 4 − (h − z)4 next steps 4[A] 3 d 4 practicalities 1 omitted models = − A(ρg)3 (hx )3 H 4 − (h − z)4 2 Numerical modelling ﬁnally get to the SIA 3 Ed Bueler ice ﬂow: a superﬁcial view ´h shallow ice approximation (SIA) • recall that the ﬂux of ice is q = b u(z) dz heat analogy & • so we integrate the expression for u(z): numerics analogy ﬁnite differences ˆ h 1 exact solutions q = − A(ρg)3 (hx )3 H 4 − (h − z)4 dz shallow ice 2 b sheets solving the SIA 1 1 is it right? = − A(ρg)3 (hx )3 H 5 − (H 5 − 0) application 2 5 the most basic shallow assumption 2 ﬂowline Stokes = − A(ρg)3 (hx )3 H 5 derivation of SIA 5 shelves and streams shallow shelf • the mass continuity equation now says: approximation (SSA) ﬂow-line ice shelf numerical solution 2 scaling up Ht = M − qx = M + A(ρg)3 (hx )3 H 5 (∗) next steps 5 x practicalities omitted models • equation (∗) matches equation (1) in the plane ﬂow case Numerical modelling have I oversold diffusivity? Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & • I have asserted that the default model for ice sheets, the SIA, numerics analogy is a diffusion ﬁnite differences exact solutions • recall the analogy: shallow ice sheets solving the SIA diffusion eqn ↔ SIA is it right? application ut = · (D u) ↔ ht = M + · (D h) the most basic shallow assumption D = D(x, y ) ↔ D = ΓH n+2 | h|n−1 ﬂowline Stokes derivation of SIA shelves and • have I oversold this diffusivity analogy? streams shallow shelf ◦ probably approximation (SSA) ﬂow-line ice shelf ◦ . . . and I’ve acknowledged there are “issues” numerical solution scaling up • but it turns out there is some robustness to my analogy next steps ◦ as follows on the next two slides practicalities omitted models Numerical modelling have I oversold diffusivity? 2 Ed Bueler DIFFUSIVE IDEA 1: rough beds have the effect of reducing ice ﬂow: a superﬁcial view diffusivity shallow ice approximation (SIA) • deﬁne the local bed topography by removing the local mean heat analogy & numerics over some range λ ≈ 10 km: analogy ﬁnite differences λ exact solutions ˜ b(x, ξ) = b(x + ξ) − b(x + ξ) dξ shallow ice −λ sheets solving the SIA • deﬁne this strange average of the local bed: is it right? application −1/n the most basic −(n+2)/n shallow assumption λ ˜ b(x, ξ) ﬂowline Stokes θ(x) = 1− dξ derivation of SIA −λ H(x) shelves and streams shallow shelf • using a multiple-scales analysis, Schoof [2003] says you will approximation (SSA) ﬂow-line ice shelf get closer to solving the Stokes equations by making these numerical solution scaling up two modiﬁcations of the SIA: next steps ◦ smooth the bed, (modelers want to do that anyway) practicalities omitted models ◦ but don’t lose track of the smoothed-away local bed roughness; use it to reduce the diffusivity: Dnew = θD where 0 ≤ θ ≤ 1 Numerical modelling have I oversold diffusivity? 3 Ed Bueler ice ﬂow: a superﬁcial view shallow ice DIFFUSIVE IDEA 2: the large-scale effect of sliding in ice approximation (SIA) streams (addressed in next section), is diffusive heat analogy & numerics • suppose that, for an ice stream modeled by the SSA equation analogy ﬁnite differences exact solutions shallow ice 2A−1/n H|ux |1/n−1 ux − C|u|m−1 u = ρgHhx sheets x solving the SIA is it right? we assume that the basal resistance term balances the application the most basic driving stress: shallow assumption ﬂowline Stokes −C|u|m−1 u = ρgHhx derivation of SIA shelves and • then the ice stream geometry evolves by a diffusion, streams shallow shelf approximation (SSA) ﬂow-line ice shelf ht = M + · (D h) where D = C H (1/m)+1 | h|(1/m)−1 numerical solution scaling up • this “outer problem” is part of a matched asymptotic next steps practicalities expansion that does a good job of tracking the grounding line omitted models in a marine ice sheet [Schoof, 2007] Numerical modelling Outline Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & numerics 1 ice ﬂow: an outsider’s superﬁcial view analogy ﬁnite differences exact solutions shallow ice 2 heat analogy & numerics sheets solving the SIA is it right? application 3 shallow ice sheets the most basic shallow assumption ﬂowline Stokes derivation of SIA 4 shelves and streams shelves and streams shallow shelf approximation (SSA) 5 next steps ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models Numerical modelling the shallow shelf approximation (SSA) stress balance Ed Bueler ice ﬂow: a superﬁcial view shallow ice is a model which applies well to large ice shelves approximation (SIA) • . . . for parts away from grounding lines heat analogy & numerics analogy • . . . and away from calving fronts ﬁnite differences exact solutions shallow ice sheets solving the SIA is it right? application the most basic shallow assumption ﬂowline Stokes derivation of SIA shelves and streams shallow shelf approximation (SSA) ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models edge of Ekström ice shelf, photo Hans Grobe, Polarstern expedition ANT-XX/2 Numerical modelling SSA stress balance 2 Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & numerics analogy ﬁnite differences exact solutions SSA also applies reasonably well to shallow ice ice streams sheets solving the SIA • . . . with modest bed topography is it right? application • . . . and weak bed strength the most basic shallow assumption ﬂowline Stokes derivation of SIA shelves and streams shallow shelf approximation (SSA) ﬂow-line ice shelf numerical solution scaling up comment: energy conservation issues—both ice temperature and basal next steps melt rate—are a major part of ice stream ﬂow modeling [Raymond, practicalities omitted models 2000], but they are not addressed in these lectures Numerical modelling what is, and is not, an ice stream? Ed Bueler ice ﬂow: a superﬁcial view • ice streams have fast ﬂow (50 to shallow ice approximation (SIA) > 1000 m a−1 ) by sliding, with heat analogy & ◦ low bed resistance and numerics analogy ◦ a critical role of liquid water at ﬁnite differences exact solutions bed [Clarke, 2005] shallow ice • “outlet glaciers” may have the sheets solving the SIA same properties, but they have is it right? application substantial vertical shear and the most basic shallow assumption proportionally less sliding ﬂowline Stokes derivation of SIA • and they have lower aspect ratio shelves and streams than “true” ice streams shallow shelf approximation (SSA) • and some outlet glaciers are really Cross sections of Jakobshavns Isbrae (a) and Whillans Ice Stream ﬂow-line ice shelf numerical solution fast, e.g. Jakobshavn Isbrae (b). Plotted without vertical scaling up exaggeration in order to better • thus: few simplifying assumptions illustrate the difference between the next steps two types. (Figure 1 in [Truffer and Echelmeyer, practicalities are possible for outlet glaciers 2003]) omitted models • not clear you should use SSA for outlet glaciers Numerical modelling stream-to-shelf ﬂow line: notation again Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & numerics analogy ﬁnite differences exact solutions shallow ice sheets solving the SIA is it right? application the most basic shallow assumption ﬂowline Stokes derivation of SIA shelves and streams shallow shelf ﬁgure modiﬁed from [Schoof, 2007] approximation (SSA) ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models Numerical modelling the SSA equation Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & • we will only work in plane ﬂow case numerics analogy • we will numerically solve only the stress balance equation ﬁnite differences exact solutions which determines the velocity in an ice shelf or ice stream: shallow ice sheets solving the SIA is it right? 2A−1/n H|ux |1/n−1 ux − C|u|m−1 u − ρgHhx = 0 (7) x application the most basic shallow assumption ﬂowline Stokes ◦ derived originally by Morland [1987], MacAyeal [1989] derivation of SIA ◦ references: Weis and others [1999], Schoof [2006; 2007] shelves and streams shallow shelf • the blue term inside parentheses is the vertically-integrated approximation (SSA) ﬂow-line ice shelf “longitudinal” or “membrane” stress numerical solution scaling up • how to solve this equation numerically? next steps practicalities • how to think about it? omitted models Numerical modelling the full monty, with a grounding line Ed Bueler here is a full ﬂow line context: ice ﬂow: a superﬁcial view shallow ice u(0) = (given) approximation (SIA) heat analogy & numerics 2A−1/n H|ux |1/n−1 ux − C|u|m−1 u − ρgHhx = 0 x analogy ﬁnite differences Ht = M − (uH)x on 0 < x < xg exact solutions ρH ≥ ρw (−b) shallow ice sheets h =H +b solving the SIA is it right? h, u, ux continuous at x = xg application the most basic shallow assumption 2A−1/n H|ux |1/n−1 ux − ρg(1 − ρ/ρw )HHx = 0 ﬂowline Stokes x derivation of SIA on xg < x < xc Ht = M − (uH)x shelves and streams ρH < ρw (−b) shallow shelf approximation (SSA) 2A−1/n H|ux |1/n−1 ux − 2 ρ(1 − ρ/ρw )gH 2 = 0 1 ﬂow-line ice shelf at x = xc numerical solution (a condition describing movement of xc ) scaling up next steps practicalities • this is the default model in the Marine Ice Sheet Model omitted models Intercomparison Project [MISMIP; Schoof and others, 2008] • effectively an open problem: what is best numerical treatment of grounding line movement in above model? Numerical modelling exact thickness and velocity for steady ice shelf Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & • we need a more limited goal! numerics analogy • ﬁrst goal: solve the steady state ice shelf in the isothermal, ﬁnite differences exact solutions plane ﬂow, and Glen law case shallow ice sheets • for this nice situation, there is a nice result: the thickness and solving the SIA velocity in the ice shelf can be completely determined5 from: is it right? application 1 ice thickness Hg at the grounding line, the most basic shallow assumption 2 ice velocity ug at the grounding line, ﬂowline Stokes derivation of SIA 3 an integrable (e.g. constant) surface balance M(x) shelves and • we do this by-hand on the next slide streams shallow shelf approximation (SSA) • we will use the by-hand result to: ﬂow-line ice shelf numerical solution ◦ understand the SSA, and scaling up ◦ verify a numerical code next steps practicalities omitted models 5 e.g. [van der Veen, 1985] Numerical modelling exact thickness and velocity for steady ice shelf 2 Ed Bueler ice ﬂow: a 1 using ﬂotation criterion h = (1 − ρ/ρw )H, and because of no bed superﬁcial view resistance, equation (7) says “derivative of something = 0”: shallow ice approximation (SIA) 1 heat analogy & 2A−1/n H|ux |1/n−1 ux − ρg(1 − ρ/ρw )H 2 =0 (8) numerics x 2 x analogy ﬁnite differences 2 use the calving-front condition to integrate (8): exact solutions 1 2A−1/n H|ux |1/n−1 ux − ρg(1 − ρ/ρw )H 2 = 0 shallow ice sheets (9) solving the SIA 2 is it right? application 3 the steady (Ht = 0) mass continuity equation can be integrated; the most basic shallow assumption here M = M0 =constant but some other M(x) are allowed: ﬂowline Stokes derivation of SIA uH = M0 (x − xg ) + ug Hg , (10) shelves and streams 4 multiply (9) by u/H, replace uH from (10), assume positive strain shallow shelf approximation (SSA) rate (ux > 0), compute nth power ﬂow-line ice shelf numerical solution 5 get u n ux = Cs (M0 (x − xg ) + ug Hg )n , where Cs is known scaling up 6 integrate the last result to ﬁnd u(x), then return to (10) to ﬁnd H(x): next steps Cs n+1 u(x)n+1 = ug + (M0 (x − xg ) + ug Hg )n+1 − (ug Hg )n+1 , practicalities omitted models M0 M0 (x − xg ) + ug Hg H(x) = u(x) Numerical modelling exact thickness and velocity for steady ice shelf 3 Ed Bueler ice ﬂow: a superﬁcial view • Matlab/Octave code exactshelf.m uses shelf length L = 200 km, shallow ice approximation (SIA) and Hg = 500, ug = 50 m a−1 heat analogy & numerics • computes this shelf geometry and velocity: analogy ﬁnite differences exact solutions shallow ice sheets solving the SIA is it right? application the most basic shallow assumption ﬂowline Stokes derivation of SIA shelves and streams exact velocity; u = 303.8539 at x=L shallow shelf 350 approximation (SSA) 300 velocity (m a-1) ﬂow-line ice shelf numerical solution 250 scaling up 200 next steps 150 practicalities 100 omitted models 50 0 50 100 150 200 x (km) Numerical modelling numerically solving the SSA Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & • the stress balance is a two-point boundary value problem numerics analogy (BVP) for the velocity u(x): ﬁnite differences 2A−1/n H|ux |1/n−1 ux − C|u|m−1 u − ρgHhx = 0, exact solutions shallow ice x sheets solving the SIA left b.c.: u(xg ) = ug , is it right? application right b.c.: 2A−1/n H|ux |1/n−1 ux − 1 ρ(1 − ρ/ρw )gH 2 = 0 at x = xc 2 the most basic shallow assumption ﬂowline Stokes • here we will prescribe the ice geometry, thus both the derivation of SIA thickness H(x) and the driving stress ρgHhx , and we will ﬁnd shelves and streams the velocity u(x) shallow shelf approximation (SSA) • I’ll describe the numerical method for a ice shelf or ice stream ﬂow-line ice shelf numerical solution scaling up • . . . but then I’ll give a code only for an ice shelf, so C = 0 and next steps h = (1 − ρ/ρw )H practicalities omitted models Numerical modelling numerically solving the SSA stress balance 2 Ed Bueler ice ﬂow: a superﬁcial view shallow ice • like most such nonlinear equations, iteration is needed to approximation (SIA) solve this one heat analogy & numerics • red term ν = A−1/n |ux |1/n−1 is the “effective viscosity”: ¯ analogy ﬁnite differences exact solutions shallow ice 2A−1/n |ux |1/n−1 Hux − C|u|m−1 u − ρgHhx = 0 sheets x solving the SIA is it right? • let application the most basic W (ux ) = 2¯H = 2A−1/n |ux |1/n−1 H ν shallow assumption ﬂowline Stokes derivation of SIA • idea: use old effective viscosity to get new velocity solution, shelves and and repeat until the solution is not changing too much and/or streams shallow shelf the differential equation is nearly-satisﬁed approximation (SSA) ﬂow-line ice shelf • idea as algorithm: from u (k −1) , ﬁnd next values u (k ) by numerical solution scaling up solving the linear problem next steps practicalities (k −1) (k ) omitted models W (ux )ux − C|u (k −1) |m−1 u (k ) − ρgHhx = 0 x Numerical modelling numerically solving the SSA stress balance 3 Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) where do you get an initial guess u (0) (x) for the velocity? heat analogy & numerics analogy • for ﬂoating ice, a velocity comes from assuming a uniform ﬁnite differences exact solutions strain rate: shallow ice sheets u (0) = γ(x − xg ) + ug solving the SIA is it right? application the most basic ◦ for example: γ is the value of ux found from the calving front shallow assumption ﬂowline Stokes stress imbalance and ug is the known velocity at the grounding derivation of SIA line shelves and streams • for grounded ice, a velocity comes from assuming the ice is shallow shelf approximation (SSA) held by basal resistance only: ﬂow-line ice shelf numerical solution 1/m scaling up u (0) = −C −1 ρgHhx next steps practicalities omitted models Numerical modelling solving the “inner” linear problem Ed Bueler ice ﬂow: a superﬁcial view shallow ice • reset x interval to be 0 < x < L, instead of xg < x < xc approximation (SIA) heat analogy & • abstract the problem to a linear BVP: numerics analogy ﬁnite differences exact solutions (W (x) ux )x − α(x) u = β(x) shallow ice sheets with boundary conditions solving the SIA is it right? application u(0) = V , ux (L) = γ the most basic shallow assumption ﬂowline Stokes derivation of SIA • W (x), α(x), β(x) are all known functions in this abstract shelves and problem streams shallow shelf approximation (SSA) • for the nonlinear SSA equation, both W (x) and α(x) will ﬂow-line ice shelf come from the previous iteration numerical solution scaling up • W (x) is needed on the staggered grid, for O(∆x 2 ) accuracy, next steps practicalities as with time-dependent diffusion problem earlier omitted models • α(x), β(x) are needed on regular grid Numerical modelling solving the “inner” linear problem 2 Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) • indices j = 1, 2, . . . , J + 1 heat analogy & numerics analogy • equal-spaced grid x1 , x2 , . . . , xJ+1 , where x1 = 0 and xJ+1 = L ﬁnite differences exact solutions • an approximation of the abstracted problem is: shallow ice sheets Wj+1/2 (uj+1 − uj ) − Wj−1/2 (uj − uj−1 ) ∗ solving the SIA − αj uj = βj is it right? ∆x 2 application the most basic shallow assumption • u1 = V given ﬂowline Stokes derivation of SIA • for right-hand boundary condition “ux (L) = γ”: shelves and streams ◦ introduce notional point xJ+2 shallow shelf ◦ approximation (SSA) uJ+2 − uJ ﬂow-line ice shelf =γ numerical solution 2∆x scaling up next steps ◦ using equation ∗ in j = J + 1 case, eliminate uJ+2 variable practicalities “by-hand” [Morton and Mayers, 2005] omitted models Numerical modelling solving the “inner” linear problem 3 Ed Bueler ice ﬂow: a superﬁcial view • so abstract linear system has matrix-vector form “ Ax = b ”: shallow ice approximation (SIA) V 1 u1 heat analogy & 2 numerics W3/2 A22 W5/2 u2 β2 ∆x u3 β3 ∆x 2 analogy W5/2 A33 ﬁnite differences . = . exact solutions .. .. . . . . . . shallow ice βJ ∆x 2 WJ−1/2 AJJ WJ+1/2 u J sheets solving the SIA AJ+1,J AJ+1,J+1 uJ+1 bJ+1 is it right? application • with diagonal entries (j = 2, 3, . . . , J) the most basic shallow assumption ﬂowline Stokes 2 derivation of SIA Ajj = −(Wj−1/2 + Wj+1/2 + αj ∆x ) shelves and streams • with special cases in last equation: shallow shelf approximation (SSA) ﬂow-line ice shelf AJ+1,J = WJ+1/2 + WJ+3/2 numerical solution scaling up 2 AJ+1,J+1 = −(WJ+1/2 + WJ+3/2 + αJ+1 ∆x ) next steps practicalities 2 omitted models bJ+1 = −2γ∆xWJ+3/2 + βJ+1 ∆x • this is a tridiagonal system Numerical modelling solving the “inner” linear problem 4 Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) function u = flowline(L,J,gamma,W,alpha,beta,V0) heat analogy & numerics dx = L / J; analogy rhs = dx^2 * beta(:); ﬁnite differences rhs(1) = V0; exact solutions rhs(J+1) = rhs(J+1) - 2 * gamma * dx * W(J+1); shallow ice sheets A = sparse(J+1,J+1); solving the SIA is it right? A(1,1) = 1.0; application for j=2:J the most basic A(j,j-1:j+1) = [ W(j-1), -(W(j-1) + W(j) + alpha(j) * dx^2), W(j) ]; shallow assumption end ﬂowline Stokes A(J+1,J) = W(J) + W(J+1); derivation of SIA A(J+1,J+1) = - (W(J) + W(J+1) + alpha(J+1) * dx^2); shelves and streams shallow shelf scale = full(max(abs(A),[],2)); approximation (SSA) for j=1:J+1, A(j,:) = A(j,:) ./ scale(j); end ﬂow-line ice shelf rhs = rhs ./ scale; numerical solution scaling up u = A \ rhs; next steps practicalities http://www.dms.uaf.edu/~bueler/karthaus/mfiles/flowline.m omitted models Numerical modelling solving the “inner” linear problem 5 Ed Bueler ice ﬂow: a • before proceeding to solve nonlinear SSA problem, we can superﬁcial view shallow ice test the “abstracted” code flowline.m approximation (SIA) • . . . by “manufacturing” solutions as in code heat analogy & numerics testflowline.m (not shown) analogy ﬁnite differences • results: works pretty well! converges at optimal rate exact solutions O(∆x 2 ) shallow ice sheets -3 solving the SIA 10 is it right? application the most basic shallow assumption 10-4 ﬂowline Stokes maximum error derivation of SIA shelves and 10-5 streams shallow shelf approximation (SSA) ﬂow-line ice shelf 10-6 numerical solution scaling up next steps 10-7 practicalities convergence rate O(dx2.00065) omitted models 10-8 10-4 10-3 10-2 10-1 dx Numerical modelling numerical solution to SSA Ed Bueler ice ﬂow: a superﬁcial view function [u,u0] = ssaflowline(p,J,H,b,ug,initchoice) shallow ice approximation (SIA) if nargin ~= 6, error(’exactly 6 input arguments required’), end heat analogy & dx = p.L / J; numerics x = (0:dx:p.L)’; analogy xstag = (dx/2:dx:p.L+dx/2)’; ﬁnite differences exact solutions alpha = p.C * ones(size(x)); h = H + b; shallow ice hx = regslope(dx,h); sheets beta = p.rho * p.g * H .* hx; solving the SIA gamma = ( 0.25 * p.A^(1/p.n) * (1 - p.rho/p.rhow) *... is it right? p.rho * p.g * H(end) )^p.n; application u0 = ssainit(p,x,beta,gamma,initchoice); the most basic u = u0; shallow assumption ﬂowline Stokes Hstag = stagav(H); derivation of SIA tol = 1.0e-14; maxdiff = Inf; shelves and W = zeros(J+1,1); streams while maxdiff > tol shallow shelf uxstag = stagslope(dx,u); approximation (SSA) W(1:J) = 2 * p.A^(-1/p.n) * Hstag .* (abs(uxstag)).^((1/p.n)-1); ﬂow-line ice shelf W(J+1) = W(J); numerical solution scaling up unew = flowline(p.L,J,gamma,W,alpha,beta,ug); maxdiff = max(abs(unew-u)); next steps u = unew; practicalities end omitted models http://www.dms.uaf.edu/~bueler/karthaus/mfiles/ssaflowline.m Numerical modelling numerical thickness and velocity for steady ice shelf Ed Bueler ice ﬂow: a • ﬁrst: let’s numerically solve the problem for which we know superﬁcial view shallow ice approximation (SIA) exact answer heat analogy & • here use very coarse 10 km grid numerics analogy -1 ﬁnite differences case with dx = 10 km has (max error) = 3.044 m a exact solutions 350 initial guess shallow ice exact solution sheets numerical solution 300 solving the SIA is it right? application 250 the most basic shallow assumption ﬂowline Stokes velocity (m a ) -1 derivation of SIA 200 shelves and streams 150 shallow shelf approximation (SSA) ﬂow-line ice shelf 100 numerical solution scaling up next steps 50 practicalities omitted models 0 0 50000 100000 150000 200000 x (km) Numerical modelling did we make any mistakes? Ed Bueler ice ﬂow: a • this does convergence analysis of ssaflowline.m: superﬁcial view J = [50 100 200 400 800 1600 3200]; dxkm = 200.0 ./ J; shallow ice approximation (SIA) for j=1:length(J) [av,maxerr(j)] = testshelf(J(j)); % calls ssaflowline.m heat analogy & numerics end analogy loglog(dxkm,maxerr,’0-’,’markersize’,16) ﬁnite differences exact solutions • goal is to see velocity error shrink by rate built-into our ﬁnite shallow ice differences, namely at rate O(∆x 2 ) sheets solving the SIA • by this standard, we made no mistakes! is it right? application 0 10 the most basic shallow assumption ﬂowline Stokes derivation of SIA 10-1 shelves and maximum error (m a-1) streams shallow shelf approximation (SSA) 10-2 ﬂow-line ice shelf numerical solution scaling up next steps 10-3 practicalities convergence rate O(dx-2.00846) omitted models 10-4 10-1 100 101 dx (km) Numerical modelling realistic ice shelf modeling Ed Bueler ice ﬂow: a • ﬂow lines (1D) are never very realistic . . . superﬁcial view shallow ice • “diagnostic” (= ﬁxed geometry) ice shelf modeling in two horizontal approximation (SIA) variables (2D) has been quite successful heat analogy & numerics • solve an elliptic PDE boundary value problem analogy ﬁnite differences • velocity measurements are available for validation exact solutions ◦ Ross ice shelf example based on [RIGGS; Bentley, 1974] data below: model (red shallow ice arrows) vs observations (black arrows) sheets ◦ this is PISM but many models can do it [MacAyeal and others, 1996] solving the SIA is it right? ◦ something of a mystery why regular grid methods do well! application the most basic shallow assumption ﬂowline Stokes derivation of SIA shelves and streams shallow shelf approximation (SSA) ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models Numerical modelling numerical solution of SSA: a summary Ed Bueler ice ﬂow: a superﬁcial view shallow ice • SSA and Blatter [1995] stress balance equations each approximation (SIA) ◦ determine horizontal velocity from geometry and b.c.s heat analogy & numerics ◦ are nonlinear problems: require iterative linearization analogy ◦ each “inner” linear problem has a sparse pattern ﬁnite differences exact solutions ∗ tridiagonal in SSA ﬂow-line (1D) case here shallow ice ∗ much less structured in SSA 2D, Blatter 2D, and Blatter 3D cases sheets solving the SIA • the Stokes stress balance is much harder because is it right? application incompressibility must be handled explicitly the most basic shallow assumption • a wise numerical modeler would give the sparse linear ﬂowline Stokes derivation of SIA problem to a numerical linear algebra software package: shelves and ◦ M ATLAB/O CTAVE, LAPACK, PETSc, Trilinos, . . . streams shallow shelf ◦ generally: modularize your code and test the parts separately! approximation (SSA) ﬂow-line ice shelf • comment on mass continuity: numerical solution scaling up ¯ ◦ the mass continuity PDE “Ht = M − · (uH)” is not very next steps ¯ diffusive if the source of u is SSA or Blatter, practicalities omitted models ◦ and there is not much theory to guide how to solve such a generic transport problem Numerical modelling scaling-up the numerical solution of stress balances Ed Bueler ice ﬂow: a superﬁcial view shallow ice • all stress balance equations used for ice ﬂow modeling are approximation (SIA) nonlinear equations for velocity u in the form heat analogy & numerics analogy ﬁnite differences F(u) = 0 exact solutions shallow ice sheets solving the SIA • the “residual” F involves various partial derivatives, nonlinear is it right? operations, and additional ﬁelds (e.g. thickness) application the most basic shallow assumption • SSA and Blatter stress balance equations can be written in ﬂowline Stokes derivation of SIA form shelves and A(u) u = g streams shallow shelf approximation (SSA) ◦ where A is a nonlinear function of u, ﬂow-line ice shelf numerical solution ◦ and A(u) acts linearly on u scaling up • we have already treated the SSA in the second way next steps practicalities • both forms are nonlinear, elliptic PDEs omitted models • Newton [1669] knew a way of solving nonlinear equations! Numerical modelling scaling-up the numerical solution of stress balances 2 Ed Bueler ice ﬂow: a • the iteration scheme we have used for SSA is “Picard”: superﬁcial view shallow ice approximation (SIA) A(u(k −1) ) u(k ) = g heat analogy & numerics analogy ◦ this iteration converges linearly, if it converges: u(k +1) − u(k ) ≤ γ u(k ) − u(k −1) , where 0 < γ < 1 ﬁnite differences exact solutions shallow ice ◦ it can be robust, but it is not fast sheets solving the SIA • Newton’s method (e.g. [Press and others, 1992]) also solves is it right? application F(u) = 0 iteratively, but by linearization: the most basic shallow assumption ﬂowline Stokes J(u(k −1) ) w = −F(u(k −1) ), derivation of SIA shelves and streams u(k ) = u(k −1) + w shallow shelf approximation (SSA) ﬂow-line ice shelf where numerical solution ∂F(u)i scaling up J(u)ij = is the Jacobian of F next steps ∂uj practicalities omitted models ◦ this iteration converges quadratically, if it converges: u(k +1) − u(k ) ≤ C u(k ) − u(k −1) 2 ◦ much faster! Numerical modelling scaling-up the numerical solution of stress balances 3 Ed Bueler ice ﬂow: a • big computers still take a long time to solve SSA equations superﬁcial view shallow ice when grids are ∼ 2 km for Greenland approximation (SIA) ◦ Antarctica is 10× the area, and that much worse! heat analogy & numerics analogy • there are various numerical paradigms for discretizing (ﬁnite ﬁnite differences exact solutions difference, ﬁnite element, spectral, . . . ), but choosing among shallow ice these is not the major scalability issue sheets solving the SIA • effective scalability of 3D stress balances (e.g. Blatter & is it right? application Stokes) at high spatial resolution, on large parallel machines the most basic requires: shallow assumption ﬂowline Stokes ◦ the “inner” linear solve must be iterative–not by Gaussian derivation of SIA elimination—and must be preconditioned shelves and streams ◦ the “outer” nonlinear iteration must offer faster-than-linear shallow shelf approximation (SSA) convergence ﬂow-line ice shelf ◦ because the exact Jacobian is too hard to evaluate—certainly numerical solution scaling up so with coupling to other physics (e.g. thermo-)—the method next steps must allow approximated Jacobians practicalities omitted models • this is the promise of “preconditioned matrix-free Newton-Krylov” (e.g. [Knoll and Keyes, 2004]) as a paradigm for big, nasty nonlinear problems Numerical modelling scaling-up the numerical solution of stress balances 4 Ed Bueler ice ﬂow: a • the PETSc [Balay and others, Level of superﬁcial view Abstraction Application Codes shallow ice 2010] library has evolved into a approximation (SIA) Newton-Krylov toolkit → TS heat analogy & (Time Stepping) numerics • my online materials turn SNES (Nonlinear Equations Solvers) analogy ﬁnite differences ssaflowline.m into a C KSP PC (Krylov Subspace Methods) exact solutions PETSc example (Preconditioners) shallow ice ssaflowline.c: sheets Matrices Vectors Index Sets solving the SIA ◦ parallel over an arbitrary is it right? BLAS MPI application number of processors the most basic shallow assumption ◦ uses Newton method (or 104 M=10 ﬂowline Stokes Picard) M=100 derivation of SIA 102 M=103 shelves and ◦ try matrix-free, different M=104 M=105 100 streams preconditioners, . . . shallow shelf approximation (SSA) ◦ result from 10-2 residual norm ﬂow-line ice shelf ssaflowline.c is at 10-4 numerical solution scaling up lower-right: 10-6 next steps ∗ evidence of super-linear practicalities convergence from Newton 10-8 omitted models method ∗ on coarsest grids, 10 or 10-10 100 points, get rapid 10-12 convergence to 10−13 of 0 1 2 3 4 5 6 Newton iteration the initial residual norm Numerical modelling Outline Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & numerics 1 ice ﬂow: an outsider’s superﬁcial view analogy ﬁnite differences exact solutions shallow ice 2 heat analogy & numerics sheets solving the SIA is it right? application 3 shallow ice sheets the most basic shallow assumption ﬂowline Stokes derivation of SIA 4 shelves and streams shelves and streams shallow shelf approximation (SSA) 5 next steps ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models Numerical modelling what technical skills are needed for Ed Bueler numerical ice sheet modeling? ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & you need: numerics analogy • comfort in a technical computing environment, usually Unix, ﬁnite differences exact solutions in which you need to know: shallow ice sheets ◦ an editor, solving the SIA ◦ a compiled language (Fortran or C), is it right? application ◦ a scripting/prototyping language (Matlab, Python, etc.), and the most basic shallow assumption ◦ a version control system (Subversion, git, etc.) ﬂowline Stokes derivation of SIA • willingness to read math, numerical analysis, computer shelves and science books streams shallow shelf approximation (SSA) • . . . and willingness to ignor some of the advice found there ﬂow-line ice shelf numerical solution • know some tools for NetCDF ﬁles scaling up next steps • get exposed to parallel computing practicalities omitted models • but, at the end of the day: physics Numerical modelling what technical skills are needed? 2 Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) • an important skill is to not re-invent the wheel heat analogy & numerics analogy • never re-invent the wheel for basic numerics: ﬁnite differences exact solutions ◦ M ATLAB, Comsol, PETSc, libmesh, Elmer, Trilinos, Triangle, shallow ice etc. handle fundamental numerical linear algebra, mesh sheets generation, ﬁnite element assembly and solve, etc. tasks; don’t solving the SIA is it right? even try to compete! application the most basic ◦ . . . except to write throw-away codes to help you understand shallow assumption numerical ideas ﬂowline Stokes derivation of SIA • sometimes there is no need to re-invent ice sheet modeling: shelves and streams ◦ open source SIA-based comprehensive models: GLIMMER, shallow shelf approximation (SSA) SICOPOLIS ﬂow-line ice shelf ◦ open source hybrid/higher-order/Stokes models: PISM, numerical solution scaling up Elmer-ice, CISM next steps practicalities • glacier and ice sheet modeling is young! there is much to do! omitted models Numerical modelling omitted models 1: thermomechanical coupling Ed Bueler 1e-23 ice ﬂow: a superﬁcial view • the softness of ice “A” is not shallow ice approximation (SIA) constant! 1e-24 heat analogy & ◦ A = A(T ) varies by more than A(T) 1e-25 numerics analogy 103 in the 50 ◦ C ice temperature Paterson-Budd form for A(T) ﬁnite differences range relevant to Antarctic ice 1e-26 exact solutions ◦ dissipation of ﬂow energy is shallow ice sheets signiﬁcant in conservation of 1e-27 -50 -40 -30 -20 T (degrees C) -10 0 solving the SIA is it right? energy equation application ◦ . . . and there is a “new” ﬂuid the most basic shallow assumption instability from feedback between ﬂowline Stokes derivation of SIA the above two items (lower ﬁgure; shelves and EISMINT II [Payne and others, 2000]) streams shallow shelf • ice temperature is part of the “long approximation (SSA) ﬂow-line ice shelf memory” of the ice sheets for past numerical solution scaling up climate next steps • conservation of energy is just as practicalities omitted models important for fast scale dynamics (involving sliding and hydrology) as it is for “long memory” questions Numerical modelling omitted models 2: surface mass balance Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & numerics analogy ﬁnite differences • ice sheet models need modeled surface mass balance exact solutions shallow ice ◦ over paleoglacial time scales sheets solving the SIA ◦ and when modeling response to future changes is it right? application • models include: the most basic shallow assumption ◦ positive degree-day (PDD) and other index models ﬂowline Stokes ◦ energy-balance models derivation of SIA shelves and • very signiﬁcant source of uncertainty for Greenland ice sheet streams shallow shelf dynamics! approximation (SSA) ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models Numerical modelling omitted models 3: solid Earth deformation Ed Bueler ice ﬂow: a 1 superﬁcial view • ice has almost 3 of the density of the viscous hot rock in the shallow ice approximation (SIA) Earth’s mantle, heat analogy & numerics • so 1000 m of ice will depress the Earth’s crust almost 300 m analogy ﬁnite differences if allowed enough time exact solutions • this changes bed topography and thus ice ﬂow, so Earth shallow ice sheets deformation is important to modeling ice ﬂow solving the SIA is it right? • read?: application the most basic ◦ [Peltier, 1998] = survey of observations and processes shallow assumption ﬂowline Stokes ◦ [Greve, 2001] = comparison of practical models derivation of SIA shelves and streams shallow shelf approximation (SSA) ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models Numerical modelling omitted models 4: numerical Stokes equations Ed Bueler ice ﬂow: a superﬁcial view • the Stokes model itselfcan be solved numerically shallow ice approximation (SIA) ◦ no shallow assumptions! heat analogy & ◦ but many more equations and unknowns than SIA, SSA, numerics analogy hybrids, Blatter, . . . ﬁnite differences exact solutions • requires explicit accounting for shallow ice sheets ◦ incompressibility = a constraint on the ﬂow solving the SIA ◦ pressure = a Lagrange multiplier for that constraint is it right? application • very tough scalability issues: can you afford the loss of the most basic shallow assumption resolution? loss of long-time runs? ﬂowline Stokes derivation of SIA • all of the success so far is at a smaller scale (e.g. below) shelves and streams shallow shelf approximation (SSA) ﬂow-line ice shelf numerical solution Figure 7 in [Maxwell and others, 2008]: scaling up Athabasca Glacier: (a) modeled velocity next steps contour lines (m a−1 ); (b) contour lines practicalities derived from measurements [Raymond, omitted models 1971] Numerical modelling omitted models 5: “higher-order” schemes Ed Bueler ice ﬂow: a superﬁcial view • both the SIA and the SSA are derived by small-parameter shallow ice approximation (SIA) arguments from the Stokes equations, so . . . heat analogy & numerics • is there a common shallow antecedent model? analogy ﬁnite differences • Schoof and Hindmarsh [2009] answer: exact solutions shallow ice ◦ yes, the Blatter [1995] model is intermediate between the sheets Stokes stress balance and both the SIA and SSA solving the SIA is it right? application Stokes the most basic shallow assumption ﬂowline Stokes derivation of SIA other higher-order? shelves and streams shallow shelf approximation (SSA) ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models Numerical modelling POSTSCRIPT on free boundary value problems Ed Bueler ice ﬂow: a superﬁcial view • ice sheet/shelf modeling means free boundary problems shallow ice approximation (SIA) • Hutter [1999] identiﬁes some below heat analogy & numerics analogy ﬁnite differences exact solutions shallow ice sheets solving the SIA is it right? application the most basic shallow assumption ﬂowline Stokes derivation of SIA shelves and streams shallow shelf approximation (SSA) ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models Numerical modelling what is a “free boundary”? Ed Bueler • a free boundary for a PDE is an unknown location at which ice ﬂow: a superﬁcial view there is a boundary condition shallow ice approximation (SIA) ◦ the location of the free boundary must be found at the same heat analogy & time as one solves the PDE problem numerics analogy ◦ in addition to the information which would already be present at ﬁnite differences a ﬁxed-location boundary, there must be enough additional exact solutions information at a free boundary to determine its location shallow ice sheets ◦ all free boundary problems are nonlinear, regardless of the solving the SIA is it right? linearity of the PDE application • classic example: consider an elastic membrane attached to a the most basic shallow assumption ﬂowline Stokes wire frame and stretched over an obstacle: derivation of SIA shelves and streams constraint: shallow shelf approximation (SSA) u≥ψ ﬂow-line ice shelf numerical solution in locations where u > ψ, solve: scaling up 2 next steps u=0 practicalities omitted models where is the free boundary, and what facts about u are true there? Numerical modelling free boundary value problem 1: polythermal ice Ed Bueler ice ﬂow: a • by volume, majority of ice sheet is cold (T < 0 ◦ C) superﬁcial view shallow ice • . . . but there is some ice which is temperate, where T = 0 ◦ C and approximation (SIA) heat analogy & there is a positive liquid fraction within the ice matrix numerics analogy • . . . so ice sheets are polythermal ﬁnite differences exact solutions • boundary between cold and temperate ice is “CTS” (= shallow ice cold-temperate transition surface): sheets solving the SIA ◦ must be found, as free boundary, when solving conservation of is it right? application energy equation (“Stefan probem”) the most basic shallow assumption ◦ can be tracked explicitly [Greve, 1997] ﬂowline Stokes ◦ or treated as a level surface of the enthalpy variable derivation of SIA [Aschwanden and Blatter, 2009] shelves and streams shallow shelf approximation (SSA) ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models Canadian-type Scandinavian-type temperate cold Numerical modelling free boundary value problem 2: for ice streaming Ed Bueler ice ﬂow: a superﬁcial view • Schoof’s [2006] insight, for diagnostic case shallow ice approximation (SIA) well-posed free boundary problem heat analogy & SSA + (plastic till) = numerics analogy for location and velocity of sliding ﬁnite differences exact solutions • “plastic till” means the basal strength (resistance) is given by shallow ice sheets a yield stress τc : τb = τc vb /|vb | solving the SIA is it right? • Schoof’s scheme is a whole ice sheet form of MacAyeal’s application the most basic [1989] individual ice stream models shallow assumption ﬂowline Stokes derivation of SIA shelves and streams shallow shelf approximation (SSA) ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models Numerical modelling free boundary problem 3: for grounded ice sheet margin Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) • side-by-side comparison, classical elastic membrane heat analogy & problem versus steady ice sheet problem numerics analogy ﬁnite differences constraint: constraint: exact solutions shallow ice u≥ψ h≥b ⇐⇒ H≥0 sheets solving the SIA is it right? where u > ψ, solve: where h > b, solve steady SIA: application the most basic 2 shallow assumption u=0 0=M+ · ΓH n+2 | h|n−1 h ﬂowline Stokes derivation of SIA shelves and streams shallow shelf approximation (SSA) ﬂow-line ice shelf numerical solution scaling up next steps practicalities omitted models Numerical modelling for SIA model, an illuminating transformation Ed Bueler ice ﬂow: a superﬁcial view • there are degenerate diffusions in other ﬁelds shallow ice approximation (SIA) • (e.g. porous media where there can be wet and dry parts) heat analogy & numerics • . . . whence transformations useful to ice sheet modeling analogy ﬁnite differences exact solutions • speciﬁcally one by Raviart [1970] (n = 3 case): η = H 8/3 shallow ice • the ﬂat bed SIA simpliﬁes: sheets solving the SIA is it right? application Ht = M + · ΓH 5 | H|2 H the most basic shallow assumption becomes ˜ ﬂowline Stokes derivation of SIA → η 3/8 =M+ · (Γ| η|2 η) t shelves and streams shallow shelf approximation (SSA) ◦ new diffusivity (red) depends on | η| only, not η ﬂow-line ice shelf ◦ the margin shape is not singular in new variable! numerical solution scaling up ∗ η is globally continuous next steps ∗ z = η(t, x, y ) is tangent to bed at margin practicalities omitted models ◦ numerical methods applied to the new equation make errors which are as expected for “good” free boundary problems [Bueler and others, 2005] Numerical modelling free boundary problems: so what? Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & numerics analogy • why does it matter that many glaciological problems are free ﬁnite differences exact solutions boundary problems in their mathematical form? shallow ice • the location of the free boundary may be the glaciological sheets solving the SIA question is it right? application • free boundaries are always locations of loss of smoothness the most basic shallow assumption relative to ﬁxed boundary solutions ﬂowline Stokes derivation of SIA ◦ numerical error may be dominated by errors near these free shelves and boundaries, streams shallow shelf ◦ hard to know whether model results at free boundaries are approximation (SSA) poor because of numerical problems or because of missing ﬂow-line ice shelf numerical solution physical processes scaling up next steps practicalities omitted models Numerical modelling The End Ed Bueler ice ﬂow: a superﬁcial view shallow ice approximation (SIA) heat analogy & numerics analogy ﬁnite differences exact solutions shallow ice sheets solving the SIA is it right? application the most basic shallow assumption ﬂowline Stokes derivation of SIA shelves and streams shallow shelf approximation (SSA) ﬂow-line ice shelf numerical solution scaling up Thanks for bearing with me! next steps practicalities omitted models