Snow model implementation

.github/workflows/ci.yml

@@ -99,6 +99,31 @@ jobs:
         env:
           TEST_GROUP: "advection"
 
+  energy_conservation:
+    name: Energy Conservation - Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}
+    runs-on: ${{ matrix.os }}
+    timeout-minutes: 120
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1.12'
+        os:
+          - ubuntu-latest
+        arch:
+          - x64
+    steps:
+      - uses: actions/checkout@v5
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: julia-actions/cache@v2
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
+        env:
+          TEST_GROUP: "energy_conservation"
+
   timestepping:
     name: Time Stepping - Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}
     runs-on: ${{ matrix.os }}
@@ -173,3 +198,28 @@ jobs:
       - uses: julia-actions/julia-runtest@v1
         env:
           TEST_GROUP: "enthalpy"
+
+  snow:
+    name: Snow Thermodynamics - Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}
+    runs-on: ${{ matrix.os }}
+    timeout-minutes: 120
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1.12'
+        os:
+          - ubuntu-latest
+        arch:
+          - x64
+    steps:
+      - uses: actions/checkout@v5
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: julia-actions/cache@v2
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
+        env:
+          TEST_GROUP: "snow"

Project.toml

@@ -1,7 +1,7 @@
 name = "ClimaSeaIce"
 uuid = "6ba0ff68-24e6-4315-936c-2e99227c95a4"
 authors = ["Climate Modeling Alliance and contributors"]
-version = "0.4.7"
+version = "0.4.8"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -16,7 +16,7 @@ SeawaterPolynomials = "d496a93d-167e-4197-9f49-d3af4ff8fe40"
 Adapt = "3, 4"
 JLD2 = "0.6.2"
 KernelAbstractions = "0.9"
-Oceananigans = "0.106, 0.107"
+Oceananigans = "0.106, 0.107, 0.108"
 RootSolvers = "0.3, 0.4, 1.0"
 Roots = "2"
 SeawaterPolynomials = "0.3.4"

docs/src/physics/thermodynamics.md

