%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.4.1, 1 node, max. 1 thread]
Info    : Started on Thu Nov  7 14:21:11 2019
Info    : Reading 'ice-shelf.geo'...
Info    : Done reading 'ice-shelf.geo'
Info    : Meshing 1D...
Info    : Meshing curve 1 (Circle)
Info    : Meshing curve 2 (Circle)
Info    : Done meshing 1D (0.000516 s)
Info    : Meshing 2D...
Info    : Meshing surface 4 (Plane, Delaunay)
Info    : Done meshing 2D (0.07548 s)
Info    : 3793 vertices 7584 elements
Info    : Writing 'ice-shelf.msh'...
Info    : Done writing 'ice-shelf.msh'
Info    : Stopped on Thu Nov  7 14:21:11 2019

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

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 0x7fcc97244750>
fig, axes = icepack.plot.subplots()
contours = icepack.plot.tricontourf(u0, 40, axes=axes)
<matplotlib.colorbar.Colorbar at 0x7fcc97104ad0>

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 -18:math:{}^circC. 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.

T = 250.0
num_timesteps = 125
dt = T / 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: 35.96795261913698 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.166799173439834 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. If you want to experiment with this example a bit more, you can increase the final time \(T\) and then compare the fluxes.


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 through either plotting or more general post-processing. Here we’ve shown how to use these functions as an input to an ice shelf flow model. There are a lot of interesting experiments you can do with synthetic ice shelves. For example, you can approximate the effect of seasonality by making the accumulation rate and inflow thickness a sinusoidal function of time. You can then add things like tidal effects and see how they alter the overall ice flow. In the next demo, we’ll show how to use these functions for simulating a real ice shelf using observational data.