tbombproof.rst - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tbombproof.rst (27714B)
---
1 .. include:: ../global.txt
2
3 .. only:: html
4
5 .. When building PDFs this will be included already.
6
7 .. include:: ../math-definitions.txt
8
9 .. _sec-bombproof:
10
11 BOMBPROOF, PISM's numerical scheme for conservation of energy
12 =============================================================
13
14
15 Introduction
16 ------------
17
18 One of the essential goals for any thermomechanically-coupled numerical ice sheet model is
19 a completely bombproof numerical scheme for the advection-conduction-reaction problem for
20 the conservation of energy within the ice. "Bombproof" means being as stable as possible
21 in as many realistic modeling contexts as possible. PISM's scheme is observed to be highly
22 robust in practice, but it is also provably stable in a significant range of
23 circumstances. The scheme is special to shallow ice sheets because it makes specific
24 tradeoff choices with respect to vertical velocity. It is generic and low order in how it
25 treats horizontal velocity. In this page we state the scheme, prove its stability
26 properties, and address the basal boundary condition.
27
28 The scheme is conditionally-stable. The length of the time step is limited only by the
29 maximum magnitude of the horizontal velocities within the ice, i.e. horizontal CFL. This
30 condition for stability is included in the PISM adaptive time-stepping technique.
31
32 Accuracy is necessarily a second goal. Our shallow scheme has truncation error `O(\Delta
33 z^2)` in many circumstances, though it reverts to lower order when it "detects trouble" in
34 the form of large vertical velocities. Overall the scheme has order `O(\Delta t,\Delta
35 x,\Delta y, \Delta z^2)` in circumstances where the vertical ice flow velocity is small
36 enough, relative to conductivity, and otherwise reverts to
37 `O(\Delta t,\Delta x,\Delta y, \Delta z^1)`.
38
39 The conservation of energy problem for the ice is in terms of an enthalpy field
40 :cite:`AschwandenBuelerKhroulevBlatter`. The current scheme supercedes the cold-ice,
41 temperature-based scheme described in the appendices of :cite:`BBL` and in
42 :cite:`BBssasliding`. Compared to a cold-ice scheme, the enthalpy formulation does a
43 better job of conserving energy, has a more-physical model for basal melt rate and
44 drainage, and can model polythermal ice with any CTS topology
45 :cite:`AschwandenBuelerKhroulevBlatter`. The finite difference implementation of the
46 enthalpy method is robust and avoids the CFL condition on vertical advection which was
47 present in the older, cold-ice scheme.
48
49 The bedrock thermal problem is solved by splitting the timestep into an update of the
50 bedrock temperature field, assuming the ice base is as constant temperature, and an update
51 of the ice enthalpy field, by the BOMBPROOF scheme here, assuming the upward heat flux
52 from the bedrock layer is constant during the timestep. For more on the implementation,
53 see the ``BedThermalUnit`` class.
54
55 The region in which the conservation of energy equation needs to be solved changes over
56 time. This is an essential complicating factor in ice sheet modeling. Also relevant is
57 that the velocity field has a complicated provenance as it comes from different stress
58 balance equations chosen at runtime. These stress balances, especially with transitions in
59 flow type, for instance at grounding lines, are incompletely understood when
60 thermomechanically-coupled. (The ``ShallowStressBalance`` instance owned by
61 ``IceModel`` could be the SIA, the SSA, a hybrid of these, or other stress balances in
62 future PISM versions.) We will therefore not make, in proving stability, assumptions about
63 the regularity of the velocity field in space or time other than boundedness.
64
65 Nor do we want the numerical scheme for advection to need any information about the
66 velocity except its value at the beginning of the time step. Thus the conservation of
67 energy timestep is assumed to be split from the mass continuity time step and its
68 associated stress balance solve. We have not considered implementing a scheme which
69 requires the Jacobian of the velocity field with respect to changes in enthalpy, for
70 example. At very least such a fully-implicit scheme would require blind iteration (e.g.
71 with no guarantee of convergence of the iteration). The scheme we propose involves no such
72 iteration.
73
74
75 Conservation of energy in a shallow ice sheet
76 ---------------------------------------------
77
78 In an enthalpy formulation :cite:`AschwandenBuelerKhroulevBlatter` (and references therein),
79 the ice sheet is regarded as a mixture of two phases, solid and liquid, so that both cold
80 and temperate ice with liquid ice matrix can be modeled. The specific enthalpy field of
81 the ice mixture is denoted `E(x,y,z,t)` and has units `J / kg`. (Within
82 the PISM documentation the symbol `H` is used for ice thickness so we use `E` for enthalpy
83 here and in the PISM source code versus "H" in :cite:`AschwandenBuelerKhroulevBlatter`.)
84 The conservation of energy equation is
85
86 .. math::
87 :label: basicEnergy
88
89 \rho \frac{dE}{dt} = -\nabla \cdot \mathbf{q} + Q,
90
91 where `\rho` is the mixture density. The mixture density is assumed to be the same as ice
92 density even if there is a nonzero liquid fraction, and the mixture is assumed to be
93 incompressible :cite:`AschwandenBuelerKhroulevBlatter`. The left and right sides of
94 equation :eq:`basicEnergy`, and thus the quantity `Q`, have units
95 `J\,\text{s}^{-1}\,\text{m}^{-3} = \text{W}\,\text{m}^{-3}`.
96
97 Neglecting the dependence of conductivity and heat capacity on temperature
98 :cite:`AschwandenBuelerKhroulevBlatter`, the heat flux in cold ice and temperate ice is
99
100 .. math::
101 :label: heatflux
102
103 \mathbf{q} =
104 \begin{cases}
105 - \frac{k_i}{c_i} \nabla E, & \text{cold ice}, \\
106 - K_0 \nabla E, & \text{temperate ice},\\
107 \end{cases}
108
109 where `k_i,c_i,K_0` are constant :cite:`AschwandenBuelerKhroulevBlatter`. The nonzero flux
110 in the temperate ice case, may be conceptualized as a regularization of the "real"
111 equation, or as a flux of latent heat carried by liquid water. Also, `dE/dt` stands for
112 the material derivative of the enthalpy field, so the expanded form of :eq:`basicEnergy`
113 is
114
115 .. math::
116
117 \rho \left(\frac{\partial E}{\partial t} + \mathbf{U}\cdot \nabla E\right)
118 = \nabla \cdot \left(\left\{\begin{matrix} \frac{k_i}{c_i} \\ K_0 \end{matrix}\right\}
119 \nabla E \right) + Q,
120
121 where `\mathbf{U}` is the three dimensional velocity, thus advection is included.
122
123 The additive quantity `Q` is the dissipation (strain-rate) heating,
124
125 .. math::
126
127 Q = \sum_{i,j=1}^3 D_{ij} \tau_{ij}
128
129 where `D_{ij}` is the strain rate tensor and `\tau_{ij}` is the deviatoric stress tensor.
130 Reference :cite:`BBssasliding` addresses how this term is computed in PISM, according to
131 the shallow stress balance approximations; see
132 ``StressBalance::compute_volumetric_strain_heating()``. (`Q` is called `\Sigma` in
133 :cite:`BBL`, :cite:`BBssasliding` and in many places in the source code.)
134
135 Friction from sliding also is a source of heating. It has units of `W / m^2 = J / (m^2 s)`, that
136 is, the same units as the heat flux `\mathbf{q}` above. In formulas we write
137
138 .. math::
139
140 F_b = - \tau_b \cdot \mathbf{u}_b,
141
142 where `\tau_b` is the basal shear stress and `\mathbf{u}_b` is the basal sliding velocity;
143 the basal shear stress is oppositely-directed to the basal velocity. For example, in the
144 plastic case `\tau_b = - \tau_c \mathbf{u}_b / |\mathbf{u}_b|` where `\tau_c` is a
145 positive scalar, the yield stress. See method
146 ``StressBalance::basal_frictional_heating()``. The friction heating is concentrated at
147 `z=0`, and it enters into the basal boundary condition and melt rate calculation,
148 addressed in section :ref:`sec-melt` below.
149
150 We use a shallow approximation of equation :eq:`basicEnergy` which lacks horizontal
151 conduction terms :cite:`Fowler`. For the initial analysis of the core BOMBPROOF
152 scheme, we specialize to cold ice. Within cold ice, the coefficient in the heat
153 flux is constant, so
154
155 .. math::
156
157 \nabla \cdot \mathbf{q} = - \frac{k_i}{c_i} \frac{\partial^2 E}{\partial z^2}.
158
159 Therefore the equation we initially analyze is
160
161 .. math::
162 :label: basicShallow
163
164 \rho_i \left(\frac{\partial E}{\partial t} + \mathbf{U}\cdot \nabla E\right) = \frac{k_i}{c_i} \frac{\partial^2 E}{\partial z^2} + Q,
165
166 We focus the analysis on the direction in which the enthalpy has largest derivative,
167 namely with respect to the vertical coordinate `z`. Rewriting equation
168 :eq:`basicShallow` to emphasize the vertical terms we have
169
170 .. math::
171 :label: vertProblem
172
173 \rho_i \left(\frac{\partial E}{\partial t} + w \frac{\partial E}{\partial z}\right)
174 = \frac{k_i}{c_i} \frac{\partial^2 E}{\partial z^2} + \Phi
175
176 where
177
178 .. math::
179
180 \Phi = Q - \rho_i \left(u \frac{\partial E}{\partial x}
181 + v \frac{\partial E}{\partial y}\right)
182
183 We assume that the surface enthalpy `E_s(t,x,y)` (K) and the geothermal flux `G(t,x,y)`
184 (`W / m^2`) at `z=0` are given. (The latter is the output of the ``energy::BedThermalUnit``
185 object, and it may come from an evolving temperature field within the upper crust, the
186 bedrock layer. If a surface temperature is given then it will be converted to enthalpy by
187 the ``EnthalpyConverter`` class.) The boundary conditions to problem :eq:`vertProblem`
188 are, therefore,
189
190 .. math::
191 :label: columnbcs
192
193 E(t,x,y,z=H) &=E_s(t,x,y),\\
194 -\frac{k_i}{c_i} \frac{\partial E}{\partial z}\Big|_{z=0} &= G.
195
196 For a temperate ice base, including any ice base below which there is liquid water,
197 the lower boundary condition is more interesting. It is addressed below in section
198 :ref:`sec-melt`.
199
200 The core BOMBPROOF scheme
201 -------------------------
202
203 For the discussion of the numerical scheme below, let `E_{ijk}^n` be our approximation to
204 the exact enthalpy `E` at the grid point with coordinates `(x_i,y_j,z_k)` at time `t_n`.
205 When `i,j` are uninteresting we suppress them and write `E_k^n`, and we will use similar
206 notation for numerical approximations to the other quantities. We put the horizontal
207 advection terms in the source term `\Phi` because we treat them explicitly, evaluating at
208 time `t_n`. (Implicit or semi-implicit treatment of horizontal advection would require a
209 coupled system distributed across processors, a difficulty which is currently avoided.)
210
211 The scheme we use for horizontal advection is explicit first-order upwinding. There is a
212 CFL condition for the scheme to be stable, in the absence of conduction, based on the
213 magnitude of the horizontal velocity components. To state the upwind scheme itself, let
214
215 .. math::
216
217 \Up{f_{\bullet}}{\alpha} =
218 \begin{cases}
219 f_i-f_{i-1}, & \alpha \ge 0, \\
220 f_{i+1}-f_i, & \alpha < 0.
221 \end{cases}
222
223 The approximate horizontal advection terms, and thus the approximation to the whole
224 term `\Phi`, are
225
226 .. math::
227
228 \Phi_{ijk}^n = \Sigma_{ijk}^n - \rho_i
229 \left( u_{ijk}^n\,\frac{\Up{E_{\bullet jk}^n}{u_{ijk}^n}}{\Delta x}
230 + v_{ijk}^n\,\frac{\Up{E_{i\bullet k}^n}{v_{ijk}^n}}{\Delta y} \right).
231
232 The CFL stability condition for this part of the scheme is
233
234
235 .. math::
236 :label: CFL
237
238 \Delta t \,\left( \left|\frac{u_{ijk}^n}{\Delta x}\right|
239 + \left|\frac{v_{ijk}^n}{\Delta y}\right| \right) \le 1.
240
241 The routine ``max_timestep_cfl_3d()`` computes the maximum of velocity
242 magnitudes. This produces a time step restriction based on the above CFL condition. Then
243 ``IceModel::max_timestep()`` implements adaptive time-stepping based on this and
244 other stability criteria.
245
246 In the analysis below we assume an equally-spaced grid `z_0,\dots,z_{M_z}`
247 with `\Delta z = z_{k+1} - z_k`. In fact PISM has a remapping scheme in each
248 column, wherein the enthalpy in a column of ice is stored on an unequally-spaced
249 vertical grid, but is mapped to a fine, equally-spaced grid for the conservation
250 of energy computation described here. (Similar structure applies to the age
251 computation. See classes ``EnthalpyModel`` and ``AgeModel``.)
252
253 The `z` derivative terms in :eq:`vertProblem` will be approximated implicitly. Let
254 `\lambda` be in the interval `0 \le \lambda \le 1`. Suppressing indices `i,j`, the
255 approximation to :eq:`vertProblem` is
256
257
258 .. math::
259 :label: bombone
260
261 \rho_i &\left(
262 \frac{E_k^{n+1} - E_k^n}{\Delta t}
263 + \lambda w_k^n \frac{E_{k+1}^{n+1} - E_{k-1}^{n+1}}{2 \Delta z}
264 + (1-\lambda) w_k^{n} \frac{\Up{E_{\bullet}^{n+1}}{w_k^{n}}}{\Delta z} \right) \\
265 &= \frac{k_i}{c_i}\, \frac{E_{k+1}^{n+1} - 2 E_{k}^{n+1} + E_{k-1}^{n+1}}{\Delta z^2} + \Phi_k^n.
266
267 Equation :eq:`bombone`, along with a determination of `\lambda` by
268 :eq:`lambdachoice` below, is the scheme BOMBPROOF. It includes two approximations
269 of vertical advection, implicit centered difference (`\lambda = 1`) and
270 implicit first-order upwinding (`\lambda=0`). They are combined using
271 nonnegative coefficients which sum to one, a convex combination. The centered
272 formula has higher accuracy,
273
274 .. math::
275
276 w_k^n \frac{E_{k+1}^{n+1} - E_{k-1}^{n+1}}{2 \Delta z}
277 = w \frac{\partial E}{\partial z} + O(\Delta t,\Delta z^2),
278
279 while the first order upwind formula has lower accuracy,
280
281 .. math::
282
283 w_k^{n} \frac{\Up{E_{\bullet}^{n+1}}{w_k^{n}}}{\Delta z}
284 = w \frac{\partial E}{\partial z} + O(\Delta t,\Delta z).
285
286 Thus we prefer to use the centered formula when possible, but we apply (implicit)
287 upwinding when it is needed for its added stability benefits.
288
289 We now rewrite :eq:`bombone` for computational purposes as one of a system of equations
290 for the unknowns `\{E_k^{n+1}\}`. In this system the coefficients will be
291 scaled so that the diagonal entries of the matrix have limit one as
292 `\Delta t\to 0`. Let
293
294 .. math::
295
296 \nu &= \frac{\Delta t}{\Delta z},\\
297 R &= \frac{k_i \Delta t}{\rho_i c_i \Delta z^2}.
298
299 Now multiply equation :eq:`bombone` by `\Delta t`, divide it by `\rho_i`,
300 and rearrange:
301
302 .. math::
303 :label: bombtwo
304
305 \left(-R - \nu w_k^n \uppair{1-\lambda/2}{\lambda/2}\right) E_{k-1}^{n+1}\\ +
306 \left(1 + 2 R + \nu w_k^n (1-\lambda) \uppair{+1}{-1}\right) E_k^{n+1}\\ +
307 \left(-R + \nu w_k^n \uppair{\lambda/2}{1-\lambda/2} \right) E_{k+1}^{n+1}\\
308 = E_k^n + \Delta t \rho_i^{-1}\Phi_k^n
309
310 Here `\uppair{a}{b} = a` when `w_k^n\ge 0` and `\uppair{a}{b} = b` when `w_k^n < 0`.
311
312 Equation :eq:`bombtwo` has coefficients which are scaled to have no units. It is
313 ready to be put in the system managed by ``enthSystemCtx``.
314
315 One way of stating the stability of first-order upwinding is to say it satisfies
316 a "maximum principle" :cite:`MortonMayers`. An example of a maximum principle
317 for this kind of finite difference scheme is that if `U_{k-1}^n,U_k^n,U_{k+1}^n`
318 are adjacent gridded values of some abstract quantity at time step `t_n`, and
319 if the next value satisfies the scheme
320
321 .. math::
322 :label: abstractexplicit
323
324 U_k^{n+1} = C_{-1} U_{k-1}^n + C_0 U_k^n + C_{+1} U_{k+1}^n
325
326 for *nonnegative* coefficients `C_i` summing to one, `C_{-1} + C_0 + C_{+1} = 1`, then it
327 follows by the triangle inequality that
328
329 .. math::
330
331 \min\{|U_{k-1}^n|, |U_k^n|, |U_{k+1}^n|\}
332 \le |U_k^{n+1}| \le \max\{|U_{k-1}^n|, |U_k^n|, |U_{k+1}^n|\}.
333
334 Thus a "wiggle" cannot appear in `\{U_k^{n+1}\}` if previous values `\{U_k^n\}` were
335 smoother. The proof below shows the corresponding "wiggle-free" property for scheme :eq:`bombtwo`.
336
337 However, the pure implicit centered difference scheme (`\lambda=1`), namely
338
339 .. math::
340 :label: centered
341
342 &\left(-R - \nu w_k^n/2\right) E_{k-1}^{n+1} + \left(1 + 2 R\right) E_k^{n+1} \\
343 &+ \left(-R + \nu w_k^n/2\right) E_{k+1}^{n+1}
344 = E_k^n + \Delta t \rho_i^{-1}\Phi_k^n
345
346 is *less stable* than implicit first-order upwinding. It is less stable in the same sense
347 that Crank-Nicolson is a less stable scheme than backwards Euler for the simplest heat
348 equation `u_t = u_{xx}` :cite:`MortonMayers`. In fact, although oscillatory modes cannot
349 grow exponentially under equation :eq:`centered`, those modes *can* appear when none are
350 present already, even in the homogeneous case `\Phi_k^n=0`.
351
352 Stability properties of the BOMBPROOF scheme
353 --------------------------------------------
354
355 We want to be precise about the phrase "unconditionally stable" for BOMBPROOF.
356 To do so we consider somewhat simplified cases which are amenable to analysis, and
357 we prove two stability properties. These stability properties identify the
358 precise advantages of BOMBPROOF.
359
360 Theorem (stating the stability properties).
361 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
362
363 Assume, for the precise but limited assertion of this theorem, that the surface
364 temperature `T_s` and the geothermal flux `G` are constant in time. Assume also that the
365 entire source function `\Phi` is identically zero (but see comments below). Fix an
366 equally-spaced vertical grid `z_0=0 < z_1 < \dots < z_N=H`, so that the upper grid point
367 coincides with the surface of the ice. With these assumptions, if
368
369 .. math::
370 :label: lambdachoice
371
372 \lambda = \min
373 \left\{
374 1, \min_{k=0,\dots,N} \left\{ \frac{2 k_i}{|w_k^n| \rho_i c_i \Delta z} \right\}
375 \right\},
376
377 reset at each time step `n`, then scheme :eq:`bombone`, :eq:`bombtwo` is
378 unconditionally-stable in the following two senses:
379
380 1. A maximum principle applies without further assumptions.
381
382 2. Suppose we freeze the coefficients of the problem to have constant values in time and
383 space. (Concretely, we assume that `\lambda` is chosen independently of the time step
384 `n`, and that `\Delta t` is the same for each time step. We assume constant vertical
385 velocity `w_k^n=w_0`. We also consider a spatially-periodic or unbounded version of our
386 problem, with no boundary conditions.) Then a von Neumann analysis of the constant
387 coefficient problem yields a growth factor less than one for all modes on the grid.
388
389 Remarks
390 ^^^^^^^
391
392 The phrases *maximum principle* and *von Neumann analysis* will be precisely illustrated
393 in the following proof. Both approaches are in :cite:`MortonMayers`. There is additional
394 information on the von Neumann analysis of implicit finite difference methods for
395 advection in :cite:`Strikwerda`.
396
397 These statements also apply in case `k_i=0`, in which case :eq:`lambdachoice` implies
398 `\lambda=0`, and the method reduces to implicit first-order upwinding. (Implicit
399 first-order upwinding has properties 1 and 2 :cite:`Strikwerda`.) The case `k_i=0` is
400 relevant because it applies to the least-transport model of temperate ice in which there
401 is zero enthalpy conduction. (One reasonable model for temperate ice is to assume no
402 transport of the liquid fraction, whether diffusive transport or otherwise, and to ignore
403 conduction along the temperature gradient, because the gradient is only from
404 pressure-melting temperature differences.)
405
406 Proof of 1
407 ^^^^^^^^^^
408
409 In the case considered for the maximum principle, with `\Phi_k^n=0`,
410 we can rewrite :eq:`bombtwo` as
411
412 .. math::
413 :label: formax
414
415 &\left(1 + 2 R + \nu w_k^n (1-\lambda) \uppair{+1}{-1}\right) E_k^{n+1} \\
416 &= E_k^n + \left(R + \nu w_k^n \uppair{1-\lambda/2}{\lambda/2}\right) E_{k-1}^{n+1}
417 + \left(R - \nu w_k^n \uppair{\lambda/2}{1-\lambda/2}\right) E_{k+1}^{n+1}.
418
419 We claim that with choice :eq:`lambdachoice` for `0 \le \lambda \le 1`, all
420 coefficients in :eq:`formax` are nonnegative. At one extreme, in
421 the upwinding case (`\lambda=0`), all the coefficients are nonnegative. Otherwise, note that
422 `\nu w_k^n (1-\lambda) \uppair{+1}{-1}` is nonnegative for any valid value
423 of `\lambda` and for any value of `w_k^n`, noting the meaning of the `\uppair{+1}{-1}`
424 symbol. Thus the coefficient on the left is always nonnegative. The coefficient of
425 `E_{k-1}^{n+1}` is clearly nonnegative for any valid value of `\lambda` if `w_k^n \ge 0`.
426 The coefficient of `E_{k+1}^{n+1}` is clearly nonnegative for any valid value of `\lambda` if
427 `w_k^n \le 0`.
428
429 Therefore the only concerns are for the coefficient of `E_{k-1}^{n+1}` when `w_k^n\le 0` and the
430 coefficient of `E_{k+1}^{n+1}` when `w_k^n\ge 0`. But if `\lambda` is smaller than
431 `2k_i/(|w_k^n| \rho_i c_i \Delta z)` then
432
433 .. math::
434
435 R - \nu |w_k^n| (\lambda/2) = \frac{k_i \Delta t}{\rho_i c_i \Delta z^2} - \frac{\Delta t |w_k^n|}{\Delta z} \frac{\lambda}{2} \ge \frac{k_i \Delta t}{\rho_i c_i \Delta z^2}
436 - \frac{\Delta t |w_k^n|}{\Delta z} \frac{k_i}{|w_k^n| \rho_i c_i \Delta z} = 0.
437
438 Thus all the coefficients in :eq:`formax` are nonnegative. On the other hand, in equation :eq:`formax`, all coefficients on the right side sum to
439
440 .. math::
441
442 1+2R+\nu w_k^n \uppair{1-\lambda}{-1+\lambda} = 1+2R+\nu w_k^n (1-\lambda) \uppair{+1}{-1},
443
444 which is exactly the coefficient on the left side of :eq:`formax`. It follows that
445
446 .. math::
447
448 E_k^{n+1} = a_k E_k^n + b_k E_{k-1}^{n+1} + c_k E_{k+1}^{n+1}
449
450 where `a_k,b_k,c_k` are positive and `a_k+b_k+c_k=1`. Thus a maximum principle applies
451 :cite:`MortonMayers`. **END OF PROOF OF 1.**
452
453 Proof of 2
454 ^^^^^^^^^^
455
456 As a von Neumann analysis is much more restrictive than the analysis above, we will be brief.
457 Let's assume the velocity is downward, `w_0<0`; the other case is similar. Equation
458 :eq:`bombtwo` becomes
459
460 .. math::
461 :label: prevon
462
463 &\left(-R - \nu w_0 (\lambda/2)\right) E_{k-1}^{n+1}
464 + \left(1 + 2 R - \nu w_0 (1-\lambda)\right) E_k^{n+1} \\
465 & + \left(-R + \nu w_0 (1-\lambda/2) \right) E_{k+1}^{n+1} = E_k^n.
466
467 The heart of the von Neumann analysis is the substitution of a growing or decaying
468 (in time index `n`) oscillatory mode on the grid of spatial wave number `\mu`:
469
470 .. math::
471
472 E_k^n = \sigma^n e^{i\mu\,(k\Delta z)}.
473
474 Here `k\Delta z = z_k` is a grid point. Such a mode is a solution to :eq:`prevon` if and
475 only if
476
477 .. math::
478
479 \sigma\Big[ &(-R - \nu w_0(\lambda/2)) e^{-i\mu\Delta z}
480 + (1 + 2 R - \nu w_0 (1-\lambda)) \\
481 &+ (-R + \nu w_0 (\lambda/2)) e^{+i\mu\Delta z}
482 + \nu w_0 (1-\lambda) e^{+i\mu\Delta z} \Big] = 1.
483
484 This equation reduces by standard manipulations to
485
486 .. math::
487
488 \sigma = \frac{1}{1 + \left(4 R - 2 \nu w_0 (1-\lambda)\right)\cos^2(\mu \Delta z/2)
489 + i\,\nu w_0 (1-\lambda/2)\sin(\mu\Delta z)}.
490
491 Note `4 R - 2 \nu w_0 (1-\lambda) \ge 0` without restrictions on
492 numerical parameters `\Delta t`, `\Delta z`, because `w_0<0` in the
493 case under consideration. Therefore
494
495 .. math::
496
497 |\sigma|^2 = \frac{1}{\left[1 + \left(4 R - 2 \nu w_0 (1-\lambda)\right)
498 \cos^2(\mu \Delta z/2)\right]^2
499 + \left[\nu w_0 (1-\lambda/2)\sin(\mu\Delta z)\right]^2}.
500
501 This positive number is less than one, so `|\sigma| < 1`. It follows that all
502 modes decay exponentially.
503 **END OF PROOF OF 2.**
504
505 Remark about our von Neumann stability analysis
506 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
507
508 The constant `\lambda` is carefully chosen in :eq:`lambdachoice` so that the maximum
509 principle 1 applies. On the other hand, both the implicit first-order upwind and the
510 implicit centered difference formulas have unconditional stability in the von Neumann
511 sense. The proof of case 2 above is thus a formality, merely showing that a convex
512 combination of unconditionally stable (von Neumann sense) schemes is still unconditionally
513 stable in the same sense.
514
515 Convergence: a consequence of the maximum principle
516 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
517
518 If we define the pointwise numerical error `e_k^n = E_k^n - E(t_n,x_i,y_j,z_k)`,
519 where `E(\dots)` is the unknown exact solution (exact enthalpy field) :cite:`MortonMayers`,
520 then :eq:`formax` implies an equality of the form
521
522 .. math::
523
524 A e_k^{n+1} = e_k^n + B_- e_{k-1}^{n+1} + B_+ e_{k+1}^{n+1} + \Delta t\, \tau_k^n
525
526 where `\tau_k^n` is the truncation error of the scheme and `A,B_\pm` are nonnegative
527 coefficients, which need no detail for now other than to note that `1 + B_- + B_+ = A`.
528 Letting `{\bar e}^n = \max_k |e_k^n|` we have, because of the positivity of coefficients,
529
530 .. math::
531 :label: prebound
532
533 A |e_k^{n+1}| \le {\bar e}^n + \left(B_- + B_+\right){\bar e}^{n+1} + \Delta t\,\bar\tau^n
534
535 for all `k`, where `\bar\tau^n = \max_k |\tau_k^n|`. Now let `k` be the index for
536 which `|e_k^{n+1}| = {\bar e}^{n+1}`. For that `k` we can replace `|e_k^{n+1}|` in
537 equation :eq:`prebound` with `{\bar e}^{n+1}`. Subtracting the same quantity from
538 each side of the resulting inequality gives
539
540 .. math::
541
542 {\bar e}^{n+1} \le {\bar e}^n + \Delta t\,\bar\tau^n,
543
544 It follows that `\bar e^n \le C \Delta t`, for some finite `C`, if `\bar e^0 = 0`
545 :cite:`MortonMayers`. Thus a maximum principle for BOMBPROOF implies convergence
546 in the standard way :cite:`MortonMayers`. This convergence proof has the same
547 assumptions as case 1 in the theorem, and thus it only *suggests* convergence
548 in any broad range of glaciologically-interesting cases.
549
550
551 Remark on nonzero source term
552 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
553
554 Now recall we assumed in Theorem 1 that the entire "source" `\Phi_k^n` was identically zero.
555 Of course this is not realistic. What we understand is provable, however, is that if a
556 numerical scheme for a linear advection/conduction equation
557
558 .. math::
559
560 u_t + A u_x = B u_{xx}
561
562 is stable in the general sense of numerical schemes for partial differential equations
563 (e.g. as defined in subsection 5.5 of :cite:`MortonMayers`) then the same scheme is stable
564 in the same general sense when applied to the equation with (linear) lower order terms:
565
566 .. math::
567
568 u_t + A u_x = B u_{xx} + C u + D.
569
570 A precise statement of this general fact is hard to find in the literature, to put it
571 mildly, but theorem 2.2.3 of :cite:`Strikwerda` is one interesting case (`B=0` and `D=0`).
572 But even the form we state with linear term (`C u + D`) is not adequate to the job because
573 of the strongly-nonlinear dependence of `\Phi` on the temperature `T` :cite:`BBL`.
574
575 Nonetheless the maximum principle is a highly-desirable form of stability because we can
576 exclude "wiggles" from the finite difference approximations of the conductive and
577 advective terms, even if the complete physics, with strain heating in particular, is not
578 yet shown to be non-explosive. Because the complete physics includes the appearance of the
579 famous "spokes" of EISMINT II, for example, a maximum principle cannot apply too
580 literally. Indeed there is an underlying fluid instability :cite:`BBL`, one that means the
581 solution of the continuum equations can include growing "wiggles" which are fluid features
582 (though not at the grid-based spatial frequency of the usual numerical wiggles). Recall
583 that, because we use first-order upwinding on the horizontal advection terms, we can
584 expect maximum principle-type stability behavior of the whole three-dimensional scheme.
585
586 .. _sec-melt:
587
588 Temperate basal boundary condition, and computing the basal melt rate
589 ---------------------------------------------------------------------
590
591 At the bottom of grounded ice, a certain amount of heat comes out of the earth and either
592 enters the ice through conduction or melts the base of the ice. On the one hand, see the
593 documentation for ``BedThermalUnit`` for the model of how much comes out of the
594 earth. On the other hand, :cite:`AschwandenBuelerKhroulevBlatter` includes a careful
595 analysis of the subglacial layer equation and the corresponding boundary conditions and
596 basal melt rate calculation, and the reader should consult that reference.
597
598 Regarding the floating case
599 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
600
601 The shelf base temperature `T_{sb}` and the melt rate `M` are supplied by ``OceanModel``.
602 Note that we make the possibly-peculiar physical choice that the shelf base temperature is
603 used as the temperature at the *top of the bedrock*, which is actually the bottom of the
604 ocean. This choice means that there should be no abrupt changes in top-of-bedrock heat
605 flux as the grounding line moves. This choice also means that the conservation of energy
606 code does not need to know about the bedrock topography or the elevation of sea level. (In
607 the future ``OceanModel`` could have a ``subshelf_bed_temperature()`` method.)