stokes-ice-tutorial

purpose

Provide an easy-to-modify, well-understood, and well-documented finite element solver for glacier modeling, for planar and 3D glacier geometries, without shallowness assumptions in the equations.

This is a tutorial, not a library, but it can be used as a starting point for other projects.

model equations

The power-law Stokes equations describe the dynamics of ice in a glacier or ice sheet as a gravity-driven, viscous, shear-thinning flow:

$$\begin{align*} -\nabla \cdot \tau + \nabla p &= \rho_i \mathbf{g} & &\text{\textsf{stress balance}} \\ \nabla \cdot \mathbf{u} &= 0 & &\text{\textsf{incompressibility}} \\ \tau &= B_n |D\mathbf{u}|^{(1/n)-1} D\mathbf{u} & &\text{\textsf{Glen-Steinemann-Nye flow law}} \end{align*}$$

The free surface of a glacier is governed by an additional equation, coupled to the above, which has the surface mass balance as a source term:

$$\begin{align*} \frac{\partial s}{\partial t} - \mathbf{u}|_s \cdot \mathbf{n}_s &= a \hspace{10mm} \text{\textsf{free-surface kinematics}} \end{align*}$$

This repository contains a practical tutorial on solving these coupled partial differential equations numerically, by using the finite element method.

The Python programs here are relatively-short and only solve idealized problems, but small modification to sources and data will make them real glacier models. We model 2D and 3D land-based glaciers with moving margins. The emphasis is on modern and robust solver techniques, but not directly on scalable solvers.

The codes are documented by slides.pdf.

prerequisites

We use the Firedrake library, which calls PETSc to solve the equations. See the Firedrake install directions to install both libraries. Note that each time you start Firedrake you will need to activate the Python virtual environment:

source venv-firedrake/bin/activate

You will need Gmsh to generate certain meshes (but only for stage1/ and stage2/). For all stages Paraview is recommended to visualize the results.

tutorial in stages

To do this tutorial, read slides.pdf and follow the stages. Each stage is essentially self-contained, with increasing sophistication.

The first three stages use fixed geometry and address how to solve the power-law Stokes equations for the dynamics of ice.

• In stage1/ and stage2/, mesh-generation and Stokes solution are separate actions. Mesh generation uses Gmsh on .geo geometry-outline files, which generates .msh mesh files.
• stage3/ combines meshing and solving into a single program by using extruded meshes; mesh management is thus entirely through Firedrake. (The base mesh could be read from a .msh if desired, by small code modifications.)
• stage3/ also allows choice of 2D (planar) or full 3D glacier geometry at runtime.

By contrast, stage4/ solves the full coupled system above, though only in 2D (planar). That is, we combine the power-law Stokes dynamics model with the free-surface equation.

• The coupled system is solved by time-stepping for a moving-margin case (by default); the PDE system includes a variational inequality.
• Several advanced and recent stabilization techniques are applied in this stage.
• The glacier could evolve in response to a climate (surface mass balance) model; this would require small code modifications.

All stages can be run in parallel.

known limitations

• We do not use any observational data from, and thus we do not actually model, any real glaciers.
• We do not model sliding or floating ice.
• The stages all currently use direct solvers, which limits scalability, especially in 3D.
• The time-stepping in stage4/, though apparently very well-behaved under the applied stabilization and adaptive techniques, remains CFL-limited.
• In the moving-margin cases for stage4/, there is a mass-conservation error committed at the free boundary.
• In stage3/ and stage4/, reading a .msh for the base mesh is easy but it requires code modifications.
• In stage4/, adding surface mass balance data, or a lapse-rate model, is easy but it requires code modifications.

version log

• v2.4 Return links to extra slides, better content/text.
• v2.2|3 Better text/citations.
• v2.1 Switch pressure elements.
• v2.0 Well-behaved glacial evolution in stage4/ using Stockholm stabilizations.
• v1.0 The first version of this tutorial was for glacier dynamics only, for fixed geometry.

other solvers based on Firedrake

Here are two open-source glacier solvers using Firedrake:

• Icepack by Dan Shapero and others is a general-purpose glacier and ice-sheet modeling framework based on Firedrake.

• I wrote an earlier open-source Stokes solver for glaciology, in the py/stokes/ directory in my McCarthy materials. It uses the same Firedrake/PETSc/Gmsh/Paraview stack.