%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from numpy import pi as π
import firedrake
from firedrake import as_vector, inner, ds
import icepack, icepack.plot, icepack.models

Synthetic ice shelf

In this demo, we’ll show how to simulate the evolution of a floating ice shelf. The example that we’ll show is an example of a model spin-up. In a spin-up experiment, the objective is to obtain a physically plausible state of some system from an initial guess by simulating its evolution for a long time. For example, it’s common to spin up climate models from a very rough initial guess for several thousand years to obtain a representative state of the atmosphere and oceans for the epoch you’re interested in.

In our case, we’ll use an idealized ice shelf geometry and a rough guess for the initial ice thickness. First we’ll solve the diagnostic equation for the velocity througout the ice shelf. We’ll then solve the prognostic equation to obtain the ice thickness at a slightly later time. By assuming a given rate of ice accumulation and melting, we can successively solve the prognostic and diagnostic equations until the system is close to a steady state. In our case, that will take about 250 years of model time and (at least on my laptop) only a few minutes of computer time.

The whole simulation can be divided into three parts:

  • Define the shape of the ice shelf and get a triangulation of the interior.

  • Define the initial guess for the ice thickness and velocity, and set a value of the rheology and accumulation rate.

  • Set the total simulation time and the number of time steps, and then iteratively update the ice thickness and velocity at each time step.

This is a pretty common workflow for a predictive model. In later demos we’ll see some variations on this procedure – incorporating real data, comparing different physics, and so forth – but the essential idea is the same throughout.


First we need to define the geometry of the ice shelf. Although we’ll be solving a synthetic problem, we’ll try to use a more-or-less realistic shape. The inflow boundary will be a circular arc centered at the origin with a radius of 200 km. The outflow boundary will be another circular arc with a much higher radius.

In the first demo, we used one of firedrake’s built-in functions to create the geometry. For more complicated shapes you’ll need to use a mesh generator, a program that turns a description of the boundary of a spatial domain into a triangulation of the interior. Two of the more popular 2D mesh generators are gmsh and Triangle. In this case we’ll use gmsh because we can create the input file entirely in Python through the package pygmsh.

We’ll first define the mesh radius and the spacing for the mesh cells.

import pygmsh

R = 200e3
δx = 5e3

Next, we’ll create an empty geometry and add some points to it. The points x1, x2 will represent the left and right endpoints of the mesh. The points center1, center2 will not actually be contained in the mesh, but rather define the centers of the two circular arcs. Finally, we’ll add the circular arcs.

geometry = pygmsh.built_in.Geometry()

x1 = geometry.add_point([-R, 0, 0], lcar=δx)
x2 = geometry.add_point([+R, 0, 0], lcar=δx)

center1 = geometry.add_point([0, 0, 0,], lcar=δx)
center2 = geometry.add_point([0, -4 * R, 0], lcar=δx)

arcs = [geometry.add_circle_arc(x1, center1, x2),
        geometry.add_circle_arc(x2, center2, x1)]

Now that we’ve added the geometric elements of our domain, we also need to tell gmsh about the topology, i.e. how all the arcs are connected to each other and how they’re oriented. The physical lines and surfaces are added so that gmsh will tag each geometric entity with a number that we can use to set different boundary conditions.

line_loop = geometry.add_line_loop(arcs)
plane_surface = geometry.add_plane_surface(line_loop)

physical_lines = [geometry.add_physical(arc) for arc in arcs]
physical_surface = geometry.add_physical(plane_surface)

This completes the definition of the input geometry. The method get_code of the geometry object returns the string describing it in the syntax that gmsh expects. We’ll write this string out to a file with the extension .geo.

with open('ice-shelf.geo', 'w') as geo_file:

Next, we’ll call gmsh from the command line on the input file we just created. The mesh generator will read the description of the domain boundary, triangulate the interior of the domain, and output a file with the extension .msh. Other mesh generators have different input and output formats, but the procedure is roughly the same.

In a jupyter notebook, you can use an exclamation mark followed by a command to execute this command at the shell rather than in Python. We’ll call gmsh from the command line with the following arguments:

  • -2: generate a 2D mesh as opposed to 3D

  • -format msh2: specify the storage format of the output file

  • -o ice-shelf.msh: name of the output file

  • ice-shelf.geo: the input data

The shell command (without the exclamation mark) is what you would use if you were working directly from the command line rather than in a notebook.

