$$$$

Blatter stress balance solver: technical details

Notation

$u$	$x$-component of ice velocity
$v$	$y$-component of ice velocity
$\mathit{u}$	2D vector $(u,v)$
$\mathit{n}$	outward-pointing normal vector at domain boundaries
$\eta$	ice viscosity (see (102))
$\rho$	ice density
${\rho }_w$	water (usually sea water) density
$g$	gravitational acceleration
$z_{\text{sea level}}$	sea level elevation
$H$	ice thickness
$s$	ice surface elevation
$B$	ice hardness
$A$	ice softness
${\epsilon }_0$	regularization parameter
$\beta$	basal resistance coefficient
$n$	Glen exponent

Introduction

This implementation is based on the PETSc example snes/tutorials/ex48.c (see [56]) and is inspired by [175], [73], [75], [72], and [176].

Define

$$\begin{array}{r}\begin{array}{rl}{\dot{\mathit{ϵ}}}_1& =(2u_x+v_y,\frac{1}{2}(u_y+v_x),\frac{1}{2}{u}_z),\\ {\dot{\mathit{ϵ}}}_2& =(\frac{1}{2}(u_y+v_x),u_x+2v_y,\frac{1}{2}{v}_z),\\ \dot{\mathit{ϵ}}& =({\dot{\mathit{ϵ}}}_1,{\dot{\mathit{ϵ}}}_2).\end{array}\end{array}$$

Using this notation, the Blatter-Pattyn (BP) stress balance equations are

(101)$$\begin{array}{r}\begin{array}{rl}-\mathrm{\nabla }\cdot (2\eta {\dot{\mathit{ϵ}}}_1)+\rho g{s}_x& =0,\\ -\mathrm{\nabla }\cdot (2\eta {\dot{\mathit{ϵ}}}_2)+\rho g{s}_y& =0.\end{array}\end{array}$$

Here “$\mathrm{\nabla }\cdot$” is the three-dimensional divergence operator and the regularized ice viscosity $\eta$ is defined by

(102)$$\begin{array}{r}\begin{array}{rl}\eta & =\frac{B}{2}{({\gamma }_{\text{BP}}+\frac{{\epsilon }_0}{2})}^{\frac{1-n}{2n}},\\ {\gamma }_{\text{BP}}& =u_x^2+v_y^2+u_x{v}_y+\frac{1}{4}(u_y+v_x{)}^2+\frac{1}{4}{u}_z^2+\frac{1}{4}{v}_z^2,\end{array}\end{array}$$

where ${\gamma }_{\text{BP}}$ is the Blatter-Pattyn approximation of the second invariant of the strain rate tensor:

(103)$$\begin{array}{r}\begin{array}{rl}{\dot{ϵ}}_{ij}& =\frac{1}{2}(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i})\\ \gamma & =\frac{1}{2}(trace({\dot{ϵ}}^2)-trace(\dot{ϵ}{)}^2).\end{array}\end{array}$$

It is also assumed that

• ice is incompressible and $trace(\dot{ϵ})=0$, and

• $\frac{\partial w}{\partial x}\ll \frac{\partial u}{\partial x}$ and $\frac{\partial w}{\partial z}\ll \frac{\partial u}{\partial y}$.

The BP approximation of the second invariant ${\gamma }_{\text{BP}}$ is obtained by omitting $w_x$ and $w_y$ and expressing $w_z$ using incompressibility.

Note

There are at least three equivalent expressions for ice viscosity in the literature: the form in (102) appears in [56] while [175] and [73] define the effective strain rate ${\dot{ϵ}}_e$:

$${\dot{ϵ}}_e^2={\dot{ϵ}}_{xx}^2+{\dot{ϵ}}_{yy}^2+{\dot{ϵ}}_{xx}{\dot{ϵ}}_{yy}+{\dot{ϵ}}_{xy}^2+{\dot{ϵ}}_{xz}^2+{\dot{ϵ}}_{yz}^2.$$

Meanwhile [177] have

$$\begin{array}{r}\begin{array}{rl}{\dot{ϵ}}_{\text{BP}}^2& ={\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2+{(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y})}^2\\ & +\frac{1}{2}{(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x})}^2+\frac{1}{2}{\left(\frac{\partial u}{\partial z}\right)}^2+\frac{1}{2}{\left(\frac{\partial v}{\partial z}\right)}^2.\end{array}\end{array}$$

Weak form

Note

Recall the product rule

$$\begin{array}{r}\begin{array}{rl}\mathrm{\nabla }\cdot (fX)& =\mathrm{\nabla }f\cdot X+f\mathrm{\nabla }\cdot X,\text{ or}\\ f\mathrm{\nabla }\cdot X& =\mathrm{\nabla }\cdot (fX)-\mathrm{\nabla }f\cdot X\end{array}\end{array}$$

and the divergence theorem:

$${\int }_{\mathrm{\Omega }}\mathrm{\nabla }\cdot (fX)={\int }_{\partial \mathrm{\Omega }}(fX)\cdot \mathit{n}ds.$$

We omit discussions of function spaces; see [175], [176] and other references mentioned above for details.

To obtain the weak form, we multiply both equations in (101) by a scalar test function $\psi$ and integrate over the domain $\mathrm{\Omega }$. For the first equation, we get

