Snow model implementation

[simone-silvestri]

Snow model implementation

This PR implements a layered snow model on top of the 0-layer slab ice model.

Algorithm

The snow layer sits on top of the ice as an independent thermodynamic layer. The ice model is completely unaware of the snow — only the boundary conditions change. The coupling is one-way: snow determines the ice top boundary condition.

Thermodynamic time step:

1. Solve for the snow surface temperature T_u via flux balance (atmospheric fluxes vs combined conductive flux), constrained by the snow melting point (0°C)
2. Compute the snow–ice interface temperature T_si analytically and pass it to the ice as a PrescribedTemperature top boundary condition
3. Ice does its thing (unchanged from the bare-ice case)
4. Snow thickness evolves: accumulation from precipitation, surface melt, and snow-ice formation (flooding)

Math

Combined conductive flux (thermal resistors in series):

F_c = (T_b - T_u) / (h_s / k_s + h_i / k_i)

where h_s, k_s ≈ 0.31 W/(m·K) are the snow thickness and conductivity, and h_i, k_i ≈ 2 W/(m·K) are the ice thickness and conductivity. When h_s = 0, this reduces to the bare-ice conductive flux k_i (T_b - T_u) / h_i.

Interface temperature (from linear temperature profile through resistors):

T_si = T_b + (T_u - T_b) · R_i / (R_s + R_i)

where R_i = h_i / k_i and R_s = h_s / k_s. When h_s = 0, R_s = 0 and T_si = T_u — the ice sees its own surface temperature, identical to the no-snow case.

Surface melt energy partition: When T_u hits the melting point, the excess energy δQ = Q_external - Q_conductive first melts snow (up to its capacity ρ_s · L_s · h_s / Δt), and any remainder melts ice from the top.

Snow-ice formation (flooding): When the snow load pushes ice below the waterline (negative freeboard), seawater floods the snow:

h_f = h_i · (1 - ρ_i / ρ_w) - h_s · ρ_s / ρ_w < 0

The flooding conversion is energy-conserving.

Snow thickness evolution:

dh_s/dt = P_s / ρ_s - G_melt - G_flooding + (h_s / ℵ) · dℵ/dt

where P_s is precipitation rate (kg/m²/s), and the last term conserves snow volume when ice concentration ℵ changes.

"No snow" case

The formulation degrades gracefully when h_s = 0:

• F_c reduces to bare-ice conductive flux
• T_si = T_u (ice sees its own surface temperature)
• Snow accumulates only where ice exists (ℵ > 0)
• Snow is removed entirely when ice disappears

Design

• Snow is a passive tracer for momentum: advected with ice velocities but does not enter the rheology
• One additional field: snow_thickness (shares concentration ℵ with sea ice)
• The IceSnowConductiveFlux type handles the combined thermal resistance
• When snow_thermodynamics is provided to SeaIceModel, the constructor automatically wires the layered coupling

Examples

The thermodynamic examples include both "bare" (ice only) and "snow" modes showing that snow insulation delays both melting and freezing.

[Copilot]

Pull request overview

This PR implements a snow model on top of the existing 0-layer slab ice model in ClimaSeaIce. The snow layer is thermodynamically coupled to the ice using a resistors-in-series approach for conductive flux, while remaining transparent to the ice model via prescribed interface temperature. The PR also renames SlabSeaIceThermodynamics to SlabThermodynamics (with a backward-compatibility alias), generalizes PhaseTransitions field names from ice-specific to generic (ice_density → density, ice_heat_capacity → heat_capacity), and fixes the latent_heat formula to use solid density instead of liquid density, matching the documented equation.

Changes:

• Adds optional snow layer (snow_thermodynamics, snow_precipitation) to SeaIceModel, with layered thermodynamic coupling, advection, flooding, and accumulation.
• Renames SlabSeaIceThermodynamics → SlabThermodynamics and generalizes PhaseTransitions fields; fixes latent_heat formula to match documentation.
• Updates all examples to compare bare-ice vs. snow-covered simulations and adds comprehensive tests for snow physics and energy conservation.

Reviewed changes

