diff --git a/Makefile b/Makefile
index 38b66b9..fff2d3b 100644
--- a/Makefile
+++ b/Makefile
@@ -6,7 +6,7 @@ dual-problems.pdf: dual-problems.tex dual-problems.bib
 	cd demos/singularity && make primal.pdf dual.pdf && cd ../../
 	cd demos/convergence-tests && make results.pdf && cd ../../
 	cd demos/slab && make tables/table_alpha0.50.txt && cd ../../
-	cd demos/gibbous-ice-shelf && make steady-state.pdf calved-thickness.pdf volumes.pdf && cd ../../
+	cd demos/gibbous-ice-shelf && make gibbous.pdf volumes.pdf && cd ../../
 	pdflatex $<
 	bibtex $(basename $<)
 	pdflatex $<
diff --git a/demos/gibbous-ice-shelf/Makefile b/demos/gibbous-ice-shelf/Makefile
index 84f70c7..379b6e3 100644
--- a/demos/gibbous-ice-shelf/Makefile
+++ b/demos/gibbous-ice-shelf/Makefile
@@ -1,4 +1,4 @@
-steady-state.pdf calved-thickness.pdf volumes.pdf &: calving.h5 make_plots.py
+gibbous.pdf volumes.pdf residuals.pdf &: calving.h5 residuals.json make_plots.py
 	python make_plots.py --calving-freq 24.0
 
 calving.h5: steady-state-fine.h5 gibbous.py
@@ -7,5 +7,8 @@ calving.h5: steady-state-fine.h5 gibbous.py
 steady-state-fine.h5: steady-state-coarse.h5 gibbous.py
 	python gibbous.py --resolution 2e3 --final-time 400.0 --num-steps 400 --input $< --output $@
 
+residuals.json: steady-state-coarse.h5 primal_dual_calving.py
+	python primal_dual_calving.py
+
 steady-state-coarse.h5: gibbous.py
 	python gibbous.py --resolution 5e3 --final-time 400.0 --num-steps 200 --output $@