!gmsh -2 -format msh2 -o ice-shelf.msh ice-shelf.geo
Info    : Running 'gmsh -2 -format msh2 -o ice-shelf.msh ice-shelf.geo' [Gmsh 4.5.6, 1 node, max. 1 thread]
Info    : Started on Mon Apr 27 08:16:19 2020
Info    : Reading 'ice-shelf.geo'...
Info    : Done reading 'ice-shelf.geo'
Info    : Meshing 1D...
Info    : [  0 %] Meshing curve 1 (Circle)
Info    : [ 50 %] Meshing curve 2 (Circle)
Info    : Done meshing 1D (0.002499 s)
Info    : Meshing 2D...
Info    : Meshing surface 4 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (0.15669 s)
Info    : 3373 nodes 6744 elements
Info    : Writing 'ice-shelf.msh'...
Info    : Done writing 'ice-shelf.msh'
Info    : Stopped on Mon Apr 27 08:16:20 2020

The terminal output from gmsh gives us some diagnostics like how many vertices and triangles it contains. This is also where gmsh will report if something went wrong – a syntax error in the .geo file, a degenerate input geometry, and so forth.

To load the mesh we pass the filename to the function firedrake.Mesh. This function will determine the mesh format based on the file extension; it works for meshes that were generated by Triangle or other programs as well.

mesh = firedrake.Mesh('ice-shelf.msh')

The colors in the figure below show how gmsh tagged the calving terminus with ID 2 and the inflow boundary as 1. This is exactly analogous to how firedrake adds tags for each side of the square geometry that we used in the previous demo. These numeric tags help us define Dirichlet (inflow) and Neumann (terminus) boundary conditions where they apply.

fig, axes = icepack.plot.subplots()
icepack.plot.triplot(mesh, axes=axes)

In the demos for real glaciers that follow, we use all of the same tools. The main difference is that the boundary arcs are drawn by hand in a geographic information system, rather than defined programatically in Python. In the repository glacier-meshes I’ve included shapefiles of the outlines of several glaciers and a program to automate the translation of a shapefile into a .geo file using pygmsh. This will be used in the demo for the Larsen Ice Shelf.

Input data

To mimic the state of a real ice shelf, we’ll pick a few angles along the inflow boundary that represent the centerlines of the ice streams that feed the shelf. We’ll then define velocity and thickness profiles along this inflow boundary. We don’t have a great idea from the outset of what the steady state of the ice shelf is; it doesn’t have a simple analytical expression in terms of the coordinates. Instead, we’ll pick a somewhat arbitrary initial profile and evolve it towards steady state.

Many ice shelves (Larsen, Ross, etc.) have several streams feeding them. Our synthetic glacier will be fed by four streams. We’ll define the inlets by picking the angles around the inflow boundary where each inlet comes in from and the width in radians. You can re-run this notebook and change the values or the number of streams to whatever you like.

inlet_angles = π * np.array([-3/4, -1/2, -1/3, -1/6])
inlet_widths = π * np.array([1/8, 1/12, 1/24, 1/12])

Next, we’ll come up with some rather arbitrary and un-physical input data. The basic idea is to make the thickness slope down as you go towards the calving terminus and away from the centerline of an inlet. Likewise the ice speed goes up as you go towards the calving terminus. In order to make this big nasty algebraic expression, we’ll create a list of the perturbation thickness and velocity for each inlet, and combine them all together at the end.

x = firedrake.SpatialCoordinate(mesh)

u_in = 300
h_in = 350
hb = 100
dh, du = 400, 250

hs, us = [], []
for θ, ϕ in zip(inlet_angles, inlet_widths):
    x0 = R * as_vector((np.cos(θ), np.sin(θ)))
    v = -as_vector((np.cos(θ), np.sin(θ)))
    L = inner(x - x0, v)
    W = x - x0 - L * v
    Rn = 2 * ϕ / π * R
    q = firedrake.max_value(1 - (W / Rn)**2, 0)
    hs.append(hb + q * ((h_in - hb) - dh * L /R))
    us.append(firedrake.exp(-4 * (W/R)**2) * (u_in + du * L / R) * v)

To combine the expressions for the thickness and velocity of each inlet into expressions for the whole ice shelf, we’ll take the maximum thickness at any point, and the sum of the velocities.