(104)$$\begin{array}{r}\begin{array}{rl}{\int }_{\mathrm{\Omega }}\psi (-\mathrm{\nabla }\cdot (2\eta {\dot{\mathit{ϵ}}}_1)+\rho g{s}_x)& =0,\\ -{\int }_{\mathrm{\Omega }}\psi \mathrm{\nabla }\cdot \left(2\eta {\dot{\mathit{ϵ}}}_1\right)+{\int }_{\mathrm{\Omega }}\psi \rho g{s}_x& =0,\\ -{\int }_{\mathrm{\Omega }}(\mathrm{\nabla }\cdot \left(\psi 2\eta {\dot{\mathit{ϵ}}}_1\right)-\mathrm{\nabla }\psi \cdot 2\eta {\dot{\mathit{ϵ}}}_1)+{\int }_{\mathrm{\Omega }}\psi \rho g{s}_x& =0,\\ -{\int }_{\mathrm{\Omega }}\mathrm{\nabla }\cdot \left(\psi 2\eta {\dot{\mathit{ϵ}}}_1\right)+{\int }_{\mathrm{\Omega }}\mathrm{\nabla }\psi \cdot 2\eta {\dot{\mathit{ϵ}}}_1+{\int }_{\mathrm{\Omega }}\psi \rho g{s}_x& =0,\\ -{\int }_{\partial \mathrm{\Omega }}\left(\psi 2\eta {\dot{\mathit{ϵ}}}_1\right)\cdot \mathit{n}ds+{\int }_{\mathrm{\Omega }}\mathrm{\nabla }\psi \cdot 2\eta {\dot{\mathit{ϵ}}}_1+{\int }_{\mathrm{\Omega }}\psi \rho g{s}_x& =0.\end{array}\end{array}$$

Similarly, multiplying the second equation by $\psi$ and integrating by parts yields

(105)$$-{\int }_{\partial \mathrm{\Omega }}\left(\psi 2\eta {\dot{\mathit{ϵ}}}_2\right)\cdot \mathit{n}ds+{\int }_{\mathrm{\Omega }}\mathrm{\nabla }\psi \cdot 2\eta {\dot{\mathit{ϵ}}}_2+{\int }_{\mathrm{\Omega }}\psi \rho g{s}_y=0.$$

We combine these and say that the weak form of (101) is

(106)$$-{\int }_{\partial \mathrm{\Omega }}\left(\psi 2\eta \dot{\mathit{ϵ}}\right)\cdot \mathit{n}ds+{\int }_{\mathrm{\Omega }}\mathrm{\nabla }\psi \cdot 2\eta \dot{\mathit{ϵ}}+{\int }_{\mathrm{\Omega }}\psi \rho g\mathrm{\nabla }s=0,$$

where $\mathrm{\nabla }s=(s_x,s_y)$.

The first term corresponds to Neumann and Robin boundary conditions; it vanishes for “natural” BC $\left(2\eta \dot{\mathit{ϵ}}\right)\cdot \mathit{n}=0$. In the basal and lateral cases this stress is nonzero. The third one corresponds to the gravitational driving stress and is replaced by a compensatory term in verification tests that use manufactured solutions.

Boundary conditions

The domain boundary consists of three disjoint parts:

1. The interface between ice and the underlying bed or (if the ice is floating) water.

   At this interface PISM uses Robin BC implementing basal sliding.

2. The interface between ice and the air[1] above it.

   The integral over this part of the boundary is zero because we assume that natural (“no stress”) boundary conditions apply. (We ignore atmospheric pressure.)

   Vertical “cliffs” at grounded margins are interpreted as approximations of the very steep, but not vertical, upper surface. Following this interpretation we use natural boundary conditions at grounded lateral margins.

3. The vertical interface between ice and air (above sea level) or water (below sea level) at marine ice margins.

   At this interface PISM uses Neumann BC corresponding to the difference between the cryostatic pressure of the ice on one side and the hydrostatic pressure of water on the other side of the interface.

In addition to this we support Dirichlet boundary conditions for verification and to de-couple unknowns at ice-free locations.

Note

Our implementation supports Dirichlet boundary conditions, but this feature is not exposed to the rest of PISM.

Unlike [56] we do not support “no-slip” BC at the base. This allows us to avoid Jacobian scaling tricks they needed to achieve good multigrid performance.

For each node that belongs to the Dirichlet boundary we assemble “trivial” equations

$$\begin{array}{r}\begin{array}{rl}u& =u_0,\\ v& =v_0,\end{array}\end{array}$$

where $(u_0,v_0)$ are given.

The implementation avoids adding contributions from adjacent elements to residual and Jacobian entries corresponding to Dirichlet locations. Prescribed velocity values are substituted into equations that depend on them (see section 3.2 of [56] for details).

Basal boundary

The boundary condition corresponding to basal sliding is

(107)$$2\eta \dot{\mathit{ϵ}}\cdot \mathit{n}=-\beta \mathit{u}.$$

In the weak form (106) this corresponds to replacing the first term:

(108)$$-{\int }_{\text{base}}(\psi 2\eta \dot{\mathit{ϵ}})\cdot \mathit{n}ds={\int }_{\text{base}}\beta \mathit{u}ds.$$

Here $\beta$ has the same meaning as in (18).

Where ice is grounded $\beta$ is determined as described in Controlling basal strength. It is assumed to be zero (corresponding to no drag) elsewhere.

The grounding line (if present) divides bottom faces of some elements into grounded parts that experience drag and floating parts that do not. This implementation uses a low-order quadrature with many equally-spaced points (a Newton-Cotes quadrature) to integrate over the part of the basal boundary containing the grounding line. Here $\beta$ is computed at each quadrature point, depending on whether the ice is grounded or floating at its location. This is similar to the SEP3 parameterization described in [178].

Note

• It may be a good idea to implement scaling of the $\beta$ coefficient according to the fraction of an element face that is grounded, similar to SEP1 (equation 7 in [178]). An approximation of the grounded fraction for an element column can be computed using existing code in PISM.

• In general the bottom face of an element is not planar and not parallel to the plane $z=0$. This means that integrals over the basal boundary should be evaluated using parameterizations of element faces (i.e. as surface integrals) and not using 2D FEM machinery.