diff --git a/demos/gibbous-ice-shelf/benchmark.py b/demos/gibbous-ice-shelf/benchmark.py
new file mode 100644
index 0000000..5167070
--- /dev/null
+++ b/demos/gibbous-ice-shelf/benchmark.py
@@ -0,0 +1,157 @@
+import argparse
+import numpy as np
+import tqdm
+import firedrake
+from firedrake import (
+    derivative,
+    NonlinearVariationalProblem,
+    NonlinearVariationalSolver,
+)
+import icepack
+import irksome
+from icepack2 import model
+from icepack2.constants import glen_flow_law as n
+import gibbous_inputs
+
+parser = argparse.ArgumentParser()
+parser.add_argument("--resolution", type=float, default=5e3)
+parser.add_argument("--final-time", type=float, default=400.0)
+parser.add_argument("--num-steps", type=int, default=200)
+parser.add_argument("--form", choices=["primal", "dual"])
+args = parser.parse_args()
+
+# Generate and load the mesh
+R = 200e3
+δx = args.resolution
+
+mesh = gibbous_inputs.make_mesh(R, δx)
+
+print(mesh.num_vertices(), mesh.num_cells())
+
+h_expr, u_expr = gibbous_inputs.make_initial_data(mesh, R)
+
+# Create some function spaces and some fields
+cg = firedrake.FiniteElement("CG", "triangle", 1)
+dg0 = firedrake.FiniteElement("DG", "triangle", 0)
+dg1 = firedrake.FiniteElement("DG", "triangle", 1)
+Q = firedrake.FunctionSpace(mesh, dg1)
+V = firedrake.VectorFunctionSpace(mesh, cg)
+
+h0 = firedrake.Function(Q).interpolate(h_expr)
+u0 = firedrake.Function(V).interpolate(u_expr)
+
+h = h0.copy(deepcopy=True)
+
+# Create some physical constants. Rather than specify the fluidity factor `A`
+# in Glen's flow law, we instead define it in terms of a stress scale `τ_c` and
+# a strain rate scale `ε_c` as
+#
+#     A = ε_c / τ_c ** n
+#
+# where `n` is the Glen flow law exponent. This way, we can do a continuation-
+# type method in `n` while preserving dimensional correctness.
+ε_c = firedrake.Constant(0.01)
+τ_c = firedrake.Constant(0.1)
+
+
+if args.form == "primal":
+    u = u0.copy(deepcopy=True)
+    solver = icepack.solvers.FlowSolver(icepack.models.IceShelf(), dirichlet_ids=[1])
+    A = firedrake.Constant(ε_c / τ_c ** n)
+    solver.diagnostic_solve(velocity=u, thickness=h, fluidity=A)
+    diagnostic_solver = solver._diagnostic_solver
+elif args.form == "dual":
+    Σ = firedrake.TensorFunctionSpace(mesh, dg0, symmetry=True)
+    Z = V * Σ
+    z = firedrake.Function(Z)
+
+    # Set up the diagnostic problem
+    u, M = firedrake.split(z)
+    fields = {
+        "velocity": u,
+        "membrane_stress": M,
+        "thickness": h,
+    }
+
+    fns = [model.viscous_power, model.ice_shelf_momentum_balance]
+
+    rheology = {
+        "flow_law_exponent": n,
+        "flow_law_coefficient": ε_c / τ_c ** n,
+    }
+
+    L = sum(fn(**fields, **rheology) for fn in fns)
+    F = derivative(L, z)
+
+    # Make an initial guess by solving a Picard linearization
+    linear_rheology = {
+        "flow_law_exponent": 1,
+        "flow_law_coefficient": ε_c / τ_c,
+    }
+
+    L_1 = sum(fn(**fields, **linear_rheology) for fn in fns)
+    F_1 = derivative(L_1, z)
+
+    qdegree = n + 2
+    bc = firedrake.DirichletBC(Z.sub(0), u0, (1,))
+    problem_params = {
+        #"form_compiler_parameters": {"quadrature_degree": qdegree},
+        "bcs": bc,
+    }
+    solver_params = {
+        "solver_parameters": {
+            "snes_type": "newtonls",
+            "ksp_type": "gmres",
+            "pc_type": "lu",
+            "pc_factor_mat_solver_type": "mumps",
+            "snes_rtol": 1e-2,
+        },
+    }
+    firedrake.solve(F_1 == 0, z, **problem_params, **solver_params)
+
+    # Create an approximate linearization
+    h_min = firedrake.Constant(1.0)
+    rfields = {
+        "velocity": u,
+        "membrane_stress": M,
+        "thickness": firedrake.max_value(h_min, h),
+    }
+
+    L_r = sum(fn(**rfields, **rheology) for fn in fns)
+    F_r = derivative(L_r, z)
+    J_r = derivative(F_r, z)
+
+    # Set up the diagnostic problem and solver
+    diagnostic_problem = NonlinearVariationalProblem(F, z, J=J_r, **problem_params)
+    diagnostic_solver = NonlinearVariationalSolver(diagnostic_problem, **solver_params)
+    diagnostic_solver.solve()
+else:
+    raise ValueError(f"--form must be either `primal` or `dual`, got {args.form}")
+
+# Set up the prognostic problem and solver
+prognostic_problem = model.mass_balance(
+    thickness=h,
+    velocity=u,
+    accumulation=firedrake.Constant(0.0),
+    thickness_inflow=h0,
+    test_function=firedrake.TestFunction(Q),
+)
+dt = firedrake.Constant(args.final_time / args.num_steps)
+t = firedrake.Constant(0.0)
+method = irksome.BackwardEuler()
+prognostic_params = {
+    "solver_parameters": {
+        "snes_type": "ksponly",
+        "ksp_type": "gmres",
+        "pc_type": "bjacobi",
+    },
+}
+prognostic_solver = irksome.TimeStepper(
+    prognostic_problem, method, t, dt, h, **prognostic_params
+)
+
+# Run the simulation
+for step in tqdm.trange(args.num_steps):
+    prognostic_solver.advance()
+    h.interpolate(firedrake.max_value(0, h))
+    diagnostic_solver.solve()
diff --git a/demos/gibbous-ice-shelf/gibbous.py b/demos/gibbous-ice-shelf/gibbous.py
index 32ec8be..b9a59b5 100644
--- a/demos/gibbous-ice-shelf/gibbous.py
+++ b/demos/gibbous-ice-shelf/gibbous.py
@@ -1,9 +1,6 @@
 import argparse
