(0.5.0) Share phase_transitions and bulk densities at the model level

[simone-silvestri]

A (Very substantial!!) refactor. Closes #133.

Summary

Reorganises the thermodynamic state of SeaIceModel so that shared
physical parameters live at the model level rather than being duplicated
inside each SlabThermodynamics. Along the way, corrects a sign error
in the Stefan correction and deletes the sea_ice_model.jl rebuild
hack.

Key moves:

• One PhaseTransitions on SeaIceModel (shared microscopic
  pure-ice ↔ water parameters), replacing one per slab.
• Two new top-level SeaIceModel fields, ice_density and
  snow_density, carrying the bulk porous-medium densities.
• SlabThermodynamics slimmed down to the state that is genuinely
  per-layer (surface T, heat BCs, internal-flux coefficient,
  concentration evolution).
• latent_heat rewritten as a per-mass function (with the correct
  sign); tendency kernels multiply externally by the appropriate bulk
  density.
• Removal of the sea_ice_model.jl:174–204 rebuild hack that was
  forcing ice_thermodynamics.top_heat_boundary_condition = PrescribedTemperature(0) and assembling the combined snow+ice flux
  on the snow side at construction time.
• New documented extension point (flux_kernel) plus passthrough for
  pre-assembled FluxFunction as internal_heat_flux.

Problem → solution

P1. Duplicated liquid properties

SlabThermodynamics used to own a PhaseTransitions, which holds
liquid_density, liquid_heat_capacity, reference_temperature, and
liquidus. These describe the same liquid water regardless of whether
it emerges from melting snow or sea ice. Two copies made it
representationally possible for the two slabs to disagree.

Fix. A single PhaseTransitions lives on SeaIceModel as
phase_transitions. Both ice_thermodynamics and
snow_thermodynamics read from it.

P2. Misplaced bulk density inside PhaseTransitions (#133)

Per @glwagner in #133: PhaseTransitions.density and
PhaseTransitions.heat_capacity refer to the microscopic pure-ice
crystal — the same substance in both the snow and sea-ice porous
media. snow_slab_thermodynamics used to override density = 330 in
the snow's PhaseTransitions, which is a bulk snow density (porous
medium, including trapped air).

Fix. PhaseTransitions is now documented and used as a purely
microscopic descriptor. The two bulk densities live as top-level
SeaIceModel fields:

• ice_density (default 900 kg/m³, bulk sea ice including brine)
• snow_density (default 330 kg/m³, bulk snow including trapped air)

Both are wrapped as Field{Center, Center, Nothing} by the
constructor (scalars become ConstantField), matching the old
model.ice_density treatment.

P3. Two roles for ρ in the latent-heat expression

The old latent_heat(pt, T) returned the volumetric form ρ × ℒ(T),
with ρ sourced from PhaseTransitions. The ρ played two distinct
physical roles simultaneously:

1. Mass-budget multiplier — converts per-mass enthalpy to
   per-volume of the melting medium. Must be the bulk density.
2. Stefan-correction factor inside
   ρℒ(T) = ρℒ₀ + (ρ_ℓ c_ℓ − ρ c)(T − T₀). Must be the
   microscopic pure-ice density.

A single ρ can play both correctly only when ρ_bulk = ρ_pure.
When the snow helper overrode density = 330, the ρ × c term
became physically incorrect.

Fix. latent_heat(pt, T) returns a per-mass Stefan-corrected
value using only microscopic properties:

$$ℒ(T) = ℒ₀ + (ρ_ℓ c_ℓ / ρ − c)(T − T₀)$$

Call sites multiply externally by ρ_ice_bulk or ρ_snow_bulk.
Matches CICE and CLM/Noah-MP conventions. (A sign error that appeared
during an earlier iteration of the refactor was detected by the energy
conservation tests — now closes to < 1e-15 across all combinations.)

P4. Rebuild hack when both layers are present

The pre-refactor SeaIceModel constructor rebuilt both thermodynamics
when snow was supplied:

• ice_thermodynamics got PrescribedTemperature(0) as its top BC so
  its surface solve would be skipped.
• snow_thermodynamics got a combined IceSnowConductiveFlux(ks, ki)
  as its internal flux.

This coupled the two slabs at construction time and made the ice's
BC depend on whether snow existed elsewhere on the model.

Fix. The combined flux is constructed inside the layered
kernel from the two conductivities. The ice's melt/freeze tendency is
factored into a pure helper ice_melt_freeze_tendency(…, Tui, …) that
takes the ice-top temperature as an argument and performs no surface
solve. The layered kernel solves for T_us, computes T_si
analytically from the resistance ratio, and calls the helper with
T_si directly. No BC rewriting.