h_expr = firedrake.Constant(hb)
for h in hs:
    h_expr = firedrake.max_value(h, h_expr)

u_expr = sum(us)

These are merely algebraic expressions. To start modeling we need to interpolate these expressions to some function spaces defined over the mesh.

Q = firedrake.FunctionSpace(mesh, family='CG', degree=2)
V = firedrake.VectorFunctionSpace(mesh, family='CG', degree=2)

h0 = firedrake.interpolate(h_expr, Q)
u0 = firedrake.interpolate(u_expr, V)

The algebra might be a little opaque, but the plots of the initial data should be more illuminating.

fig, axes = icepack.plot.subplots()
contours = icepack.plot.tricontourf(h0, 40, axes=axes)
<matplotlib.colorbar.Colorbar at 0x7fc4e99de690>
fig, axes = icepack.plot.subplots()
contours = icepack.plot.tricontourf(u0, 40, axes=axes)
<matplotlib.colorbar.Colorbar at 0x7fc4e9986ad0>

As the simulation progresses, the ice streams, as represented by areas with higher thickness, will spread and grow.


To actually start solving things, we’ll make an IceShelf object that represents the physical model we’re using. Each model type has two methods, diagnostic_solve and prognostic_solve, that solve the PDEs for the ice velocity and thickness.

ice_shelf = icepack.models.IceShelf()
Help on method diagnostic_solve in module icepack.models.ice_shelf:

diagnostic_solve(u0, h, dirichlet_ids, tol=1e-06, **kwargs) method of icepack.models.ice_shelf.IceShelf instance
    Solve for the ice velocity from the thickness

    u0 : firedrake.Function
        Initial guess for the ice velocity; the Dirichlet boundaries
        are taken from u0
    h : firedrake.Function
        Ice thickness
    dirichlet_ids : list of int
        list of integer IDs denoting the parts of the boundary where
        Dirichlet conditions should be applied
    tol : float
        dimensionless tolerance for when to terminate Newton's method

    u : firedrake.Function
        Ice velocity

    Other parameters
        All other keyword arguments will be passed on to the
        viscosity and gravity functions that were set when this
        model object was initialized

For now we’ll use a fluidity that assumes a spatially constant temperature of -18C. The code itself assumes that all temperatures are in Kelvin. In later demos, we’ll use a spatially variable fluidity coefficient and show how to infer it from observations.

T = firedrake.Constant(255.15)
A = firedrake.interpolate(icepack.rate_factor(T), Q)

Finally, we need to prescribe the boundary conditions for the velocity and the convergence tolerance for the nonlinear solver. Since the boundary conditions and convergence tolerance are usually the same throughout a simulation, we’ll pack them into a dictionary call opts and reuse it throughout.

To compute the velocity field, we call the diagnostic solve method of the ice shelf object. This function takes in the initial guess u0 for the velocity, the thickness, the fluidity, and all the extra options we packed into opts.

opts = {'dirichlet_ids': [1], 'tol': 1e-12}
h = h0.copy(deepcopy=True)
u = ice_shelf.diagnostic_solve(u0=u0, h=h, A=A, **opts)

Here we passed all the arguments to diagnostic_solve as keywords. This convention will be used throughout all the demos. Passing arguments by keyword is slightly more verbose than passing them by position, but it saves you the trouble of remembering what order everything goes in. On top of that, many fields are only passed in as keyword arguments. For example, the fluidity parameter is not a positional argument of the diagnostic solve routine. The reason for this choice, as we’ll see later, is that it’s much easier to swap out components of the model physics for your own customized versions.

The following plot shows streamlines of the velocity field. These kind of plots are useful for showing areas of convergence and divergence of the flow field.

fig, axes = icepack.plot.subplots()
streamlines = icepack.plot.streamplot(u, precision=1e3, density=2e3, axes=axes)
fig.colorbar(streamlines, label='meters/year')

To project the state of the ice shelf forward in time, we’ll use the prognostic solve method. The prognostic solver updates the ice thickness forward by a given timestep given the accumulation rate and velocity. We then update the velocity using the diagnostic solver at each timestep. The following code runs the model forward for several years until the ice shelf is roughly in steady state.

final_time = 400.
num_timesteps = 200
dt = final_time / num_timesteps
a = firedrake.Constant(0.0)