-import subprocess
 import numpy as np
-from numpy import pi as π
 import tqdm
-import pygmsh
 import firedrake
 from firedrake import (
     inner,
@@ -17,6 +14,7 @@
 import irksome
 from icepack2 import model
 from icepack2.constants import glen_flow_law as n
+import gibbous_inputs
 
 parser = argparse.ArgumentParser()
 parser.add_argument("--input")
@@ -31,60 +29,8 @@
 R = 200e3
 δx = args.resolution
 
-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),
-]
-
-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)
-
-with open("ice-shelf.geo", "w") as geo_file:
-    geo_file.write(geometry.get_code())
-
-command = "gmsh -2 -format msh2 -v 0 -o ice-shelf.msh ice-shelf.geo"
-subprocess.run(command.split())
-
-mesh = firedrake.Mesh("ice-shelf.msh")
-
-# Generate the initial data
-inlet_angles = π * np.array([-3 / 4, -1 / 2, -1 / 3, -1 / 6])
-inlet_widths = π * np.array([1 / 8, 1 / 12, 1 / 24, 1 / 12])
-
-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 * firedrake.as_vector((np.cos(θ), np.sin(θ)))
-    v = -firedrake.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)
-
-h_expr = firedrake.Constant(hb)
-for h in hs:
-    h_expr = firedrake.max_value(h, h_expr)
-
-u_expr = sum(us)
+mesh = gibbous_inputs.make_mesh(R, δx)
+h_expr, u_expr = gibbous_inputs.make_initial_data(mesh, R)
 
 # Create some function spaces and some fields
 cg = firedrake.FiniteElement("CG", "triangle", 1)
@@ -95,8 +41,8 @@
 Σ = firedrake.TensorFunctionSpace(mesh, dg0, symmetry=True)
 Z = V * Σ
 
-h0 = firedrake.interpolate(h_expr, Q)
-u0 = firedrake.interpolate(u_expr, V)
+h0 = firedrake.Function(Q).interpolate(h_expr)
+u0 = firedrake.Function(V).interpolate(u_expr)
 
 h = h0.copy(deepcopy=True)
 z = firedrake.Function(Z)
@@ -218,6 +164,7 @@
 
 # Set up a calving mask
 R = firedrake.Constant(60e3)
+x = firedrake.SpatialCoordinate(mesh)
 y = firedrake.Constant((0.0, R))
 mask = firedrake.conditional(inner(x - y, x - y) < R**2, 0.0, 1.0)
 
diff --git a/demos/gibbous-ice-shelf/gibbous_inputs.py b/demos/gibbous-ice-shelf/gibbous_inputs.py
new file mode 100644
index 0000000..54a318e
--- /dev/null
+++ b/demos/gibbous-ice-shelf/gibbous_inputs.py
@@ -0,0 +1,64 @@
+import subprocess
+import numpy as np
+from numpy import pi as π
+import pygmsh
+import firedrake
+from firedrake import inner
+
+
+def make_mesh(R, δx):
+    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),
+    ]
+
+    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)
+
+    with open("ice-shelf.geo", "w") as geo_file:
+        geo_file.write(geometry.get_code())
+
+    command = "gmsh -2 -format msh2 -v 0 -o ice-shelf.msh ice-shelf.geo"
+    subprocess.run(command.split())
+
+    return firedrake.Mesh("ice-shelf.msh")
+
+
+def make_initial_data(mesh, R):
+    inlet_angles = π * np.array([-3 / 4, -1 / 2, -1 / 3, -1 / 6])
+    inlet_widths = π * np.array([1 / 8, 1 / 12, 1 / 24, 1 / 12])
+
+    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 * firedrake.as_vector((np.cos(θ), np.sin(θ)))
+        v = -firedrake.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)
+
+    h_expr = firedrake.Constant(hb)
+    for h in hs:
+        h_expr = firedrake.max_value(h, h_expr)
+
+    u_expr = sum(us)
+    return h_expr, u_expr
diff --git a/demos/gibbous-ice-shelf/make_plots.py b/demos/gibbous-ice-shelf/make_plots.py
index ab90692..54375ab 100644
--- a/demos/gibbous-ice-shelf/make_plots.py
+++ b/demos/gibbous-ice-shelf/make_plots.py
@@ -1,4 +1,5 @@
 import argparse