Copilot reviewed 24 out of 24 changed files in this pull request and generated 3 comments.

Show a summary per file

File	Description
src/SeaIceThermodynamics/SeaIceThermodynamics.jl	Renames PhaseTransitions fields, fixes latent_heat formula
src/SeaIceThermodynamics/slab_sea_ice_thermodynamics.jl	Renames struct, adds SlabSnowThermodynamics constructor, removes thermodynamic_tendency field
src/SeaIceThermodynamics/slab_heat_and_tracer_fluxes.jl	Adds IceSnowConductiveFlux, ice_snow_conductive_flux, and interface_temperature
src/SeaIceThermodynamics/slab_thermodynamics_tendencies.jl	Variable renames (cosmetic)
src/SeaIceThermodynamics/thermodynamic_time_step.jl	Adds layered snow+ice kernel, ice_volume_update, snow accumulation, flooding
src/sea_ice_model.jl	Adds snow_thickness, snow_thermodynamics, snow_precipitation to model; wires layered coupling
src/sea_ice_fe_step.jl	Forward Euler snow stepping support
src/sea_ice_rk_substep.jl	Runge-Kutta snow stepping support
src/tracer_tendency_kernel_functions.jl	Snow advection tendency via div_Uℵh
src/ClimaSeaIce.jl	Exports new types
test/test_snow_thermodynamics.jl	New: comprehensive snow tests
test/test_energy_conservation.jl	New: energy conservation tests with/without snow
test/test_time_stepping.jl	Adds snow to time-stepping test matrix
test/test_sea_ice_advection.jl	Snow advection test
test/test_checkpointing.jl	Snow checkpointing test
test/runtests.jl	Registers new test groups
.github/workflows/ci.yml	Adds energy_conservation CI job
examples/freezing_bucket.jl	Bare vs snow comparison
examples/freezing_of_a_lake.jl	Bare vs snow comparison with energy budgets
examples/melting_in_spring.jl	Bare vs snow comparison
examples/arctic_basin_seasonal_cycle.jl	Bare vs snow with albedo dependence
examples/perpetual_night.jl	Rename to SlabThermodynamics
docs/src/timestepping.md	Updates show output for snow
docs/src/physics/thermodynamics.md	Adds snow thermodynamics documentation

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

You can also share your feedback on Copilot code review. Take the survey.

[glwagner]

I think the main issue with this PR is that it does not implement an interface for snow models. Instead it hard code a snow model.

This can work if we take the approach that we only want to use a 1-layer snow model with a 0 layer sea ice thermodynamics. If we take that stance then a more general interface for snow modeling would be tackled in the future.

However, I think more generally what we need are interfaces for both thermodynamics and snow. In the long run, we probably want to use a snow model that comes from a land modeling package like ClimaLand or Terrarium.

[simone-silvestri]

I tried to keep the implementation as general as possible such that

• sea ice is agnostic to snow, so we can use any other snow model we want
• thermodynamics of both sea ice and snow can be changed via the thermodynamics keyword argument.
  Definitely we have to figure out how much work would it be to implement another thermodynamic, but I hope if we have a clear understanding of the equations we should have a rather generalizable framework (maybe we can revisit it a bit probably it needs some more work).

but yeah, this PR hinges on the fact that the possibility of a snow model is embedded directly inside the sea ice model. I think the difference with what you are proposing is to have a completely separate snow model? Something like

snow = SnowModel(grid)
sea_ice = SeaIceMode(grid)

# then the coupling done through another model like EarthSystemModel

Am I understanding it correctly?

[simone-silvestri]

We could probably think at different thermodynamic formulations to see if they fit within this framework easily?

[glwagner]

I understand this to be a Claude-ism who loves backwards compatibility, but I don't think we should maintain this

[glwagner]

ah, aftrer reading through the PR I understand this is also wrong. You can define a constructor

sea_ice_slab_thermodynamics

[glwagner]

for slab models, don't snow and sea ice have identical thermodynamics structure? Why do we need a new type for snow?

[glwagner]

