diff --git a/docs/Project.toml b/docs/Project.toml
index b29fb873..1dda9289 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -18,6 +18,7 @@ JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
+NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
diff --git a/docs/make.jl b/docs/make.jl
index e8d7e32b..6762c18c 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -44,7 +44,7 @@ makedocs(
   pagesonly = true,
   linkcheck = true,
   linkcheck_ignore = [r"agupubs\.onlinelibrary\.wiley\.com", # This gives a 403 Forbidden
-                      r"Decapodes\.jl/dev"], # 404, probably due to bad self-rerference
+                      r"Decapodes\.jl/dev"], # 404, probably due to bad self-reference
   pages     = Any[
     "Decapodes.jl" => "index.md",
     "Overview" => "overview/overview.md",
@@ -57,6 +57,7 @@ makedocs(
     "Glacial Flow" => "ice_dynamics/ice_dynamics.md",
     "Grigoriev Ice Cap" => "grigoriev/grigoriev.md", # Requires ice_dynamics
     "Budyko-Sellers-Halfar" => "bsh/budyko_sellers_halfar.md", # Requires ice_dynamics
+    "Halfar-EBM-Water" => "ebm_melt/ebm_melt.md",
     "Halfar-NS" => "halmo/halmo.md", # Requires grigoriev
     "NHS" => "nhs/nhs_lite.md",
     "Pipe Flow" => "poiseuille/poiseuille.md",
diff --git a/docs/src/ebm_melt/ebm_melt.md b/docs/src/ebm_melt/ebm_melt.md
new file mode 100644
index 00000000..520a582f
--- /dev/null
+++ b/docs/src/ebm_melt/ebm_melt.md
@@ -0,0 +1,334 @@
+# Haflar-EBM-Water
+
+This docs page demonsrates a composition of the Halfar model of ice dynamics with the Budyko-Sellers energy-balance model (EBM) defining temperature dynamics. Surface temperature affects the rate at which ice diffuses, and melts ice into a water density term. This water density then diffuses across the domain.
+
+We execute these dynamics on real sea ice thickness data provided by the Polar Science Center.
+
+``` @example DEC
+# Import Dependencies 
+
+# AlgebraicJulia Dependencies
+using Catlab
+using Catlab.Graphics
+using CombinatorialSpaces
+using DiagrammaticEquations
+using Decapodes
+
+# External Dependencies
+using ComponentArrays
+using CoordRefSystems
+using CairoMakie
+using LinearAlgebra
+using MLStyle
+using NetCDF
+using OrdinaryDiffEq
+Point3D = Point3{Float64}
+```
+
+First, we load a mesh. We will execute these dynamics on the sphere:
+
+``` @example DEC
+# Load a mesh
+s_plots = loadmesh(Icosphere(7));
+s = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s_plots);
+subdivide_duals!(s, Barycenter());
+wireframe(s_plots)
+```
+
+## Load data
+
+The data provided by the Polar Science Center is given as a NetCDF file. Ice thickness is a matrix with the same dimensions as a matrix provided Latitude and Longitude of the associated point on the Earth's surface. We need to convert between polar and Cartesian coordinates to use this data on our mesh.
+
+``` @example DEC
+ice_thickness_file = "piomas20c.heff.1901.2010.v1.0.nc"
+run(`curl -o $ice_thickness_file https://pscfiles.apl.uw.edu/axel/piomas20c/v1.0/monthly/piomas20c.heff.1901.2010.v1.0.nc`)
+
+# Use ncinfo(ice_thickness_file) to interactively get information on variables.
+# Sea ice thickness ("sit") has dimensions of [y, x, time].
+# y,x index into "Latitude" and "Longitude" variables.
+# Time is in units of days since 1901-01-01.
+lat = ncread(ice_thickness_file, "Latitude")
+lon = ncread(ice_thickness_file, "Longitude")
+sit = ncread(ice_thickness_file, "sit")
+
+# Convert latitude from [90, -90] to [0, 180] for convenience.
+lat .= -lat .+ 90
+
+# Convert mesh points from Cartesian to spherical coordinates.
+p_sph = map(point(s)) do p
+  p = convert(Spherical, Cartesian(p...))
+  [rad2deg(p.θ).val, rad2deg(p.ϕ).val]
+end
+
+# Note: You can instead use an algebraic parameterization, rather than nearest-neighbor interpolation.
+# Note: You can alternatively set a value to 0.0 if the distance to the nearest-neighbor is greater than some threshold.
+sit_sph_idxs = map(p_sph) do p
+  argmin(map(i -> sqrt((lat[i] - p[1])^2 + (lon[i] - p[2])^2), eachindex(lat)))
+end
+
+sit_sph = map(sit_sph_idxs, p_sph) do i, p
+  ((p[1] > maximum(lat)) || isnan(sit[i])) ? 0.0f0 : sit[i]
+end
+
+f = Figure()
+ax = LScene(f[1,1], scenekw=(lights=[],))
+msh = mesh!(ax, s_plots, color=sit_sph)
+Colorbar(f[1,2], msh)
+f
+```
+
+## Define the model
+
+Here, the Halfar model of ice dynamics is recalled, as well as the Budyko-Sellers EBM. These these models are composed individually. They are then coupled together via `warming` and `melting` components.
+
+``` @example DEC
+halfar_eq2 = @decapode begin
+  (h, melt)::Form0
+  Γ::Form1
+  n::Constant
+
+  ∂ₜ(h)  == ∘(⋆, d, ⋆)(Γ * d(h) * avg₀₁(mag(♯(d(h)))^(n-1)) * avg₀₁(h^(n+2))) - melt
+end
+
+glens_law = @decapode begin
+  (A,Γ)::Form1
+  (ρ,g,n)::Constant
+  
+  Γ == (2/(n+2))*A*(ρ*g)^n
+end
+
+ice_dynamics_composition_diagram = @relation () begin
+  dynamics(Γ,n)
+  stress(Γ,n)
+end
+
+ice_dynamics = apex(oapply(ice_dynamics_composition_diagram,
+  [Open(halfar_eq2, [:Γ,:n]),
+   Open(glens_law, [:Γ,:n])]))
+
+energy_balance = @decapode begin
+  (Tₛ, ASR, OLR, HT)::Form0
+  C::Constant
+
+  ∂ₜ(Tₛ) == (ASR - OLR + HT) ./ C
+end
+
+absorbed_shortwave_radiation = @decapode begin
+  (Q, ASR)::Form0
+  α::Constant
+
+  ASR == (1 .- α) .* Q
+end
+
+outgoing_longwave_radiation = @decapode begin
+  (Tₛ, OLR)::Form0
+  (A,B)::Constant
+
+  OLR == A .+ (B .* Tₛ)
+end
+
+heat_transfer = @decapode begin
+  (HT, Tₛ)::Form0
+  (D,cosϕᵖ,cosϕᵈ)::Constant
+
+  HT == (D ./ cosϕᵖ) .* ⋆(d(cosϕᵈ .* ⋆(d(Tₛ))))
+end
+
+insolation = @decapode begin
+  Q::Form0
+  cosϕᵖ::Constant
+
+  Q == 450 * cosϕᵖ
+end
+
+budyko_sellers_composition_diagram = @relation () begin
+  energy(Tₛ, ASR, OLR, HT)
+  absorbed_radiation(Q, ASR)
+  outgoing_radiation(Tₛ, OLR)
+  diffusion(Tₛ, HT, cosϕᵖ)
+  insolation(Q, cosϕᵖ)
+end
+
+budyko_sellers = apex(oapply(budyko_sellers_composition_diagram,
+  [Open(energy_balance, [:Tₛ, :ASR, :OLR, :HT]),
+   Open(absorbed_shortwave_radiation, [:Q, :ASR]),
+   Open(outgoing_longwave_radiation, [:Tₛ, :OLR]),
+   Open(heat_transfer, [:Tₛ, :HT, :cosϕᵖ]),
+   Open(insolation, [:Q, :cosϕᵖ])]))
+nothing # hide
+```
+
+``` @example DEC
+warming = @decapode begin
+  Tₛ::Form0
+  A::Form1
+
+  A == avg₀₁(5.8282*10^(-0.236 * Tₛ)*1.01e-19)
+end
+
+melting = @decapode begin
+  (Tₛ, h, melt, water)::Form0
+  Dₕ₂ₒ::Constant
+
+  melt == (Tₛ - 15)*1e-16*h
+  ∂ₜ(water) == melt + Dₕ₂ₒ*Δ(water)
+end
+
+budyko_sellers_halfar_water_composition_diagram = @relation () begin
+  budyko_sellers(Tₛ)
+
+  warming(A, Tₛ)
+
+  melting(Tₛ, h, melt)
+
+  halfar(A, h, melt)
+end
+```
+
+``` @example DEC
+budyko_sellers_halfar_water = apex(oapply(budyko_sellers_halfar_water_composition_diagram,
+  [Open(budyko_sellers, [:Tₛ]),
+   Open(warming, [:A, :Tₛ]),
+   Open(melting, [:Tₛ, :h, :melt]),
+   Open(ice_dynamics, [:stress_A, :dynamics_h, :dynamics_melt])]))
+```
+
+## Define initial conditions
+
+The initial data must be specified for state variables, as well as constants and parameters.
+
+``` @example DEC
+# This is a primal 0-form, with values at vertices.
+cosϕᵖ = map(x -> cos(x[1]), point(s))
+# This is a dual 0-form, with values at edge centers.
+cosϕᵈ = map(edges(s)) do e
+  (cos(point(s, src(s, e))[1]) + cos(point(s, tgt(s, e))[1])) / 2
+end
+
+α₀ = 0.354
+α₂ = 0.25
+α = map(point(s)) do ϕ
+  α₀ + α₂*((1/2)*(3*ϕ[1]^2 - 1))
+end
+A = 210
+B = 2
+f = 0.70
+ρ = 1025
+cw = 4186
+H = 70
+C = map(point(s)) do ϕ
+  f * ρ * cw * H
+end
+D = 0.6
+
+# Isothermal initial conditions:
+Tₛ₀ = map(point(s)) do ϕ
+  15.0
+end
+
+water = map(point(s)) do _
+  0.0
+end
+
+Dₕ₂ₒ = 1e-16
+
+n = 3
+halfar_ρ = 910
+g = 9.8
+
+h₀ = sit_sph
+# Store these values to be passed to the solver.
+u₀ = ComponentArray(
+  Tₛ = Tₛ₀,
+  h = h₀,
+  melting_water = water)
+
+constants_and_parameters = (
+  budyko_sellers_absorbed_radiation_α = α,
+  budyko_sellers_outgoing_radiation_A = A,
+  budyko_sellers_outgoing_radiation_B = B,
+  budyko_sellers_energy_C = C,
+  budyko_sellers_diffusion_D = D,
+  budyko_sellers_cosϕᵖ = cosϕᵖ,
+  budyko_sellers_diffusion_cosϕᵈ = cosϕᵈ,
+  halfar_n = n,
+  halfar_stress_ρ = halfar_ρ,
+  halfar_stress_g = g,
+  melting_Dₕ₂ₒ = Dₕ₂ₒ)
+
+# Define how symbols map to Julia functions
+
+function generate(sd, my_symbol; hodge=GeometricHodge())
+  op = @match my_symbol begin
+    :♯ => begin
+      sharp_mat = ♯_mat(sd, AltPPSharp())
+      x -> sharp_mat * x
+    end
+    :mag => x -> begin
+      norm.(x)
+    end
+    :^ => (x,y) -> x .^ y
+    :* => (x,y) -> x .* y
+    x => error("Unmatched operator $my_symbol")
+  end
+  return (args...) -> op(args...)
+end
+```
+
+## Generate and run simulation 
+
+The composed model is generated and executed for 100 years. The model is run twice to demonstrate the speed of the model after the simulation code is precompiled by the first run.
+
+We will save the final ice thickness data in a .jld2 file, an HDF5-compatible file format. We will also save the latitude and longitude of the points on the sphere.
+
+``` @example DEC
+sim = eval(gensim(budyko_sellers_halfar_water))
+fₘ = sim(s, generate)
+
+tₑ = 100.0
+
+@info("Precompiling Solver")
+prob = ODEProblem(fₘ, u₀, (0, 1e-4), constants_and_parameters)
+soln = solve(prob, Tsit5())
+soln.retcode != :Unstable || error("Solver was not stable")
+
+@info("Solving")
+prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters)
+soln = solve(prob, Tsit5())
+@show soln.retcode
+@info("Done")
+
+save("ice.jld2",
+  Dict("lat" => map(x -> x[1], p_sph), "lon" => map(x -> x[2], p_sph), "ice" => soln(tₑ).h))
+
+(extrema(soln(0.0).h), extrema(soln(tₑ).h))
+```
+
+## Visualize 
+
+Let's visualize the initial conditions for ice height and the ice height after 100 years.
+
+### Initial ice height
+
+``` @example DEC
+f = Figure()
+ax = LScene(f[1,1], scenekw=(lights=[],))
+msh = mesh!(ax, s_plots, color=soln.u[end].h)
+Colorbar(f[1,2], msh)
+f
+```
+
+### Ice height after 100 years
+
+``` @example DEC
+f = Figure()
+ax = LScene(f[1,1], scenekw=(lights=[],))
+msh = mesh!(ax, s_plots, color=soln.u[end].h)
+Colorbar(f[1,2], msh)
+f
+```
+
+```@example DEC
+run(`rm $ice_thickness_file`) # hide
+```
+