+import json
 import numpy as np
 import matplotlib.pyplot as plt
 from matplotlib_scalebar.scalebar import ScaleBar
@@ -17,10 +18,10 @@
 with firedrake.CheckpointFile("calving.h5", "r") as chk:
     timesteps_3 = chk.h5pyfile["timesteps"][:]
     mesh = chk.load_mesh()
-    hs_3 = [
-        chk.load_function(mesh, "thickness", idx)
-        for idx in range(len(timesteps_3))
-    ]
+    num_steps = len(timesteps_3)
+    hs = [chk.load_function(mesh, "thickness", idx) for idx in range(num_steps)]
+    us = [chk.load_function(mesh, "velocity", idx) for idx in range(num_steps)]
+    Ms = [chk.load_function(mesh, "membrane_stress", idx) for idx in range(num_steps)]
 
 # Get the results from the initial coarse spin-up and project them into the
 # function space we used for the calving experiment
@@ -46,50 +47,63 @@
     u_steady = chk.load_function(mesh_2, "velocity", len(timesteps_2) - 1)
 
 # Make a plot showing the steady state at the end of the final spin-up
-fig, axes = plt.subplots(nrows=1, ncols=3, sharex=True, sharey=True)
-for ax in axes:
+fig, axes = plt.subplots(nrows=2, ncols=3, sharex=True, sharey=True, layout="compressed")
+for ax in axes.flatten():
     ax.set_aspect("equal")
     ax.get_xaxis().set_visible(False)
     ax.get_yaxis().set_visible(False)
 
-axes[0].set_title("Thickness")
-colors = tripcolor(h_steady, axes=axes[0])
-fig.colorbar(colors, label="meters", orientation="horizontal", pad=0.04, ax=axes[0])
+axes[0, 0].set_title("Thickness")
+kw = {"xy": (0.02, 0.08), "xycoords": "axes fraction"}
+axes[0, 0].annotate("a)", **kw)
+colors = tripcolor(h_steady, vmin=0.0, axes=axes[0, 0])
+fig.colorbar(colors, label="meters", orientation="horizontal", pad=0.04, ax=axes[0, 0])
 
-axes[1].set_title("Velocity")
-colors = streamplot(u_steady, resolution=10e3, axes=axes[1])
-fig.colorbar(colors, label="meters/year", orientation="horizontal", pad=0.04, ax=axes[1])
+axes[0, 1].set_title("Velocity")
+axes[0, 1].annotate("b)", **kw)
+colors = streamplot(u_steady, resolution=10e3, seed=1729, axes=axes[0, 1])
+fig.colorbar(colors, label="meters/year", orientation="horizontal", pad=0.04, ax=axes[0, 1])
 
-axes[2].set_title("Membrane stress")
+axes[0, 2].set_title("Membrane stress")
+axes[0, 2].annotate("c)", **kw)
 elt = firedrake.FiniteElement("DG", "triangle", M_steady.ufl_element().degree())
-S = firedrake.FunctionSpace(M_steady.ufl_domain(), elt)
-m = firedrake.Function(S).interpolate(1e3 * sqrt(inner(M_steady, M_steady)))
-colors = tripcolor(m, axes=axes[2])
-fig.colorbar(colors, label="kPa", orientation="horizontal", pad=0.04, ax=axes[2])
+S = firedrake.FunctionSpace(mesh, elt)
+m = firedrake.Function(S).interpolate(1e3 * sqrt(inner(Ms[0], Ms[0])))
+colors = tripcolor(m, axes=axes[0, 2])
+fig.colorbar(colors, label="kPa", orientation="horizontal", pad=0.04, ax=axes[0, 2])
 
-scalebar = ScaleBar(1, units="m", length_fraction=0.25, location="lower right")
-axes[0].add_artist(scalebar)
+axes[1, 0].set_title("Calved thickness")
+axes[1, 0].annotate("d)", **kw)
+index = (timesteps_3 <= args.calving_freq).argmin() + 1
+colors = tripcolor(hs[index], vmin=0.0, axes=axes[1, 0])
 
