diff --git a/examples/arctic_basin_seasonal_cycle.jl b/examples/arctic_basin_seasonal_cycle.jl
index 7f343f66..ffbef5c3 100644
--- a/examples/arctic_basin_seasonal_cycle.jl
+++ b/examples/arctic_basin_seasonal_cycle.jl
@@ -1,59 +1,90 @@
-#####
-##### The purpose of this script is to reproduce some results from Semtner (1976)
-#####
+# # Arctic basin seasonal cycle example
+#
+# This example reproduces results from Semtner (1976), which simulates the seasonal
+# cycle of sea ice in the Arctic basin. The model uses climatological forcing data
+# from Fletcher (1965) for shortwave radiation, longwave radiation, sensible heat
+# flux, and latent heat flux. This example demonstrates how to:
+#
+#   * use time-varying climatological forcing data,
+#   * set up seasonal cycles with `FieldTimeSeries` and cyclical indexing,
+#   * combine multiple heat flux components,
+#   * simulate multi-year seasonal cycles.
+#
+# ## Install dependencies
+#
+# First let's make sure we have all required packages installed.
+#
+# ```julia
+# using Pkg
+# pkg"add Oceananigans, ClimaSeaIce, CairoMakie"
+# ```
+#
+# ## The physical domain
+#
+# We use a single grid point to model a horizontally uniform ice column:
 
 using Oceananigans
 using Oceananigans.Units
 using Oceananigans.Units: Time
-using Oceananigans.OutputReaders
 using ClimaSeaIce
+using CairoMakie
 
-# Generate a 0D grid for a single column slab model 
 grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat))
 
-# Forcings (Semtner 1976, table 1), originally tabulated by Fletcher (1965)
-# Note: these are in kcal, which was deprecated in the ninth General Conference on Weights and
-# Measures in 1948. We convert these to Joules (and then to Watts) below.
+# ## Climatological forcing data
+#
+# The forcing data comes from Semtner (1976, table 1), originally tabulated by
+# Fletcher (1965). Note: these are in kcal, which was deprecated in the ninth
+# General Conference on Weights and Measures in 1948. We convert these to Joules
+# (and then to Watts) below.
+#
 # Month:        Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep   Oct    Nov    Dec
-tabulated_shortwave = - [   0,      0,    1.9,    9.9,   17.7,   19.2,   13.6,    9.0,    3.7,   0.4,      0,      0] .* 1e4 # to convert to kcal / m²
-tabulated_longwave  = - [10.4,   10.3,   10.3,   11.6,   15.1,   18.0,   19.1,   18.7,   16.5,  13.9,   11.2,   10.9] .* 1e4 # to convert to kcal / m²
-tabulated_sensible  = - [1.18,   0.76,   0.72,   0.29,  -0.45,  -0.39,  -0.30,  -0.40,  -0.17,   0.1,   0.56,   0.79] .* 1e4 # to convert to kcal / m²
-tabulated_latent    = - [   0,  -0.02,  -0.03,   -0.09, -0.46,  -0.70,  -0.64,  -0.66,  -0.39,  -0.19, -0.01,  -0.01] .* 1e4 # to convert to kcal / m²
 
-# Pretend every month is just 30 days
+tabulated_shortwave = - [   0,      0,    1.9,    9.9,   17.7,   19.2,   13.6,    9.0,    3.7,   0.4,      0,      0] .* 1e4 # kcal m⁻²
+tabulated_longwave  = - [10.4,   10.3,   10.3,   11.6,   15.1,   18.0,   19.1,   18.7,   16.5,  13.9,   11.2,   10.9] .* 1e4 # kcal m⁻²
+tabulated_sensible  = - [1.18,   0.76,   0.72,   0.29,  -0.45,  -0.39,  -0.30,  -0.40,  -0.17,   0.1,   0.56,   0.79] .* 1e4 # kcal m⁻²
+tabulated_latent    = - [   0,  -0.02,  -0.03,   -0.09, -0.46,  -0.70,  -0.64,  -0.66,  -0.39,  -0.19, -0.01,  -0.01] .* 1e4 # kcal m⁻²
+
+# We assume every month has 30 days and set up the time array:
+
 Nmonths = 12
 month_days = 30
 year_days = month_days * Nmonths
 times_days = 15:30:(year_days - 15)
 times = times_days .* day # times in seconds
 
-# Sea ice emissivity
+# Sea ice emissivity:
+
 ϵ = 1
 
-# Convert fluxes to the right units
+# Convert fluxes to the right units (W m⁻²):
+
 kcal_to_joules = 4184
-tabulated_shortwave .*= kcal_to_joules / (month_days * days) 
+tabulated_shortwave .*= kcal_to_joules / (month_days * days)
 tabulated_longwave  .*= kcal_to_joules / (month_days * days) .* ϵ
 tabulated_sensible  .*= kcal_to_joules / (month_days * days)
 tabulated_latent    .*= kcal_to_joules / (month_days * days)
 
-using CairoMakie
+# Let's visualize the forcing data:
 
 fig = Figure()
-ax = Axis(fig[1, 1])
-lines!(ax, times, tabulated_shortwave)
-lines!(ax, times, tabulated_longwave)
-lines!(ax, times, tabulated_sensible)
-lines!(ax, times, tabulated_latent)
-
-display(fig)
-
-# Make them into a FieldTimeSeries for better manipulation
-
-Rs = FieldTimeSeries{Nothing, Nothing, Nothing}(grid, times; time_indexing = Cyclical())
-Rl = FieldTimeSeries{Nothing, Nothing, Nothing}(grid, times; time_indexing = Cyclical())
-Qs = FieldTimeSeries{Nothing, Nothing, Nothing}(grid, times; time_indexing = Cyclical())
-Ql = FieldTimeSeries{Nothing, Nothing, Nothing}(grid, times; time_indexing = Cyclical())
+ax = Axis(fig[1, 1], xlabel="Time (days)", ylabel="Heat flux (W m⁻²)")
+lines!(ax, times ./ day, tabulated_shortwave, label="Shortwave")
+lines!(ax, times ./ day, tabulated_longwave, label="Longwave")
+lines!(ax, times ./ day, tabulated_sensible, label="Sensible")
+lines!(ax, times ./ day, tabulated_latent, label="Latent")
+axislegend(ax)
+current_figure() #hide
+
+# ## Creating time-varying flux functions
+#
+# We create `FieldTimeSeries` objects with cyclical indexing so the forcing
+# repeats annually:
+
+Rs = FieldTimeSeries{Nothing, Nothing, Nothing}(grid, times; time_indexing = Oceananigans.OutputReaders.Cyclical())
+Rl = FieldTimeSeries{Nothing, Nothing, Nothing}(grid, times; time_indexing = Oceananigans.OutputReaders.Cyclical())
+Qs = FieldTimeSeries{Nothing, Nothing, Nothing}(grid, times; time_indexing = Oceananigans.OutputReaders.Cyclical())
+Ql = FieldTimeSeries{Nothing, Nothing, Nothing}(grid, times; time_indexing = Oceananigans.OutputReaders.Cyclical())
 
 for (i, time) in enumerate(times)
     set!(Rs[i], tabulated_shortwave[i:i])
@@ -62,20 +93,30 @@ for (i, time) in enumerate(times)
     set!(Ql[i], tabulated_latent[i:i])
 end
 
+# We define flux functions that linearly interpolate the time series:
+
 @inline function linearly_interpolate_flux(i, j, grid, Tₛ, clock, model_fields, flux)
     t  = Time(clock.time)
     return flux[i, j, 1, t]
 end
 
+# For shortwave radiation, we account for surface albedo:
+
 @inline function linearly_interpolate_solar_flux(i, j, grid, Tₛ, clock, model_fields, flux)
     Q = linearly_interpolate_flux(i, j, grid, Tₛ, clock, model_fields, flux)
-    α = ifelse(Tₛ < -0.1, 0.75, 0.64)
+    α = ifelse(Tₛ < -0.1, 0.75, 0.64) # albedo depends on surface temperature
     return Q * (1 - α)
 end
 
-# NOTE: Semtner (1976) uses a wrong value for the Stefan-Boltzmann constant (roughly 2% higher) to match
-# the results of Maykut & Unterstainer (196(5)). We use the same value here for comparison purposes.
-σ = 5.67e-8 * 1.02 # Wrong value!!
+# !!! note "Slightly wrong value for Stefan-Boltzmann constant"
+#
+#     Semtner (1976) uses a wrong value for the Stefan-Boltzmann constant
+#     (roughly 2% higher) to match the results of Maykut & Untersteiner (1965).
+#     We use here the same value here for comparison purposes.
+
+σ = 5.67e-8 * 1.02 # Wrong value, but used for comparison!
+
+# We create flux functions for each component:
 
 Q_shortwave = FluxFunction(linearly_interpolate_solar_flux, parameters=Rs)
 Q_longwave  = FluxFunction(linearly_interpolate_flux,       parameters=Rl)
@@ -83,14 +124,27 @@ Q_sensible  = FluxFunction(linearly_interpolate_flux,       parameters=Qs)
 Q_latent    = FluxFunction(linearly_interpolate_flux,       parameters=Ql)
 Q_emission  = RadiativeEmission(emissivity=ϵ, stefan_boltzmann_constant=σ)
 
+# We combine all top heat fluxes:
+
 top_heat_flux = (Q_shortwave, Q_longwave, Q_sensible, Q_latent, Q_emission)
 
+# ## Building the model
+#
+# We assemble the model and initialize it with 30 cm of ice at full concentration:
+
 model = SeaIceModel(grid; top_heat_flux)
-set!(model, h=0.3, ℵ=1) # We start from 300cm of ice and full concentration
+set!(model, h=0.3, ℵ=1) # Start from 30 cm of ice and full concentration
+
+# ## Running the simulation
+#
+# We run the simulation for 30 years with an 8-hour time step:
+
+simulation = Simulation(model, Δt=8hours, stop_time=30 * 360days)
 
-simulation = Simulation(model, Δt=8hours, stop_time= 30 * 360days)
+# ## Collecting data
+#
+# We set up a callback to accumulate time series data:
 
-# Accumulate data
 series = []
 
 function accumulate_timeseries(sim)
@@ -106,7 +160,9 @@ simulation.callbacks[:save] = Callback(accumulate_timeseries)
 
 run!(simulation)
 
-# Extract and visualize data
+# ## Visualizing the results
+#
+# We extract the time series data and visualize the evolution:
 
 t = [datum[1] for datum in series]
 h = [datum[2] for datum in series]
