BOMBPROOF, PISM’s numerical scheme for conservation of energy

Introduction

One of the essential goals for any thermomechanically-coupled numerical ice sheet model is a completely bombproof numerical scheme for the advection-conduction-reaction problem for the conservation of energy within the ice. “Bombproof” means being as stable as possible in as many realistic modeling contexts as possible. PISM’s scheme is observed to be highly robust in practice, but it is also provably stable in a significant range of circumstances. The scheme is special to shallow ice sheets because it makes specific tradeoff choices with respect to vertical velocity. It is generic and low order in how it treats horizontal velocity. In this page we state the scheme, prove its stability properties, and address the basal boundary condition.

The scheme is conditionally-stable. The length of the time step is limited only by the maximum magnitude of the horizontal velocities within the ice, i.e. horizontal CFL. This condition for stability is included in the PISM adaptive time-stepping technique.

Accuracy is necessarily a second goal. Our shallow scheme has truncation error $O(\mathrm{\Delta }{z}^2)$ in many circumstances, though it reverts to lower order when it “detects trouble” in the form of large vertical velocities. Overall the scheme has order $O(\mathrm{\Delta }t,\mathrm{\Delta }x,\mathrm{\Delta }y,\mathrm{\Delta }{z}^2)$ in circumstances where the vertical ice flow velocity is small enough, relative to conductivity, and otherwise reverts to $O(\mathrm{\Delta }t,\mathrm{\Delta }x,\mathrm{\Delta }y,\mathrm{\Delta }{z}^1)$.

The conservation of energy problem for the ice is in terms of an enthalpy field [22]. The current scheme supercedes the cold-ice, temperature-based scheme described in the appendices of [16] and in [10]. Compared to a cold-ice scheme, the enthalpy formulation does a better job of conserving energy, has a more-physical model for basal melt rate and drainage, and can model polythermal ice with any CTS topology [22]. The finite difference implementation of the enthalpy method is robust and avoids the CFL condition on vertical advection which was present in the older, cold-ice scheme.

The bedrock thermal problem is solved by splitting the timestep into an update of the bedrock temperature field, assuming the ice base is as constant temperature, and an update of the ice enthalpy field, by the BOMBPROOF scheme here, assuming the upward heat flux from the bedrock layer is constant during the timestep. For more on the implementation, see the BedThermalUnit class.

The region in which the conservation of energy equation needs to be solved changes over time. This is an essential complicating factor in ice sheet modeling. Also relevant is that the velocity field has a complicated provenance as it comes from different stress balance equations chosen at runtime. These stress balances, especially with transitions in flow type, for instance at grounding lines, are incompletely understood when thermomechanically-coupled. (The ShallowStressBalance instance owned by IceModel could be the SIA, the SSA, a hybrid of these, or other stress balances in future PISM versions.) We will therefore not make, in proving stability, assumptions about the regularity of the velocity field in space or time other than boundedness.

Nor do we want the numerical scheme for advection to need any information about the velocity except its value at the beginning of the time step. Thus the conservation of energy timestep is assumed to be split from the mass continuity time step and its associated stress balance solve. We have not considered implementing a scheme which requires the Jacobian of the velocity field with respect to changes in enthalpy, for example. At very least such a fully-implicit scheme would require blind iteration (e.g. with no guarantee of convergence of the iteration). The scheme we propose involves no such iteration.

Conservation of energy in a shallow ice sheet

In an enthalpy formulation [22] (and references therein), the ice sheet is regarded as a mixture of two phases, solid and liquid, so that both cold and temperate ice with liquid ice matrix can be modeled. The specific enthalpy field of the ice mixture is denoted $E(x,y,z,t)$ and has units $J/kg$. (Within the PISM documentation the symbol $H$ is used for ice thickness so we use $E$ for enthalpy here and in the PISM source code versus “H” in [22].) The conservation of energy equation is

(82)$$\rho \frac{dE}{dt}=-\mathrm{\nabla }\cdot \mathbf{q}+Q,$$

where $\rho$ is the mixture density. The mixture density is assumed to be the same as ice density even if there is a nonzero liquid fraction, and the mixture is assumed to be incompressible [22]. The left and right sides of equation (82), and thus the quantity $Q$, have units $J{\text{s}}^{-1}{\text{m}}^{-3}=\text{W}{\text{m}}^{-3}$.