Marine margins

We assume that marine ice margins consist of vertical cliffs, i.e. the outward-pointing normal vector has the form $\mathit{n}=(\cdot ,\cdot ,0)$.

The Neumann boundary condition at marine margins is

(109)$$\begin{array}{r}\begin{array}{rl}2\eta \dot{\mathit{ϵ}}\cdot \mathit{n}& =p_{\text{ice}}-p_{\text{water}},\text{where}\\ p_{\text{ice}}& =\rho g(s-z),\\ p_{\text{water}}& ={\rho }_{w}gmax(z_{\text{sea level}}-z,0).\end{array}\end{array}$$

In other words, this boundary condition corresponds to the difference between the cryostatic pressure of the ice on one side of the interface and the hydrostatic pressure of water on the other. The atmospheric pressure is ignored. Equation (109) is a generalization of equation 18 in [175].

Just like in the implementation of the basal boundary condition near grounding lines, we use a low order quadrature with many equally-spaced points to approximate integrals over element faces intersected by the sea level. This should improve the quality of the approximation of this boundary condition (note that the right hand side of (109) is continuous but not continuously differentiable).

Note

We need to evaluate the importance of the quadrature choice described above. Does using depth-dependent BC matter? (We could simplify the code and use depth-averaged BC it if it does not.) Should we use lots of quadrature points?

Note

CISM 2.1 [73] uses a depth-averaged lateral boundary condition without justification.

The lateral BC described in [175] is equivalent to the one described here, but the implementation in Albany/FELIX (and therefore in MALI) matches the one in CISM.

Implementations in CISM and MALI do not use this boundary condition at grounded margins because doing so appears to produce over-estimates of the ice speed near grounded ice margins. Our experiments show the same behavior.

Solver implementation

Discretization

To create a non-linear algebraic system of equations approximating (106), we create a hexahedral mesh on the domain $\mathrm{\Omega }$ and use $Q_1$ Galerkin finite elements.

Let ${\varphi }_j$ be the scalar trial function associated with the node $j$, then the FE approximation of the solution $\mathit{u}$ has the form

(110)$$\begin{array}{r}\begin{array}{rl}u& =\sum _j{\varphi }_j{u}_j,\\ v& =\sum _j{\varphi }_j{v}_j.\end{array}\end{array}$$

Then the problem is

Find $u_j$, $v_j$ ($j=1,\dots ,N$) so that

$$-{\int }_{\partial \mathrm{\Omega }}\left({\psi }_{i}2\eta \dot{\mathit{ϵ}}\right)\cdot \mathit{n}ds+{\int }_{\mathrm{\Omega }}\mathrm{\nabla }{\psi }_i\cdot 2\eta \dot{\mathit{ϵ}}+{\int }_{\mathrm{\Omega }}{\psi }_i\rho g\mathrm{\nabla }s=0$$

holds for all $i=1,\dots ,N$, where $N$ is the number of nodes in the mesh, subject to the boundary conditions.

As in section 3 of [56], we write the discretization of (106) as an algebraic system $F(U)=0$ with Jacobian $J(U)$ and solve this nonlinear system using Newton iterations requiring approximations of $\delta U$ in

(111)$$J(U)\delta U=-F(U),$$

where

(112)$$\begin{array}{r}\begin{array}{rl}F(U)& =F^1(U)+F^2(U)+F^3(U),\\ F^1& =-{\int }_{\partial \mathrm{\Omega }}(\psi 2\eta \dot{\mathit{ϵ}})\cdot \mathit{n}ds,\\ F^2& ={\int }_{\mathrm{\Omega }}\mathrm{\nabla }\psi \cdot 2\eta \dot{\mathit{ϵ}},\\ F^3& ={\int }_{\mathrm{\Omega }}\psi \rho g\mathrm{\nabla }s.\end{array}\end{array}$$

(compare to (106)) and $J(U)$ is a square sparse matrix containing one row per node in the mesh and at most 54 non-zero entries per row (there are $2$ unknowns per node and each node belongs to at most 8 elements forming a $3\cdot 3\cdot 3=27$ node “neighborhood”).

Residual evaluation

The residual evaluation is performed in the usual manner, by iterating over all the elements containing ice (see Model domain and mesh structure) and adding element contributions to the “global” residual vector.

The residual itself can be broken up into the following parts:

1. The basal boundary term implementing the basal drag (part of $F^1$)

2. The “main” part ($F^2$)

3. The source term corresponding to the driving stress ($F^3$)

4. The top boundary part (zero in actual simulations because we use the natural BC at the top boundary; can be non-zero in verification tests, part of $F^1$)

5. The marine boundary part implementing stress (Neumann) BC at the calving front (part of $F^1$).

6. Residual at Dirichlet nodes.

This decomposition makes it possible to use source terms dictated by the choice of a manufactured solution while keeping (and testing) the rest of the code.

Note

We integrate over the whole domain ($\mathrm{\Omega }$) below (see (113)) for simplicity. In actuality each integral is over the intersection of supports of test and trial functions appearing in the integrand.

Main residual contribution

(113)$$\begin{array}{r}\begin{array}{rl}{F}_{i,u}^2& ={\int }_{\mathrm{\Omega }}\mathrm{\nabla }{\psi }_i\cdot 2\eta {\dot{\mathit{ϵ}}}_1\\ & ={\int }_{\mathrm{\Omega }}\eta (\frac{\partial {\psi }_i}{\partial x}(4u_x+2v_y)+\frac{\partial {\psi }_i}{\partial y}(u_y+v_x)+\frac{\partial {\psi }_i}{\partial z}{u}_z),\\ F_{i,v}^2& ={\int }_{\mathrm{\Omega }}\mathrm{\nabla }{\psi }_i\cdot 2\eta {\dot{\mathit{ϵ}}}_2\\ & ={\int }_{\mathrm{\Omega }}\eta (\frac{\partial {\psi }_i}{\partial y}(2u_x+4v_y)+\frac{\partial {\psi }_i}{\partial x}(u_y+v_x)+\frac{\partial {\psi }_i}{\partial z}{v}_z).\end{array}\end{array}$$