-fig.savefig("steady-state.pdf", bbox_inches="tight")
+axes[1, 1].set_title("Speed change")
+axes[1, 1].annotate("e)", **kw)
+V = us[0].function_space()
+δu = firedrake.Function(V).interpolate(us[index] - us[0])
+colors = firedrake.tripcolor(δu, axes=axes[1, 1])
+fig.colorbar(colors, label="meters/year", orientation="horizontal", pad=0.04, ax=axes[1, 1])
+
+axes[1, 2].set_title("Stress change")
+axes[1, 2].annotate("f)", **kw)
+ΔM = Ms[index] - Ms[0]
+δM = firedrake.Function(S).project(1e3 * sqrt(inner(ΔM, ΔM)))
+colors = firedrake.tripcolor(δM, axes=axes[1, 2])
+fig.colorbar(colors, label="kPa", orientation="horizontal", pad=0.04, ax=axes[1, 2])
+
+scalebar = ScaleBar(1, units="m", length_fraction=0.25, location="lower right")
+axes[0, 0].add_artist(scalebar)
+fig.get_layout_engine().set(h_pad=8/72)
+fig.savefig("gibbous.pdf", bbox_inches="tight")
 
-# Make a plot showing the thickness immediately after calving
-fig, axes = plt.subplots()
-axes.set_aspect("equal")
-axes.get_xaxis().set_visible(False)
-axes.set_ylabel("northing (meters)")
-axes.ticklabel_format(style="sci", axis="both", scilimits=(0, 0))
-index = (timesteps_3 <= args.calving_freq).argmin() + 1
-colors = tripcolor(hs_3[index], vmin=0.0, axes=axes)
-fig.colorbar(colors, label="meters", orientation="horizontal", pad=0.04)
-fig.savefig("calved-thickness.pdf", bbox_inches="tight")
 
 # Make a plot showing the volumes through the entire experiment
 t1 = timesteps_1[-1]
 t2 = timesteps_1[-1] + timesteps_2[-1]
 t3 = t2 + timesteps_3[-1]
 timesteps = np.concatenate((timesteps_1, timesteps_2 + t1, timesteps_3 + t2))
-hs = hs_1 + hs_2 + hs_3
-volumes = np.array([firedrake.assemble(h * dx) / 1e9 for h in hs])
+hs_all = hs_1 + hs_2 + hs
+volumes = np.array([firedrake.assemble(h * dx) / 1e9 for h in hs_all])
 
 fig, ax = plt.subplots(figsize=(6, 3.2))
 ax.set_xlabel("Time (years)")
@@ -109,3 +123,15 @@
 ax.annotate("fine\nspin-up", xy=(t1, tlevel), xytext=((t1 + t2) / 2, tlevel), **kwargs)
 ax.annotate("calving\ncycle", xy=(t2, tlevel), xytext=((t2 + t3) / 2, tlevel), **kwargs)
 fig.savefig("volumes.pdf", bbox_inches="tight")