Neglecting the dependence of conductivity and heat capacity on temperature [22], the heat flux in cold ice and temperate ice is

(83)$$\begin{array}{r}\mathbf{q}=\{\begin{array}{ll}-\frac{k_i}{c_i}\mathrm{\nabla }E,& \text{cold ice},\\ -K_0\mathrm{\nabla }E,& \text{temperate ice},\end{array}\end{array}$$

where $k_i,c_i,K_0$ are constant [22]. The nonzero flux in the temperate ice case, may be conceptualized as a regularization of the “real” equation, or as a flux of latent heat carried by liquid water. Also, $dE/dt$ stands for the material derivative of the enthalpy field, so the expanded form of (82) is

$$\begin{array}{r}\rho (\frac{\partial E}{\partial t}+\mathbf{U}\cdot \mathrm{\nabla }E)=\mathrm{\nabla }\cdot \left(\left\{\begin{array}{c}\frac{k_i}{c_i}\\ K_0\end{array}\right\}\mathrm{\nabla }E\right)+Q,\end{array}$$

where $\mathbf{U}$ is the three dimensional velocity, thus advection is included.

The additive quantity $Q$ is the dissipation (strain-rate) heating,

$$Q=\sum _{i,j=1}^3{D}_{ij}{\tau }_{ij}$$

where $D_{ij}$ is the strain rate tensor and ${\tau }_{ij}$ is the deviatoric stress tensor. Reference [10] addresses how this term is computed in PISM, according to the shallow stress balance approximations; see StressBalance::compute_volumetric_strain_heating(). ($Q$ is called $\mathrm{\Sigma }$ in [16], [10] and in many places in the source code.)

Friction from sliding also is a source of heating. It has units of $W/m^2=J/(m^{2}s)$, that is, the same units as the heat flux $\mathbf{q}$ above. In formulas we write

$$F_b=-{\tau }_b\cdot {\mathbf{u}}_b,$$

where ${\tau }_b$ is the basal shear stress and ${\mathbf{u}}_b$ is the basal sliding velocity; the basal shear stress is oppositely-directed to the basal velocity. For example, in the plastic case ${\tau }_b=-{\tau }_c{\mathbf{u}}_b/|{\mathbf{u}}_b|$ where ${\tau }_c$ is a positive scalar, the yield stress. See method StressBalance::basal_frictional_heating(). The friction heating is concentrated at $z=0$, and it enters into the basal boundary condition and melt rate calculation, addressed in section Temperate basal boundary condition, and computing the basal melt rate below.

We use a shallow approximation of equation (82) which lacks horizontal conduction terms [38]. For the initial analysis of the core BOMBPROOF scheme, we specialize to cold ice. Within cold ice, the coefficient in the heat flux is constant, so

$$\mathrm{\nabla }\cdot \mathbf{q}=-\frac{k_i}{c_i}\frac{{\partial }^{2}E}{\partial z^2}.$$

Therefore the equation we initially analyze is

(84)$${\rho }_i(\frac{\partial E}{\partial t}+\mathbf{U}\cdot \mathrm{\nabla }E)=\frac{k_i}{c_i}\frac{{\partial }^{2}E}{\partial z^2}+Q,$$

We focus the analysis on the direction in which the enthalpy has largest derivative, namely with respect to the vertical coordinate $z$. Rewriting equation (84) to emphasize the vertical terms we have

(85)$${\rho }_i(\frac{\partial E}{\partial t}+w\frac{\partial E}{\partial z})=\frac{k_i}{c_i}\frac{{\partial }^{2}E}{\partial z^2}+\mathrm{\Phi }$$

where

$$\mathrm{\Phi }=Q-{\rho }_i(u\frac{\partial E}{\partial x}+v\frac{\partial E}{\partial y})$$

We assume that the surface enthalpy $E_s(t,x,y)$ (K) and the geothermal flux $G(t,x,y)$ ($W/m^2$) at $z=0$ are given. (The latter is the output of the energy::BedThermalUnit object, and it may come from an evolving temperature field within the upper crust, the bedrock layer. If a surface temperature is given then it will be converted to enthalpy by the EnthalpyConverter class.) The boundary conditions to problem (85) are, therefore,

(86)$$\begin{array}{rl}E(t,x,y,z=H)& =E_s(t,x,y),\\ -\frac{k_i}{c_i}\frac{\partial E}{\partial z}{|}_{z=0}& =G.\end{array}$$

For a temperate ice base, including any ice base below which there is liquid water, the lower boundary condition is more interesting. It is addressed below in section Temperate basal boundary condition, and computing the basal melt rate.

The core BOMBPROOF scheme

For the discussion of the numerical scheme below, let $E_{ijk}^n$ be our approximation to the exact enthalpy $E$ at the grid point with coordinates $(x_i,y_j,z_k)$ at time $t_n$. When $i,j$ are uninteresting we suppress them and write $E_k^n$, and we will use similar notation for numerical approximations to the other quantities. We put the horizontal advection terms in the source term $\mathrm{\Phi }$ because we treat them explicitly, evaluating at time $t_n$. (Implicit or semi-implicit treatment of horizontal advection would require a coupled system distributed across processors, a difficulty which is currently avoided.)

The scheme we use for horizontal advection is explicit first-order upwinding. There is a CFL condition for the scheme to be stable, in the absence of conduction, based on the magnitude of the horizontal velocity components. To state the upwind scheme itself, let

$$\begin{array}{r}\mathrm{Up}\left(f_{\bullet }|\alpha \right)=\{\begin{array}{ll}{f}_i-f_{i-1},& \alpha \ge 0,\\ f_{i+1}-f_i,& \alpha <0.\end{array}\end{array}$$

The approximate horizontal advection terms, and thus the approximation to the whole term $\mathrm{\Phi }$, are

$${\mathrm{\Phi }}_{ijk}^n={\mathrm{\Sigma }}_{ijk}^n-{\rho }_i(u_{ijk}^n\frac{\mathrm{Up}\left(E_{\bullet jk}^n|u_{ijk}^n\right)}{\mathrm{\Delta }x}+v_{ijk}^n\frac{\mathrm{Up}\left(E_{i\bullet k}^n|v_{ijk}^n\right)}{\mathrm{\Delta }y}).$$

The CFL stability condition for this part of the scheme is

(87)$$\mathrm{\Delta }t(\left|\frac{u_{ijk}^n}{\mathrm{\Delta }x}\right|+\left|\frac{v_{ijk}^n}{\mathrm{\Delta }y}\right|)\le 1.$$

The routine max_timestep_cfl_3d() computes the maximum of velocity magnitudes. This produces a time step restriction based on the above CFL condition. Then IceModel::max_timestep() implements adaptive time-stepping based on this and other stability criteria.

In the analysis below we assume an equally-spaced grid $z_0,\dots ,z_{M_z}$ with $\mathrm{\Delta }z=z_{k+1}-z_k$. In fact PISM has a remapping scheme in each column, wherein the enthalpy in a column of ice is stored on an unequally-spaced vertical grid, but is mapped to a fine, equally-spaced grid for the conservation of energy computation described here. (Similar structure applies to the age computation. See classes EnthalpyModel and AgeModel.)

The $z$ derivative terms in (85) will be approximated implicitly. Let $\lambda$ be in the interval $0\le \lambda \le 1$. Suppressing indices $i,j$, the approximation to (85) is

(88)$$\begin{array}{rl}{\rho }_i& (\frac{E_k^{n+1}-E_k^n}{\mathrm{\Delta }t}+\lambda w_k^n\frac{E_{k+1}^{n+1}-E_{k-1}^{n+1}}{2\mathrm{\Delta }z}+(1-\lambda )w_k^n\frac{\mathrm{Up}\left(E_{\bullet }^{n+1}|w_k^n\right)}{\mathrm{\Delta }z})\\ & =\frac{k_i}{c_i}\frac{E_{k+1}^{n+1}-2E_k^{n+1}+E_{k-1}^{n+1}}{\mathrm{\Delta }{z}^2}+{\mathrm{\Phi }}_k^n.\end{array}$$

Equation (88), along with a determination of $\lambda$ by (92) below, is the scheme BOMBPROOF. It includes two approximations of vertical advection, implicit centered difference ($\lambda =1$) and implicit first-order upwinding ($\lambda =0$). They are combined using nonnegative coefficients which sum to one, a convex combination. The centered formula has higher accuracy,

$$w_k^n\frac{E_{k+1}^{n+1}-E_{k-1}^{n+1}}{2\mathrm{\Delta }z}=w\frac{\partial E}{\partial z}+O(\mathrm{\Delta }t,\mathrm{\Delta }{z}^2),$$

while the first order upwind formula has lower accuracy,

$$w_k^n\frac{\mathrm{Up}\left(E_{\bullet }^{n+1}|w_k^n\right)}{\mathrm{\Delta }z}=w\frac{\partial E}{\partial z}+O(\mathrm{\Delta }t,\mathrm{\Delta }z).$$

Thus we prefer to use the centered formula when possible, but we apply (implicit) upwinding when it is needed for its added stability benefits.

We now rewrite (88) for computational purposes as one of a system of equations for the unknowns $\{E_k^{n+1}\}$. In this system the coefficients will be scaled so that the diagonal entries of the matrix have limit one as $\mathrm{\Delta }t\to 0$. Let

$$\begin{array}{rl}\nu & =\frac{\mathrm{\Delta }t}{\mathrm{\Delta }z},\\ R& =\frac{k_i\mathrm{\Delta }t}{{\rho }_i{c}_i\mathrm{\Delta }{z}^2}.\end{array}$$

Now multiply equation (88) by $\mathrm{\Delta }t$, divide it by ${\rho }_i$, and rearrange:

(89)$$\begin{array}{r}(-R-\nu w_k^n\left\{\begin{array}{c}1-\lambda /2\\ \lambda /2\end{array}\right\})E_{k-1}^{n+1}\\ +(1+2R+\nu w_k^n(1-\lambda )\left\{\begin{array}{c}+1\\ -1\end{array}\right\})E_k^{n+1}\\ +(-R+\nu w_k^n\left\{\begin{array}{c}\lambda /2\\ 1-\lambda /2\end{array}\right\})E_{k+1}^{n+1}\\ =E_k^n+\mathrm{\Delta }t{\rho }_i^{-1}{\mathrm{\Phi }}_k^n\end{array}$$

Here $\left\{\begin{array}{c}a\\ b\end{array}\right\}=a$ when $w_k^n\ge 0$ and $\left\{\begin{array}{c}a\\ b\end{array}\right\}=b$ when $w_k^n<0$.

Equation (89) has coefficients which are scaled to have no units. It is ready to be put in the system managed by enthSystemCtx.

One way of stating the stability of first-order upwinding is to say it satisfies a “maximum principle” [74]. An example of a maximum principle for this kind of finite difference scheme is that if $U_{k-1}^n,U_k^n,U_{k+1}^n$ are adjacent gridded values of some abstract quantity at time step $t_n$, and if the next value satisfies the scheme

(90)$$U_k^{n+1}=C_{-1}{U}_{k-1}^n+C_0{U}_k^n+C_{+1}{U}_{k+1}^n$$

for nonnegative coefficients $C_i$ summing to one, $C_{-1}+C_0+C_{+1}=1$, then it follows by the triangle inequality that

$$min\{|U_{k-1}^n|,|U_k^n|,|U_{k+1}^n|\}\le |U_k^{n+1}|\le max\{|U_{k-1}^n|,|U_k^n|,|U_{k+1}^n|\}.$$

Thus a “wiggle” cannot appear in $\{U_k^{n+1}\}$ if previous values $\{U_k^n\}$ were smoother. The proof below shows the corresponding “wiggle-free” property for scheme (89).

However, the pure implicit centered difference scheme ($\lambda =1$), namely

(91)$$\begin{array}{rl}& (-R-\nu w_k^n/2)E_{k-1}^{n+1}+(1+2R)E_k^{n+1}\\ & +(-R+\nu w_k^n/2)E_{k+1}^{n+1}=E_k^n+\mathrm{\Delta }t{\rho }_i^{-1}{\mathrm{\Phi }}_k^n\end{array}$$

is less stable than implicit first-order upwinding. It is less stable in the same sense that Crank-Nicolson is a less stable scheme than backwards Euler for the simplest heat equation $u_t=u_{xx}$ [74]. In fact, although oscillatory modes cannot grow exponentially under equation (91), those modes can appear when none are present already, even in the homogeneous case ${\mathrm{\Phi }}_k^n=0$.

Stability properties of the BOMBPROOF scheme

We want to be precise about the phrase “unconditionally stable” for BOMBPROOF. To do so we consider somewhat simplified cases which are amenable to analysis, and we prove two stability properties. These stability properties identify the precise advantages of BOMBPROOF.

Theorem (stating the stability properties).

Assume, for the precise but limited assertion of this theorem, that the surface temperature $T_s$ and the geothermal flux $G$ are constant in time. Assume also that the entire source function $\mathrm{\Phi }$ is identically zero (but see comments below). Fix an equally-spaced vertical grid $z_0=0<z_1<\cdots <z_N=H$, so that the upper grid point coincides with the surface of the ice. With these assumptions, if

(92)$$\lambda =min\{1,\underset{k=0,\dots ,N}{min}\left\{\frac{2k_i}{|w_k^n|{\rho }_i{c}_i\mathrm{\Delta }z}\right\}\},$$

reset at each time step $n$, then scheme (88), (89) is unconditionally-stable in the following two senses:

1. A maximum principle applies without further assumptions.

2. Suppose we freeze the coefficients of the problem to have constant values in time and space. (Concretely, we assume that $\lambda$ is chosen independently of the time step $n$, and that $\mathrm{\Delta }t$ is the same for each time step. We assume constant vertical velocity $w_k^n=w_0$. We also consider a spatially-periodic or unbounded version of our problem, with no boundary conditions.) Then a von Neumann analysis of the constant coefficient problem yields a growth factor less than one for all modes on the grid.

Remarks

The phrases maximum principle and von Neumann analysis will be precisely illustrated in the following proof. Both approaches are in [74]. There is additional information on the von Neumann analysis of implicit finite difference methods for advection in [172].

These statements also apply in case $k_i=0$, in which case (92) implies $\lambda =0$, and the method reduces to implicit first-order upwinding. (Implicit first-order upwinding has properties 1 and 2 [172].) The case $k_i=0$ is relevant because it applies to the least-transport model of temperate ice in which there is zero enthalpy conduction. (One reasonable model for temperate ice is to assume no transport of the liquid fraction, whether diffusive transport or otherwise, and to ignore conduction along the temperature gradient, because the gradient is only from pressure-melting temperature differences.)

Proof of 1

In the case considered for the maximum principle, with ${\mathrm{\Phi }}_k^n=0$, we can rewrite (89) as

(93)$$\begin{array}{rl}& (1+2R+\nu w_k^n(1-\lambda )\left\{\begin{array}{c}+1\\ -1\end{array}\right\})E_k^{n+1}\\ & =E_k^n+(R+\nu w_k^n\left\{\begin{array}{c}1-\lambda /2\\ \lambda /2\end{array}\right\})E_{k-1}^{n+1}+(R-\nu w_k^n\left\{\begin{array}{c}\lambda /2\\ 1-\lambda /2\end{array}\right\})E_{k+1}^{n+1}.\end{array}$$

We claim that with choice (92) for $0\le \lambda \le 1$, all coefficients in (93) are nonnegative. At one extreme, in the upwinding case ($\lambda =0$), all the coefficients are nonnegative. Otherwise, note that $\nu w_k^n(1-\lambda )\left\{\begin{array}{c}+1\\ -1\end{array}\right\}$ is nonnegative for any valid value of $\lambda$ and for any value of $w_k^n$, noting the meaning of the $\left\{\begin{array}{c}+1\\ -1\end{array}\right\}$ symbol. Thus the coefficient on the left is always nonnegative. The coefficient of $E_{k-1}^{n+1}$ is clearly nonnegative for any valid value of $\lambda$ if $w_k^n\ge 0$. The coefficient of $E_{k+1}^{n+1}$ is clearly nonnegative for any valid value of $\lambda$ if $w_k^n\le 0$.

Therefore the only concerns are for the coefficient of $E_{k-1}^{n+1}$ when $w_k^n\le 0$ and the coefficient of $E_{k+1}^{n+1}$ when $w_k^n\ge 0$. But if $\lambda$ is smaller than $2k_i/(|w_k^n|{\rho }_i{c}_i\mathrm{\Delta }z)$ then

$$R-\nu |w_k^n|(\lambda /2)=\frac{k_i\mathrm{\Delta }t}{{\rho }_i{c}_i\mathrm{\Delta }{z}^2}-\frac{\mathrm{\Delta }t|w_k^n|}{\mathrm{\Delta }z}\frac{\lambda }{2}\ge \frac{k_i\mathrm{\Delta }t}{{\rho }_i{c}_i\mathrm{\Delta }{z}^2}-\frac{\mathrm{\Delta }t|w_k^n|}{\mathrm{\Delta }z}\frac{k_i}{|w_k^n|{\rho }_i{c}_i\mathrm{\Delta }z}=0.$$

Thus all the coefficients in (93) are nonnegative. On the other hand, in equation (93), all coefficients on the right side sum to

$$1+2R+\nu w_k^n\left\{\begin{array}{c}1-\lambda \\ -1+\lambda \end{array}\right\}=1+2R+\nu w_k^n(1-\lambda )\left\{\begin{array}{c}+1\\ -1\end{array}\right\},$$

which is exactly the coefficient on the left side of (93). It follows that

$$E_k^{n+1}=a_k{E}_k^n+b_k{E}_{k-1}^{n+1}+c_k{E}_{k+1}^{n+1}$$

where $a_k,b_k,c_k$ are positive and $a_k+b_k+c_k=1$. Thus a maximum principle applies [74]. END OF PROOF OF 1.

Proof of 2

As a von Neumann analysis is much more restrictive than the analysis above, we will be brief. Let’s assume the velocity is downward, $w_0<0$; the other case is similar. Equation (89) becomes

(94)$$\begin{array}{rl}& (-R-\nu w_0(\lambda /2))E_{k-1}^{n+1}+(1+2R-\nu w_0(1-\lambda ))E_k^{n+1}\\ & +(-R+\nu w_0(1-\lambda /2))E_{k+1}^{n+1}=E_k^n.\end{array}$$

The heart of the von Neumann analysis is the substitution of a growing or decaying (in time index $n$) oscillatory mode on the grid of spatial wave number $\mu$:

$$E_k^n={\sigma }^n{e}^{i\mu (k\mathrm{\Delta }z)}.$$

Here $k\mathrm{\Delta }z=z_k$ is a grid point. Such a mode is a solution to (94) if and only if

$$\begin{array}{rl}\sigma [& (-R-\nu w_0(\lambda /2))e^{-i\mu \mathrm{\Delta }z}+(1+2R-\nu w_0(1-\lambda ))\\ & +(-R+\nu w_0(\lambda /2))e^{+i\mu \mathrm{\Delta }z}+\nu w_0(1-\lambda )e^{+i\mu \mathrm{\Delta }z}]=1.\end{array}$$

This equation reduces by standard manipulations to

$$\sigma =\frac{1}{1+(4R-2\nu w_0(1-\lambda )){\cos }^2(\mu \mathrm{\Delta }z/2)+i\nu w_0(1-\lambda /2)\sin (\mu \mathrm{\Delta }z)}.$$

Note $4R-2\nu w_0(1-\lambda )\ge 0$ without restrictions on numerical parameters $\mathrm{\Delta }t$, $\mathrm{\Delta }z$, because $w_0<0$ in the case under consideration. Therefore

$$|\sigma {|}^2=\frac{1}{{[1+(4R-2\nu w_0(1-\lambda )){\cos }^2(\mu \mathrm{\Delta }z/2)]}^2+{[\nu w_0(1-\lambda /2)\sin (\mu \mathrm{\Delta }z)]}^2}.$$

This positive number is less than one, so $|\sigma |<1$. It follows that all modes decay exponentially. END OF PROOF OF 2.

Remark about our von Neumann stability analysis

The constant $\lambda$ is carefully chosen in (92) so that the maximum principle 1 applies. On the other hand, both the implicit first-order upwind and the implicit centered difference formulas have unconditional stability in the von Neumann sense. The proof of case 2 above is thus a formality, merely showing that a convex combination of unconditionally stable (von Neumann sense) schemes is still unconditionally stable in the same sense.

Convergence: a consequence of the maximum principle

If we define the pointwise numerical error $e_k^n=E_k^n-E(t_n,x_i,y_j,z_k)$, where $E(\dots )$ is the unknown exact solution (exact enthalpy field) [74], then (93) implies an equality of the form

$$A{e}_k^{n+1}=e_k^n+B_{-}{e}_{k-1}^{n+1}+B_{+}{e}_{k+1}^{n+1}+\mathrm{\Delta }t{\tau }_k^n$$

where ${\tau }_k^n$ is the truncation error of the scheme and $A,B_{\pm }$ are nonnegative coefficients, which need no detail for now other than to note that $1+B_{-}+B_{+}=A$. Letting ${\overline{e}}^n=\underset{k}{max}|e_k^n|$ we have, because of the positivity of coefficients,

(95)$$A|e_k^{n+1}|\le {\overline{e}}^n+(B_{-}+B_{+}){\overline{e}}^{n+1}+\mathrm{\Delta }t{\overline{\tau }}^n$$

for all $k$, where ${\overline{\tau }}^n=\underset{k}{max}|{\tau }_k^n|$. Now let $k$ be the index for which $|e_k^{n+1}|={\overline{e}}^{n+1}$. For that $k$ we can replace $|e_k^{n+1}|$ in equation (95) with ${\overline{e}}^{n+1}$. Subtracting the same quantity from each side of the resulting inequality gives

$${\overline{e}}^{n+1}\le {\overline{e}}^n+\mathrm{\Delta }t{\overline{\tau }}^n,$$

It follows that ${\overline{e}}^n\le C\mathrm{\Delta }t$, for some finite $C$, if ${\overline{e}}^0=0$ [74]. Thus a maximum principle for BOMBPROOF implies convergence in the standard way [74]. This convergence proof has the same assumptions as case 1 in the theorem, and thus it only suggests convergence in any broad range of glaciologically-interesting cases.

Remark on nonzero source term

Now recall we assumed in Theorem 1 that the entire “source” ${\mathrm{\Phi }}_k^n$ was identically zero. Of course this is not realistic. What we understand is provable, however, is that if a numerical scheme for a linear advection/conduction equation

$$u_t+A{u}_x=B{u}_{xx}$$

is stable in the general sense of numerical schemes for partial differential equations (e.g. as defined in subsection 5.5 of [74]) then the same scheme is stable in the same general sense when applied to the equation with (linear) lower order terms:

$$u_t+A{u}_x=B{u}_{xx}+Cu+D.$$

A precise statement of this general fact is hard to find in the literature, to put it mildly, but theorem 2.2.3 of [172] is one interesting case ($B=0$ and $D=0$). But even the form we state with linear term ($Cu+D$) is not adequate to the job because of the strongly-nonlinear dependence of $\mathrm{\Phi }$ on the temperature $T$ [16].

Nonetheless the maximum principle is a highly-desirable form of stability because we can exclude “wiggles” from the finite difference approximations of the conductive and advective terms, even if the complete physics, with strain heating in particular, is not yet shown to be non-explosive. Because the complete physics includes the appearance of the famous “spokes” of EISMINT II, for example, a maximum principle cannot apply too literally. Indeed there is an underlying fluid instability [16], one that means the solution of the continuum equations can include growing “wiggles” which are fluid features (though not at the grid-based spatial frequency of the usual numerical wiggles). Recall that, because we use first-order upwinding on the horizontal advection terms, we can expect maximum principle-type stability behavior of the whole three-dimensional scheme.

Temperate basal boundary condition, and computing the basal melt rate

At the bottom of grounded ice, a certain amount of heat comes out of the earth and either enters the ice through conduction or melts the base of the ice. On the one hand, see the documentation for BedThermalUnit for the model of how much comes out of the earth. On the other hand, [22] includes a careful analysis of the subglacial layer equation and the corresponding boundary conditions and basal melt rate calculation, and the reader should consult that reference.

Regarding the floating case

The shelf base temperature $T_{sb}$ and the melt rate $M$ are supplied by OceanModel. Note that we make the possibly-peculiar physical choice that the shelf base temperature is used as the temperature at the top of the bedrock, which is actually the bottom of the ocean. This choice means that there should be no abrupt changes in top-of-bedrock heat flux as the grounding line moves. This choice also means that the conservation of energy code does not need to know about the bedrock topography or the elevation of sea level. (In the future OceanModel could have a subshelf_bed_temperature() method.)

---

Previous

Up

Next