for step in range(num_timesteps + 1):
    h = ice_shelf.prognostic_solve(dt, h0=h, a=a, u=u, h_inflow=h0)
    u = ice_shelf.diagnostic_solve(u0=u, h=h, A=A, **opts)

    print('.' if step % 10 == 0 else '', end='', flush=True)

The main loop to run the simulation consists of alternating calls to the prognostic and diagnostic solve methods of the model object. We’ll see this pattern repeated in all of the demos that follow. In some cases we’ll add in extra physics, like updating the temperature or damage field, but the core idea is the same.

We’ve added a bit of feedback to the loop by printing out a period every 10 iterations. These kinds of diagnostics are helpful to know that the simulation is actually running and that it hasn’t frozen. But you can also put in whatever extra code you want here. For example, you might want to make plots of the thickness and velocity, print out some physical quantity like the total flux of ice out of the calving front, or accumulate the fields into a list so that you can analyze the entire time series later.

To wrap things up, we’ll make a plot of the final ice thickness and velocity. The initial thickness profile of each ice stream, which flattened out in the middle of the shelf, has extended all the way to the terminus.

fig, axes = icepack.plot.subplots()
contours = icepack.plot.tricontourf(h, 40, axes=axes)
fig.colorbar(contours, label='meters')

A stream plot of the ice velocity shows that the areas of greatest flow divergence have moved around relative to the initial velocity.

fig, axes = icepack.plot.subplots()
streamlines = icepack.plot.streamplot(u, precision=1e3, density=2e3, axes=axes)
fig.colorbar(streamlines, label='meters/year')

A useful quantity to know is how much ice is moving through the system. In the previous demo, we showed how to symbolically create functionals of the solution. We can use this to calculate the total ice flux through the calving terminus.

ν = firedrake.FacetNormal(mesh)
flux = h * inner(u, ν) * ds(2)

The extra argument to ds indicates that we want to integrate over just the calving terminus rather than over the entire domain boundary. The variable flux is a symbolic representation of this quantity in terms of the input fields; it isn’t a number but rather a description of how to compute a number. The function firedrake.assemble does the actual work of evaluating the integral.

print('Flux: {} km^3 / year'.format(firedrake.assemble(flux) / 1e9))
Flux: 38.11244102070711 km^3 / year

To check the degree to which the system really is in steady state, we can compute the flux along the inflow boundary rather than along the terminus. This should be equal to the flux out of the domain since we’ve set the accumulation rate to be zero.

influx = -h * inner(u, ν) * ds(1)
print('Influx: {} km^3 / year'.format(firedrake.assemble(influx) / 1e9))
Influx: 38.14274214232329 km^3 / year

The influx and outflux are reasonably close, with the influx slightly exceeding the outflux. The ice shelf will thicken, increasing the flux out of the terminus, until the two values equilibrate.


In the last demo, we saw how to create a mesh, define a finite element space, interpolate functions to that finite element space, and analyze the results. Here we’ve shown how to use these functions as an input to an ice shelf flow model. If you want to experiment with this example a bit more, you can run this notebook yourself and:

  • change the curvature of the inflow boundary or the calving terminus

  • change the stream widths, thicknesses, or velocities

  • add a sinusoidal variation in the accumulation rate to approximate seasonal cycles

In the next demo, we’ll show how to use these functions for simulating a real ice shelf in the Antarctic Peninsula using observational data from satellites. As a bonus epilogue to this demo, we’ll show how to also include the physics of damage and crevassing.


The simulation we just ran uses the prognostic and diagnostic model physics of ice shelf flow. These models successively update the ice thickness and velocity. We assumed that the ice rheology was constant in space in order to simplify things, but this is never true for real ice shelves. Here we’ll show how to add a model for how crevasses and rifts affect glacier flow.

As ice becomes more damaged, the same driving stress must be transferred through a smaller total cross-sectional area of the ice. We can incorporate the effect of crevasses and rifts into large-scale models of glacier velocity by increasing the fluidity of the ice column. The extent of the damage is quantified in a scalar field \(D\) which takes values between 0 (undamaged) and 1 (completely damaged). Following Albrecht and Levermann, we define the damaged rate factor \(A_D\) to be

\[A_D = (1 - D)^{-n}A(T)\]