Here $F_{i,u}$ if the contribution to the $u$-component of the residual at the $i$-th node, etc.

Driving stress contribution

(114)$$\begin{array}{r}\begin{array}{rl}{F}_{i,u}^3& ={\int }_{\mathrm{\Omega }}{\psi }_i\rho g{s}_x,\\ F_{i,v}^3& ={\int }_{\mathrm{\Omega }}{\psi }_i\rho g{s}_y.\end{array}\end{array}$$

Basal contribution

(115)$$\begin{array}{r}\begin{array}{rl}{F}_{i,u}^1& ={\int }_{\text{base}}{\psi }_i\beta u,\\ F_{i,v}^1& ={\int }_{\text{base}}{\psi }_i\beta v.\end{array}\end{array}$$

Marine boundary contribution

(116)$$\begin{array}{r}\begin{array}{rl}{F}_{i,u}^1& =-{\int }_{\partial \mathrm{\Omega }}{\psi }_i(p_{\text{ice}}-p_{\text{water}}){\mathit{n}}_x,\\ F_{i,v}^1& =-{\int }_{\partial \mathrm{\Omega }}{\psi }_i(p_{\text{ice}}-p_{\text{water}}){\mathit{n}}_y,\end{array}\end{array}$$

where

$$\begin{array}{r}\begin{array}{rl}{p}_{\text{ice}}& =\rho g(s-z),\\ p_{\text{water}}& ={\rho }_{w}gmax(z_{\text{sea level}}-z,0).\end{array}\end{array}$$

Residual at Dirichlet BC locations

(117)$$\begin{array}{r}\begin{array}{rl}{F}_{i,u}& =u-u_0,\\ F_{i,v}& =v-v_0,\end{array}\end{array}$$

where $(u_0,v_0)$ is the prescribed velocity.

Jacobian evaluation

We use an analytical (as opposed to approximated using finite differences or automatic differentiation) Jacobian computed using formulas listed below.

Here we focus on the derivation of the Jacobian contribution corresponding to the “main” part of the residual ($F^2$, see (112) and (113)). The only other non-trivial contribution comes from the basal boundary condition.

This Jacobian contribution has four parts:

$$\begin{array}{r}\begin{array}{rl}{K}_{i,j,u,u}& =\frac{\partial F_{i,u}^2}{\partial u_j},\\ K_{i,j,u,v}& =\frac{\partial F_{i,u}^2}{\partial v_j},\\ K_{i,j,v,u}& =\frac{\partial F_{i,v}^2}{\partial u_j},\\ K_{i,j,v,v}& =\frac{\partial F_{i,v}^2}{\partial v_j},\end{array}\end{array}$$

with $F_{\cdot ,\cdot }^2$ defined by (113).

Let $G_{i,k}=\mathrm{\nabla }{\psi }_i\cdot 2{\dot{\mathit{ϵ}}}_k$, then

$$F_{i,u}^2={\int }_{\mathrm{\Omega }}\eta G_{i,1}$$

and (using the product rule) we get

(118)$$\begin{array}{r}\begin{array}{rl}{K}_{i,j,u,u}& ={\int }_{\mathrm{\Omega }}\eta \frac{\partial G_{i,1}}{\partial u_j}+\frac{\partial \eta }{\partial u_j}{G}_{i,1},\\ K_{i,j,u,v}& ={\int }_{\mathrm{\Omega }}\eta \frac{\partial G_{i,1}}{\partial v_j}+\frac{\partial \eta }{\partial v_j}{G}_{i,1},\\ K_{i,j,v,u}& ={\int }_{\mathrm{\Omega }}\eta \frac{\partial G_{i,2}}{\partial u_j}+\frac{\partial \eta }{\partial u_j}{G}_{i,2},\\ K_{i,j,v,v}& ={\int }_{\mathrm{\Omega }}\eta \frac{\partial G_{i,2}}{\partial v_j}+\frac{\partial \eta }{\partial v_j}{G}_{i,2}.\end{array}\end{array}$$

The derivatives of $\eta$ are computed using the chain rule:

$$\begin{array}{r}\begin{array}{rl}\frac{\partial \eta }{\partial u_j}& =\frac{\partial \eta }{\partial \gamma }\frac{\partial \gamma }{\partial u_j},\\ \frac{\partial \eta }{\partial v_j}& =\frac{\partial \eta }{\partial \gamma }\frac{\partial \gamma }{\partial v_j}.\end{array}\end{array}$$

Taking derivatives of $\eta$ and $\gamma$ (102) gives

(119)$$\begin{array}{r}\begin{array}{rl}\frac{\partial \eta }{\partial \gamma }& =\frac{B}{2}\cdot \frac{1-n}{2n}\cdot {(\gamma +\frac{{\epsilon }_0}{2})}^{\frac{1-n}{2n}-1}\\ & =\eta \cdot \frac{1-n}{2n}\cdot {(\gamma +\frac{{\epsilon }_0}{2})}^{-1}.\end{array}\end{array}$$

and

