From a1de817d373e2d558f6b774dde68e49f945c1c42 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 23 Apr 2026 18:46:14 +0200 Subject: [PATCH 1/7] add a convervation --- examples/coupled_conservation.jl | 344 +++++++++++++++++++++ src/Oceans/assemble_net_ocean_fluxes.jl | 18 +- src/SeaIces/assemble_net_sea_ice_fluxes.jl | 4 +- test/test_conservation.jl | 251 +++++++++++++++ 4 files changed, 609 insertions(+), 8 deletions(-) create mode 100644 examples/coupled_conservation.jl create mode 100644 test/test_conservation.jl diff --git a/examples/coupled_conservation.jl b/examples/coupled_conservation.jl new file mode 100644 index 000000000..ac1c87f0e --- /dev/null +++ b/examples/coupled_conservation.jl @@ -0,0 +1,344 @@ +# # Coupled energy conservation +# +# In this example, we run a minimal-physics `OceanSeaIceModel` through a freeze-then-melt cycle and verify +# that the coupled energy budget closes. The setup is a single ocean column with thermodynamics-only sea ice, a snow +# layer on top of the ice, and a uniform prescribed atmosphere. We drive two 30-day phases: a cold phase with +# light snowfall that grows the ice, and a warm phase with rainfall that melts it. +# +# The invariant we check at the end of the run is +# +# ``` +# ΔE = ∫ Q dt +# ``` +# +# where `E = Hₒ + Eis` is the ocean heat content plus the ice+snow stored latent energy, and `Q` is the atmospheric +# heat flux into the coupled system. Closure to machine precision requires that every internal flux cancels exactly +# between the components. +# +# ## Install dependencies +# +# ```julia +# using Pkg +# pkg"add Oceananigans, NumericalEarth, ClimaSeaIce, CairoMakie" +# ``` + +using Oceananigans +using Oceananigans.Units +using Oceananigans.Fields: interior +using Oceananigans.Operators: Azᶜᶜᶜ + +using ClimaSeaIce +using ClimaSeaIce.SeaIceThermodynamics: latent_heat + +using NumericalEarth +using NumericalEarth.Atmospheres: PrescribedAtmosphere +using NumericalEarth.Diagnostics: atmosphere_ocean_heat_flux, frazil_heat_flux +using NumericalEarth.EarthSystemModels: OceanSeaIceModel, Radiation, update_state! +using NumericalEarth.Oceans: ocean_simulation + +using CairoMakie +using Printf + +# ## Constant latent heat for diagnostic closure +# +# `ClimaSeaIce`'s slab mass balance uses a temperature-dependent latent heat, `ℒ(T) = ℒ₀ + (ρℓ cℓ / ρi − cᵢ)(T − T₀)`, +# with `ℰu = ρi · ℒ(T_u)` at the top interface and `ℰb = ρi · ℒ(T_b)` at the bottom. A single state-based +# `Eis = − ℵ · ρi · ℒ · h · Az` cannot simultaneously match freeze at `T_b` and top-melt at 0 ᵒC: the 4.7 kJ/kg gap accounts +# for a ~1% residual scaling with top-melt mass. To isolate coupler-side bookkeeping from this intrinsic `ℒ(T)` mismatch we +# locally override `latent_heat` to the constant `pt.reference_latent_heat`. This is a diagnostic choice for the present +# example and does not modify upstream. + +@inline function ClimaSeaIce.SeaIceThermodynamics.latent_heat( + pt::ClimaSeaIce.SeaIceThermodynamics.PhaseTransitions, T) + return pt.reference_latent_heat +end + +# ## Grid, ocean, sea ice, and atmosphere +# +# We build a single ocean column 100 m deep with 10 vertical levels on a `(Flat, Flat, Bounded)` `RectilinearGrid`. +# The ocean is started just above freezing at `S = 34`, with advection and Coriolis turned off and +# `CATKEVerticalDiffusivity` providing vertical mixing. + +arch = CPU() + +grid = RectilinearGrid(arch; + size = 10, + z = (-100, 0), + topology = (Flat, Flat, Bounded)) + +ocean = ocean_simulation(grid; + momentum_advection = nothing, + tracer_advection = nothing, + closure = CATKEVerticalDiffusivity(), + coriolis = nothing, + bottom_drag_coefficient = 0) + +Tᵢ = -1.5 # ᵒC, just above freezing at S = 34 +Sᵢ = 34.0 # psu +set!(ocean.model, T = Tᵢ, S = Sᵢ) + +# Sea ice has thermodynamics only and starts with `h = 1 m`, `ℵ = 1`, and a 10 cm snow layer. + +sea_ice = sea_ice_simulation(grid, ocean; + dynamics = nothing, + advection = nothing) + +set!(sea_ice.model, h = 1, ℵ = 1, hs = 0.10) + +# The atmosphere is prescribed and spatially uniform. It lives on its own scalar grid and we overwrite +# the `parent` array in place at the start of each phase. + +atmosphere_grid = RectilinearGrid(arch; size=(), topology=(Flat, Flat, Flat)) +atmosphere = PrescribedAtmosphere(atmosphere_grid, [0.0, 1e9]) +radiation = Radiation() + +coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) + +# ## Helpers +# +# `set_atmosphere!` fills every `FieldTimeSeries` of the prescribed atmosphere with scalar constants +# so the atmospheric forcing is spatio-temporally uniform. + +function set_atmosphere!(atmosphere, T, q, u, v, p, ℐꜜˢʷ, ℐꜜˡʷ, Jᶜ, Jˢⁿ) + fill!(parent(atmosphere.tracers.T), T ) + fill!(parent(atmosphere.tracers.q), q ) + fill!(parent(atmosphere.velocities.u), u ) + fill!(parent(atmosphere.velocities.v), v ) + fill!(parent(atmosphere.pressure), p ) + fill!(parent(atmosphere.downwelling_radiation.shortwave), ℐꜜˢʷ) + fill!(parent(atmosphere.downwelling_radiation.longwave), ℐꜜˡʷ) + fill!(parent(atmosphere.freshwater_flux.rain), Jᶜ ) + fill!(parent(atmosphere.freshwater_flux.snow), Jˢⁿ ) + return nothing +end + +# Physical constants read directly from the coupled model so the budget diagnostics stay consistent with +# the constants the model itself is using. `Az` is the horizontal cell area (unity for `(Flat, Flat, Bounded)`). + +ρi = coupled_model.sea_ice.model.sea_ice_density[1, 1, 1] +ρs = coupled_model.sea_ice.model.snow_density[1, 1, 1] +ℒ₀ = ClimaSeaIce.SeaIceThermodynamics.latent_heat(coupled_model.sea_ice.model.phase_transitions, 0.0) +ρᵒᶜ = coupled_model.interfaces.ocean_properties.reference_density +cᵒᶜ = coupled_model.interfaces.ocean_properties.heat_capacity + +Az = Azᶜᶜᶜ(1, 1, 1, grid) + +# Volume integral of the ocean temperature, built once and re-used every step via `compute!`. The underlying `Integral` operation +# keeps a reference to the live `T` field, so each `compute!` re-evaluates over the current state. + +∫T = Field(Integral(ocean.model.tracers.T)) + +# `column_state` returns a snapshot of the diagnostic quantities we track: ice and snow geometry, the ice+snow stored latent energy, +# and the ocean heat content. + +function column_state(coupled_model) + h = first(interior(coupled_model.sea_ice.model.ice_thickness)) + ℵ = first(interior(coupled_model.sea_ice.model.ice_concentration)) + hs = first(interior(coupled_model.sea_ice.model.snow_thickness)) + + Eis = -ℵ * (ρi * ℒ₀ * h + ρs * ℒ₀ * hs) * Az + Hₒ = ρᵒᶜ * cᵒᶜ * first(compute!(∫T)) + + return (; h, ℵ, hs, Eis, Hₒ) +end + +# `net_top_heat_flux` returns the atmospheric energy input to the coupled (ice + ocean) system in Watts: +# `Q_atm = − (ΣQt + ΣQao) · Az`, where `ΣQt` is the sea-ice top heat flux per cell and `ΣQao` is the per-cell +# atmosphere-to-ocean flux over the open-water fraction. The ocean-side piece comes from `atmosphere_ocean_heat_flux`, +# which subtracts the frazil and interface contributions internally so it never picks up a spurious ocean / ice exchange term. + +function net_top_heat_flux(coupled_model) + ΣQt = first(interior(coupled_model.interfaces.net_fluxes.sea_ice.top.heat)) + ΣQao = first(interior(atmosphere_ocean_heat_flux(coupled_model))) + return -(ΣQt + ΣQao) * Az +end + +# ## Running the freeze-melt cycle +# +# We run two 30-day phases at `Δt = 10 min`: a cold phase with light snowfall that grows the ice, then a warm phase with +# rain and strong radiation that melts it back. The run is driven by a single `Simulation` spanning the full cycle. +# Two callbacks do the bookkeeping: `budget_callback` records state at every step, and `phase_switch_callback` swaps the +# atmosphere from freeze to melt at `t = Δτ`. + +Δt = 20minutes +Δτ = 40days + +simulation = Simulation(coupled_model; Δt, stop_time = 2Δτ) + +freeze_phase = (T = 253.15, + q = 1.0e-4, + u = 2.0, + v = 0.0, + p = 101325.0, + ℐꜜˢʷ = 50.0, + ℐꜜˡʷ = 180.0, + Jᶜ = 0.0, + Jˢⁿ = 1.0e-5) + +melt_phase = (T = 278.15, + q = 5.0e-3, + u = 2.0, + v = 0.0, + p = 101325.0, + ℐꜜˢʷ = 250.0, + ℐꜜˡʷ = 320.0, + Jᶜ = 5.0e-6, + Jˢⁿ = 0.0) + +# We keep a history of the budget-relevant quantities at every time step. + +history = (t = Float64[], + phase = Int[], + h = Float64[], + ℵ = Float64[], + hs = Float64[], + Eis = Float64[], + Hₒ = Float64[], + Q = Float64[], + 𝒬ᶠʳᶻ = Float64[]) + +function record!(history, coupled_model, phase_id, Q) + st = column_state(coupled_model) + 𝒬f = first(interior(frazil_heat_flux(coupled_model))) + push!(history.t, coupled_model.clock.time) + push!(history.phase, phase_id) + push!(history.h, st.h) + push!(history.ℵ, st.ℵ) + push!(history.hs, st.hs) + push!(history.Eis, st.Eis) + push!(history.Hₒ, st.Hₒ) + push!(history.Q, Q) + push!(history.𝒬ᶠʳᶻ, 𝒬f) + return nothing +end + +# The `budget_callback` reads the current `phase_ctx` — a small mutable box holding the current phase id +# and the snowfall enthalpy `Qᵖ` to add to the net atmospheric flux. `phase_switch_callback` is the one that +# updates this context at the phase boundary. + +phase_ctx = Ref((; phase_id = 1, Qᵖ = - freeze_phase.Jˢⁿ * ℒ₀ * Az)) + +function budget_callback(simulation) + ctx = phase_ctx[] + Q = net_top_heat_flux(simulation.model) + ctx.Qᵖ + record!(history, simulation.model, ctx.phase_id, Q) + return nothing +end + +# At `t = Δτ` the atmosphere swaps from freeze to melt. `update_state!` would zero the pending frazil flux +# (the ocean is at `Tₘ`, no supercooling), stranding the latent energy already deposited into the ocean by +# the last freeze step's frazil mutation. We preserve `𝒬ᶠʳᶻ` across the refresh and add it back into the +# sea-ice bottom heat flux that the slab will read. The callback also overwrites the just-recorded `Q` entry +# with the melt-phase starting flux, which is the flux that will drive the next step under rectangle-at-start +# integration. +# +# Oceananigans fires every scheduled callback once at initialization to sync its schedule, so we guard against +# the `t = 0` fire — we only want to switch at the actual phase boundary. + +function phase_switch_callback(simulation) + simulation.model.clock.time < Δτ && return nothing + + set_atmosphere!(atmosphere, melt_phase.T, melt_phase.q, melt_phase.u, melt_phase.v, + melt_phase.p, melt_phase.ℐꜜˢʷ, melt_phase.ℐꜜˡʷ, melt_phase.Jᶜ, melt_phase.Jˢⁿ) + + Qᵖ = - melt_phase.Jˢⁿ * ℒ₀ * Az + 𝒬ᶠʳᶻ = simulation.model.interfaces.sea_ice_ocean_interface.fluxes.frazil_heat + ΣQb = simulation.model.interfaces.net_fluxes.sea_ice.bottom.heat + 𝒬⁻ = first(interior(𝒬ᶠʳᶻ)) # pending frazil + update_state!(simulation.model) + interior(𝒬ᶠʳᶻ, 1, 1, 1) .= 𝒬⁻ + interior(ΣQb, 1, 1, 1) .+= 𝒬⁻ + + phase_ctx[] = (; phase_id = 2, Qᵖ) + history.Q[end] = net_top_heat_flux(simulation.model) + Qᵖ + return nothing +end + +add_callback!(simulation, budget_callback, IterationInterval(1)) +add_callback!(simulation, phase_switch_callback, SpecifiedTimes([Δτ])) + +# Initialize the coupler for the freeze atmosphere. Oceananigans' `run!` will re-run `update_state!` at init +# and then fire `budget_callback` once, which serves as the `t = 0` seed entry for the history. + +set_atmosphere!(atmosphere, freeze_phase.T, freeze_phase.q, freeze_phase.u, freeze_phase.v, + freeze_phase.p, freeze_phase.ℐꜜˢʷ, freeze_phase.ℐꜜˡʷ, freeze_phase.Jᶜ, freeze_phase.Jˢⁿ) + +@info "Running 60-day coupled freeze/melt cycle…" +run!(simulation) + +# ## Budget analysis +# +# Cumulative atmospheric input uses rectangle-at-start integration to match what the coupler actually applied +# during each step. + +t = history.t +τ = t ./ day # days axis + +∫Q = similar(t) +∫Q[1] = 0.0 +for n in 2:length(t) + ∫Q[n] = ∫Q[n-1] + history.Q[n-1] * (t[n] - t[n-1]) +end + +# The frazil mass gain is deposited by `compute_sea_ice_ocean_fluxes!` at the end of step `n` +# (mutating ocean `T` and writing `𝒬ᶠʳᶻ`) but the corresponding ice mass gain is consumed only during +# step `n + 1`. At a diagnostic snapshot the ocean shows the warming while the ice has not yet grown. +# We anticipate this one-step pending quantity by adding `𝒬ᶠʳᶻ(n) · Δt⁺ · Az` to `Eis(n)` so the energy +# budget closure is not polluted by bookkeeping lag. + +Δt⁺ = similar(t) +for n in 1:(length(t) - 1) + Δt⁺[n] = t[n+1] - t[n] +end +Δt⁺[end] = Δt⁺[end-1] +δE = history.𝒬ᶠʳᶻ .* Δt⁺ .* Az + +Ẽᵢₛ = history.Eis .+ δE +ΔE = (Ẽᵢₛ .+ history.Hₒ) .- (Ẽᵢₛ[1] + history.Hₒ[1]) +R = ΔE .- ∫Q + +# ## Visualizing the budget +# +# The plot shows the two stored components (ocean heat content and ice+snow stored latent energy), the cumulative +# match between `ΔE` and `∫Q dt`, and the residual on both absolute and relative log scales. + +set_theme!(Theme(fontsize=16, linewidth=2)) + +fig = Figure(size=(1100, 900)) +axE1 = Axis(fig[1, 1], ylabel = "Ocean heat content (J)", title = "Ocean heat content") +axE2 = Axis(fig[2, 1], ylabel = "Ice + snow stored E (J)", title = "Ice + snow stored latent energy") +axE3 = Axis(fig[3, 1], ylabel = "Cumulative (J)", title = "ΔE vs. ∫Q dt") +axE4 = Axis(fig[4, 1], ylabel = "Residual (J)", title = "Energy residual = ΔE − ∫Q dt") +axE5 = Axis(fig[5, 1], ylabel = "log₁₀|rel residual|", title = "Relative energy residual", xlabel = "Time (days)", ) + +lines!(axE1, τ, history.Hₒ, color = :royalblue) +lines!(axE2, τ, history.Eis, color = :orange) +lines!(axE3, τ, ΔE, label = "ΔE", color = :black) +lines!(axE3, τ, ∫Q, label = "∫Q dt", color = :crimson, linestyle = :dash) +lines!(axE4, τ, R, color = :seagreen) +ε = log10.(abs.(R ./ max(maximum(abs.(ΔE)), 1))) +lines!(axE5, τ[2:end], ε[2:end], color = :seagreen) +axislegend(axE3, position = :lt) +hlines!(axE4, [0]; color = :gray, linestyle = :dot) +for ax in (axE1, axE2, axE3, axE4, axE5) + vlines!(ax, [Δτ / day]; color = :gray, linestyle = :dot, linewidth = 1) +end + +save("coupled_conservation_energy.png", fig) +nothing #hide + +# ![](coupled_conservation_energy.png) + +# ## Per-phase summary + +nᶠ = findlast(p -> p == 1, history.phase) + +ΔEᶠ = Ẽᵢₛ[nᶠ] + history.Hₒ[nᶠ] - Ẽᵢₛ[1] - history.Hₒ[1] +ΔEᵐ = Ẽᵢₛ[end] + history.Hₒ[end] - Ẽᵢₛ[nᶠ] - history.Hₒ[nᶠ] +∫Qᶠ = ∫Q[nᶠ] +∫Qᵐ = ∫Q[end] - ∫Q[nᶠ] + +@printf(" freeze: ΔE = %+.3e J ∫Q = %+.3e J residual = %+.2e (%.1e rel)\n", ΔEᶠ, ∫Qᶠ, ΔEᶠ - ∫Qᶠ, abs(ΔEᶠ - ∫Qᶠ) / max(abs(ΔEᶠ), 1)) +@printf(" melt : ΔE = %+.3e J ∫Q = %+.3e J residual = %+.2e (%.1e rel)\n", ΔEᵐ, ∫Qᵐ, ΔEᵐ - ∫Qᵐ, abs(ΔEᵐ - ∫Qᵐ) / max(abs(ΔEᵐ), 1)) +@printf(" full-cycle relative residual: %.1e\n", abs(R[end]) / max(maximum(abs.(ΔE)), 1)) +nothing #hide diff --git a/src/Oceans/assemble_net_ocean_fluxes.jl b/src/Oceans/assemble_net_ocean_fluxes.jl index a1b1d0dcd..21378696b 100644 --- a/src/Oceans/assemble_net_ocean_fluxes.jl +++ b/src/Oceans/assemble_net_ocean_fluxes.jl @@ -44,6 +44,7 @@ function update_net_ocean_fluxes!(coupled_model, ocean_model, grid) ℐꜜˡʷ = atmosphere_fields.ℐꜜˡʷ.data) freshwater_flux = atmosphere_fields.Jᶜ.data + snowfall_flux = atmosphere_fields.Jˢⁿ.data ice_concentration = sea_ice_concentration(sea_ice) ocean_surface_salinity = EarthSystemModels.ocean_surface_salinity(ocean_model) @@ -66,6 +67,7 @@ function update_net_ocean_fluxes!(coupled_model, ocean_model, grid) ice_concentration, downwelling_radiation, freshwater_flux, + snowfall_flux, atmos_ocean_properties, ocean_properties) @@ -83,6 +85,7 @@ end sea_ice_concentration, downwelling_radiation, freshwater_flux, + snowfall_flux, atmos_ocean_properties, ocean_properties) @@ -101,6 +104,7 @@ end Tₛ = convert_to_kelvin(ocean_properties.temperature_units, Tₛ) Jᶜ = freshwater_flux[i, j, 1] # Prescribed freshwater (condensate) flux + Jˢⁿ = snowfall_flux[i, j, 1] # Prescribed snowfall component of Jᶜ ℐꜜˢʷ = downwelling_radiation.ℐꜜˢʷ[i, j, 1] # Downwelling shortwave radiation ℐꜜˡʷ = downwelling_radiation.ℐꜜˡʷ[i, j, 1] # Downwelling longwave radiation 𝒬ᵀ = get_possibly_zero_flux(atmos_ocean_fluxes, :sensible_heat)[i, j, 1] # sensible or "conductive" heat flux @@ -136,13 +140,13 @@ end # Convert from a mass flux to a volume flux (aka velocity) # by dividing with the ocean reference density. # Also switch the sign, for some reason we are given freshwater flux as positive down. + # Rain falls on both ice and ocean but is routed to the ocean always (melt-pond + # effects are absent); snow accumulates on ice and only reaches the ocean + # through the open-water fraction (1-ℵ). Evaporation is restricted to the + # open-water fraction (turbulent fluxes are per-ice elsewhere). ρᵒᶜ⁻¹ = 1 / ocean_properties.reference_density - ΣFao = - Jᶜ * ρᵒᶜ⁻¹ - - # Add the contribution from the turbulent water vapor flux, which has - # a different sign convention as the prescribed water mass fluxes (positive upwards) - Jᵛᵒᶜ = Jᵛ * ρᵒᶜ⁻¹ - ΣFao += Jᵛᵒᶜ + Jʳⁿ = Jᶜ - Jˢⁿ # rain mass flux + ΣFao = - (Jʳⁿ + (1 - ℵᵢ) * Jˢⁿ) * ρᵒᶜ⁻¹ + (1 - ℵᵢ) * Jᵛ * ρᵒᶜ⁻¹ # Compute fluxes for u, v, T, and S from momentum, heat, and freshwater fluxes τˣ = net_ocean_fluxes.u @@ -173,6 +177,6 @@ end # Tracer fluxes Jᵀ[i, j, 1] = ifelse(inactive, zero(grid), Jᵀao + Jᵀio) # Jᵀao is already multiplied by the sea ice concentration - Jˢ[i, j, 1] = ifelse(inactive, zero(grid), (1 - ℵᵢ) * Jˢao + Jˢio) + Jˢ[i, j, 1] = ifelse(inactive, zero(grid), Jˢao + Jˢio) end end diff --git a/src/SeaIces/assemble_net_sea_ice_fluxes.jl b/src/SeaIces/assemble_net_sea_ice_fluxes.jl index daf8f00e3..f9b3af822 100644 --- a/src/SeaIces/assemble_net_sea_ice_fluxes.jl +++ b/src/SeaIces/assemble_net_sea_ice_fluxes.jl @@ -97,7 +97,9 @@ end ℐₜˢʷ = transmitted_shortwave_radiation(i, j, kᴺ, grid, time, α, ℐꜜˢʷ) ℐₐˡʷ = absorbed_longwave_radiation(i, j, kᴺ, grid, time, ϵ, ℐꜜˡʷ) - ΣQt = (ℐₜˢʷ + ℐₐˡʷ + ℐꜛˡʷ + 𝒬ᵀ + 𝒬ᵛ) * (ℵi > 0) # If ℵi == 0 there is no heat flux from the top! + # Atmosphere–sea-ice fluxes are computed at the ice surface + # (per unit ice area), hence we multiply by ℵ + ΣQt = (ℐₜˢʷ + ℐₐˡʷ + ℐꜛˡʷ + 𝒬ᵀ + 𝒬ᵛ) * ℵi ΣQb = 𝒬ᶠʳᶻ + 𝒬ⁱⁿᵗ # Mask fluxes over land for convenience diff --git a/test/test_conservation.jl b/test/test_conservation.jl new file mode 100644 index 000000000..150dd4ba0 --- /dev/null +++ b/test/test_conservation.jl @@ -0,0 +1,251 @@ +include("runtests_setup.jl") + +using Oceananigans.Units +using Oceananigans.Fields: interior +using Oceananigans.Operators: Azᶜᶜᶜ, volume +using Oceananigans.AbstractOperations: KernelFunctionOperation +using Oceananigans.TimeSteppers: time_step! + +using ClimaSeaIce +using ClimaSeaIce.SeaIceThermodynamics: latent_heat + +using NumericalEarth.Atmospheres: PrescribedAtmosphere +using NumericalEarth.Diagnostics: atmosphere_ocean_heat_flux, frazil_heat_flux +using NumericalEarth.EarthSystemModels: OceanSeaIceModel, Radiation, update_state! +using NumericalEarth.Oceans: ocean_simulation +using NumericalEarth.SeaIces: sea_ice_simulation + +# Diagnostic-side override: ClimaSeaIce's slab mass balance uses a +# T-dependent latent heat. A single state-based +# `Eᵢₛ = −ℵ · ρᵢ · ℒ · h · Az` cannot match both freeze at the ice-base +# temperature and top-melt at 0 ᵒC under `ℒ(T)`. Overriding `latent_heat` +# to the constant `pt.reference_latent_heat` isolates coupler-side +# bookkeeping errors from this intrinsic mismatch. The override is local to +# this test's process; it does not touch any package source. + +@inline function ClimaSeaIce.SeaIceThermodynamics.latent_heat( + pt::ClimaSeaIce.SeaIceThermodynamics.PhaseTransitions, T) + return pt.reference_latent_heat +end + +@testset "Coupled energy conservation" begin + + arch = CPU() + + grid = RectilinearGrid(arch; + size = 10, + z = (-100, 0), + topology = (Flat, Flat, Bounded)) + + ocean = ocean_simulation(grid; + momentum_advection = nothing, + tracer_advection = nothing, + closure = CATKEVerticalDiffusivity(), + coriolis = nothing, + bottom_drag_coefficient = 0) + + Tᵢ = -1.5 + Sᵢ = 34.0 + set!(ocean.model, T = Tᵢ, S = Sᵢ) + + sea_ice = sea_ice_simulation(grid, ocean; + dynamics = nothing, + advection = nothing, + ice_salinity = 0) + + set!(sea_ice.model, h = 1.0, ℵ = 1.0, hs = 0.10) + + atmosphere_grid = RectilinearGrid(arch; size=(), topology=(Flat, Flat, Flat)) + atmosphere = PrescribedAtmosphere(atmosphere_grid, [0.0, 1e9]) + radiation = Radiation() + + coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) + + function set_atmosphere!(atmosphere; + T_air, q_air, u_air, v_air, p_air, + SW_down, LW_down, rain_flux, snow_flux) + for (fts, value) in ((atmosphere.tracers.T, T_air), + (atmosphere.tracers.q, q_air), + (atmosphere.velocities.u, u_air), + (atmosphere.velocities.v, v_air), + (atmosphere.pressure, p_air), + (atmosphere.downwelling_radiation.shortwave, SW_down), + (atmosphere.downwelling_radiation.longwave, LW_down), + (atmosphere.freshwater_flux.rain, rain_flux), + (atmosphere.freshwater_flux.snow, snow_flux)) + fill!(parent(fts), value) + end + return nothing + end + + ρᵢ = coupled_model.sea_ice.model.sea_ice_density[1, 1, 1] + ρₛ = coupled_model.sea_ice.model.snow_density[1, 1, 1] + ℒ₀ = ClimaSeaIce.SeaIceThermodynamics.latent_heat(coupled_model.sea_ice.model.phase_transitions, 0.0) + ρᵒᶜ = coupled_model.interfaces.ocean_properties.reference_density + cᵒᶜ = coupled_model.interfaces.ocean_properties.heat_capacity + + Az = Azᶜᶜᶜ(1, 1, 1, grid) + Vᶜ = sum(KernelFunctionOperation{Center, Center, Center}(volume, grid, Center(), Center(), Center())) # for the virtual-salt FW diagnostic + + # Built once and re-used every step via `compute!`; the operations keep + # a reference to the live `T`/`S` fields so the reduction tracks them. + ∫T = Field(Integral(ocean.model.tracers.T)) + mean_S = Field(Average(ocean.model.tracers.S)) + + function column_state(coupled_model) + h = first(interior(coupled_model.sea_ice.model.ice_thickness)) + ℵ = first(interior(coupled_model.sea_ice.model.ice_concentration)) + hs = first(interior(coupled_model.sea_ice.model.snow_thickness)) + + Mᵢₛ = (ρᵢ * h + ρₛ * hs) * ℵ * Az + Eᵢₛ = - ℵ * (ρᵢ * ℒ₀ * h + ρₛ * ℒ₀ * hs) * Az + Hₒ = ρᵒᶜ * cᵒᶜ * first(compute!(∫T)) + Sₒ = first(compute!(mean_S)) + + return (; h, ℵ, hs, Mᵢₛ, Eᵢₛ, Hₒ, Sₒ) + end + + function net_top_heat_flux(coupled_model) + ΣQt = first(interior(coupled_model.interfaces.net_fluxes.sea_ice.top.heat)) + ΣQao = first(interior(atmosphere_ocean_heat_flux(coupled_model))) + return -(ΣQt + ΣQao) * Az + end + + ocean_virtual_freshwater(Sₒ, Sᵣ) = - ρᵒᶜ * Vᶜ * (Sₒ - Sᵣ) / Sᵣ + + # Short phases so the test stays under a few seconds of wall time. + Δt = 10minutes + Δτ = 2days # phase duration + Nsteps = Int(Δτ ÷ Δt) + + freeze_phase = (T_air = 253.15, + q_air = 1.0e-4, + u_air = 2.0, v_air = 0.0, + p_air = 101325.0, + SW_down = 50.0, + LW_down = 180.0, + rain_flux = 0.0, + snow_flux = 1.0e-5) + + melt_phase = (T_air = 278.15, + q_air = 5.0e-3, + u_air = 2.0, v_air = 0.0, + p_air = 101325.0, + SW_down = 250.0, + LW_down = 320.0, + rain_flux = 5.0e-6, + snow_flux = 0.0) + + history = (t = Float64[], + phase = Int[], + h = Float64[], ℵ = Float64[], hs = Float64[], + Mᵢₛ = Float64[], Eᵢₛ = Float64[], + Hₒ = Float64[], Sₒ = Float64[], + Ṁ = Float64[], Q = Float64[], + 𝒬ᶠʳᶻ = Float64[]) + + function record!(history, coupled_model, phase_id, Ṁ, Q) + st = column_state(coupled_model) + 𝒬f = first(interior(frazil_heat_flux(coupled_model))) + push!(history.t, coupled_model.clock.time) + push!(history.phase, phase_id) + push!(history.h, st.h) + push!(history.ℵ, st.ℵ) + push!(history.hs, st.hs) + push!(history.Mᵢₛ, st.Mᵢₛ) + push!(history.Eᵢₛ, st.Eᵢₛ) + push!(history.Hₒ, st.Hₒ) + push!(history.Sₒ, st.Sₒ) + push!(history.Ṁ, Ṁ) + push!(history.Q, Q) + push!(history.𝒬ᶠʳᶻ, 𝒬f) + return nothing + end + + record!(history, coupled_model, 0, 0.0, 0.0) + + # The phase-boundary `update_state!` would zero the pending frazil flux + # written at the end of the previous phase's final step, stranding the + # latent energy that was already deposited into the ocean by that same + # frazil mutation. We preserve `𝒬ᶠʳᶻ` across the refresh and add it + # back into the sea-ice bottom heat flux that the slab will read. + function run_phase!(coupled_model, spec, phase_id; Nsteps, Δt, history, atmosphere) + set_atmosphere!(atmosphere; + T_air = spec.T_air, q_air = spec.q_air, + u_air = spec.u_air, v_air = spec.v_air, + p_air = spec.p_air, + SW_down = spec.SW_down, LW_down = spec.LW_down, + rain_flux = spec.rain_flux, snow_flux = spec.snow_flux) + + Ṁ = (spec.rain_flux + spec.snow_flux) * Az # freshwater input + Qᵖ = - spec.snow_flux * ℒ₀ * Az # snowfall enthalpy + + 𝒬ᶠʳᶻ = coupled_model.interfaces.sea_ice_ocean_interface.fluxes.frazil_heat + ΣQb = coupled_model.interfaces.net_fluxes.sea_ice.bottom.heat + 𝒬⁻ = first(interior(𝒬ᶠʳᶻ)) # pending frazil + update_state!(coupled_model) + 𝒬ᶠʳᶻ[1, 1, 1] = 𝒬⁻ + ΣQb[1, 1, 1] += 𝒬⁻ + + history.Q[end] = net_top_heat_flux(coupled_model) + Qᵖ + history.Ṁ[end] = Ṁ + + for _ in 1:Nsteps + time_step!(coupled_model, Δt) + Q = net_top_heat_flux(coupled_model) + Qᵖ + record!(history, coupled_model, phase_id, Ṁ, Q) + end + end + + run_phase!(coupled_model, freeze_phase, 1; Nsteps, Δt, history, atmosphere) + run_phase!(coupled_model, melt_phase, 2; Nsteps, Δt, history, atmosphere) + + t = history.t + + # --- Energy budget --- + + ∫Q = similar(t); ∫Q[1] = 0.0 + for n in 2:length(t) + ∫Q[n] = ∫Q[n-1] + history.Q[n-1] * (t[n] - t[n-1]) + end + + Δt⁺ = similar(t) + for n in 1:(length(t) - 1) + Δt⁺[n] = t[n+1] - t[n] + end + Δt⁺[end] = Δt⁺[end-1] + δE = history.𝒬ᶠʳᶻ .* Δt⁺ .* Az + + Ẽᵢₛ = history.Eᵢₛ .+ δE + ΔE = (Ẽᵢₛ .+ history.Hₒ) .- (Ẽᵢₛ[1] + history.Hₒ[1]) + Rₑ = ΔE .- ∫Q + + εₑ = abs(Rₑ[end]) / max(maximum(abs.(ΔE)), 1) + @test εₑ < 1e-10 + + # --- Freshwater budget --- + # + # TODO: the coupled freshwater budget does not currently close. The + # numbers below are computed for reference / plotting purposes but we + # do not assert on them; fix the underlying bookkeeping first, then + # add `@test εₘ < 1e-10`. + + ∫Ṁ = similar(t); ∫Ṁ[1] = 0.0 + for n in 2:length(t) + ∫Ṁ[n] = ∫Ṁ[n-1] + history.Ṁ[n-1] * (t[n] - t[n-1]) + end + + Sᵣ = history.Sₒ[1] + Mᶠʷ = @. ocean_virtual_freshwater(history.Sₒ, Sᵣ) + ΔM = (history.Mᵢₛ .+ Mᶠʷ) .- (history.Mᵢₛ[1] + Mᶠʷ[1]) + Rₘ = ΔM .- ∫Ṁ + εₘ = abs(Rₘ[end]) / max(maximum(abs.(ΔM)), 1) # not tested (broken) + + # --- Sanity checks on the physics --- + + nᶠ = findlast(p -> p == 1, history.phase) + + @test history.h[nᶠ] > history.h[1] # ice grew during freeze + @test history.h[end] < history.h[nᶠ] # ice shrank during melt + @test history.ℵ[nᶠ] ≈ 1.0 # stays fully ice-covered through freeze +end From 2a5b49c4f22f5f9adf49a0c5cfec9ed043efc1a7 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 23 Apr 2026 18:46:24 +0200 Subject: [PATCH 2/7] add the example to make --- docs/make.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/make.jl b/docs/make.jl index f3f332f8d..3203efa48 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -30,6 +30,7 @@ mkpath(OUTPUT_DIR) # on pushes to `main`/tags, or when the `build all examples` label is added to a PR. examples = [ Example("Single-column ocean simulation", "single_column_os_papa_simulation", true), + Example("Coupled energy and freshwater conservation", "coupled_conservation", true), Example("One-degree ocean--sea ice simulation", "one_degree_simulation", false), Example("Near-global ocean simulation", "near_global_ocean_simulation", false), Example("Global climate simulation", "global_climate_simulation", false), From 88ed1810b5148e0f46865763557d604d0415aee9 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 30 Apr 2026 17:46:17 +1000 Subject: [PATCH 3/7] Apply suggestions from code review Co-authored-by: Navid C. Constantinou --- examples/coupled_conservation.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/examples/coupled_conservation.jl b/examples/coupled_conservation.jl index ac1c87f0e..34b8e9707 100644 --- a/examples/coupled_conservation.jl +++ b/examples/coupled_conservation.jl @@ -48,10 +48,8 @@ using Printf # locally override `latent_heat` to the constant `pt.reference_latent_heat`. This is a diagnostic choice for the present # example and does not modify upstream. -@inline function ClimaSeaIce.SeaIceThermodynamics.latent_heat( - pt::ClimaSeaIce.SeaIceThermodynamics.PhaseTransitions, T) - return pt.reference_latent_heat -end +@inline ClimaSeaIce.SeaIceThermodynamics.latent_heat(pt::ClimaSeaIce.SeaIceThermodynamics.PhaseTransitions, T) = + pt.reference_latent_heat # ## Grid, ocean, sea ice, and atmosphere # From e115c3bdb72ef26a05ae111232b3540ab728e8b6 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 30 Apr 2026 17:52:05 +1000 Subject: [PATCH 4/7] Refactor comments for consistency and clarity Updated comments to use consistent multiplication syntax and improved clarity. --- examples/coupled_conservation.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/examples/coupled_conservation.jl b/examples/coupled_conservation.jl index 34b8e9707..d4a0694e3 100644 --- a/examples/coupled_conservation.jl +++ b/examples/coupled_conservation.jl @@ -24,7 +24,6 @@ using Oceananigans using Oceananigans.Units -using Oceananigans.Fields: interior using Oceananigans.Operators: Azᶜᶜᶜ using ClimaSeaIce @@ -41,9 +40,9 @@ using Printf # ## Constant latent heat for diagnostic closure # -# `ClimaSeaIce`'s slab mass balance uses a temperature-dependent latent heat, `ℒ(T) = ℒ₀ + (ρℓ cℓ / ρi − cᵢ)(T − T₀)`, -# with `ℰu = ρi · ℒ(T_u)` at the top interface and `ℰb = ρi · ℒ(T_b)` at the bottom. A single state-based -# `Eis = − ℵ · ρi · ℒ · h · Az` cannot simultaneously match freeze at `T_b` and top-melt at 0 ᵒC: the 4.7 kJ/kg gap accounts +# `ClimaSeaIce`'s slab mass balance uses a temperature-dependent latent heat, `ℒ(T) = ℒ₀ + (ρℓ * cℓ / ρi − cᵢ)(T − T₀)`, +# with `ℰu = ρi * ℒ(T_u)` at the top interface and `ℰb = ρi * ℒ(T_b)` at the bottom. A single state-based +# `Eis = − ℵ * ρi * ℒ * h * Az` cannot simultaneously match freeze at `T_b` and top-melt at 0 ᵒC: the 4.7 kJ/kg gap accounts # for a ~1% residual scaling with top-melt mass. To isolate coupler-side bookkeeping from this intrinsic `ℒ(T)` mismatch we # locally override `latent_heat` to the constant `pt.reference_latent_heat`. This is a diagnostic choice for the present # example and does not modify upstream. @@ -141,7 +140,7 @@ function column_state(coupled_model) end # `net_top_heat_flux` returns the atmospheric energy input to the coupled (ice + ocean) system in Watts: -# `Q_atm = − (ΣQt + ΣQao) · Az`, where `ΣQt` is the sea-ice top heat flux per cell and `ΣQao` is the per-cell +# `Q_atm = − (ΣQt + ΣQao) * Az`, where `ΣQt` is the sea-ice top heat flux per cell and `ΣQao` is the per-cell # atmosphere-to-ocean flux over the open-water fraction. The ocean-side piece comes from `atmosphere_ocean_heat_flux`, # which subtracts the frazil and interface contributions internally so it never picks up a spurious ocean / ice exchange term. @@ -281,7 +280,7 @@ end # The frazil mass gain is deposited by `compute_sea_ice_ocean_fluxes!` at the end of step `n` # (mutating ocean `T` and writing `𝒬ᶠʳᶻ`) but the corresponding ice mass gain is consumed only during # step `n + 1`. At a diagnostic snapshot the ocean shows the warming while the ice has not yet grown. -# We anticipate this one-step pending quantity by adding `𝒬ᶠʳᶻ(n) · Δt⁺ · Az` to `Eis(n)` so the energy +# We anticipate this one-step pending quantity by adding `𝒬ᶠʳᶻ(n) * Δt⁺ * Az` to `Eis(n)` so the energy # budget closure is not polluted by bookkeeping lag. Δt⁺ = similar(t) From 43042a86529c3f7616f9d49dbd2047d307777d84 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 30 Apr 2026 18:00:06 +1000 Subject: [PATCH 5/7] Refactor variable assignments and improve formatting --- test/test_conservation.jl | 37 +++++++++++++++++++++---------------- 1 file changed, 21 insertions(+), 16 deletions(-) diff --git a/test/test_conservation.jl b/test/test_conservation.jl index 150dd4ba0..0fa1dcbe8 100644 --- a/test/test_conservation.jl +++ b/test/test_conservation.jl @@ -89,7 +89,7 @@ end # Built once and re-used every step via `compute!`; the operations keep # a reference to the live `T`/`S` fields so the reduction tracks them. - ∫T = Field(Integral(ocean.model.tracers.T)) + ∫T = Field(Integral(ocean.model.tracers.T)) mean_S = Field(Average(ocean.model.tracers.S)) function column_state(coupled_model) @@ -115,7 +115,7 @@ end # Short phases so the test stays under a few seconds of wall time. Δt = 10minutes - Δτ = 2days # phase duration + Δτ = 2days # phase duration Nsteps = Int(Δτ ÷ Δt) freeze_phase = (T_air = 253.15, @@ -136,12 +136,17 @@ end rain_flux = 5.0e-6, snow_flux = 0.0) - history = (t = Float64[], + history = (t = Float64[], phase = Int[], - h = Float64[], ℵ = Float64[], hs = Float64[], - Mᵢₛ = Float64[], Eᵢₛ = Float64[], - Hₒ = Float64[], Sₒ = Float64[], - Ṁ = Float64[], Q = Float64[], + h = Float64[], + ℵ = Float64[], + hs = Float64[], + Mᵢₛ = Float64[], + Eᵢₛ = Float64[], + Hₒ = Float64[], + Sₒ = Float64[], + Ṁ = Float64[], + Q = Float64[], 𝒬ᶠʳᶻ = Float64[]) function record!(history, coupled_model, phase_id, Ṁ, Q) @@ -184,8 +189,8 @@ end ΣQb = coupled_model.interfaces.net_fluxes.sea_ice.bottom.heat 𝒬⁻ = first(interior(𝒬ᶠʳᶻ)) # pending frazil update_state!(coupled_model) - 𝒬ᶠʳᶻ[1, 1, 1] = 𝒬⁻ - ΣQb[1, 1, 1] += 𝒬⁻ + 𝒬ᶠʳᶻ[1, 1, 1] = 𝒬⁻ + ΣQb[1, 1, 1] += 𝒬⁻ history.Q[end] = net_top_heat_flux(coupled_model) + Qᵖ history.Ṁ[end] = Ṁ @@ -236,16 +241,16 @@ end end Sᵣ = history.Sₒ[1] - Mᶠʷ = @. ocean_virtual_freshwater(history.Sₒ, Sᵣ) - ΔM = (history.Mᵢₛ .+ Mᶠʷ) .- (history.Mᵢₛ[1] + Mᶠʷ[1]) - Rₘ = ΔM .- ∫Ṁ - εₘ = abs(Rₘ[end]) / max(maximum(abs.(ΔM)), 1) # not tested (broken) + Mᶠʷ = @. ocean_virtual_freshwater(history.Sₒ, Sᵣ) + ΔM = (history.Mᵢₛ .+ Mᶠʷ) .- (history.Mᵢₛ[1] + Mᶠʷ[1]) + Rₘ = ΔM .- ∫Ṁ + εₘ = abs(Rₘ[end]) / max(maximum(abs.(ΔM)), 1) # not tested (broken) # --- Sanity checks on the physics --- nᶠ = findlast(p -> p == 1, history.phase) - @test history.h[nᶠ] > history.h[1] # ice grew during freeze - @test history.h[end] < history.h[nᶠ] # ice shrank during melt - @test history.ℵ[nᶠ] ≈ 1.0 # stays fully ice-covered through freeze + @test history.h[nᶠ] > history.h[1] # ice grew during freeze + @test history.h[end] < history.h[nᶠ] # ice shrank during melt + @test history.ℵ[nᶠ] ≈ 1.0 # stays fully ice-covered through freeze end From 79aa4644e663ac3bbd881cf7cf3370ba806588f9 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 30 Apr 2026 18:03:46 +1000 Subject: [PATCH 6/7] =?UTF-8?q?Refactor=20integral=20initialization=20for?= =?UTF-8?q?=20Q=20and=20=E1=B9=80?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- test/test_conservation.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/test/test_conservation.jl b/test/test_conservation.jl index 0fa1dcbe8..d7794c683 100644 --- a/test/test_conservation.jl +++ b/test/test_conservation.jl @@ -209,7 +209,8 @@ end # --- Energy budget --- - ∫Q = similar(t); ∫Q[1] = 0.0 + ∫Q = similar(t) + ∫Q[1] = 0.0 for n in 2:length(t) ∫Q[n] = ∫Q[n-1] + history.Q[n-1] * (t[n] - t[n-1]) end @@ -222,8 +223,8 @@ end δE = history.𝒬ᶠʳᶻ .* Δt⁺ .* Az Ẽᵢₛ = history.Eᵢₛ .+ δE - ΔE = (Ẽᵢₛ .+ history.Hₒ) .- (Ẽᵢₛ[1] + history.Hₒ[1]) - Rₑ = ΔE .- ∫Q + ΔE = (Ẽᵢₛ .+ history.Hₒ) .- (Ẽᵢₛ[1] + history.Hₒ[1]) + Rₑ = ΔE .- ∫Q εₑ = abs(Rₑ[end]) / max(maximum(abs.(ΔE)), 1) @test εₑ < 1e-10 @@ -235,7 +236,8 @@ end # do not assert on them; fix the underlying bookkeeping first, then # add `@test εₘ < 1e-10`. - ∫Ṁ = similar(t); ∫Ṁ[1] = 0.0 + ∫Ṁ = similar(t) + ∫Ṁ[1] = 0.0 for n in 2:length(t) ∫Ṁ[n] = ∫Ṁ[n-1] + history.Ṁ[n-1] * (t[n] - t[n-1]) end From f34aa844b39b6f33707b2654c8f3e7207baa0add Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 30 Apr 2026 18:04:11 +1000 Subject: [PATCH 7/7] Fix comment formatting in test_conservation.jl --- test/test_conservation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_conservation.jl b/test/test_conservation.jl index d7794c683..85b10de4f 100644 --- a/test/test_conservation.jl +++ b/test/test_conservation.jl @@ -17,7 +17,7 @@ using NumericalEarth.SeaIces: sea_ice_simulation # Diagnostic-side override: ClimaSeaIce's slab mass balance uses a # T-dependent latent heat. A single state-based -# `Eᵢₛ = −ℵ · ρᵢ · ℒ · h · Az` cannot match both freeze at the ice-base +# `Eᵢₛ = −ℵ * ρᵢ * ℒ * h * Az` cannot match both freeze at the ice-base # temperature and top-melt at 0 ᵒC under `ℒ(T)`. Overriding `latent_heat` # to the constant `pt.reference_latent_heat` isolates coupler-side # bookkeeping errors from this intrinsic mismatch. The override is local to