@@ -121,16 +177,20 @@ fig = Figure(size=(1000, 1200))
 axT = Axis(fig[1, 1], xlabel="Time (days)", ylabel="Top temperature (ᵒC)")
 axh = Axis(fig[2, 1], xlabel="Time (days)", ylabel="Ice thickness (m)")
 axℵ = Axis(fig[3, 1], xlabel="Time (days)", ylabel="Ice concentration (-)")
-axQ = Axis(fig[4, 1], xlabel="Time (days)", ylabel="Top heat flux (W/m²)")
+axQ = Axis(fig[4, 1], xlabel="Time (days)", ylabel="Top heat flux (W m⁻²)")
 
-lines!(axT, t / day, T)
-lines!(axh, t / day, h)
-lines!(axℵ, t / day, ℵ)
-lines!(axQ, t / day, Q)
+lines!(axT, t ./ day, T)
+lines!(axh, t ./ day, h)
+lines!(axℵ, t ./ day, ℵ)
+lines!(axQ, t ./ day, Q)
 
-display(fig)
+current_figure() #hide
 
 save("ice_timeseries.png", fig)
 nothing # hide
 
 # ![](ice_timeseries.png)
+#
+# The results show the seasonal cycle of sea ice, with ice growing in winter and
+# melting in summer. After several years, the model reaches a quasi-equilibrium
+# seasonal cycle.
diff --git a/examples/diffusive_ice_column_model.jl b/examples/diffusive_ice_column_model.jl
index 93a298b5..d0a70f52 100644
--- a/examples/diffusive_ice_column_model.jl
+++ b/examples/diffusive_ice_column_model.jl
@@ -1,3 +1,28 @@
+# # Diffusive ice column model example
+#
+# This example demonstrates the use of `EnthalpyMethodSeaIceModel`, which resolves
+# the vertical structure of sea ice using an enthalpy-based formulation. The model
+# includes internal heat conduction and tracks the phase transitions between ice and
+# water within the ice column. This example demonstrates how to:
+#
+#   * set up an `EnthalpyMethodSeaIceModel` with molecular diffusivity,
+#   * prescribe time-varying boundary conditions,
+#   * compute ice thickness from the temperature profile,
+#   * visualize the evolution of temperature, enthalpy, porosity, and diffusivity.
+#
+# ## Install dependencies
+#
+# First let's make sure we have all required packages installed.
+#
+# ```julia
+# using Pkg
+# pkg"add Oceananigans, ClimaSeaIce, CairoMakie"
+# ```
+#
+# ## The physical domain
+#
+# We build a one-dimensional vertical grid with 20 grid points and 10 cm resolution:
+
 using ClimaSeaIce.EnthalpyMethodSeaIceModels: EnthalpyMethodSeaIceModel, MolecularDiffusivity
 using Oceananigans
 using Oceananigans.Units
@@ -5,33 +30,39 @@ using Oceananigans.Operators: Δzᶜᶜᶠ
 using Oceananigans.Grids: znode
 using CairoMakie
 
-#####
-##### Set up a EnthalpyMethodSeaIceModel
-#####
-
-# Build a grid with 10 cm resolution
 grid = RectilinearGrid(size=20, z=(-1, 0), topology=(Flat, Flat, Bounded))
 
-# Set up a simple problem and build the ice model
+# ## Model configuration
+#
+# We set up a simple problem with molecular diffusivity. The diffusivity differs
+# for ice and water phases:
+
 closure = MolecularDiffusivity(grid, κ_ice=1e-5, κ_water=1e-6)
 
-#####
-##### Create temperature boundary conditions
-#####
+# ## Temperature boundary conditions
+#
+# We define time-varying boundary conditions for the top (air-ice interface) and
+# bottom (ice-ocean interface) of the ice column:
+
+initial_air_ice_temperature = -5 # °C
+top_T_amplitude = 5 # °C
+top_T_slope = -0.5 / day # °C s⁻¹
+
+# Information about ocean cooling:
 
-initial_air_ice_temperature = -5
-top_T_amplitude = 5
-top_T_slope = -0.5 / day # ᵒC s⁻¹
+initial_ice_ocean_temperature = 1.1 # °C
+bottom_T_slope = -0.1 / day # °C s⁻¹
 
-# Information about ocean cooling
-initial_ice_ocean_temperature = 1.1
-bottom_T_slope = -0.1 / day # ᵒC s⁻¹
+# The top temperature oscillates with a daily cycle and has a long-term cooling trend:
 
-# Calculate BCs
 air_ice_temperature(x, y, t) = top_T_slope * t + top_T_amplitude * sin(2π*t/day) + initial_air_ice_temperature
+
+# The bottom temperature cools linearly:
+
 ice_ocean_temperature(x, y, t) = bottom_T_slope * t + initial_ice_ocean_temperature
 
-# Plot boundary condition functions
+# Let's visualize the boundary condition functions:
+
 dt = 10minutes
 tf = 10days
 t = 0:dt:tf
@@ -42,32 +73,39 @@ ax = Axis(fig[1, 1], title="Boundary Conditions", xlabel="Time (s)", ylabel="Tem
 lines!(ax, t, air_ice_temperature.(0, 0, t), label="Air-ice surface temperature")
 lines!(ax, t, ice_ocean_temperature.(0, 0, t), label="Ice-ocean temperature")
 axislegend(ax)
-     
-display(fig)
+
+current_figure() #hide
+
+# We create the boundary conditions:
 
 top_T_bc = ValueBoundaryCondition(air_ice_temperature)
 bottom_T_bc = ValueBoundaryCondition(ice_ocean_temperature)
 T_bcs = FieldBoundaryConditions(top=top_T_bc, bottom=bottom_T_bc)
 
+# ## Building the model
+#
+# We assemble the model with the grid, closure, and boundary conditions:
+
 model = EnthalpyMethodSeaIceModel(; grid, closure, boundary_conditions=(; T=T_bcs))
 
-# Initialize and run
-set!(model, T=initial_ice_ocean_temperature)
+# We initialize the temperature field:
 
-H = model.state.H
+set!(model, T=initial_ice_ocean_temperature)
 
-# We're using explicit time stepping. The CFL condition is
+# ## Time stepping
 #
-#   Δt < Δz² / κ ≈ 0.1² / 1e-6 ≈ 1e4,
+# We're using explicit time stepping. The CFL number is ``κ Δt / Δz²`` and thus
+# to ensure that remains below, e.g., 0.1 we choose our timestep ``Δt`` accordingly.
 
 κ = 1e-5
-Δz = Δzᶜᶜᶠ(1, 1, 1, grid)
-Δt = 0.1 * Δz^2 / κ
+Δz_min = minimum_zspacing(grid)
+Δt = 0.1 * Δz_min^2 / κ
+
 simulation = Simulation(model; Δt)
 
-#####
-##### Set up diagnostics
-#####
+# ## Diagnostics
+#
+# We set up diagnostics to compute ice thickness and collect profiles:
 
 const c = Center()
 
@@ -97,8 +135,10 @@ end
 
 simulation.callbacks[:thickness] = Callback(compute_ice_thickness, IterationInterval(1))
 
-tt = [] 
-Tt = [] 
+# We also collect profiles of temperature, enthalpy, porosity, and diffusivity:
+
+tt = []
+Tt = []
 Ht = []
 ϕt = []
 κt = []
@@ -111,7 +151,7 @@ function grab_profiles!(sim)
     κ = sim.model.closure.κ
 
     # Extract interior data (excluding halos)
-    Ti = interior(T, 1, 1, :)  
+    Ti = interior(T, 1, 1, :)
     Hi = interior(H, 1, 1, :)
     ϕi = interior(ϕ, 1, 1, :)
     κi = interior(κ, 1, 1, :)
@@ -126,23 +166,25 @@ function grab_profiles!(sim)
     return nothing
 end
 
-simulation.callbacks[:grabber] = Callback(grab_profiles!, TimeInterval(1hour)) #SpecifiedTimes(10minutes, 30minutes, 1hour))
+simulation.callbacks[:grabber] = Callback(grab_profiles!, TimeInterval(1hour))
 
-#####
-##### Run the simulation
-#####
+# ## Running the simulation
+#
+# We run the simulation for 10 days:
 
 simulation.stop_time = 10days
 run!(simulation)
 
-# Make a plot
+# ## Visualizing the results
+#
+# We create an interactive visualization with sliders to explore the evolution:
 
 fig = Figure()
 
 axT = Axis(fig[1, 1], xlabel="Temperature (ᵒC)", ylabel="z (m)")
 axH = Axis(fig[1, 2], xlabel="Enthalpy (J m⁻³)", ylabel="z (m)")
 axϕ = Axis(fig[1, 3], xlabel="Porosity", ylabel="z (m)")
-axκ = Axis(fig[1, 4], xlabel="Diffusivity", ylabel="z (m)")
+axκ = Axis(fig[1, 4], xlabel="Diffusivity (m² s⁻¹)", ylabel="z (m)")
 axh = Axis(fig[2, 1:4], xlabel="Time (hours)", ylabel="Ice thickness (m)")
 axq = Axis(fig[3, 1:4], xlabel="Time (hours)", ylabel="Surface temperatures (ᵒC)")
 
@@ -156,14 +198,6 @@ n = slider.value
 
 z = znodes(model.state.T)
 
-#=
-# TODO: calculate analytical solution
-ℒ = model.fusion_enthalpy
-ΔT = ocean_temperature - atmosphere_temperature 
-c = ice_heat_capacity 
-f(λ) = λ * exp(λ^2) * erf(λ) - St / sqrt(π)
-=#
-
 tn = @lift tt[$n]
 tnh = @lift tt[$n] / hour
 Tn = @lift Tt[$n]
@@ -176,7 +210,6 @@ scatterlines!(axH, Hn, z; label)
 scatterlines!(axϕ, ϕn, z; label)
 scatterlines!(axκ, κn, z; label)
 
-lines!(axh, th ./ hour, ht)
 lines!(axh, th ./ hour, ht)
 
 hn = @lift begin
@@ -192,5 +225,8 @@ vlines!(axq, tnh)
 axislegend(axq, position=:lb)
 axislegend(axT, position=:lb)
 
-display(fig)
+current_figure() #hide
 
+# The visualization shows the evolution of the temperature profile, enthalpy, porosity
+# (liquid fraction), and thermal diffusivity within the ice column. The ice thickness
+# evolves as the column freezes and melts in response to the boundary conditions.
diff --git a/examples/freezing_bucket.jl b/examples/freezing_bucket.jl
index 692645c0..47a26217 100644
--- a/examples/freezing_bucket.jl
+++ b/examples/freezing_bucket.jl
@@ -1,130 +1,150 @@
-# # A freezing bucket
+# # Freezing bucket example
 #
-# A common laboratory experiment freezes an insultated bucket of water
-# from the top down, using a metal lid to keep the top of the bucket
-# at some constant, very cold temperature. In this example, we simulate such
-# a scenario using the `SeaIceModel`. Here, the bucket is perfectly insulated
-# and infinitely deep, like many buckets are: if the `Simulation` is run for longer,
-# the ice will keep freezing, and freezing, and will never run out of water.
-# Also, the water in the infinite bucket is (somehow) all at the same temperature,
-# in equilibrium with the ice-water interface (and therefore fixed at the melting
-# temperature). Yep, this kind of thing happens _all the time_.
+# A common laboratory experiment freezes an insulated bucket of water from the top
+# down, using a metal lid to keep the top of the bucket at some constant, very cold
+# temperature. In this example, we simulate such a scenario using the `SeaIceModel`.
+# Here, the bucket is perfectly insulated and infinitely deep: if the `Simulation`
+# is run for longer, the ice will keep freezing, and freezing, and will never run
+# out of water. Also, the water in the infinite bucket is (somehow) all at the same
+# temperature, in equilibrium with the ice-water interface (and therefore fixed at
+# the melting temperature). This example demonstrates how to:
 #
-# We start by `using Oceananigans` to bring in functions for building grids
-# and `Simulation`s and the like.
+#   * use `SlabSeaIceThermodynamics` with prescribed boundary conditions,
+#   * configure internal heat conduction,
+#   * use `FluxFunction` for frazil ice formation,
+#   * collect and visualize time series data.
+#
+# ## Install dependencies
+#
+# First let's make sure we have all required packages installed.
+#
+# ```julia
+# using Pkg
+# pkg"add Oceananigans, ClimaSeaIce, CairoMakie"
+# ```
+#
+# ## Using `ClimaSeaIce.jl`
+#
+# Write
 
 using Oceananigans
 using Oceananigans.Units
-
-# Next we `using ClimaSeaIce` to get some ice-specific names.
-
 using ClimaSeaIce
+using CairoMakie
 
-# # An infinitely deep bucket with a single grid point
+# to load ClimaSeaIce functions and objects into our script.
+#
+# ## An infinitely deep bucket with a single grid point
 #
-# Perhaps surprisingly, we need just one grid point
-# to model an possibly infinitely thick slab of ice with `SlabSeaIceModel`.
-# We would only need more than 1 grid point if our boundary conditions
-# vary in the horizontal direction.
+# Perhaps surprisingly, we need just one grid point to model a possibly infinitely
+# thick slab of ice with `SlabSeaIceModel`. We would only need more than 1 grid point
+# if our boundary conditions vary in the horizontal direction:
 
 grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat))
 