(120)$$\begin{array}{r}\begin{array}{rl}\frac{\partial \gamma }{\partial u_j}& =2u_x\frac{\partial u_x}{\partial u_j}+2v_y\frac{\partial v_y}{\partial u_j}+\frac{\partial u_x}{\partial u_j}{v}_y+u_x\frac{\partial v_y}{\partial u_j}+\frac{1}{4}\cdot 2(u_y+v_x)(\frac{\partial u_y}{\partial u_j}+\frac{\partial v_x}{\partial u_j})\\ & +\frac{1}{4}\cdot 2u_z\frac{\partial u_z}{\partial u_j}+\frac{1}{4}\cdot 2v_z\frac{\partial v_z}{\partial u_j}\\ & =2u_x\frac{\partial {\phi }_j}{\partial x}+v_y\frac{\partial {\phi }_j}{\partial x}+\frac{1}{2}\frac{\partial {\phi }_j}{\partial y}(u_y+v_x)+\frac{1}{2}{u}_z\frac{\partial {\phi }_j}{\partial z},\\ \frac{\partial \gamma }{\partial v_j}& =2u_x\frac{\partial u_x}{\partial v_j}+2v_y\frac{\partial v_y}{\partial v_j}+\frac{\partial u_x}{\partial v_j}{v}_y+u_x\frac{\partial v_y}{\partial v_j}+\frac{1}{4}\cdot 2(u_y+v_x)(\frac{\partial u_y}{\partial v_j}+\frac{\partial v_x}{\partial v_j})\\ & +\frac{1}{4}\cdot 2u_z\frac{\partial u_z}{\partial v_j}+\frac{1}{4}\cdot 2v_z\frac{\partial v_z}{\partial v_j}\\ & =2v_y\frac{\partial {\phi }_j}{\partial y}+u_x\frac{\partial {\phi }_j}{\partial y}+\frac{1}{2}\frac{\partial {\phi }_j}{\partial x}(u_y+v_x)+\frac{1}{2}{v}_z\frac{\partial {\phi }_j}{\partial z}.\end{array}\end{array}$$

The derivatives of $G_{\cdot ,\cdot }$ are

(121)$$\begin{array}{r}\begin{array}{rl}\frac{\partial G_{i,1}}{\partial u_j}& =4\frac{\partial {\psi }_i}{\partial x}\frac{\partial {\phi }_j}{\partial x}+\frac{\partial {\psi }_i}{\partial y}\frac{\partial {\phi }_j}{\partial y}+\frac{\partial {\psi }_i}{\partial z}\frac{\partial {\phi }_j}{\partial z},\\ \frac{\partial G_{i,1}}{\partial v_j}& =2\frac{\partial {\psi }_i}{\partial x}\frac{\partial {\phi }_j}{\partial y}+\frac{\partial {\psi }_i}{\partial y}\frac{\partial {\phi }_j}{\partial x},\\ \frac{\partial G_{i,2}}{\partial u_j}& =2\frac{\partial {\psi }_i}{\partial y}\frac{\partial {\phi }_j}{\partial x}+\frac{\partial {\psi }_i}{\partial x}\frac{\partial {\phi }_j}{\partial y},\\ \frac{\partial G_{i,2}}{\partial v_j}& =4\frac{\partial {\psi }_i}{\partial y}\frac{\partial {\phi }_j}{\partial y}+\frac{\partial {\psi }_i}{\partial x}\frac{\partial {\phi }_j}{\partial x}+\frac{\partial {\psi }_i}{\partial z}\frac{\partial {\phi }_j}{\partial z}.\end{array}\end{array}$$

To compute $\frac{\partial \gamma }{\partial u_j}$ (120) and $\frac{\partial G_{\cdot ,\cdot }}{\partial \cdot }$ (121) we use FE basis expansions of $u$ and $v$ (110), which imply:

$$\begin{array}{r}\begin{array}{rl}\frac{\partial u}{\partial u_j}& ={\varphi }_j,\\ \frac{\partial u_x}{\partial u_j}& =\frac{\partial {\varphi }_j}{\partial x},\end{array}\end{array}$$

and so on.

Basal contribution to the Jacobian

(122)$$\begin{array}{r}\begin{array}{rl}{K}_{i,j,u,u}& ={\int }_{\text{base}}{\psi }_i{\varphi }_j(\beta +\frac{\partial \beta }{\partial \alpha }{u}^2),\\ K_{i,j,u,v}& ={\int }_{\text{base}}{\psi }_i{\varphi }_j\frac{\partial \beta }{\partial \alpha }uv,\\ K_{i,j,v,u}& ={\int }_{\text{base}}{\psi }_i{\varphi }_j\frac{\partial \beta }{\partial \alpha }vu,\\ K_{i,j,v,v}& ={\int }_{\text{base}}{\psi }_i{\varphi }_j(\beta +\frac{\partial \beta }{\partial \alpha }{v}^2).\end{array}\end{array}$$

Here $\frac{\partial \beta }{\partial \alpha }$ is the derivative of $\beta$ with respect to $\alpha =\frac{1}{2}|\mathit{u}{|}^2=\frac{1}{2}(u^2+v^2)$ (one of the outputs of PISM’s basal resistance parameterizations).

Note that

$$\begin{array}{r}\begin{array}{rl}\frac{\partial \alpha }{\partial u_j}& =\frac{1}{2}\cdot (2u)\cdot \frac{\partial u}{\partial u_j}\\ & =u{\varphi }_j,\\ \frac{\partial \alpha }{\partial v_j}& =v{\varphi }_j.\end{array}\end{array}$$

Recall (see (115)) that the basal sliding contribution to the residual is equal to

$$\begin{array}{r}\begin{array}{rl}{F}_u^1& ={\int }_{\text{base}}{\psi }_i\beta u,\\ F_v^1& ={\int }_{\text{base}}{\psi }_i\beta v,\end{array}\end{array}$$

and so the basal contribution to the Jacobian consists of partial derivatives of these with respect to $u_j$ and $v_j$.

For example,