Suggested change

	function SlabSnowThermodynamics(grid;
	function snow_slab_thermodynamics(grid;

after combing through the PR looking for the type SnowSlabThermodynamics I finally realized that this is not a type... and you are reusing SlabThermodynamics for snow and sea ice (which is great, provided we use the right names)

[glwagner]

unfortunately I don't appear to be able to suggest changes and also comment :-( so I have to "request changes

[glwagner]

We could probably think at different thermodynamic formulations to see if they fit within this framework easily?

Your top level description should mention prominently that you have generalized SlabThermodynamics to encompass snow and sea ice... this is rather important...

Now that I understand the PR better I think it is a good addition.

We can reserve a more general infrastructure for the cases where both snow and sea ice have a prognostic temperature/energy (eg going beyond mere "slab" models). For the slab models, simplicity is useful which I think motivates non-generalizability.

For the record, what I meant is that there would be functions that could dispatch on the type of snow model being used. However, the slab model is an especially simple case. We would probably need a case with prognostic temperature/energy to motivate the implementation of a more generic interface (and I suspect things get a lot tricker as well, due to the need for an implicit solve that spans all layers).

[glwagner]

Doesn't this have implications for the skin temperature solver? Ie we need to solve through both layers and the solution includes both the sea-ice/snow interface temperature as well as the snow/air interface temperature. Does this come into play in the examples as well?

There is a big diff for freezing_of_a_lake but I don't understand how snow is incorporated into it. For each example can you please describe what is going on with the snow? I think we should leave some examples unchanged, ie ice only cases. For the melting cases, it is interesting to show the transition from snow-covered to ice-only to ice-free. Can you past some figures in the PR?

[simone-silvestri]

For the record, what I meant is that there would be functions that could dispatch on the type of snow model being used. However, the slab model is an especially simple case. We would probably need a case with prognostic temperature/energy to motivate the implementation of a more generic interface (and I suspect things get a lot tricker as well, due to the need for an implicit solve that spans all layers).

Indeed, for these simple slab models the implicit solve is rather simple. At the moment, the slab snow is aware of the sea ice and not the other way around, so basically the snow solves for the top temperature with the knowledge of the sea ice state and then computes the sea-ice snow interface temperature analytically which passes to the sea ice as a PrescribedTemperature top boundary condition. I wonder how complicated this will become when we have a multi-layer sea ice.

[glwagner]

At the moment, the slab snow is aware of the sea ice and not the other way around, so basically the snow solves for the top temperature with the knowledge of the sea ice state and then computes the sea-ice snow interface temperature analytically which passes to the sea ice as a PrescribedTemperature top boundary condition

What happens when there is no snow.

[simone-silvestri]

for slab / slab, the snow computes the top surface temperature using the resistors in series approach

where Tb is the bottom temperature of the sea ice and Tu is the top temperature of the snow.
Then the interface temperature is analytical

where $$R_i = h_i / k_i$$ and $$R_s = h_s / k_s$$

[simone-silvestri]

I am not sure if the same approach could be used also with multi-layer sea-ice where snow sees only the first sea ice layer from the top or if we need to do an implicit solve over the whole column to find the top temperature

[glwagner]

Thank you. It's clear how this captures the "no snow" case identically as if the snow were not included at all.

Please put the algorithm + math in the top-level of this PR and in the docs!

[simone-silvestri]

@glwagner I should have addressed all the comments

[glwagner]

dont we need categories for snow and ice? confused what is going on here.

[glwagner]

semantics is still a little messy, but overall good

confused about PhaseTransitions though. I think we need a modular design for thermodynamic properties, with PhaseTransitions there is repeated info it seems. Eg if we have slab snow and slab ice, then the properties of liquid water are repeated. They can even be inconsistent. So that's a bad design.

We can improve in the future. It might be worth looking at now Breeze does this.

One thing to note; the thermodynamic properties should ultimately represent microscopic properties (eg microscopic density, etc). There is also a bulk density (for snow) which depends on porosity and would dynamically evolve.

The design should be enlightened to all of these facts.

[simone-silvestri]

Sounds good. I ll merge this and start a PR to generalize the phase transitions. Maybe we want only one thermodynamics field with thermodynamics.ice and thermodynamics.snow and thermodynamics.phase_transitions