-# Next, we build our model of an ice slab freezing into a bucket.
-# We start by defining a constant internal `ConductiveFlux` with 
-# ice_conductivity 
+# ## Model configuration
+#
+# We build our model of an ice slab freezing into a bucket. We start by defining
+# a constant internal `ConductiveFlux` with ice conductivity:
 
-conductivity = 2 # kg m s⁻³ K⁻¹
+conductivity = 2 # W m⁻¹ K⁻¹
 internal_heat_flux = ConductiveFlux(; conductivity)
 
-# Note that other units besides Celsius _can_ be used, but that requires
-# setting model.phase_transitions` with appropriate parameters.
-# We set the ice heat capacity and density as well,
+# Note that other units besides Celsius _can_ be used, but that requires setting
+# `model.phase_transitions` with appropriate parameters. We set the ice heat capacity
+# and density as well:
 
 ice_heat_capacity = 2100 # J kg⁻¹ K⁻¹
 ice_density = 900 # kg m⁻³
 phase_transitions = PhaseTransitions(; ice_heat_capacity, ice_density)
 
-# We set the top ice temperature,
+# We set the top ice temperature:
 
-top_temperature = -10 # ᵒC
-top_heat_boundary_condition = PrescribedTemperature(-10)
+top_temperature = -10 # °C
+top_heat_boundary_condition = PrescribedTemperature(top_temperature)
 
-# Construct the ice_thermodynamics of sea ice, for this we use a simple
-# slab sea ice representation of ice_thermodynamics
+# ## Constructing the thermodynamics
+#
+# We construct the ice thermodynamics using a simple slab sea ice representation:
 
 ice_thermodynamics = SlabSeaIceThermodynamics(grid;
-                                          internal_heat_flux,
-                                          phase_transitions,
-                                          top_heat_boundary_condition)
+                                              internal_heat_flux,
+                                              phase_transitions,
+                                              top_heat_boundary_condition)
 
-# We also prescribe a frazil ice heat flux that stops when the ice has reached a concentration of 1.
-# This heat flux represents the initial ice formation from a liquid bucket.
+# ## Frazil ice formation
+#
+# We also prescribe a frazil ice heat flux that stops when the ice has reached a
+# concentration of 1. This heat flux represents the initial ice formation from a
+# liquid bucket:
 
 @inline frazil_ice_formation(i, j, grid, Tuᵢ, clock, fields) = - (1 - fields.ℵ[i, j, 1]) # W m⁻²
 
 bottom_heat_flux = FluxFunction(frazil_ice_formation)
 
-# Then we assemble it all into a model.
+# ## Building the model
+#
+# Then we assemble it all into a model:
 
 model = SeaIceModel(grid; ice_thermodynamics, bottom_heat_flux)
 
-# Note that the default bottom heat boundary condition for `SlabSeaIceThermodynamics` is
-# `IceWaterThermalEquilibrium` with freshwater. That's what we want!
+# Note that the default bottom heat boundary condition for `SlabSeaIceThermodynamics`
+# is `IceWaterThermalEquilibrium` with freshwater. That's what we want!
 
 model.ice_thermodynamics.heat_boundary_conditions.bottom
 
-# Ok, we're ready to freeze the bucket for 10 straight days.
-# The ice will start forming suddenly due to the frazil ice heat flux and then eventually
-# grow more slowly.
+# ## Running a simulation
+#
+# We're ready to freeze the bucket for 10 straight days. The ice will start forming
+# suddenly due to the frazil ice heat flux and then eventually grow more slowly:
 
 simulation = Simulation(model, Δt=10minute, stop_time=10days)
 
-# # Collecting data and running the simulation
+# ## Collecting data
 #
-# Before simulating the freezing bucket, we set up a `Callback` to create
-# a timeseries of the ice thickness saved at every time step.
+# Before simulating the freezing bucket, we set up a `Callback` to create a time
+# series of the ice thickness saved at every time step:
 
-## Container to hold the data
 timeseries = []
 
-## Callback function to collect the data from the `sim`ulation
 function accumulate_timeseries(sim)
     h = sim.model.ice_thickness
     ℵ = sim.model.ice_concentration
     push!(timeseries, (time(sim), first(h), first(ℵ)))
 end
 
-## Add the callback to `simulation`
 simulation.callbacks[:save] = Callback(accumulate_timeseries)
 
-# Now we're ready to hit the Big Red Button (it should run pretty quick):
+# Now we're ready to run the simulation:
+
 run!(simulation)
 
-# # Visualizing the result
+# ## Visualizing the results
 #
-# It'd be a shame to run such a "cool" simulation without looking at the
-# results. We'll visualize it with Makie.
-
-using CairoMakie
+# It'd be a shame to run such a "cool" simulation without looking at the results.
+# We'll visualize it with Makie:
 
 # `timeseries` is a `Vector` of `Tuple`. So we have to do a bit of processing
 # to build `Vector`s of time `t` and thickness `h`. It's not much work though:
+
 t = [datum[1] for datum in timeseries]
 h = [datum[2] for datum in timeseries]
 ℵ = [datum[3] for datum in timeseries]
 V = h .* ℵ
 
 # Just for fun, we also compute the velocity of the ice-water interface:
+
 dVdt = @. (h[2:end] .* ℵ[2:end] - h[1:end-1] .* ℵ[1:end-1]) / simulation.Δt
 
 # All that's left, really, is to put those `lines!` in an `Axis`:
+
 set_theme!(Theme(fontsize=24, linewidth=4))
 
 fig = Figure(size=(1600, 700))
 
 axh = Axis(fig[1, 1], xlabel="Time (days)", ylabel="Ice thickness (cm)")
 axℵ = Axis(fig[1, 2], xlabel="Time (days)", ylabel="Ice concentration (-)")
-axV = Axis(fig[1, 3], xlabel="Ice Volume (cm)", ylabel="Freezing rate (μm s⁻¹)")
+axV = Axis(fig[1, 3], xlabel="Ice volume (cm)", ylabel="Freezing rate (μm s⁻¹)")
 
 lines!(axh, t ./ day, 1e2 .* h)
 lines!(axℵ, t ./ day, ℵ)
@@ -134,7 +154,6 @@ save("freezing_bucket.png", fig)
 nothing # hide
 
 # ![](freezing_bucket.png)
-
+#
 # If you want more ice, you can increase `simulation.stop_time` and
 # `run!(simulation)` again (or just re-run the whole script).
-
diff --git a/examples/freezing_of_a_lake.jl b/examples/freezing_of_a_lake.jl
index 51a50ed2..9973bf8b 100644
--- a/examples/freezing_of_a_lake.jl
+++ b/examples/freezing_of_a_lake.jl
@@ -1,12 +1,28 @@
-# # Freezing of a lake
+# # Freezing of a lake example
 #
 # In this example we simulate the freezing of a lake in winter. The lake is
-# represented by four points that start at 1ᵒC and are cooled down by an atmosphere with 
-# temperatures of -20, -10, -1, and -0.1ᵒC. The lake is 10 m deep and not subjected to 
-# radiative transfer (this lake is covered in a wind tunnel under which we blow some cold air).
+# represented by four grid points that start at 1°C and are cooled down by an
+# atmosphere with temperatures of -20, -10, -5, and 0°C. The lake is 10 m deep
+# and not subjected to radiative transfer (this lake is covered in a wind tunnel
+# under which we blow some cold air). This example demonstrates how to:
 #