$$\begin{array}{r}\begin{array}{rl}{K}_{i,j,u,u}& =\frac{\partial \left({\int }_{\text{base}}{\psi }_i\beta u\right)}{\partial u_j}\\ & ={\int }_{\text{base}}{\psi }_i\frac{\partial \left(\beta u\right)}{\partial u_j}\\ & ={\int }_{\text{base}}{\psi }_i(\beta \frac{\partial u}{\partial u_j}+\frac{\partial \beta }{\partial u_j}u)\\ & ={\int }_{\text{base}}{\psi }_i(\beta {\varphi }_j+\frac{\partial \beta }{\partial \alpha }\cdot \frac{\partial \alpha }{\partial u_j}u)\\ & ={\int }_{\text{base}}{\psi }_i(\beta {\varphi }_j+\frac{\partial \beta }{\partial \alpha }{u}^2{\varphi }_j)\\ & ={\int }_{\text{base}}{\psi }_i{\varphi }_j(\beta +\frac{\partial \beta }{\partial \alpha }{u}^2).\end{array}\end{array}$$

Jacobian at Dirichlet locations

The Jacobian at Dirichlet locations is set to $1$. Together with the residual at Dirichlet locations this completes the assembly of trivial equations at these locations.

More specifically (due to the interleaved ordering of the unknowns: $u_j$, $v_j$, $u_{j+1}$, $v_{v+1}$, …) this requires setting the corresponding block of the Jacobian to the $2\times 2$ identity matrix.

Note

It may be interesting to see if varying the scaling of Jacobian entries at Dirichlet nodes affects the condition number of the Jacobian matrix. See the variable scaling in the code and set

-bp_ksp_view_singularvalues

to see if it has an effect.

Iceberg elimination

As described in [72], isolated patches of ice with low basal resistance and patches connected only via a single node (a “hinge”) are problematic because the system (101) determines ice velocity up to rigid rotations and translations.

This is not a new issue: both FD and FEM solvers of the SSA stress balance require some form of iceberg elimination. We use a connected component labeling algorithm to identify patches of floating ice. In the FEM context this requires inspecting elements: two elements are considered connected if they share a boundary.

Note

We could improve this mechanism by implementing a version of the method described in [72]: instead of removing floating ice not connected to grounded ice (PISM’s approach) they

1. Identify all the connected patches of ice.

2. Remove patches of ice which have no mesh nodes with

   $$\beta >{\beta }_{\text{threshold}}.$$

Preconditioning

Back in 2013 Brown et al [56] showed that multigrid can act as an effective preconditioner for BP stress balance solvers. However, all test cases considered in that work use periodic boundary conditions and therefore avoid all the issues related to the moving ice margin present in complete ice sheet models. Our tests suggests that a naive implementation assembling trivial equations at “ice free” nodes combined with standard geometric multigrid (coarsening in all 3 directions) is not likely to succeed and we need a different approach.

So, we use semi-coarsening in the vertical direction even though Brown et al state that “semi-coarsening is unattractive”. One of the arguments against semi-coarsening is the larger number of multigrid levels needed: semi-coarsening gives a smaller reduction in the number of unknowns from one level to the next (factor of 2 instead of 8 in the full multigrid approach). Our tests show that “aggressive” semi-coarsening (i.e. using coarsening factors larger than 2 and as high as 8) appears to be effective, allowing one to achieve similar reductions in the number of unknowns from one level to the next in a multigrid hierarchy.

The second argument against semi-coarsening is deeper: spatial variations in the sliding parameter $\beta$ may lead to the kind of anisotropy that cannot be addressed by coarsening in the vertical direction (see chapter 7 in [179] for a discussion). Still, we are encouraged by results published by Tuminaro et al [72], who used a similar mesh structure and an approach equivalent to using geometric multigrid with semi-coarsening in the vertical direction for the finer part of the hierarchy and algebraic multigrid for coarser levels. [2]

Inspired by [72], we use geometric multigrid to build a mesh hierarchy with the coarsest level containing a small number (2 or 3) of vertical levels combined with algebraic multigrid as a preconditioner for the solver on the coarsest level. (Semi-coarsening in the vertical direction cannot reduce the number of unknowns in the $x$ and $y$ directions and the coarsest problem is likely to be too large for the redundant direct solver, which is PETSc’s default.)

In addition to this, we follow [56] in ordering unknowns so that columns are contiguous (and $u$ and $v$ are interleaved), allowing ILU factorization to compute a good approximation of ice velocities in areas where SIA is applicable.

Vertical grid sizes compatible with coarsening

Ideally, the coarsest mesh in the hierarchy should have 2 nodes in the $z$ direction, i.e. be one element thick. If $N$ is the coarsening factor, the total number of vertical levels ($M_z$) has to have the form

(123)$$M_z=A\cdot N^k+1$$

for some positive integer $A$ (ideally $A=1$), so that the mesh hierarchy containing $k$ levels will include levels with

$$A+1,A\cdot N+1,A\cdot N^2+1,\dots ,A\cdot N^k+1$$

nodes in the $z$ direction.

This means that for a given stress_balance­.blatter­.coarsening_factor and number of multigrid levels -bp_pc_mg_levels k the value of $M_z$ cannot be chosen at random.

Controlling using PETSc options

The PETSc SNES object solving the BP system uses the bp_ prefix for all command-line options.

To choose a preconditioner, set

-bp_pc_type XXX

where XXX is the name of a preconditioner.

Run

pism -stress_balance blatter [other options] -help | grep "-bp"

to get the list of all PETSc options controlling this solver.

To use a geometric multigrid preconditioner with $N$ levels, set

-bp_pc_type mg -bp_pc_mg_levels N

An “aggressive” (i.e. greater than 2) coarsening factor may work well. Use stress_balance­.blatter­.coarsening_factor to set it.