+
+
+# Make a plot showing the number of Newton iterations required during the
+# calving phase of the experiment
+with open("residuals.json", "r") as residuals_file:
+    residuals = json.load(residuals_file)
+
+fig, ax = plt.subplots(figsize=(6.4, 3.2))
+ax.bar(list(range(len(residuals))), [len(r) for r in residuals])
+ax.set_xlabel("Timestep (years)")
+ax.set_ylabel("Iterations")
+fig.savefig("residuals.pdf", bbox_inches="tight")
diff --git a/demos/gibbous-ice-shelf/primal_dual_calving.py b/demos/gibbous-ice-shelf/primal_dual_calving.py
new file mode 100644
index 0000000..25d9844
--- /dev/null
+++ b/demos/gibbous-ice-shelf/primal_dual_calving.py
@@ -0,0 +1,118 @@
+import json
+import numpy as np
+from numpy import pi as π
+import tqdm
+import firedrake
+from firedrake import inner, derivative
+import irksome
+from icepack2 import model, solvers
+from icepack2.constants import glen_flow_law as n
+
+# Load the input data
+with firedrake.CheckpointFile("steady-state-coarse.h5", "r") as chk:
+    mesh = chk.load_mesh()
+    idx = len(chk.h5pyfile["timesteps"]) - 1
+    h = chk.load_function(mesh, "thickness", idx=idx)
+    h0 = h.copy(deepcopy=True)
+    u = chk.load_function(mesh, "velocity", idx=idx)
+    M = chk.load_function(mesh, "membrane_stress", idx=idx)
+
+Q = h.function_space()
+V = u.function_space()
+Σ = M.function_space()
+
+Z = V * Σ
+z = firedrake.Function(Z)
+z.sub(0).assign(u)
+z.sub(1).assign(M);
+
+# Set up the Lagrangian
+ε_c = firedrake.Constant(0.01)
+τ_c = firedrake.Constant(0.1)
+
+u, M = firedrake.split(z)
+fields = {
+    "velocity": u,
+    "membrane_stress": M,
+    "thickness": h,
+}
+
+fns = [model.viscous_power, model.ice_shelf_momentum_balance]
+
+rheology = {
+    "flow_law_exponent": n,
+    "flow_law_coefficient": ε_c / τ_c ** n,
+}
+
+L = sum(fn(**fields, **rheology) for fn in fns)
+F = derivative(L, z)
+
+# Set up a regularized Lagrangian
+h_min = firedrake.Constant(0.1)
+rfields = {
+    "velocity": u,
+    "membrane_stress": M,
+    "thickness": firedrake.max_value(h_min, h),
+}
+
+L_r = sum(fn(**rfields, **rheology) for fn in fns)
+F_r = derivative(L_r, z)
+H_r = derivative(F_r, z)
+
+# Set up the mass balance equation
+prognostic_problem = model.mass_balance(
+    thickness=h,
+    velocity=u,
+    accumulation=firedrake.Constant(0.0),
+    thickness_inflow=h0,
+    test_function=firedrake.TestFunction(Q),
+)
+
+final_time = 400.0
+num_steps = 400
+
+dt = firedrake.Constant(final_time / num_steps)
+t = firedrake.Constant(0.0)
+method = irksome.BackwardEuler()
+prognostic_params = {
+    "solver_parameters": {
+        "snes_type": "ksponly",
+        "ksp_type": "gmres",
+        "pc_type": "bjacobi",
+    },
+}
+prognostic_solver = irksome.TimeStepper(
+    prognostic_problem, method, t, dt, h, **prognostic_params
+)
+
+# Set up the diagnostic solver and boundary conditions and do an initial solve
+bc = firedrake.DirichletBC(Z.sub(0), firedrake.Constant((0, 0)), (1,))
+problem = solvers.ConstrainedOptimizationProblem(L, z, H=H_r, bcs=bc)
+diagnostic_solver = solvers.NewtonSolver(problem, tolerance=1e-4)
+
+residuals = [list(diagnostic_solver.solve())]
+
+# Create the calving mask -- this describes how we'll remove ice
+radius = firedrake.Constant(60e3)
+x = firedrake.SpatialCoordinate(mesh)
+y = firedrake.Constant((0.0, radius))
+mask = firedrake.conditional(inner(x - y, x - y) < radius**2, 0.0, 1.0)
+
+# Run the simulation
+calving_frequency = 24.0
+time_since_calving = 0.0
+
+for step in tqdm.trange(num_steps):
+    prognostic_solver.advance()
+
+    if time_since_calving > calving_frequency:
+        h.interpolate(mask * h)
+        time_since_calving = 0.0
+    time_since_calving += float(dt)
+    h.interpolate(firedrake.max_value(0, h))
+
+    residuals.append(list(diagnostic_solver.solve()))
+
+# Save the results to disk
+with open("residuals.json", "w") as residuals_file:
+    json.dump(residuals, residuals_file)
diff --git a/dual-problems.tex b/dual-problems.tex
index 324f7bc..a49c030 100644
--- a/dual-problems.tex
+++ b/dual-problems.tex
@@ -242,15 +242,16 @@ \subsection{Primal action principles}
 In other words, being the Euler-Lagrange equations for some action functional is a special property of only a restricted class of PDEs.
 There is no rote procedure to determine what the action functional might be, but in attempting to construct one, the main mathematical hurdle to overcome is computing the anti-derivatives of certain terms in the differential equation.
 
-The constitutive relation \eqref{eq:constitutive-relation} for $M$ can be expressed as the derivative of a certain scalar quantity:
+The constitutive relation \eqref{eq:constitutive-relation} for $M$ can be expressed as the derivative of a certain scalar quantity.
+In index notation,
 \begin{equation}
