diff --git a/demos/larsen/simulate.py b/demos/larsen/simulate.py
index cf895c1..0e747e1 100644
--- a/demos/larsen/simulate.py
+++ b/demos/larsen/simulate.py
@@ -1,4 +1,5 @@
 import argparse
+import subprocess
 import numpy as np
 import geojson
 import xarray
@@ -11,7 +12,7 @@
 parser = argparse.ArgumentParser()
 parser.add_argument("--input", default="larsen-extrapolated.h5")
 parser.add_argument("--degree", type=int, default=1)
-parser.add_argument("--output")
+parser.add_argument("--output", default="larsen-simulation.h5")
 args = parser.parse_args()
 
 # Read in the starting data
@@ -132,22 +133,45 @@
 h_problem = firedrake.NonlinearVariationalProblem(G, h)
 h_solver = firedrake.NonlinearVariationalSolver(h_problem)
 
-for step in range(num_steps):
-    h_solver.solve()
-    h.interpolate(firedrake.max_value(0, h))
+# Load in the outline of Larsen in 2019
+outline_filename = icepack.datasets.fetch_outline("larsen-2019")
+with open(outline_filename, "r") as outline_file:
+    outline = geojson.load(outline_file)
+
+geometry = icepack.meshing.collection_to_geo(outline)
+with open("larsen-2019.geo", "w") as geometry_file:
+    geometry_file.write(geometry.get_code())
+
+command = "gmsh -2 -v 0 -o larsen-2019.msh larsen-2019.geo"
+subprocess.run(command.split())
+
+# Create a mask to remove ice where the calving event occurred
+mesh_2019 = firedrake.Mesh("larsen-2019.msh")
+Q_2019 = firedrake.FunctionSpace(mesh_2019, "CG", 1)
+μ = firedrake.Function(Q_2019)
+μ.assign(Constant(1.0))
+μ = firedrake.project(μ, Q)
+μ.interpolate(firedrake.max_value(0, firedrake.min_value(1, μ)))
+
+h_c = firedrake.Constant(1.0)
+with firedrake.CheckpointFile(args.output, "w") as chk:
+    chk.save_function(h, name="thickness", idx=0)
+
+    for step in range(num_steps):
+        h_solver.solve()
+        h.interpolate(firedrake.conditional(h < h_c, 0, h))
+        h_n.assign(h)
+        u_solver.solve()
+        chk.save_function(h, name="thickness", idx=step + 1)
+
+    print("IT CALVING NOW")
+    h.interpolate(μ * h)
     h_n.assign(h)
     u_solver.solve()
 
-import matplotlib.pyplot as plt
-fig, axes = plt.subplots()
-axes.set_aspect("equal")
-colors = firedrake.tripcolor(h, axes=axes)
-fig.colorbar(colors)
-plt.show()
-
-fig, axes = plt.subplots()
-axes.set_aspect("equal")
-δh = firedrake.interpolate(h - h0, Q)
-colors = firedrake.tripcolor(δh, vmin=-50, vmax=+50, cmap="RdBu_r", axes=axes)
-fig.colorbar(colors)
-plt.show()
+    for step in range(num_steps):
+        h_solver.solve()
+        h.interpolate(firedrake.conditional(h < h_c, 0, h))
+        h_n.assign(h)
+        u_solver.solve()
+        chk.save_function(h, name="thickness", idx=step + num_steps + 1)