See Vertical grid sizes compatible with coarsening for the discussion of the relationship between the number of vertical levels, number of multigrid levels, and the coarsening factor.

Set

-bp_mg_coarse_pc_type gamg

to use PETSc’s GAMG on the coarsest multigrid level.

Note

It would be interesting to compare different preconditioning options on the coarsest MG level (GAMG, Hypre BoomerAMG, …).

Additional code needed to support geometric multigrid

To support the multigrid preconditioner we need to be able re-discretize the system on the mesh provided by PETSc in our to residual and Jacobian evaluators. In general, this requires

• re-computing grid-related constants ($\mathrm{\Delta }x$, etc) using the grid (avoid using hard-wired constants, e.g. computed using the fine grid), and

• additional code to restrict gridded inputs from the fine grid mesh to coarser meshes.

This solver does not support coarsening in horizontal directions, so gridded two-dimensional inputs can be used on all multigrid levels. The grid spacing ($\mathrm{\Delta }x$, $\mathrm{\Delta }y$) remains the same as well.

To transfer the one three-dimensional gridded input field (ice hardness), we create interpolation matrices mapping from a coarse level to the next (finer) level in the hierarchy. The transpose of this matrix is used as a restriction operator.

Parameter continuation as a recovery mechanism

As in [175], we can start with a large ${\epsilon }_0$, find an approximate solution, then use it as an initial guess for the next solve with a reduced ${\epsilon }_0$.

Note

Not implemented yet.

Model domain and mesh structure

The domain is

$$\begin{array}{r}\begin{array}{rl}x& \in [x_{\text{min}},x_{\text{max}}],\\ y& \in [y_{\text{min}},y_{\text{max}}],\\ z& \in [z_{\text{min}},z_{\text{min}}+H],\end{array}\end{array}$$

where $[x_{\text{min}},x_{\text{max}}]\times [y_{\text{min}},y_{\text{max}}]$ is the “map plane” domain corresponding to the maximum ice extent, $z_{\text{min}}$ is the bottom ice surface elevation (equal to bed elevation where ice is grounded and determined using sea level, ice thickness, and the floatation criterion where floating) and $H$ is the ice thickness.

Coordinates of the mesh nodes have the form

(124)$$\begin{array}{r}\begin{array}{rl}{x}_i& =x_{\text{min}}+i\cdot \mathrm{\Delta }x,\\ y_j& =y_{\text{min}}+j\cdot \mathrm{\Delta }y,\\ z_k& =z_{\text{min}}(x_i,y_j)+k\cdot \frac{H(x_i,y_j)}{M_z-1}.\end{array}\end{array}$$

Each element’s projection onto the plane $z=0$ is a rectangle with sides $\mathrm{\Delta }x$ and $\mathrm{\Delta }y$, but the spacing between nodes in the vertical direction is not constant: each vertical column of nodes contains $M_z$ nodes with the spacing of $H/(M_z-1)$. This mesh structure is exactly the same as the one used in [56] and CISM 2.1 [73].

Supporting evolving ice extent

The ice extent changes as the model runs and the solver implementation has to allow for these changes.

To simplify the logic used to identify the interior of the ice volume and its lateral boundaries we compute the “type” of each node in the mesh. Given a threshold $H_{\text{min}}$ (see stress_balance­.ice_free_thickness_standard), we say that

• an element contains ice if ice thickness at all its nodes equals to or exceeds $H_{\text{min}}$,

• a node is interior (within the ice) if all the elements it belongs to contain ice,

• a node is exterior (outside the ice volume) if no element it belongs to contains ice,

• a node that is neither interior nor exterior is a boundary node.

  Only elements containing ice are included in the residual and Jacobian evaluation.

We prescribe zero $(0,0)$ velocity at exterior nodes by marking them as Dirichlet BC locations.

In addition to this, we need to identify vertical faces of elements at lateral boundaries.

An element face is a part of the lateral boundary if all four of its nodes are boundary nodes.

Solver inputs and outputs

The BP solver uses the following inputs:

• basal yield stress (${\tau }_c$),

• ice thickness,

• bed elevation,

• sea level elevation,

• ice enthalpy (used to compute ice hardness)

And provides the following outputs:

• $u$ and $v$ components of ice velocity at the nodes of the mesh (saved to output files to be used as an initial guess later)

• vertically-averaged $u$ and $v$ (used in the mass continuity step, i.e. to update ice geometry)

• basal frictional heating (used to couple stress balance to energy conservation)

Steps performed by the solver

1. Compute ice bottom elevation.

2. Compute floatation function $f$ ($f\le 0$ if ice is grounded, $f>0$ if floating).

3. Compute node type.

4. Compute ice hardness at the nodes of the mesh.

5. Call PETSc’s SNESSolve().

6. Extract basal velocity and compute basal frictional heating.

7. Compute vertically-averaged ice velocity.

Integration with the rest of PISM

Conservation of energy

Coupling to PISM’s energy balance models requires

• all 3 components of ice velocity on PISM’s grid, and

• strain heating.

These are computed by using piecewise-linear interpolation in the vertical direction to put $u$, $v$ on PISM’s grid, after which vertical velocity $w$ and strain heating are computed using existing code.

Testing and verification

Verification test XY

Exact solution:

$$\begin{array}{r}\begin{array}{rl}u& =\exp (x)\sin (2\pi y)\\ v& =\exp (x)\cos (2\pi y)\end{array}\end{array}$$

Domain:

$x\in [0,1]$, $y\in [0,1]$

Boundary conditions:

Dirichlet BC corresponding to the exact solution on all lateral boundaries ($x=0$, $x=1$, $y=0$, $y=1$). Natural BC at the top and bottom boundaries.