-    M = \frac{d}{d\dot\varepsilon}\left(\frac{2n}{n + 1}A^{-\frac{1}{n}}|\dot\varepsilon|_{\mathscr C}^{\frac{1}{n} + 1}\right)
+    M_{ij} = \frac{\partial}{\partial\dot\varepsilon_{ij}}\left(\frac{2n}{n + 1}A^{-\frac{1}{n}}|\dot\varepsilon|_{\mathscr C}^{\frac{1}{n} + 1}\right)
     \label{eq:M-anti-derivative}
 \end{equation}
 if we think of the strain rate tensor as an independent variable and briefly forget that it is the symmetrized velocity gradient.
 Likewise, for the sliding relation \eqref{eq:sliding-law}, we can observe that
 \begin{equation}
-    \tau = -\frac{d}{du}\left(\frac{m}{m + 1}C|u|^{\frac{1}{m} + 1}\right).
+    \tau_i = -\frac{\partial}{\partial u_i}\left(\frac{m}{m + 1}C|u|^{\frac{1}{m} + 1}\right).
     \label{eq:tau-anti-derivative}
 \end{equation}
 Using equations \eqref{eq:M-anti-derivative} and \eqref{eq:tau-anti-derivative}, one can show that the SSA momentum balance are the Euler-Lagrange equations for the following action functional \citep{dukowicz2010consistent}:
@@ -428,7 +429,11 @@ \subsection{Comparison against primal form on gibbous ice shelf}
 We first run a spin-up of the system to steady state using the coupled mass and momentum balance equations with both the primal and dual forms.
 We check that the velocity obtained by solving the dual form is within discretization error of the velocity obtained with the primal form, which offers an additional degree of verification that we are solving the equations right.
 We additionally evaluate the total wall clock time to run this experiment using both the primal and dual form.