@@ -168,17 +168,17 @@ sensible_heat(Tu, clock, fields, parameters) = parameters.coefficient * (paramet
 flux = FluxFunction(sensible_heat; parameters = (coefficient = 15.0, T_air = -10.0))
 ```
 
-## Putting it together: SlabSeaIceThermodynamics
+## Putting it together: SlabThermodynamics
 
-The [`SlabSeaIceThermodynamics`](@ref) struct combines all thermodynamic components:
+The [`SlabThermodynamics`](@ref) struct combines all thermodynamic components:
 
 ```@example thermodynamics
 using Oceananigans
-using ClimaSeaIce.SeaIceThermodynamics: SlabSeaIceThermodynamics, ConductiveFlux, MeltingConstrainedFluxBalance
+using ClimaSeaIce.SeaIceThermodynamics: SlabThermodynamics, ConductiveFlux, MeltingConstrainedFluxBalance
 
 grid = RectilinearGrid(size=(10, 10), x=(0, 1e5), y=(0, 1e5), topology=(Periodic, Periodic, Flat))
 
-thermodynamics = SlabSeaIceThermodynamics(grid;
+thermodynamics = SlabThermodynamics(grid;
     top_heat_boundary_condition = MeltingConstrainedFluxBalance(),
     internal_heat_flux = ConductiveFlux(Float64; conductivity = 2.0))
 
@@ -187,3 +187,85 @@ thermodynamics
 
 When coupled to a [`SeaIceModel`](@ref), these thermodynamics automatically compute thickness
 tendencies based on the configured heat fluxes and boundary conditions.
+
+## Snow thermodynamics
+
+ClimaSeaIce supports a **layered snow model** where a snow layer sits on top of the ice slab.
+Snow insulates the ice by adding thermal resistance in series, reducing the conductive flux
+and slowing ice growth.
+
+### Snow layer physics
+
+When snow is present, the combined conductive flux through the snow+ice slab uses
+resistors in series:
+```math
+F_c = \frac{T_b - T_u}{h_s / k_s + h_i / k_i}
+```
+where ``h_s, k_s`` are the snow thickness and conductivity, and ``h_i, k_i`` are the ice
+thickness and conductivity. Snow typically has ``k_s \approx 0.31`` W/(m·K) compared to
+``k_i \approx 2`` W/(m·K) for ice, making it an effective insulating layer.
+
+The snow surface temperature ``T_u`` is determined by balancing the atmospheric heat
+fluxes against the combined conductive flux, using the same
+[`MeltingConstrainedFluxBalance`](@ref) approach as bare ice. The melting point for snow
+is 0°C (rather than the salinity-dependent liquidus for ice).
+
+### Interface temperature
+
+The snow-ice interface temperature ``T_{si}`` is computed analytically from the surface
+temperature using the thermal resistance ratio:
+```math
+T_{si} = T_b + (T_u - T_b) \frac{R_i}{R_s + R_i}
+```
+where ``R_i = h_i / k_i`` and ``R_s = h_s / k_s``. This temperature is passed to the ice
+thermodynamics as a prescribed boundary condition, so the ice layer is completely
+unaware of the snow layer above it.
+
+### Snow processes
+
+The snow model includes three processes that modify the snow thickness each time step:
+
+1. **Surface melt**: When the surface temperature reaches the melting point, any excess
+   conductive flux first melts snow before reaching the ice. The snow absorbs energy
+   ``\rho_s L_s \Delta h_s`` before excess energy melts ice from the top.
+
+2. **Snowfall accumulation**: Snow precipitates at a prescribed rate (kg/m²/s) and
+   accumulates only where ice is present (``\aleph > 0``).
+
+3. **Snow-ice formation (flooding)**: When the snow load pushes the ice surface below
+   the waterline (negative freeboard), seawater floods the snow layer and converts it
+   to ice. The freeboard criterion is:
+   ```math
+   h_f = h_i (1 - \rho_i / \rho_w) - h_s \rho_s / \rho_w < 0
+   ```
+   The flooding step formulated like this is **energy-conserving**
+
+### Snow volume conservation
+
+Snow thickness ``h_s`` represents the thickness over the ice-covered fraction ``\aleph``.
+The total snow energy per unit area is ``\rho_s L_s h_s \aleph``, consistent with how ice
+energy is ``\mathscr{L} h \aleph``. When ice concentration changes (e.g., during
+consolidation), the snow thickness is adjusted to conserve the snow volume
+``h_s \times \aleph``.
+
+### Setting up a snow-covered model
+
+Use [`snow_slab_thermodynamics`](@ref) for convenient construction with snow defaults:
+
+```@example thermodynamics
+using ClimaSeaIce
+
+grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat))
+
+snow_thermodynamics = snow_slab_thermodynamics(grid)
+
+model = SeaIceModel(grid;
+    snow_thermodynamics,
+    snowfall = 3e-6) # kg/m²/s
+
+set!(model, h=1, ℵ=1, hs=0.1) # 10 cm snow on 1 m ice
+```
+
+The [`SeaIceModel`](@ref) constructor automatically wires the layered coupling:
+the snow's internal heat flux is replaced with the combined snow+ice conductive flux,
+and the ice's top boundary condition becomes a prescribed interface temperature.

docs/src/timestepping.md

@@ -39,6 +39,7 @@ SeaIceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
 ├── grid: 16×16×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×1 halo
 ├── timestepper: ForwardEulerTimeStepper
 ├── ice_thermodynamics: SlabThermodynamics
+├── snow_thermodynamics: Nothing
 ├── advection: Nothing
 └── external_heat_fluxes:
     ├── top: Int64
@@ -61,6 +62,7 @@ SeaIceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
 ├── grid: 16×16×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×1 halo
 ├── timestepper: SplitRungeKuttaTimeStepper(3)
 ├── ice_thermodynamics: SlabThermodynamics
+├── snow_thermodynamics: Nothing
 ├── advection: Nothing
 └── external_heat_fluxes:
     ├── top: Int64

examples/freezing_bucket.jl

@@ -9,7 +9,7 @@
 # temperature, in equilibrium with the ice-water interface (and therefore fixed at
 # the melting temperature). This example demonstrates how to:
 #
-#   * use `SlabSeaIceThermodynamics` with prescribed boundary conditions,
+#   * use `SlabThermodynamics` with prescribed boundary conditions,
 #   * configure internal heat conduction,
 #   * use `FluxFunction` for frazil ice formation,
 #   * collect and visualize time series data.
@@ -54,9 +54,9 @@ internal_heat_flux = ConductiveFlux(; conductivity)
 # `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)
+heat_capacity = 2100 # J kg⁻¹ K⁻¹
+density = 900 # kg m⁻³
+phase_transitions = PhaseTransitions(; heat_capacity, density)
 
 # We set the top ice temperature:
 
@@ -67,10 +67,10 @@ top_heat_boundary_condition = PrescribedTemperature(top_temperature)
 #
 # 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)
+ice_thermodynamics = SlabThermodynamics(grid;
+                                        internal_heat_flux,
+                                        phase_transitions,
+                                        top_heat_boundary_condition)
 
 # ## Frazil ice formation
 #
@@ -88,7 +88,7 @@ bottom_heat_flux = FluxFunction(frazil_ice_formation)
 
 model = SeaIceModel(grid; ice_thermodynamics, bottom_heat_flux)
 
-# Note that the default bottom heat boundary condition for `SlabSeaIceThermodynamics`
+# Note that the default bottom heat boundary condition for `SlabThermodynamics`
 # is `IceWaterThermalEquilibrium` with freshwater. That's what we want!
 
 model.ice_thermodynamics.heat_boundary_conditions.bottom