Comments:

This test uses a constant ice hardness, $n=3$, and has no variation in the $z$ direction. It is similar to the $x-y$ MMS verification test in [175], section 4.1 (we use Dirichlet boundary conditions along the whole lateral boundary instead of Robin conditions derived from the chosen exact solution).

The compensatory term is computed using SymPy; please see the code in src/stressbalance/blatter/verification.

Verification test XZ

Exact solution:

$$\begin{array}{r}\begin{array}{rl}u& =\frac{2A(g\rho {)}^n((s-z{)}^{n+1}-H^{n+1}){\left|s_x\right|}^{n-1}{s}_x}{n+1}-\frac{Hg\rho s_x}{\beta }\\ v& =0\end{array}\end{array}$$

Domain:

$$\begin{array}{r}\begin{array}{rl}x& \in [-L,L],\\ z& \in [s(x)-H,s(x)],\\ s(x)& =s_0-\alpha x^2.\end{array}\end{array}$$

Boundary conditions:

• Dirichlet BC corresponding to the exact solution at $x=-L$ and $x=L$.

• Basal ($z=s(x)-H$) BC is a combination (i.e. sum) of the BC in Basal boundary and a compensatory term derived using the exact solution.

• Top surface ($z=s(x)$) BC is derived using the exact solution.

Comments:

This test

• uses a constant ice hardness,

• Glen exponent $n=3$,

• has no variation (and is periodic) in the $y$ direction,

• uses a constant basal resistance coefficient $\beta$.

It is similar to the $x-z$ MMS verification test [175], section 4.2 (again, we use Dirichlet BC at lateral boundaries instead of Robin conditions stated in the paper).

See [175] and the code in src/stressbalance/blatter/verification for details.

Verification test XZ-CFBC

This setup tests the “calving front boundary condition” (see Marine margins).

Exact solution:

$$\begin{array}{r}\begin{array}{rl}u& =\frac{(\rho -{\rho }_w)gL}{2B\pi }\sin \left(\frac{\pi x}{L}\right)z,\\ v& =0\end{array}\end{array}$$

Domain:

$x\in [0,L]$, $z\in [-H,0]$

Boundary conditions:

• Dirichlet BC at $x=0$.

• Uses the BC described in Marine margins at $x=L$.

Comments:

This test uses the Glen exponent of $1$ (constant viscosity) and has no variation in the $y$ direction.

The sea level is set to $0$, overriding the floatation criterion to ensure that all the ice is submerged.

Verification test XZ-VV (van der Veen profile)

This setup tests the implementation of the basal sliding boundary condition (see Basal boundary) using the van der Veen shelf profile [180]:

(125)$$\begin{array}{r}\begin{array}{rl}H(x)& ={[\frac{4Cx}{Q_0}+H_0^{-4}]}^{-\frac{1}{4}},\\ u(x)& =\frac{Q_0}{H(x)},\\ v(x)& =0,\\ C& ={\left(\frac{\alpha g{\rho }_i}{2B}\right)}^3,\end{array}\end{array}$$

where $Q_0$ is the flux at the left boundary and $H_0$ is the corresponding ice thickness.

The surface elevation $s$ and bed elevation $b$ are defined by

$$\begin{array}{r}\begin{array}{rl}s(x)& =\alpha H(x),\\ b(x)& =(\alpha -1)H(x)\end{array}\end{array}$$

for some positive constant $\alpha$. (A free-floating shelf corresponds to $\alpha =1-{\rho }_i/{\rho }_w$).

Note

Functions $(u,v)$ in (125) solve (101) exactly in the interior of the domain

$$\begin{array}{r}\begin{array}{rl}0\le x_{\text{min}}& \le x\le x_{\text{max}},\\ b(x)& \le z\le s(x)\end{array}\end{array}$$

if the Glen exponent $n=3$. No compensatory source term is needed.

We use a Dirichlet BC at the left boundary:

$$\begin{array}{r}\begin{array}{rl}u(x_{\text{min}})& =\frac{Q_0}{H_0},\\ v(x_{\text{min}})& =0\end{array}\end{array}$$

and a stress BC at the right boundary:

(126)$$2\eta \dot{\mathit{ϵ}}\cdot {\mathit{n}}_{\text{right}}=(\alpha g{\rho }_{i}H(x),0)$$

The stress BC at the top surface is

(127)$$2\eta \dot{\mathit{ϵ}}\cdot {\mathit{n}}_{\text{top}}=(\frac{2B{C}^{\frac{4}{3}}\alpha H^6(x)}{\sqrt{C^2{\alpha }^2{H}^{10}(x)+Q_0^2}},0)$$

The boundary condition at the bottom surface has the form

$$\begin{array}{r}\begin{array}{rl}2\eta \dot{\mathit{ϵ}}\cdot {\mathit{n}}_{\text{bottom}}& =-\beta \mathit{u},\\ \beta & =\frac{2B{C}^{\frac{4}{3}}(\alpha -1)H^7(x)}{Q_0\sqrt{C^2{(\alpha -1)}^2{H}^{10}(x)+Q_0^2}}\end{array}\end{array}$$

Known issues and future work

• Eliminate “wiggles” near areas with steep surface slopes.

• Implement drag at lateral boundaries in fjords and alpine valleys.

• Implement parameter continuation.

• Couple to melange back pressure parameterizations by replacing $p_{\text{water}}$ in Marine margins.

Footnotes

[1]

Strictly speaking the top surface of the ice may be in contact with firn or snow as well as air, but these details are not relevant here.

[2]

The code developed by [72] uses the algebraic multigrid framework throughout, i.e. even for the part of the hierarchy where the mesh structure allows one to use “geometric” coarsening.

---

Previous

Up

Next