-# We start by `using Oceananigans` to bring in functions for building grids
-# and `Simulation`s and the like.
+#   * couple a simple lake model with sea ice thermodynamics,
+#   * use `FluxFunction` for complex boundary conditions,
+#   * track energy budgets.
+#
+# ## Install dependencies
+#
+# First let's make sure we have all required packages installed.
+#
+# ```julia
+# using Pkg
+# pkg"add Oceananigans, ClimaSeaIce, CairoMakie"
+# ```
+#
+# ## The physical domain
+#
+# We generate a one-dimensional grid with 4 grid cells to model different lake
+# columns subject to different atmospheric temperatures:
 
 using Oceananigans
 using Oceananigans.Units
@@ -14,23 +30,24 @@ using ClimaSeaIce
 using ClimaSeaIce.SeaIceThermodynamics: latent_heat
 using CairoMakie
 
-# Generate a 1D grid for difference ice columns subject to different solar insolation
-
 grid = RectilinearGrid(size=4, x=(0, 1), topology=(Periodic, Flat, Flat))
 
+# ## Atmospheric forcing
+#
 # The sensible heat flux from the atmosphere is represented by a `FluxFunction`.
+# We define the atmospheric parameters:
 
 atmosphere = (
-    transfer_coefficient     = 1e-3,  # Unitless
+    transfer_coefficient     = 1e-3,  # unitless
     atmosphere_density       = 1.225, # kg m⁻³
-    atmosphere_heat_capacity = 1004,  # 
-    atmosphere_temperature   = [-20, -10, -5, 0],    # ᵒC
+    atmosphere_heat_capacity = 1004,  # J kg⁻¹ K⁻¹
+    atmosphere_temperature   = [-20, -10, -5, 0], # °C
     atmosphere_wind_speed    = 5,     # m s⁻¹
     atmosphere_ice_flux      = [0.0, 0.0, 0.0, 0.0], # W m⁻²
 )
-    
-# Flux is positive (cooling by fluxing heat up away from upper surface)
-# when Tₐ < Tᵤ:
+
+# The flux is positive (cooling by fluxing heat upward away from the upper surface)
+# when the atmosphere temperature is less than the surface temperature:
 
 @inline function sensible_heat_flux(i, j, grid, Tᵤ, clock, fields, atmosphere)
     Cₛ = atmosphere.transfer_coefficient