-We expect that the dual form is slower because it has more unknowns: 3 additional degrees of freedom per triangle for the independent components of the membrane stress tensor.
+The primal problem using $CG(1)$ elements for the velocity has two unknowns for each vertex of the mesh.
+The dual problem using $CG(1) \times DG(0)$ elements for the velocity and stress has an additional three unknowns per triangle.
+The Euler formula (\#vertices - \#edges + \#triangles $\approx$ 2) implies that there are approximately twice as many triangles as there are vertices.
+Consequently there are about 4x as many degrees of freedom when solving the dual problem as there are for the primal problem.
+Assuming naively that the time to solution scales linearly with the number of unknowns, we would then expect that solving the dual problem is 4x as expensive as solving the primal problem.
 
 As a third and final phase of this experiment, we run the same simulation, but every 24 years we set the ice thickness to 0 in a prescribed region near the terminus.
 This forcing mimics the effect of a large iceberg calving event.
@@ -568,9 +573,10 @@ \subsection{Gibbous ice shelf} \label{sec:gibbous-ice-shelf}
 
 \begin{figure}[h]
     \begin{center}
-        \includegraphics[width=0.85\linewidth]{demos/gibbous-ice-shelf/steady-state.pdf}
+        \includegraphics[width=0.85\linewidth]{demos/gibbous-ice-shelf/gibbous.pdf}
     \end{center}
-    \caption{The thickness, velocity, and magnitude of the membrane stress tensor in steady state.}
+    \caption{The thickness (a), velocity (b), and magnitude of the membrane stress tensor (c) in steady state, and the thickness (d), magnitude of the velocity change (e), and magnitude of the stress change (f) immediately after the calving event.
+    We remove a semi-circular segment from the end of the shelf with a prescribed center and radius.}
     \label{fig:gibbous}
 \end{figure}
 
@@ -584,35 +590,37 @@ \subsection{Gibbous ice shelf} \label{sec:gibbous-ice-shelf}
     \label{fig:gibbous-calving-volumes}
 \end{figure}
 
-\begin{figure}[h]
-    \begin{center}
-        \includegraphics[width=0.85\linewidth]{demos/gibbous-ice-shelf/calved-thickness.pdf}
-    \end{center}
-    \caption{The thickness of the ice shelf immediately after the calving event.
-    We remove a semi-circular segment from the end of the shelf with a prescribed center and radius.
-    The nonlinear solver for the velocity and membrane stress still converges with this thickness input.}
-    \label{fig:gibbous-calving-thickness}
-\end{figure}
-
 For the spin-up phase of the experiment, we did an initial run for 400 years on a mesh with a 5km resolution, at which point the system is close to steady state.
 We then projected these fields to a finer mesh with a resolution of 2km and use them as the initial state for a further 400 years of spin-up.
-The results are shown in figure \ref{fig:gibbous} and are identical to those obtained from the primal form of the problem up to discretization error.
-\textcolor{red}{Talk about the benchmark.}
+The results are shown in figure \ref{fig:gibbous}(a)-(c) and are identical to those obtained from the primal form of the problem up to discretization error.
+
+When we used the spin-up phase of the experiment as a benchmark to measure the performance of the dual and primal solvers, we found that the dual problem required between 2.5x and 2.7x as much time.
+These results were consistent across different mesh resolutions and when run several times on multiple machines.
+Since the dual problem has 4x as many unknowns, the added cost that we found experimentally is less than what we would expect if we naively assumed that cost is proportional to the number of degrees of freedom.
 
-In the calving phase of the experiment, our solver for the dual problem still performed gracefully in ice-free areas.
+In the calving phase of the experiment, our solver for the dual problem still worked in ice-free areas.
 This feature offers the possibility of implementing physically-based calving models in a simple way.
 Figure \ref{fig:gibbous-calving-volumes} shows the evolution of the volume of ice in the shelf over the two spin-up phases and the calving phase.
 The 24-year recurrence interval is not enough time for the ice to advance back to the original edge of the computational domain.
 Using a longer interval would allow the calving terminus to advance back to its original position.
-Figure \ref{fig:gibbous-calving-thickness} shows the thickness of the ice shelf immediately after the calving event.
+Figure \ref{fig:gibbous}(d) shows the thickness of the ice shelf immediately after the calving event (e) and (f) show the magnitude of the change in speed and stress respectively.
+In particular, the stress field shows a discontinuity at the new calving front as expected.
+
+\begin{figure}[h]
+    \begin{center}
+        \includegraphics[width=0.75\linewidth]{demos/gibbous-ice-shelf/residuals.pdf}
+    \end{center}
+    \caption{The number of Newton iterations to compute the ice velocity at each step of the calving phase of the experiment.
+    Calving occurs every 24 years.}
+    \label{fig:gibbous-residuals}
+\end{figure}
 
-The number of nonlinear solver iterations to recompute the ice velocity after the calving event is much greater than after a normal timestep.
+Figure \ref{fig:gibbous-residuals} shows the number of Newton iterations needed to compute the ice velocity at each timestep.
+The number of iterations after a calving event is much greater than after a normal timestep.
 This behavior is to be expected if we run the solver to convergence because it introduces a type of shock into the system.
 As the system relaxes back, the number of iterations decreases again.
-Moreover, we had to do some manual adjustment of the convergence tolerances.
 Several different strategies can alleviate the need for manual adjustment (see the appendix).
-\textcolor{red}{Give a plot showing the number of Newton solves.
-Do a comparison to the primal form with e.g. a small minimum thickness.}
+\textcolor{red}{Do a comparison to the primal form with e.g. a small minimum thickness.}
 
 \subsection{Larsen C Ice Shelf}
 
@@ -709,6 +717,9 @@ \section{Discussion}
 The number of unknowns in the dual formulation is greater than in the primal form, thus putting more pressure on computer memory.
 The resulting linear systems are indefinite rather than positive-definite.
 Finally, since the dual form is a mixed problem, it is possible to make bad choices of finite element basis, whereas almost any basis will work for the primal form.
+We did find, however, that the increased cost of solving the dual problem was not as high as one might expect just based on counting the number of degrees of freedom.
+We used a fairly naive solution approach (direct factorization) for the linear system in each step of Newton's method in the benchmark for both the primal and dual forms.
+There may be significant room for improvement on these benchmarks through the use of more sophisticated techniques such as Schur complement preconditioners that use static condensation of the stress degrees of freedom \citep{boffi2013mixed}.
 
 There are several promising avenues of future work on this problem.
 Including the stress tensor as an unknown and the constitutive relation as an equation to be solved opens up several possibilities for modifying the physics.