diff --git a/demos/gibbous-ice-shelf/gibbous.py b/demos/gibbous-ice-shelf/gibbous.py
index 4e1d4ed..e275837 100644
--- a/demos/gibbous-ice-shelf/gibbous.py
+++ b/demos/gibbous-ice-shelf/gibbous.py
@@ -14,6 +14,7 @@
     NonlinearVariationalProblem,
     NonlinearVariationalSolver,
 )
+import irksome
 from icepack2 import model
 from icepack2.constants import glen_flow_law as n
 
@@ -87,10 +88,11 @@
 
 # Create some function spaces and some fields
 cg = firedrake.FiniteElement("CG", "triangle", 1)
-dg = firedrake.FiniteElement("DG", "triangle", 0)
-Q = firedrake.FunctionSpace(mesh, cg)
+dg0 = firedrake.FiniteElement("DG", "triangle", 0)
+dg1 = firedrake.FiniteElement("DG", "triangle", 1)
+Q = firedrake.FunctionSpace(mesh, dg1)
 V = firedrake.VectorFunctionSpace(mesh, cg)
-Σ = firedrake.TensorFunctionSpace(mesh, dg, symmetry=True)
+Σ = firedrake.TensorFunctionSpace(mesh, dg0, symmetry=True)
 Z = V * Σ
 
 h0 = firedrake.interpolate(h_expr, Q)
@@ -190,16 +192,26 @@
 diagnostic_solver.solve()
 
 # Set up the prognostic problem and solver
-h_n = h.copy(deepcopy=True)
-φ = firedrake.TestFunction(Q)
+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)
-flux_cells = ((h - h_n) / dt * φ - inner(h * u, grad(φ))) * dx
-ν = firedrake.FacetNormal(mesh)
-flux_in = h0 * firedrake.min_value(0, inner(u, ν)) * φ * ds
-flux_out = h * firedrake.max_value(0, inner(u, ν)) * φ * ds
-G = flux_cells + flux_in + flux_out
-prognostic_problem = NonlinearVariationalProblem(G, h)
-prognostic_solver = NonlinearVariationalSolver(prognostic_problem)
+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 a calving mask
 R = firedrake.Constant(60e3)
@@ -216,7 +228,7 @@
     chk.save_function(h, name="thickness", idx=0)
 
     for step in tqdm.trange(args.num_steps):
-        prognostic_solver.solve()
+        prognostic_solver.advance()
 
         if args.calving_freq != 0.0:
             if time_since_calving > args.calving_freq:
@@ -225,7 +237,6 @@
             time_since_calving += float(dt)
 
         h.interpolate(firedrake.max_value(0, h))
-        h_n.assign(h)
 
         diagnostic_solver.solve()