P5. Liquidus plumbing

The liquidus is only ever used via bottom_temperature(bc, liquidus)
for IceWaterThermalEquilibrium. An interim design stored liquidus
inside each slab's internal_heat_flux.parameters, forcing the
SeaIceModel constructor to rebind it from
model.phase_transitions.liquidus into every slab.

Fix. SlabThermodynamics stores only the raw flux coefficient
(e.g. ConductiveFlux). The FluxFunction wrapper is assembled
inside the tendency kernel, pulling liquidus from
model.phase_transitions at call time. No per-slab liquidus state.

Breaking changes to the UI

SlabThermodynamics and its helpers

PhaseTransitions is no longer stored inside SlabThermodynamics.
Construct it once on the model:

# OLD
ice_thermodynamics = SlabThermodynamics(grid;
    phase_transitions = PhaseTransitions(; liquid_density = 1020))

# NEW
model = SeaIceModel(grid;
    ice_thermodynamics = SlabThermodynamics(grid),
    phase_transitions = PhaseTransitions(; liquid_density = 1020))

snow_slab_thermodynamics no longer accepts density,
heat_capacity, or reference_latent_heat:

# OLD
snow_thermodynamics = snow_slab_thermodynamics(grid;
    density = 330,
    heat_capacity = 2090,
    reference_latent_heat = 334e3)

# NEW
model = SeaIceModel(grid;
    snow_thermodynamics = snow_slab_thermodynamics(grid),  # just conductivity
    snow_density = 330,
    phase_transitions = PhaseTransitions(; reference_latent_heat = 334e3))

SlabThermodynamics.internal_heat_flux now stores the raw flux
coefficient rather than a pre-wrapped FluxFunction:

# OLD (introspection)
model.ice_thermodynamics.internal_heat_flux.parameters.flux.conductivity

# NEW
model.ice_thermodynamics.internal_heat_flux.conductivity

SeaIceModel

New keyword arguments with safe defaults:

SeaIceModel(grid;
    phase_transitions = PhaseTransitions(eltype(grid)),  # NEW
    ice_density       = 900,                              # kg/m³, existed
    snow_density      = 330,                              # NEW
    ice_thermodynamics  = sea_ice_slab_thermodynamics(grid),
    snow_thermodynamics = nothing,
    ...)

latent_heat

Now returns per-mass enthalpy. If you were using it for diagnostics,
multiply by the bulk density of the medium:

# OLD
enthalpy_density = latent_heat(phase_transitions, T)   # J/m³

# NEW
enthalpy_density = model.ice_density[i,j,1] * latent_heat(phase_transitions, T)

The new UI in action

Bare sea ice (default):

using ClimaSeaIce, Oceananigans

grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat))
model = SeaIceModel(grid)                # uses defaults everywhere
set!(model, h=1, ℵ=1)
time_step!(model, 3600)

Layered snow + ice:

model = SeaIceModel(grid;
    ice_thermodynamics  = sea_ice_slab_thermodynamics(grid),
    snow_thermodynamics = snow_slab_thermodynamics(grid),
    ice_density         = 900,
    snow_density        = 330,
    snowfall            = 6e-5) # kg m⁻² s⁻¹
set!(model, h=1, ℵ=1, hs=0.1)
time_step!(model, 3600)

Dynamics-only (no phase physics):

dynamics = SeaIceMomentumEquation(grid; rheology = ViscousRheology(ν=1000))
model = SeaIceModel(grid;
    dynamics,
    ice_thermodynamics = nothing,
    ice_density        = 900,   # still required by momentum equations
    advection          = WENO())

Custom liquidus / phase transition parameters:

model = SeaIceModel(grid;
    phase_transitions = PhaseTransitions(;
        liquidus        = LinearLiquidus(; slope = 0.0543),
        reference_latent_heat = 335e3))

Pluggable internal flux. Four shapes accepted, documented in
docs/src/physics/thermodynamics.md:

# 1. Default ConductiveFlux
SlabThermodynamics(grid; internal_heat_flux = ConductiveFlux(conductivity=2))

# 2. Bare Function (signature (i,j,grid,Tu,clock,fields,parameters)->Q)
my_flux(i,j,grid,Tu,clock,fields,p) = ...
SlabThermodynamics(grid; internal_heat_flux = my_flux)

# 3. Pre-assembled FluxFunction (passthrough, no parameter injection)
SlabThermodynamics(grid; internal_heat_flux = FluxFunction(my_flux; parameters=(...,), ...))