@@ -41,30 +58,33 @@ atmosphere = (
     ℵ  = fields.ℵ[i, j, 1]
     Qₐ = atmosphere.atmosphere_ice_flux
 
-    Qₐ[i] =  ifelse(ℵ == 0, zero(grid), Cₛ * ρₐ * cₐ * uₐ * (Tᵤ - Tₐ))
+    Qₐ[i] = ifelse(ℵ == 0, zero(grid), Cₛ * ρₐ * cₐ * uₐ * (Tᵤ - Tₐ))
 
-    return Qₐ[i] 
+    return Qₐ[i]
 end
 
-# We also evolve a bucket freshwater lake that cools down and freezes from below
-# generating fresh sea-ice (or lake-ice in this case?).
-# We set the Δt of the lake to 10 minutes. This time step will be used to also for the sea-ice
-# model.
+# ## Lake model
+#
+# We also evolve a simple bucket freshwater lake that cools down and freezes from
+# above, generating fresh sea ice (or lake ice in this case). We set the time step
+# of the lake to 10 minutes, which will also be used for the sea ice model:
 
 lake = (
     lake_density         = 1000, # kg m⁻³
-    lake_heat_capacity   = 4000,  # 
-    lake_temperature     = [1.0, 1.0, 1.0, 1.0], # ᵒC
+    lake_heat_capacity   = 4000,  # J kg⁻¹ K⁻¹
+    lake_temperature     = [1.0, 1.0, 1.0, 1.0], # °C
     lake_depth           = 10, # m
     lake_ice_flux        = [0.0, 0.0, 0.0, 0.0], # W m⁻²
     atmosphere_lake_flux = [0.0, 0.0, 0.0, 0.0], # W m⁻²
     Δt                   = 10minutes
 )
 
+# ## Bottom heat flux function
+#
 # We build a flux function that serves three purposes:
-# 1. computing the cooling of the lake caused by the atmosphere
+# 1. Computing the cooling of the lake caused by the atmosphere
 # 2. If the lake temperature is low enough, freezing the lake from above
-# 3. and adding the heat flux to the bottom of the ice
+# 3. Adding the heat flux to the bottom of the ice
 
 @inline function advance_lake_and_frazil_flux(i, j, grid, Tuᵢ, clock, fields, parameters)
     atmos = parameters.atmosphere
@@ -81,40 +101,47 @@ lake = (
     Δ  = lake.lake_depth
     Δt = lake.Δt
     ℵ  = fields.ℵ[i, j, 1]
-    
+
     Qᵗ = lake.atmosphere_lake_flux
     Qᵇ = lake.lake_ice_flux
-    
-    Qₐ = Cₛ * ρₐ * cₐ * uₐ * (Tₐ - Tₒ[i]) * (1 - ℵ)
-
-    Tₒ[i] = Tₒ[i] + Qₐ / (ρₒ * cₒ) * Δt
 
-    Qᵢ = ρₒ * cₒ * (Tₒ[i] - 0) / Δt * Δ # W m⁻²
-    Qᵢ = min(Qᵢ, zero(Qᵢ)) # We only freeze, not melt
-
-    Tₒ[i] = ifelse(Qᵢ == 0, Tₒ[i], zero(Qᵢ))
+    Qₐ = Cₛ * ρₐ * cₐ * uₐ * (Tₐ - Tₒ[i]) * (1 - ℵ)
+    Tⁿ = Tₒ[i] + Qₐ / (ρₒ * cₒ) * Δt
+    Qᵢ = ρₒ * cₒ * (Tⁿ - 0) / Δt * Δ # W m⁻²
+    Qᵢ = min(Qᵢ, zero(Qᵢ))
 
+    Tₒ[i] = ifelse(Qᵢ == 0, Tⁿ, zero(Qᵢ))
     Qᵗ[i] = Qₐ
     Qᵇ[i] = Qᵢ
 
     return Qᵢ
 end
 
+# ## Building the model
+#
+# We assemble the model with the top and bottom heat flux functions:
+
 top_heat_flux    = FluxFunction(sensible_heat_flux; parameters=atmosphere)
 bottom_heat_flux = FluxFunction(advance_lake_and_frazil_flux; parameters=(; lake, atmosphere))
 
 model = SeaIceModel(grid;
                     ice_consolidation_thickness = 0.05, # m
-                    top_heat_flux, 
+                    top_heat_flux,
                     bottom_heat_flux)
 
-# We initialize all the columns with open water (0 thickness and 0 concentration)
+# We initialize all columns with open water (0 thickness and 0 concentration):
 
 set!(model, h=0, ℵ=0)
 
+# ## Running a simulation
+#
+# We set up a simulation that runs for 20 days:
+
 simulation = Simulation(model, Δt=lake.Δt, stop_time=20days)
 
-# The data is accumulated in a timeseries for visualization.
+# ## Collecting data
+#
+# We set up callbacks to accumulate time series data and energy budgets:
 
 timeseries = []
 
@@ -133,7 +160,8 @@ end
 
 simulation.callbacks[:save] = Callback(accumulate_timeseries)
 
-# accumulate energy
+# We also accumulate energy for budget analysis:
+
 Ei = []
 Qa = []
 Ql = []
@@ -142,19 +170,20 @@ function accumulate_energy(sim)
     h  = sim.model.ice_thickness
     ℵ  = sim.model.ice_concentration
     PT = sim.model.ice_thermodynamics.phase_transitions
-    ℰ  = latent_heat(PT, 0) # ice is at 0ᵒC
-
-    push!(Ei, deepcopy(@. - h * ℵ * ℰ))
+    ℰ  = latent_heat(PT, 0) # ice is at 0°C
+    En = - h .* ℵ .* ℰ
+    push!(Ei, deepcopy(En))
     push!(Qa, deepcopy(atmosphere.atmosphere_ice_flux))
     push!(Ql, deepcopy(lake.lake_ice_flux))
 end
 
-simulation.callbacks[:save]   = Callback(accumulate_timeseries)
 simulation.callbacks[:energy] = Callback(accumulate_energy)
 
 run!(simulation)
 
-# Extract and visualize data
+# ## Visualizing the results
+#
+# We extract the time series data for all four columns:
 
 t  = [datum[1]  for datum in timeseries]
 h1 = [datum[2]  for datum in timeseries]
@@ -174,6 +203,9 @@ L2 = [datum[15] for datum in timeseries]
 L3 = [datum[16] for datum in timeseries]
 L4 = [datum[17] for datum in timeseries]
 
+# And visualize the evolution of ice thickness, concentration, surface temperature,
+# and lake temperature:
+
 set_theme!(Theme(fontsize=18, linewidth=3))
 
 fig = Figure(size=(1000, 900))
@@ -183,32 +215,35 @@ axℵ = Axis(fig[2, 1], xlabel="Time (days)", ylabel="Ice concentration (-)")
 axh = Axis(fig[3, 1], xlabel="Time (days)", ylabel="Ice thickness (m)")
 axL = Axis(fig[4, 1], xlabel="Time (days)", ylabel="Lake temperature (ᵒC)")
 
-lines!(axT, t / day, T1)
-lines!(axℵ, t / day, ℵ1)
-lines!(axh, t / day, h1)
-lines!(axL, t / day, L1)
+lines!(axT, t ./ day, T1)
+lines!(axℵ, t ./ day, ℵ1)
+lines!(axh, t ./ day, h1)
+lines!(axL, t ./ day, L1)
 
-lines!(axT, t / day, T2)
-lines!(axℵ, t / day, ℵ2)
-lines!(axh, t / day, h2)
-lines!(axL, t / day, L2)
+lines!(axT, t ./ day, T2)
+lines!(axℵ, t ./ day, ℵ2)
+lines!(axh, t ./ day, h2)
+lines!(axL, t ./ day, L2)
 
-lines!(axT, t / day, T3)
-lines!(axℵ, t / day, ℵ3)
-lines!(axh, t / day, h3)
-lines!(axL, t / day, L3)
+lines!(axT, t ./ day, T3)
+lines!(axℵ, t ./ day, ℵ3)
+lines!(axh, t ./ day, h3)
+lines!(axL, t ./ day, L3)
 
-lines!(axT, t / day, T4)
-lines!(axℵ, t / day, ℵ4)
-lines!(axh, t / day, h4)
-lines!(axL, t / day, L4)
+lines!(axT, t ./ day, T4)
+lines!(axℵ, t ./ day, ℵ4)
+lines!(axh, t ./ day, h4)
+lines!(axL, t ./ day, L4)
 
 save("freezing_in_winter.png", fig)
 nothing # hide
 
 # ![](freezing_in_winter.png)
 
-# Extract and visualize energy
+# ## Energy budget analysis
+#
+# We also visualize the energy budget to verify conservation:
+
 Ei1 = [datum[1]  for datum in Ei]
 Qa1 = [datum[1]  for datum in Qa]
 Ql1 = [datum[1]  for datum in Ql]
@@ -224,36 +259,37 @@ Ql4 = [datum[4]  for datum in Ql]
 
 fig = Figure(size=(1000, 900))
 
-axE = Axis(fig[1, 1], xlabel="Time (days)", ylabel="Sea Ice energy (J)")
-axA = Axis(fig[2, 1], xlabel="Time (days)", ylabel="Atmosphere HF")
-axL = Axis(fig[3, 1], xlabel="Time (days)", ylabel="Lake HF")
-axB = Axis(fig[4, 1], xlabel="Time (days)", ylabel="Heat budget")
+axE = Axis(fig[1, 1], xlabel="Time (days)", ylabel="Sea ice energy (J)")
+axA = Axis(fig[2, 1], xlabel="Time (days)", ylabel="Atmosphere heat flux (W m⁻²)")
+axL = Axis(fig[3, 1], xlabel="Time (days)", ylabel="Lake heat flux (W m⁻²)")
+axB = Axis(fig[4, 1], xlabel="Time (days)", ylabel="Heat budget residual (W m⁻²)")
 
 dEi1 = (Ei1[2:end] - Ei1[1:end-1]) ./ 10minutes
 dEi2 = (Ei2[2:end] - Ei2[1:end-1]) ./ 10minutes
 dEi3 = (Ei3[2:end] - Ei3[1:end-1]) ./ 10minutes
 dEi4 = (Ei4[2:end] - Ei4[1:end-1]) ./ 10minutes
-tpE  = t[2:end] 
+tpE  = t[2:end]
 
-lines!(axE, t / day, Ei1)
-lines!(axA, t / day, Qa1)
-lines!(axL, t / day, Ql1)
+lines!(axE, t ./ day, Ei1)
+lines!(axA, t ./ day, Qa1)
+lines!(axL, t ./ day, Ql1)
 
-lines!(axE, t / day, Ei2)
-lines!(axA, t / day, Qa2)
-lines!(axL, t / day, Ql2)
+lines!(axE, t ./ day, Ei2)
+lines!(axA, t ./ day, Qa2)
+lines!(axL, t ./ day, Ql2)
 
-lines!(axE, t / day, Ei3)
-lines!(axA, t / day, Qa3)
-lines!(axL, t / day, Ql3)
+lines!(axE, t ./ day, Ei3)
+lines!(axA, t ./ day, Qa3)
+lines!(axL, t ./ day, Ql3)
 
-lines!(axE, t / day, Ei4)
-lines!(axA, t / day, Qa4)
-lines!(axL, t / day, Ql4)
+lines!(axE, t ./ day, Ei4)
+lines!(axA, t ./ day, Qa4)
+lines!(axL, t ./ day, Ql4)
 
-lines!(axB, tpE / day, dEi1 .- (.- Qa1 .+ Ql1)[2:end])
-lines!(axB, tpE / day, dEi2 .- (.- Qa2 .+ Ql2)[2:end])
-lines!(axB, tpE / day, dEi3 .- (.- Qa3 .+ Ql3)[2:end])
+# The budget residual should be close to zero
+lines!(axB, tpE ./ day, dEi1 .- (.- Qa1 .+ Ql1)[2:end])
+lines!(axB, tpE ./ day, dEi2 .- (.- Qa2 .+ Ql2)[2:end])
+lines!(axB, tpE ./ day, dEi3 .- (.- Qa3 .+ Ql3)[2:end])
 
 save("energy_budget.png", fig)
 nothing # hide
diff --git a/examples/ice_advected_by_anticyclone.jl b/examples/ice_advected_by_anticyclone.jl
index a6ba4559..9f80dbc5 100644
--- a/examples/ice_advected_by_anticyclone.jl
+++ b/examples/ice_advected_by_anticyclone.jl
@@ -1,5 +1,32 @@
-# # Sea ice advected by an atmospheric anticyclone
+# # Sea ice advected by an atmospheric anticyclone example
 #
+# This example simulates sea ice advected by an atmospheric anticyclone, based on
+# the experiment found in the paper:
+#
+# > Simulating Linear Kinematic Features in Viscous-Plastic Sea Ice Models on
+# > Quadrilateral and Triangular Grids With Different Variable Staggering
+#
+# The ice is driven by a moving anticyclonic atmospheric eddy and interacts with
+# a cyclonic ocean eddy. This example demonstrates how to:
+#
+#   * set up a two-dimensional sea ice model with time-varying atmospheric forcing,
+#   * prescribe moving atmospheric and ocean eddies,
+#   * use elasto-visco-plastic rheology with Coriolis effects,
+#   * save and visualize time series data.
+#
+# ## Install dependencies
+#
+# First let's make sure we have all required packages installed.
+#
+# ```julia
+# using Pkg
+# pkg"add Oceananigans, ClimaSeaIce, CairoMakie"
+# ```
+#
+# ## The physical domain
+#
+# We set up a two-dimensional bounded domain:
+
 using ClimaSeaIce
 using Oceananigans
 using Oceananigans.Units
@@ -8,17 +35,12 @@ using Oceananigans.BoundaryConditions
 using Printf
 using CairoMakie
 
-# The experiment found in the paper:
-# Simulating Linear Kinematic Features in Viscous-Plastic Sea Ice Models
-# on Quadrilateral and Triangular Grids With Different Variable Staggering
-
 arch = CPU()
 
 L  = 512kilometers
-𝓋ₒ = 0.01 # m / s maximum ocean speed
-𝓋ₐ = 30.0 # m / s maximum atmospheric speed modifier
+𝓋ₒ = 0.01 # m s⁻¹ maximum ocean speed
+𝓋ₐ = 30.0 # m s⁻¹ maximum atmospheric speed modifier
 
-# 2 km domain
 grid = RectilinearGrid(arch;
                        size = (128, 128),
                           x = (0, L),
@@ -26,9 +48,9 @@ grid = RectilinearGrid(arch;
                        halo = (7, 7),
                    topology = (Bounded, Bounded, Flat))
 
-#####
-##### Value boundary conditions for velocities
-#####
+# ## Boundary conditions
+#
+# We set value boundary conditions for velocities:
 
 u_bcs = FieldBoundaryConditions(north=ValueBoundaryCondition(0),
                                 south=ValueBoundaryCondition(0))
@@ -36,11 +58,10 @@ u_bcs = FieldBoundaryConditions(north=ValueBoundaryCondition(0),
 v_bcs = FieldBoundaryConditions(west=ValueBoundaryCondition(0),
                                 east=ValueBoundaryCondition(0))
 
-#####
-##### Ocean sea-ice stress
-#####
+# ## Ocean sea-ice stress
+#
+# We set up constant ocean velocities corresponding to a cyclonic eddy:
 
-# Constant ocean velocities corresponding to a cyclonic eddy
 Uₒ = XFaceField(grid)
 Vₒ = YFaceField(grid)
 
@@ -50,22 +71,23 @@ fill_halo_regions!((Uₒ, Vₒ))
 
 τₒ = SemiImplicitStress(uₑ=Uₒ, vₑ=Vₒ)
 
-####
-#### Atmosphere - sea ice stress
-####
+# ## Atmosphere-sea ice stress
+#
+# We set up atmospheric velocities corresponding to an anticyclonic eddy moving
+# northeast. The eddy center moves with time:
 
-Uₐ = XFaceField(grid)
-Vₐ = YFaceField(grid)
+@inline center(t) = 256kilometers + 51.2kilometers * t / day
+@inline radius(x, y, t) = sqrt((x - center(t))^2 + (y - center(t))^2)
+@inline speed(x, y, t) = 1 / 100 * exp(- radius(x, y, t) / 100kilometers)
 
-# Atmospheric velocities corresponding to an anticyclonic eddy moving north-east
-@inline center(t) = 256kilometers + 51.2kilometers * t / 86400
-@inline radius(x, y, t)  = sqrt((x - center(t))^2 + (y - center(t))^2)
-@inline speed(x, y, t)   = 1 / 100 * exp(- radius(x, y, t) / 100kilometers)
+@inline ua_time(x, y, t) = - 𝓋ₐ * speed(x, y, t) * (  cosd(72) * (x - center(t)) + sind(72) * (y - center(t))) / 1000
+@inline va_time(x, y, t) = - 𝓋ₐ * speed(x, y, t) * (- sind(72) * (x - center(t)) + cosd(72) * (y - center(t))) / 1000
 
-@inline ua_time(x, y, t) = - 𝓋ₐ * speed(x, y, t) * (cosd(72) * (x - center(t)) + sind(72) * (y - center(t))) / 1000
-@inline va_time(x, y, t) = - 𝓋ₐ * speed(x, y, t) * (cosd(72) * (y - center(t)) - sind(72) * (x - center(t))) / 1000
+# Initialize the stress at time t = 0:
+
+Uₐ = XFaceField(grid)
+Vₐ = YFaceField(grid)
 
-# Initialize the stress at time t = 0
 set!(Uₐ, (x, y) -> ua_time(x, y, 0))
 set!(Vₐ, (x, y) -> va_time(x, y, 0))
 
@@ -74,12 +96,10 @@ fill_halo_regions!((Uₐ, Vₐ))
 τₐu = Field(- Uₐ * sqrt(Uₐ^2 + Vₐ^2) * 1.3 * 1.2e-3)
 τₐv = Field(- Vₐ * sqrt(Uₐ^2 + Vₐ^2) * 1.3 * 1.2e-3)
 
-#####
-##### Numerical details
-#####
-
-# We use an elasto-visco-plastic rheology and WENO seventh order
-# for advection of h and ℵ
+# ## Model configuration
+#
+# We use an elasto-visco-plastic rheology and WENO seventh order for advection
+# of ice thickness and concentration. We include Coriolis effects:
 
 dynamics = SeaIceMomentumEquation(grid;
                                   top_momentum_stress = (u=τₐu, v=τₐv),
@@ -88,28 +108,30 @@ dynamics = SeaIceMomentumEquation(grid;
 
 model = SeaIceModel(grid;
                     dynamics,
-                    ice_thermodynamics = nothing, # No ice_thermodynamics here
-                    timestepper = :SplitRungeKutta3,
+                    ice_thermodynamics = nothing, # No thermodynamics here
                     advection = WENO(order=7),
                     boundary_conditions = (u=u_bcs, v=v_bcs))
 
-# We start with a concentration of ℵ = 1 and an
-# initial height field with perturbations around 0.3 m
+# ## Initial conditions
+#
+# We start with a concentration of `ℵ = 1` and an initial height field with
+# perturbations around 0.3 meters:
 
 h₀(x, y) = 0.3 + 0.005 * (sin(60 * x / 1000kilometers) + sin(30 * y / 1000kilometers))
 
 set!(model, h = h₀)
 set!(model, ℵ = 1)
 
-#####
-##### Setup the simulation
-#####
-
-# run the model for 2 days
+# ## Running the simulation
+#
+# We run the model for 2 days with a 2-minute time step:
 
 simulation = Simulation(model, Δt = 2minutes, stop_time = 2days)
 
-# Remember to evolve the wind stress field in time!
+# ## Time-varying wind stress
+#
+# We need to evolve the wind stress field in time. We set up a callback to
+# update the atmospheric stress at each time step:
 
 function compute_wind_stress(sim)
     time = sim.model.clock.time
@@ -128,6 +150,11 @@ end
 
 simulation.callbacks[:top_stress] = Callback(compute_wind_stress, IterationInterval(1))
 
+# ## Output writer
+#
+# We set up an output writer to save the ice thickness, concentration, and
+# velocity fields:
+
 h = model.ice_thickness
 ℵ = model.ice_concentration
 u, v = model.velocities
@@ -136,44 +163,23 @@ outputs = (; h, u, v, ℵ)
 
 simulation.output_writers[:sea_ice] = JLD2Writer(model, outputs;
                                                  filename = "sea_ice_advected_by_anticyclone.jld2",
-                                                 including = [:grid],
                                                  schedule = IterationInterval(5),
                                                  overwrite_existing = true)
 
-wall_time = [time_ns()]
-
-function progress(sim)
-    h = sim.model.ice_thickness
-    ℵ = sim.model.ice_concentration
-    u = sim.model.velocities.u
-    v = sim.model.velocities.v
-
-    hmax = maximum(interior(h))
-    ℵmin = minimum(interior(ℵ))
-    umax = maximum(interior(u)), maximum(interior(v))
-    step_time = 1e-9 * (time_ns() - wall_time[1])
-
-    @info @sprintf("Time: %s, Iteration %d, Δt %s, max(vel): (%.2e, %.2e), max(h): %.2f, min(ℵ): %.2f, wtime: %s \n",
-                   prettytime(sim.model.clock.time),
-                   sim.model.clock.iteration,
-                   prettytime(sim.Δt),
-                   umax..., hmax, ℵmin, prettytime(step_time))
 
-     wall_time[1] = time_ns()
-end
-
-simulation.callbacks[:progress] = Callback(progress, IterationInterval(100))
+# We can finally run the simulation
 
 run!(simulation)
 
-# Load output
+# ## Visualizing the results
+#
+# We load the saved time series and create an animation:
 
 htimeseries = FieldTimeSeries("sea_ice_advected_by_anticyclone.jld2", "h")
 utimeseries = FieldTimeSeries("sea_ice_advected_by_anticyclone.jld2", "u")
 vtimeseries = FieldTimeSeries("sea_ice_advected_by_anticyclone.jld2", "v")
 ℵtimeseries = FieldTimeSeries("sea_ice_advected_by_anticyclone.jld2", "ℵ")
 
-# Visualize!
 Nt = length(htimeseries)
 iter = Observable(1)
 
@@ -182,20 +188,25 @@ hi = @lift(interior(htimeseries[$iter], :, :, 1))
 ui = @lift(interior(utimeseries[$iter], :, :, 1))
 vi = @lift(interior(vtimeseries[$iter], :, :, 1))
 
-fig = Figure()
-ax1 = Axis(fig[1, 1], title = "sea ice thickness")
-ax2 = Axis(fig[1, 2], title = "sea ice concentration")
-ax3 = Axis(fig[2, 1], title = "zonal velocity")
-ax4 = Axis(fig[2, 2], title = "meridional velocity")
+fig = Figure(size = (1200, 1200))
+ax1 = Axis(fig[1, 1], title = "Sea ice thickness (m)")
+ax2 = Axis(fig[1, 2], title = "Sea ice concentration")
+ax3 = Axis(fig[2, 1], title = "Zonal velocity (m s⁻¹)")
+ax4 = Axis(fig[2, 2], title = "Meridional velocity (m s⁻¹)")
 
 heatmap!(ax1, hi, colormap = :magma, colorrange = (0.23, 0.37))
 heatmap!(ax2, ℵi, colormap = Reverse(:deep), colorrange = (0.9, 1))
-heatmap!(ax3, ui, colorrange = (-0.1, 0.1))
-heatmap!(ax4, vi, colorrange = (-0.1, 0.1))
+heatmap!(ax3, ui, colorrange = (-0.1, 0.1), colormap = :balance)
+heatmap!(ax4, vi, colorrange = (-0.1, 0.1), colormap = :balance)
 
 CairoMakie.record(fig, "sea_ice_advected_by_anticyclone.mp4", 1:Nt, framerate = 8) do i
     iter[] = i
+    @info "Rendering frame $i"
 end
 nothing #hide
 
 # ![](sea_ice_advected_by_anticyclone.mp4)
+#
+# The animation shows how the ice is advected and deformed by the moving
+# anticyclonic atmospheric eddy. Linear kinematic features (leads and ridges)
+# form as the ice responds to the complex stress patterns.
diff --git a/examples/ice_advected_on_coastline.jl b/examples/ice_advected_on_coastline.jl
index 3efb0b67..4047db80 100644
--- a/examples/ice_advected_on_coastline.jl
+++ b/examples/ice_advected_on_coastline.jl
@@ -1,3 +1,27 @@
+# # Ice advected on coastline example
+#
+# This example simulates a solid block of ice moving against a triangular coastline
+# in a periodic channel. The ice is driven by atmospheric winds and interacts with
+# the coastline through immersed boundaries. This example demonstrates how to:
+#
+#   * set up a two-dimensional sea ice model with immersed boundaries,
+#   * prescribe atmospheric and oceanic stresses,
+#   * use elasto-visco-plastic rheology with split-explicit time stepping,
+#   * visualize the evolution of ice thickness, concentration, and velocity.
+#
+# ## Install dependencies
+#
+# First let's make sure we have all required packages installed.
+#
+# ```julia
+# using Pkg
+# pkg"add Oceananigans, ClimaSeaIce, CairoMakie"
+# ```
+#
+# ## The physical domain
+#
+# We set up a two-dimensional periodic channel with a triangular coastline:
+
 using ClimaSeaIce
 using Oceananigans
 using Oceananigans.Units
@@ -5,8 +29,6 @@ using Oceananigans.Operators
 using Printf
 using CairoMakie
 
-# A solid block of ice moving against a triangular coastline in a periodic channel
-
 Lx = 512kilometers
 Ly = 256kilometers
 Nx = 256
@@ -16,27 +38,33 @@ y_max = Ly / 2
 
 arch = CPU()
 
-𝓋ₐ = 10.0   # m / s 
-Cᴰ = 1.2e-3 # Atmosphere - sea ice drag coefficient
-ρₐ = 1.3    # kg/m³
+# ## Grid configuration
+#
+# We create a rectilinear grid with periodic boundaries in x and bounded boundaries in y:
 
-# 2 km domain
-grid = RectilinearGrid(arch; size = (Nx, Ny), 
-                                x = (-Lx/2, Lx/2), 
-                                y = (0, Ly), 
+grid = RectilinearGrid(arch; size = (Nx, Ny),
+                                x = (-Lx/2, Lx/2),
+                                y = (0, Ly),
                              halo = (4, 4),
                          topology = (Periodic, Bounded, Flat))
 
-bottom(x, y) = ifelse(y > y_max, 0, 
+# We define a triangular coastline using an immersed boundary:
+
+bottom(x, y) = ifelse(y > y_max, 0,
                ifelse(abs(x / Lx) * Nx + y / Ly * Ny > 24, 0, 1))
 
 grid = ImmersedBoundaryGrid(grid, GridFittedBoundary(bottom))
 
-#####
-##### Setup atmospheric and oceanic forcing
-#####
+# ## Atmospheric and oceanic forcing
+#
+# We set up atmospheric wind stress. The atmosphere has a constant wind speed:
+
+𝓋ₐ = 10.0   # m s⁻¹
+Cᴰ = 1.2e-3 # atmosphere-sea ice drag coefficient
+ρₐ = 1.3    # kg m⁻³
+
+# We create a field for the atmospheric wind and compute the stress:
 
-# Atmosphere - sea ice stress
 Ua = XFaceField(grid)
 τᵤ = Field(- ρₐ * Cᴰ * Ua^2)
 set!(Ua, 𝓋ₐ)
@@ -44,53 +72,55 @@ compute!(τᵤ)
 Oceananigans.BoundaryConditions.fill_halo_regions!(τᵤ)
 τᵥ = 0.0
 
-#####
-##### Ocean stress (a zero-velocity ocean with a drag)
-#####
+# The ocean stress is represented by a semi-implicit stress with zero ocean velocity:
 
 τₒ = SemiImplicitStress()
 
-#####
-##### Numerical details
-#####
+# ## Model configuration
+#
+# We use an elasto-visco-plastic rheology and WENO seventh order for advection
+# of ice thickness and concentration:
 
-# We use an elasto-visco-plastic rheology and WENO seventh order 
-# for advection of h and ℵ
-dynamics = SeaIceMomentumEquation(grid; 
+dynamics = SeaIceMomentumEquation(grid;
                                   top_momentum_stress = (u=τᵤ, v=τᵥ),
-                                  bottom_momentum_stress = τₒ, 
+                                  bottom_momentum_stress = τₒ,
                                   rheology = ElastoViscoPlasticRheology(),
                                   solver = SplitExplicitSolver(substeps=150))
 
+# We set boundary conditions for the velocity field:
+
 u_bcs = FieldBoundaryConditions(top = nothing, bottom = nothing,
                                 north = ValueBoundaryCondition(0),
                                 south = ValueBoundaryCondition(0))
 
-#Define the model! 
-model = SeaIceModel(grid; 
+# We define the model with WENO advection and no thermodynamics:
+
+model = SeaIceModel(grid;
                     advection = WENO(order=7),
                     dynamics = dynamics,
                     boundary_conditions = (; u=u_bcs),
                     ice_thermodynamics = nothing)
 
-# We start with a concentration of ℵ = 1 everywhere
+# We initialize the model with uniform ice thickness and concentration:
+
 set!(model, h = 1)
 set!(model, ℵ = 1)
 
-#####
-##### Setup the simulation
-#####
+# ## Running the simulation
+#
+# We run the model for 2 days with a 2-minute time step:
 
-# run the model for 10 day
-simulation = Simulation(model, Δt = 2minutes, stop_time=2days) 
+simulation = Simulation(model, Δt = 2minutes, stop_time=2days)
+
+# ## Collecting data
+#
+# We set up containers to hold the time series data:
 
-# Container to hold the data
 htimeseries = []
 ℵtimeseries = []
 utimeseries = []
 vtimeseries = []
 
-## Callback function to collect the data from the `sim`ulation
 function accumulate_timeseries(sim)
     h = sim.model.ice_thickness
     ℵ = sim.model.ice_concentration
@@ -102,34 +132,14 @@ function accumulate_timeseries(sim)
     push!(vtimeseries, deepcopy(Array(interior(v))))
 end
 
-wall_time = Ref(time_ns())
-
-function progress(sim) 
-    h = sim.model.ice_thickness
-    ℵ = sim.model.ice_concentration
-    u = sim.model.velocities.u
-    v = sim.model.velocities.v
-
-    hmax = maximum(interior(h))
-    ℵmin = minimum(interior(ℵ))
-    umax = maximum(interior(u)), maximum(interior(v))
-    step_time = 1e-9 * (time_ns() - wall_time[])
-
-    @info @sprintf("Time: %s, Iteration %d, Δt %s, max(vel): (%.2e, %.2e), max(trac): %.2f, %.2f, wtime: %s \n",
-                   prettytime(sim.model.clock.time),
-                   sim.model.clock.iteration,
-                   prettytime(sim.Δt),
-                   umax..., hmax, ℵmin, prettytime(step_time))
-
-     wall_time[] = time_ns()
-end
-
-simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))
 simulation.callbacks[:save]     = Callback(accumulate_timeseries, IterationInterval(10))
 
 run!(simulation)
 
-# Visualize!
+# ## Visualizing the results
+#
+# We create an animation of the ice thickness, concentration, and velocity fields:
+
 Nt = length(htimeseries)
 iter = Observable(1)
 
@@ -139,22 +149,26 @@ ui = @lift(utimeseries[$iter][:, :, 1])
 vi = @lift(vtimeseries[$iter][:, :, 1])
 
 fig = Figure(size = (1000, 500))
-ax = Axis(fig[1, 1], title = "sea ice thickness")
-heatmap!(ax, hi, colormap = :magma,         colorrange = (0.0, 2.0))
+ax = Axis(fig[1, 1], title = "Sea ice thickness (m)")
+heatmap!(ax, hi, colormap = :magma, colorrange = (0.0, 2.0))
 
-ax = Axis(fig[1, 2], title = "sea ice concentration")
+ax = Axis(fig[1, 2], title = "Sea ice concentration")
 heatmap!(ax, ℵi, colormap = Reverse(:deep), colorrange = (0.0, 1))
 
-ax = Axis(fig[2, 1], title = "zonal velocity")
+ax = Axis(fig[2, 1], title = "Zonal velocity (m s⁻¹)")
 heatmap!(ax, ui, colorrange = (0, 0.12), colormap = :balance)
 
-ax = Axis(fig[2, 2], title = "meridional velocity")
+ax = Axis(fig[2, 2], title = "Meridional velocity (m s⁻¹)")
 heatmap!(ax, vi, colorrange = (-0.025, 0.025), colormap = :bwr)
 
 CairoMakie.record(fig, "sea_ice_advected_on_coastline.mp4", 1:Nt, framerate = 8) do i
     iter[] = i
-    @info "doing iter $i"
+    @info "Rendering frame $i"
 end
 nothing #hide
 
-# ![](sea_ice_advected_on_coastline.mp4)
\ No newline at end of file
+# ![](sea_ice_advected_on_coastline.mp4)
+#
+# The animation shows how the ice block moves and deforms as it interacts with
+# the triangular coastline. The ice accumulates and thickens near the coast due
+# to the no-slip boundary conditions.
diff --git a/examples/melting_in_spring.jl b/examples/melting_in_spring.jl
index 8f3c7673..c58c689d 100644
--- a/examples/melting_in_spring.jl
+++ b/examples/melting_in_spring.jl
@@ -1,41 +1,65 @@
-# # Melting in Spring
+# # Melting in spring example
 #
-# A simulation that mimicks the melting of (relatively thick) sea ice in the spring
-# when the sun is shining. The ice is subject to solar insolation and sensible heat
-# fluxes from the atmosphere. Different cells show how the ice melts at different rates
-# depending on the amount of solar insolation they receive.
+# This example simulates the melting of relatively thick sea ice in spring when
+# the sun is shining. The ice is subject to solar insolation and sensible heat
+# fluxes from the atmosphere. Different grid cells show how the ice melts at
+# different rates depending on the amount of solar insolation they receive.
+# This example demonstrates how to:
 #
-# We start by `using Oceananigans` to bring in functions for building grids
-# and `Simulation`s and the like.
+#   * set up a one-dimensional model with multiple grid cells,
+#   * prescribe spatially varying solar insolation,
+#   * use `FluxFunction` for parameterized heat fluxes,
+#   * visualize the evolution of multiple ice columns.
+#
+# ## Install dependencies
+#
+# First let's make sure we have all required packages installed.
+#
+# ```julia
+# using Pkg
+# pkg"add Oceananigans, ClimaSeaIce, CairoMakie"
+# ```
+#
+# ## The physical domain
+#
+# We generate a one-dimensional grid with 4 grid cells to model different ice
+# columns subject to different solar insolation:
 
 using Oceananigans
 using Oceananigans.Units
 using ClimaSeaIce
+using ClimaSeaIce.SeaIceThermodynamics.HeatBoundaryConditions: RadiativeEmission
 using CairoMakie
 
-# Generate a 1D grid for difference ice columns subject to different solar insolation
-
 grid = RectilinearGrid(size=4, x=(0, 1), topology=(Periodic, Flat, Flat))
 
-# Build a model of an ice slab that has internal conductive fluxes
-# and that emits radiation from its top surface.
+# ## Top boundary conditions
+#
+# We prescribe different solar insolation values for each grid cell, ranging
+# from -600 to -1200 W m⁻² (negative values indicate downward/into the ice):
 
 solar_insolation = [-600, -800, -1000, -1200] # W m⁻²
 solar_insolation = reshape(solar_insolation, (4, 1, 1))
+
+# The ice also emits longwave radiation from its top surface:
+
 outgoing_radiation = RadiativeEmission()
 
+# ## Sensible heat flux parameterization
+#
 # The sensible heat flux from the atmosphere is represented by a `FluxFunction`.
+# We define the parameters for the bulk formula:
 
 parameters = (
-    transfer_coefficient     = 1e-3,  # Unitless
+    transfer_coefficient     = 1e-3,  # unitless
     atmosphere_density       = 1.225, # kg m⁻³
-    atmosphere_heat_capacity = 1004,  # 
-    atmosphere_temperature   = -5,    # ᵒC
+    atmosphere_heat_capacity = 1004,  # J kg⁻¹ K⁻¹
+    atmosphere_temperature   = -5,    # °C
     atmosphere_wind_speed    = 5      # m s⁻¹
 )
 
-# Flux is positive (cooling by fluxing heat up away from upper surface)
-# when Tₐ < Tᵤ:
+# The flux is positive (cooling by fluxing heat upward away from the upper surface)
+# when the atmosphere temperature is less than the surface temperature:
 
 @inline function sensible_heat_flux(i, j, grid, Tᵤ, clock, fields, parameters)
     Cₛ = parameters.transfer_coefficient
@@ -45,23 +69,36 @@ parameters = (
     uₐ = parameters.atmosphere_wind_speed
     ℵ  = fields.ℵ[i, j, 1]
 
-    return Cₛ * ρₐ * cₐ * uₐ * (Tᵤ - Tₐ) * ℵ 
+    return Cₛ * ρₐ * cₐ * uₐ * (Tᵤ - Tₐ) * ℵ
 end
 
 aerodynamic_flux = FluxFunction(sensible_heat_flux; parameters)
+
+# We combine all top heat fluxes into a tuple:
+
 top_heat_flux = (outgoing_radiation, solar_insolation, aerodynamic_flux)
 
+# ## Building the model
+#
+# We assemble the model with an ice consolidation thickness of 0.05 m:
+
 model = SeaIceModel(grid;
                     ice_consolidation_thickness = 0.05, # m
                     top_heat_flux)
 
-# We initialize all the columns with a 1 m thick slab of ice with 100% ice concentration.
-            
+# We initialize all columns with a 1 m thick slab of ice with 100% ice concentration:
+
 set!(model, h=1, ℵ=1)
 
+# ## Running a simulation
+#
+# We set up a simulation that runs for 30 days with a 10-minute time step:
+
 simulation = Simulation(model, Δt=10minute, stop_time=30days)
 
-# The data is accumulated in a timeseries for visualization.
+# ## Collecting data
+#
+# We set up a callback to accumulate time series data for all four columns:
 
 timeseries = []
 
@@ -80,7 +117,9 @@ simulation.callbacks[:save] = Callback(accumulate_timeseries)
 
 run!(simulation)
 
-# Extract and visualize data
+# ## Visualizing the results
+#
+# We extract the time series data for all four columns:
 
 t  = [datum[1]  for datum in timeseries]
 h1 = [datum[2]  for datum in timeseries]
@@ -96,6 +135,8 @@ h4 = [datum[11] for datum in timeseries]
 ℵ4 = [datum[12] for datum in timeseries]
 T4 = [datum[13] for datum in timeseries]
 
+# And visualize the evolution of ice thickness, concentration, and surface temperature:
+
 set_theme!(Theme(fontsize=18, linewidth=3))
 
 fig = Figure(size=(1000, 800))
@@ -104,24 +145,27 @@ axT = Axis(fig[1, 1], xlabel="Time (days)", ylabel="Top temperature (ᵒC)")
 axℵ = Axis(fig[2, 1], xlabel="Time (days)", ylabel="Ice concentration (-)")
 axh = Axis(fig[3, 1], xlabel="Time (days)", ylabel="Ice thickness (m)")
 
-lines!(axT, t / day, T1)
-lines!(axℵ, t / day, ℵ1)
-lines!(axh, t / day, h1)
+lines!(axT, t ./ day, T1)
+lines!(axℵ, t ./ day, ℵ1)
+lines!(axh, t ./ day, h1)
 
-lines!(axT, t / day, T2)
-lines!(axℵ, t / day, ℵ2)
-lines!(axh, t / day, h2)
+lines!(axT, t ./ day, T2)
+lines!(axℵ, t ./ day, ℵ2)
+lines!(axh, t ./ day, h2)
 
-lines!(axT, t / day, T3)
-lines!(axℵ, t / day, ℵ3)
-lines!(axh, t / day, h3)
+lines!(axT, t ./ day, T3)
+lines!(axℵ, t ./ day, ℵ3)
+lines!(axh, t ./ day, h3)
 
-lines!(axT, t / day, T4)
-lines!(axℵ, t / day, ℵ4)
-lines!(axh, t / day, h4)
+lines!(axT, t ./ day, T4)
+lines!(axℵ, t ./ day, ℵ4)
+lines!(axh, t ./ day, h4)
 
 save("melting_in_spring.png", fig)
 nothing # hide
 
 # ![](melting_in_spring.png)
-
+#
+# The results show that ice columns receiving more solar insolation melt faster,
+# as expected. The ice concentration and thickness decrease more rapidly for
+# columns with higher incoming solar radiation.
diff --git a/examples/perpetual_night.jl b/examples/perpetual_night.jl
index 6c004670..db4bfb98 100644
--- a/examples/perpetual_night.jl
+++ b/examples/perpetual_night.jl
@@ -1,33 +1,79 @@
+# # Perpetual night example
+#
+# In this example, we simulate sea ice under perpetual night conditions, where
+# the ice surface emits longwave radiation but receives no incoming solar radiation.
+# This example demonstrates how to:
+#
+#   * use `RadiativeEmission` for outgoing longwave radiation,
+#   * set initial ice thickness,
+#   * collect and visualize time series data.
+#
+# ## Install dependencies
+#
+# First let's make sure we have all required packages installed.
+#
+# ```julia
+# using Pkg
+# pkg"add Oceananigans, ClimaSeaIce, CairoMakie"
+# ```
+#
+# ## The physical domain
+#
+# We use a single grid point to model a horizontally uniform ice slab:
+
 using Oceananigans
 using Oceananigans.Units
 using ClimaSeaIce
-using ClimaSeaIce.HeatBoundaryConditions: RadiativeEmission
+using ClimaSeaIce.SeaIceThermodynamics.HeatBoundaryConditions: RadiativeEmission
 using CairoMakie
 
-# Generate a zero-dimensional grid for a single column slab model 
 grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat))
 
-# Build a model of an ice slab that has internal conductive fluxes
-# and that emits radiation from its top surface.
-model = SlabSeaIceModel(grid; top_heat_flux=RadiativeEmission())
-set!(model, h=0.01)
+# ## Model configuration
+#
+# We build a model of an ice slab that has internal conductive heat fluxes
+# and emits longwave radiation from its top surface. The `RadiativeEmission`
+# boundary condition implements the Stefan-Boltzmann law for blackbody radiation:
+
+ice_thermodynamics = SlabSeaIceThermodynamics(grid; top_heat_boundary_condition=MeltingConstrainedFluxBalance())
+
+top_flux = Field{Nothing, Nothing, Nothing}(grid)
+interior(top_flux) .= - 200.0
+model = SeaIceModel(grid; top_heat_flux=(RadiativeEmission(), top_flux), ice_thermodynamics)
+
+# We initialize the ice with a small thickness:
+
+set!(model, h=0.01) # m
+
+# ## Running a simulation
+#
+# We set up a simulation that runs for 40 days with a 1-hour time step:
 
 simulation = Simulation(model, Δt=1hour, stop_time=40days)
 
-# Accumulate data
+# ## Collecting data
+#
+# Before running the simulation, we set up a callback to accumulate time series
+# data of the ice thickness and top surface temperature:
+
 timeseries = []
 
 function accumulate_timeseries(sim)
-    T = model.top_surface_temperature
+    T = model.ice_thermodynamics.top_surface_temperature
     h = model.ice_thickness
     push!(timeseries, (time(sim), first(h), first(T)))
 end
 
 simulation.callbacks[:save] = Callback(accumulate_timeseries)
 
+# Now we run the simulation:
+
 run!(simulation)
 
-# Extract and visualize data
+# ## Visualizing the results
+#
+# We extract the time series data and visualize the evolution of ice thickness
+# and top surface temperature:
 
 t = [datum[1] for datum in timeseries]
 h = [datum[2] for datum in timeseries]
@@ -40,8 +86,10 @@ fig = Figure(size=(1000, 800))
 axT = Axis(fig[1, 1], xlabel="Time (days)", ylabel="Top temperature (ᵒC)")
 axh = Axis(fig[2, 1], xlabel="Time (days)", ylabel="Ice thickness (m)")
 
-lines!(axT, t / day, T)
-lines!(axh, t / day, h)
+lines!(axT, t ./ day, T)
+lines!(axh, t ./ day, h)
 
-display(fig)
+current_figure() #hide
 
+# Under perpetual night conditions, the ice will cool and grow as it loses
+# heat through radiative emission at the surface.
diff --git a/examples/simple_freezing_bucket.jl b/examples/simple_freezing_bucket.jl
index 249f2b9b..9a994c6f 100644
--- a/examples/simple_freezing_bucket.jl
+++ b/examples/simple_freezing_bucket.jl
@@ -1,14 +1,60 @@
+# # Simple freezing bucket example
+#
+# This is ClimaSeaIce.jl's simplest example: a single grid point model of sea ice
+# freezing in a bucket. This example demonstrates how to:
+#
+#   * load `ClimaSeaIce.jl`,
+#   * instantiate a `SlabSeaIceModel` with prescribed boundary conditions,
+#   * configure internal heat conduction.
+#
+# ## Install dependencies
+#
+# First let's make sure we have all required packages installed.
+#
+# ```julia
+# using Pkg
+# pkg"add Oceananigans, ClimaSeaIce"
+# ```
+#
+# ## Using `ClimaSeaIce.jl`
+#
+# Write
+
 using Oceananigans
 using Oceananigans.Units
 using ClimaSeaIce
 
-top_temperature = -10 # C
-conductivity = 2 # W m² K⁻¹ .. ?
+# to load ClimaSeaIce functions and objects into our script.
+#
+# ## A single grid point model
+#
+# For a simple freezing bucket, we only need a single grid point since we're
+# modeling a horizontally uniform ice slab. We use a `Flat` topology in all
+# directions to create a zero-dimensional grid:
 
 grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat))
 
+# ## Model configuration
+#
+# We configure the model with a prescribed top temperature and internal heat
+# conduction. The top temperature is kept constant at -10°C, simulating a cold
+# metal lid on top of the bucket:
+
+top_temperature = -10 # °C
 top_heat_boundary_condition = PrescribedTemperature(top_temperature)
+
+# We specify internal heat conduction through the ice with a thermal conductivity
+# of 2 W m⁻¹ K⁻¹ (typical for sea ice):
+
+conductivity = 2 # W m⁻¹ K⁻¹
 internal_heat_flux = ConductiveFlux(; conductivity)
 
+# ## Building the model
+#
+# We assemble these components into a `SlabSeaIceModel`:
+
 model = SlabSeaIceModel(grid; top_heat_boundary_condition, internal_heat_flux)
 
+# The model is now ready to use! By default, the bottom boundary condition
+# is `IceWaterThermalEquilibrium`, which maintains the ice-water interface
+# at the melting temperature.
diff --git a/src/SeaIceThermodynamics/HeatBoundaryConditions/boundary_fluxes.jl b/src/SeaIceThermodynamics/HeatBoundaryConditions/boundary_fluxes.jl
index 6e6c5c7f..15871b26 100644
--- a/src/SeaIceThermodynamics/HeatBoundaryConditions/boundary_fluxes.jl
+++ b/src/SeaIceThermodynamics/HeatBoundaryConditions/boundary_fluxes.jl
@@ -1,4 +1,3 @@
-
 #####
 ##### Fluxes
 #####
diff --git a/validation/ice_ocean_model/cooling_then_warming_ocean.jl b/validation/ice_ocean_model/cooling_then_warming_ocean.jl
index cc7701f3..eeda503b 100644
--- a/validation/ice_ocean_model/cooling_then_warming_ocean.jl
+++ b/validation/ice_ocean_model/cooling_then_warming_ocean.jl
@@ -4,7 +4,7 @@ using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity
 using ClimaSeaIce
 using ClimaSeaIce.SeaIceThermodynamics: melting_temperature
 using SeawaterPolynomials: TEOS10EquationOfState
-using ClimaSeaIce.HeatBoundaryConditions: RadiativeEmission, IceWaterThermalEquilibrium
+using ClimaSeaIce.SeaIceThermodynamics.HeatBoundaryConditions: RadiativeEmission, IceWaterThermalEquilibrium
 using CairoMakie
 
 import Oceananigans.Simulations: time_step!, time
diff --git a/validation/ice_ocean_model/melting_baroclinicity.jl b/validation/ice_ocean_model/melting_baroclinicity.jl
index 9ed8326e..be6cfc0d 100644
--- a/validation/ice_ocean_model/melting_baroclinicity.jl
+++ b/validation/ice_ocean_model/melting_baroclinicity.jl
@@ -8,7 +8,7 @@ using SeawaterPolynomials: TEOS10EquationOfState, heat_expansion, haline_contrac
 
 using ClimaSeaIce
 using ClimaSeaIce: melting_temperature
-using ClimaSeaIce.HeatBoundaryConditions: RadiativeEmission, IceWaterThermalEquilibrium
+using ClimaSeaIce.SeaIceThermodynamics.HeatBoundaryConditions: RadiativeEmission, IceWaterThermalEquilibrium
 
 using Printf
 using CairoMakie