where \(n = 3\) is the Glen flow law exponent and \(A(T)\) is the rate factor for undamaged ice at temperature \(T\). As the damage field approaches is maximum value of 1, the fluidity goes to infinity. A large-scale calving event would occur before this limit is reached for a real glacier, but we won’t include this effect here.

The ice shelf model includes several functions for calculating the net power dissipation of the ice shelf. The velocity that solves the diagnostic equations minimizes this power dissipation functional. To incorporate damage, we need to create our own function that calculates the viscous power dissipation in the presence of damage. The default function to calculate the viscous power dissipation lives in the module icepack.models.viscosity, so all we need to do is pass the damaged fluidity \((1 - D)^{-n}\cdot A\) instead of just \(A\) itself.

from icepack.constants import glen_flow_law as n
from icepack.models.viscosity import viscosity_depth_averaged
def viscosity_damaged(u, h, A, D):
    return viscosity_depth_averaged(u, h, (1 - D)**(-n) * A)

We can then create a new ice shelf model object that will use viscosity_damaged to calculate the power dissipation instead of the default.

ice_shelf = icepack.models.IceShelf(viscosity=viscosity_damaged)

The procedure is exactly the same to modify any other part of the physics of any model. For example, you could also change how the back-pressure at the calving terminus is computed in order to account for the presence of ice melange. In later tutorials, we’ll look at alternative models for basal friction, and we’ll swap out the model physics in exactly the same way.

We’ve decided how to include damage in the model, but we also need some rule for how damage will evolve. One approach would be to use fracture mechanics to explicitly model individual crevasses, but this approach would be very computationally expensive. Instead, the continuum damage mechanics approximation posits that there is some bulk transport law for damage that looks something like this PDE:

\[\frac{\partial D}{\partial t} + u\cdot\nabla D = f_o - f_c\]

where \(f_o\) is a crevasse opening term and \(f_c\) is a crevasse closing or healing term. Crevasses open when the stress is sufficiently high, and crevasses close when the strain rate is low or compressional. The real physics comes in deciding the exact details of how these terms depend on the stress, strain rate, and other variables. In the following, we’ll use the model from Albrecht and Levermann, but there are many others.

The class icepack.models.DamageTransport includes a solver for this damage transport law. We’ll initialize the damage field to be zero everywhere and we’ll advect in zero damage from upstream. When we run the simulation forward in time, the

damage_transport = icepack.models.DamageTransport()

Δ = firedrake.FunctionSpace(mesh, family='DG', degree=1)
D = firedrake.Function(Δ)
D_inflow = firedrake.Constant(0.)

Next, we’ll run a new simulation that includes interleaved updates for ice thickness and velocity, just like before, but also an update for the damage field according to the transport law we wrote down above. The simulation will go for long enough to propagate the damage field through an entire residence time of the ice shelf.

for step in range(num_timesteps + 1):
    h = ice_shelf.prognostic_solve(dt, h0=h, a=a, u=u, h_inflow=h0)
    D = damage_transport.solve(dt, D0=D, u=u, A=A, D_inflow=D_inflow)
    u = ice_shelf.diagnostic_solve(u0=u, h=h, A=A, D=D, **opts)

    print('.' if step % 10 == 0 else '', end='', flush=True)

Finally we’ll plot the result. An interesting feature you can see in the final result is that the damage field is highest in between the ice streams. These aren’t the areas where the velocity is highest, but rather where the shear strain rate is highest.

fig, axes = icepack.plot.subplots()
colors = icepack.plot.tripcolor(D, axes=axes)
<matplotlib.colorbar.Colorbar at 0x7fc4e25581d0>

If you want to experiment with this problem further, you can try:

  • changing the parameters of the damage model, such as the critical stress at which damage grows; see DamageTransport.__init__ for a full list and the default values

  • changing the inflow velocities of each ice stream; the larger the speed contrasts between neighboring streams, the greater the damage between them

  • adding an ice rise to the domain to see how shear around this obstacle affects damage.

Damage mechanics is one of many things you might want to simulate in order to accurately capture the physics of glacier flow. In later demos, we’ll show how to incorporate heat transfer. We could also include models describing things like the buoyant meltwater plume underneath the ice shelf or the densification of firn on the shelf surface. The details are different in each case, but the essential procedure is the same: first decide how the addition of a new field will affect the flow, and then decide how this field will evolve in time. We used the continuum damage formulation of Albrecht and Levermann, but there are many other models we could have chosen.