# Blatter stress balance solver: technical details¶

## Notation¶

\(u\) |
\(x\)-component of ice velocity |

\(v\) |
\(y\)-component of ice velocity |

\(\uu\) |
2D vector \((u, v)\) |

\(\n\) |
outward-pointing normal vector at domain boundaries |

\(\eta\) |
ice viscosity (see (101)) |

\(\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 |

\(\varepsilon_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

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

Here “\(\nabla\cdot\)” is the three-dimensional divergence operator and the regularized ice viscosity \(\eta\) is defined by

where \(\gamma_{\text{BP}}\) is the Blatter-Pattyn approximation of the second invariant of the strain rate tensor:

It is also assumed that

ice is incompressible and \(\mathop{trace}(\dot \epsilon) = 0\), and

\(\diff{w}{x} \ll \diff{u}{x}\) and \(\diff{w}{z} \ll \diff{u}{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 (101) appears in [56] while [175] and [73] define the effective strain rate \(\dot\epsilon_e\):

Meanwhile [177] have

## Weak form¶

Note

Recall the product rule

and the divergence theorem:

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 (100) by a *scalar* test
function \(\psi\) and integrate over the domain \(\Omega\). For the first equation, we get

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

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

where \(\nabla s = ( s_x, s_y )\).

The first term corresponds to Neumann and Robin boundary conditions; it vanishes for “natural” BC \(\left(2\, \eta\, \E \right) \cdot \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:

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.

The interface between ice and the air1 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.

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

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

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

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

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 \(\n = (\cdot,\cdot,0)\).

The Neumann boundary condition at marine margins is

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 (108) 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
(108) 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 (105), we create a hexahedral mesh on the domain \(\Omega\) and use \(Q_1\) Galerkin finite elements.

Let \(\phi_j\) be the scalar trial function associated with the node \(j\), then the FE approximation of the solution \(\uu\) has the form

Then the problem is

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

\[{-\Ib \left(\psi_i 2\, \eta\, \E\right) \cdot \n\, ds} + {\Id \nabla \psi_i \cdot 2\, \eta\, \E} + {\Id \psi_i\, \rho\, g\, \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 (105) 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

where

(compare to (105)) 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\cdot3\cdot3 = 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:

The

*basal boundary*term implementing the basal drag (part of \(F^1\))The “main” part (\(F^2\))

The

*source term*corresponding to the driving stress (\(F^3\))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\))The

*marine boundary*part implementing stress (Neumann) BC at the calving front (part of \(F^1\)).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 (\(\Omega\)) below (see (112)) for
simplicity. In actuality each integral is over the *intersection of supports of test
and trial functions* appearing in the integrand.

#### Main residual contribution¶

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

#### Driving stress contribution¶

#### Basal contribution¶

#### Marine boundary contribution¶

where

#### Residual at Dirichlet BC locations¶

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 (111) and (112)). The only other non-trivial contribution comes from the basal boundary condition.

This Jacobian contribution has four parts:

with \(F^2_{\cdot, \cdot}\) defined by (112).

Let \(G_{i, k} = \nabla\psi_i \cdot 2\, \E_{k}\), then

and (using the product rule) we get

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

Taking derivatives of \(\eta\) and \(\gamma\) (101) gives

and

The derivatives of \(G_{\cdot, \cdot}\) are

To compute \(\diff{\gamma}{u_j}\) (119) and \(\diff{G_{\cdot, \cdot}}{\cdot}\) (120) we use FE basis expansions of \(u\) and \(v\) (109), which imply:

and so on.

#### Basal contribution to the Jacobian¶

Here \(\dbeta\) is the derivative of \(\beta\) with respect to \(\alpha = \frac12 |\uu|^2 = \frac12 \left( u^2 + v^2 \right)\) (one of the outputs of PISM’s basal resistance parameterizations).

Note that

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

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

For example,

#### 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 (100) 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

Identify all the connected patches of ice.

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

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

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

```
pismr -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 (\(\dx\), 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 (\(\dx\), \(\dy\)) 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 \(\varepsilon_0\), find an approximate solution, then use it as an initial guess for the next solve with a reduced \(\varepsilon_0\).

Note

Not implemented yet.

### Model domain and mesh structure¶

The domain is

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

Each element’s projection onto the plane \(z = 0\) is a rectangle with sides \(\dx\) and
\(\dy\), 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 \(\Hmin\) (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 \(\Hmin\),

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

boundarynodes.

### 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¶

Compute ice bottom elevation.

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

Compute node type.

Compute ice hardness at the nodes of the mesh.

Call PETSc’s

`SNESSolve()`

.Extract basal velocity and compute basal frictional heating.

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{align}\begin{aligned}u &= \exp(x) \sin(2 \pi y)\\v &= \exp(x) \cos(2 \pi y)\end{aligned}\end{align} \]
- 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{align}\begin{aligned}u &= \frac{2 A (g \rho)^{n} \left((s - z)^{n + 1} - H^{n + 1}\right) \left|s_x\right|^{n - 1} s_x}{n + 1} - \frac{H g \rho s_x}{\beta}\\v &= 0\end{aligned}\end{align} \]
- Domain
- \[ \begin{align}\begin{aligned}x &\in [-L, L],\\z &\in [s(x) - H, s(x)],\\s(x) &= s_0 - \alpha\, x^2.\end{aligned}\end{align} \]
- 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{align}\begin{aligned}u &= \frac{(\rho - \rho_w) g L}{2 B \pi} \sin\p{\frac{\pi x}L} z,\\v &= 0\end{aligned}\end{align} \]
- 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]:

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

for some positive constant \(\alpha\). (A free-floating shelf corresponds to \(\alpha = 1 - \rho_i / \rho_w\)).

Note

Functions \((u, v)\) in (124) solve (100) *exactly* in the
interior of the domain

if the Glen exponent \(n = 3\). No compensatory source term is needed.

We use a Dirichlet BC at the left boundary:

and a stress BC at the right boundary:

The stress BC at the top surface is

The boundary condition at the bottom surface has the form

## 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 |