# 4. Custom struct via flux_kernel extension (one-line pirate)
struct NonLinearConductiveFlux{T}; k0::T; α::T; end
ClimaSeaIce.SeaIceThermodynamics.flux_kernel(::NonLinearConductiveFlux) = my_nonlinear_kernel
SlabThermodynamics(grid; internal_heat_flux = NonLinearConductiveFlux(2.0, 0.01))

Pros and cons

Pros

• Correctness. Single shared PhaseTransitions makes it structurally
  impossible for snow and ice to disagree about the liquid; the old
  duplication had no guardrail.
• Physical transparency. The volumetric ρ × ℒ(T) split into
  per-mass ℒ(T) × bulk density disentangles the two roles of ρ. The
  previously-hidden sign behaviour of the Stefan correction is now
  covered by machine-precision energy-conservation tests.
• Pluggability. Swapping in a non-default snow model (or ice
  thermodynamics) is simpler because the two slabs are no longer
  welded together by the constructor hack. Custom internal fluxes are
  supported via four documented paths.
• Smaller constructor. Removing the rebuild hack deletes ~30 lines
  of state surgery from sea_ice_model.jl.
• Defaults are safe. Existing scripts that only touched
  ice_thermodynamics and snow_thermodynamics keyword args continue
  to work, because the new phase_transitions and snow_density
  default sensibly.

Cons

• Breaking change. Any user code setting density, heat_capacity,
  or reference_latent_heat on snow_slab_thermodynamics, or reading
  phase_transitions off a SlabThermodynamics, or indexing into
  .internal_heat_flux.parameters.flux, needs to be updated. Users of
  default construction are unaffected.
• Slightly larger kwarg surface on SeaIceModel. Two new kwargs
  (phase_transitions, snow_density) add API surface. Both have
  safe defaults.
• Layered kernel still assumes conductivity-based coupling. The
  layered snow+ice kernel hardcodes IceSnowConductiveFlux(ks, ki)
  built from each layer's .conductivity. A user who wants a
  non-conductive layered coupling still has to touch the layered
  kernel. A future PR can lift this via a
  combined_internal_flux(ice, snow) trait.
• SlabThermodynamics.internal_heat_flux type is looser. It now
  holds a raw flux coefficient and the tendency kernel wraps it
  inline. Users who relied on the pre-wrapped FluxFunction shape
  for introspection need to adapt.

Test plan

• test_snow_thermodynamics.jl — 20/20 (layered physics,
  flooding, accumulation, melt ordering).
• test_energy_conservation.jl — 6/6 at machine precision
  (< 1e-15) across bare/snow/freeze/melt/precipitation.
• test_sea_ice_advection.jl — 10/10 (bare, layered,
  dynamics-only).
• test_checkpointing.jl — 155/155 (all pickup modes; bare,
  ice+dynamics, snow+dynamics).
• Time-stepping smoke (3 thermo paths × 2 steppers) — 5/5. Full
  144-combo matrix in test_time_stepping.jl skipped locally for
  CI time; will be exercised on CI.

🤖 Generated with Claude Code

[glwagner]

Quick comment, for this:

Two new top-level SeaIceModel fields, ice_density and snow_density, carrying the bulk porous-medium densities.

I suggest defining "ice" as "pure ice" and then use the term "sea ice" to refer to the porous matrix of pure ice and brine. Therefore, the top-level SeaIceModel property should be sea_ice_density.

[glwagner]

Is the latent heat solution correct? For a bulk media representation, the latent heat will be the weighted sum of components, eg

$$ \mathscr{L} = \sum_i q_i \mathscr{L}_i $$

where $i$ denotes the component and q is the mass fraction of the component; ie $q_i = \rho_i / \rho$ where $\rho_i = m_i / V$; ie the "bulk component density" (mass of component divided by total volume of the region) and $\rho$ is the bulk density. This is what we do for the atmosphere, which is a four-component media: dry air, vapor, liquid, and ice.

For sea ice, we have two components, ice and brine. For snow, we have ice and air. In the case of snow, it may be a reasonable approximation to neglect the latent heat of the air. But we shouldn't do this flippantly or by prescription (also it may not even really be necessary, since it won't save any computation and leads to come confusion about the physics -- we just need to know the latent heat of saturated air?)

Also for snow we may need to have the density and latent heat of air as PhaseTransitions properties.

[simone-silvestri]

I have changed the ice_density to be sea_ice_density. For the latent heat, this is the simplified approach of CICE. We can track the latent heat problem in an issue, I think it would be rather simple to implement the more correct "multi-species" approach, and we can implement it in a fresh PR