diff --git a/Project.toml b/Project.toml index e33fb9936..6fbdea39e 100644 --- a/Project.toml +++ b/Project.toml @@ -43,6 +43,9 @@ SpeedyWeather = "9e226e20-d153-4fed-8a5b-493def4f21a9" WorldOceanAtlasTools = "04f20302-f1b9-11e8-29d9-7d841cb0a64a" XESMF = "2e0b0046-e7a1-486f-88de-807ee8ffabe5" +[sources] +ClimaSeaIce = {rev = "ss/refactor-thermodynamics", url = "https://github.com/CliMA/ClimaSeaIce.jl.git"} + [extensions] NumericalEarthBreezeExt = "Breeze" NumericalEarthCDSAPIExt = "CDSAPI" diff --git a/docs/Project.toml b/docs/Project.toml index 82347bf37..6ba01cd44 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -21,8 +21,8 @@ SpeedyWeather = "9e226e20-d153-4fed-8a5b-493def4f21a9" XESMF = "2e0b0046-e7a1-486f-88de-807ee8ffabe5" [sources] +ClimaSeaIce = {url = "https://github.com/CliMA/ClimaSeaIce.jl", rev = "ss/refactor-thermodynamics"} NumericalEarth = {path = ".."} -Breeze = {url = "https://github.com/NumericalEarth/Breeze.jl/", rev = "main"} [compat] Documenter = "1" diff --git a/docs/src/NumericalEarth.bib b/docs/src/NumericalEarth.bib index 84e4628c2..f97dcd243 100644 --- a/docs/src/NumericalEarth.bib +++ b/docs/src/NumericalEarth.bib @@ -84,9 +84,9 @@ @article{holland1999modeling doi={10.1175/1520-0485(1999)029<1787:MTIOIA>2.0.CO;2} } -@article{hieronymus2021comparison, - title={A comparison of ocean-ice flux parametrizations}, - author={Hieronymus, Magnus and Holtermann, Peter and Gr{\"a}we, Ulf and Burchard, Hans}, +@article{shi2021sensitivity, + title={Sensitivity of {Northern Hemisphere} climate to ice--ocean interface heat flux parameterizations}, + author={Shi, Xiaoxu and Notz, Dirk and Liu, Jiping and Yang, Hu and Lohmann, Gerrit}, journal={Geoscientific Model Development}, volume={14}, number={8}, diff --git a/docs/src/interface_fluxes.md b/docs/src/interface_fluxes.md index b776bc475..2d0714510 100644 --- a/docs/src/interface_fluxes.md +++ b/docs/src/interface_fluxes.md @@ -817,7 +817,7 @@ where: - ``\lambda_1, \lambda_2`` are liquidus coefficients The ratio ``R = \alpha_h / \alpha_s`` (typically around 35) reflects the different molecular diffusivities of heat and -salt, with heat diffusing faster than salt [hieronymus2021comparison](@citep). +salt, with heat diffusing faster than salt [shi2021sensitivity](@citep). ```@example interface_fluxes using NumericalEarth.EarthSystemModels: ThreeEquationHeatFlux @@ -975,4 +975,4 @@ Note: The `ComponentInterfaces` call above is illustrative; it requires fully co The implementations follow: - [holland1999modeling](@citet): foundational three-equation model for ice shelf-ocean interaction -- [hieronymus2021comparison](@citet): comparison of different ocean-ice flux parameterizations +- [shi2021sensitivity](@citet): sensitivity of Northern Hemisphere climate to ice-ocean interface heat flux parameterizations diff --git a/src/Atmospheres/interpolate_atmospheric_state.jl b/src/Atmospheres/interpolate_atmospheric_state.jl index 56cf48727..e0c8248a7 100644 --- a/src/Atmospheres/interpolate_atmospheric_state.jl +++ b/src/Atmospheres/interpolate_atmospheric_state.jl @@ -31,6 +31,7 @@ function interpolate_state!(exchanger, grid, atmosphere::PrescribedAtmosphere, c ℐꜜˡʷ = atmosphere.downwelling_radiation.longwave downwelling_radiation = (shortwave=ℐꜜˢʷ.data, longwave=ℐꜜˡʷ.data) freshwater_flux = map(ϕ -> ϕ.data, atmosphere.freshwater_flux) + snowfall_flux = haskey(atmosphere.freshwater_flux, :snow) ? atmosphere.freshwater_flux.snow.data : nothing atmosphere_pressure = atmosphere.pressure.data # Extract info for time-interpolation @@ -51,7 +52,8 @@ function interpolate_state!(exchanger, grid, atmosphere::PrescribedAtmosphere, c q = atmosphere_fields.q.data, ℐꜜˢʷ = atmosphere_fields.ℐꜜˢʷ.data, ℐꜜˡʷ = atmosphere_fields.ℐꜜˡʷ.data, - Jᶜ = atmosphere_fields.Jᶜ.data) + Jᶜ = atmosphere_fields.Jᶜ.data, + Jˢⁿ = atmosphere_fields.Jˢⁿ.data) kernel_parameters = interface_kernel_parameters(grid) @@ -74,6 +76,7 @@ function interpolate_state!(exchanger, grid, atmosphere::PrescribedAtmosphere, c atmosphere_pressure, downwelling_radiation, freshwater_flux, + snowfall_flux, atmosphere_backend, atmosphere_time_indexing) @@ -101,6 +104,7 @@ end atmos_pressure, downwelling_radiation, prescribed_freshwater_flux, + snowfall_flux, atmos_backend, atmos_time_indexing) @@ -124,9 +128,12 @@ end ℐꜜˢʷ = interp_atmos_time_series(downwelling_radiation.shortwave, atmos_args...) ℐꜜˡʷ = interp_atmos_time_series(downwelling_radiation.longwave, atmos_args...) - # Usually precipitation + # Total precipitation (rain + snow) Mh = interp_atmos_time_series(prescribed_freshwater_flux, atmos_args...) + # Snowfall only (for sea ice snow accumulation) + Ms = interp_atmos_time_series(snowfall_flux, atmos_args...) + # Convert atmosphere velocities (usually defined on a latitude-longitude grid) to # the frame of reference of the native grid kᴺ = size(exchange_grid, 3) # index of the top ocean cell @@ -141,6 +148,7 @@ end surface_atmos_state.ℐꜜˢʷ[i, j, 1] = ℐꜜˢʷ surface_atmos_state.ℐꜜˡʷ[i, j, 1] = ℐꜜˡʷ surface_atmos_state.Jᶜ[i, j, 1] = Mh + surface_atmos_state.Jˢⁿ[i, j, 1] = Ms end end @@ -148,9 +156,11 @@ end ##### Utility for interpolating tuples of fields ##### +@inline interp_atmos_time_series(::Nothing, X, time, grid, args...) = 0 + # Note: assumes loc = (c, c, nothing) (and the third location should not matter.) -@inline interp_atmos_time_series(J::AbstractArray, x_itp::FractionalIndices, t_itp, args...) = - interpolate(x_itp, t_itp, J, args...) +@inline interp_atmos_time_series(J::AbstractArray, X::FractionalIndices, time, args...) = + interpolate(X, time, J, args...) @inline interp_atmos_time_series(J::AbstractArray, X, time, grid, args...) = interpolate(X, time, J, (Center(), Center(), nothing), grid, args...) diff --git a/src/Atmospheres/prescribed_atmosphere_regridder.jl b/src/Atmospheres/prescribed_atmosphere_regridder.jl index da93d1766..2c1a8aec9 100644 --- a/src/Atmospheres/prescribed_atmosphere_regridder.jl +++ b/src/Atmospheres/prescribed_atmosphere_regridder.jl @@ -9,7 +9,8 @@ function ComponentExchanger(atmosphere::PrescribedAtmosphere, grid) q = Field{Center, Center, Nothing}(grid), ℐꜜˢʷ = Field{Center, Center, Nothing}(grid), ℐꜜˡʷ = Field{Center, Center, Nothing}(grid), - Jᶜ = Field{Center, Center, Nothing}(grid)) + Jᶜ = Field{Center, Center, Nothing}(grid), + Jˢⁿ = Field{Center, Center, Nothing}(grid)) return ComponentExchanger(state, regridder) end diff --git a/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl index 1292d6311..317b337b9 100644 --- a/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl @@ -84,7 +84,8 @@ end # Sea ice properties uˢⁱ = zero(FT) # ℑxᶜᵃᵃ(i, j, 1, grid, interior_state.u) vˢⁱ = zero(FT) # ℑyᵃᶜᵃ(i, j, 1, grid, interior_state.v) - hˢⁱ = interior_state.h[i, j, 1] + hˢⁱ = interior_state.hi[i, j, 1] + hˢⁿ = interior_state.hs[i, j, 1] hc = interior_state.hc[i, j, 1] ℵᵢ = interior_state.ℵ[i, j, 1] Tₛ = interface_temperature[i, j, 1] @@ -92,10 +93,8 @@ end end # Build thermodynamic and dynamic states in the atmosphere and interface. - # Notation: - # ⋅ 𝒰 ≡ "dynamic" state vector (thermodynamics + reference height + velocity) ℂᵃᵗ = atmosphere_properties.thermodynamics_parameters - zᵃᵗ = atmosphere_properties.surface_layer_height # elevation of atmos variables relative to interface + zᵃᵗ = atmosphere_properties.surface_layer_height local_atmosphere_state = (z = zᵃᵗ, u = uᵃᵗ, @@ -106,7 +105,7 @@ end h_bℓ = atmosphere_state.h_bℓ) downwelling_radiation = (; ℐꜜˢʷ, ℐꜜˡʷ) - local_interior_state = (u=uˢⁱ, v=vˢⁱ, T=Tᵒᶜ, S=Sᵒᶜ, h=hˢⁱ, hc=hc) + local_interior_state = (u=uˢⁱ, v=vˢⁱ, T=Tᵒᶜ, S=Sᵒᶜ, hi=hˢⁱ, hs=hˢⁿ, hc=hc) # Estimate initial interface state (FP32 compatible) u★ = convert(FT, 1f-4) diff --git a/src/EarthSystemModels/InterfaceComputations/component_interfaces.jl b/src/EarthSystemModels/InterfaceComputations/component_interfaces.jl index 65d6adbe8..fa4efa420 100644 --- a/src/EarthSystemModels/InterfaceComputations/component_interfaces.jl +++ b/src/EarthSystemModels/InterfaceComputations/component_interfaces.jl @@ -281,7 +281,14 @@ function atmosphere_sea_ice_interface(grid, temperature_formulation, velocity_formulation) - interface_temperature = sea_ice.model.ice_thermodynamics.top_surface_temperature + # When snow is present, the atmosphere interacts with the snow surface; + # otherwise with the ice top surface. + snow_thermo = sea_ice.model.snow_thermodynamics + interface_temperature = if isnothing(snow_thermo) + sea_ice.model.ice_thermodynamics.top_surface_temperature + else + snow_thermo.top_surface_temperature + end return AtmosphereInterface(fluxes, ai_flux_formulation, interface_temperature, properties) end @@ -419,7 +426,7 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; sea_ice_properties = (reference_density = sea_ice_reference_density, heat_capacity = sea_ice_heat_capacity, freshwater_density = freshwater_density, - liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus, + liquidus = sea_ice.model.phase_transitions.liquidus, temperature_units = sea_ice_temperature_units) else sea_ice_properties = nothing diff --git a/src/EarthSystemModels/InterfaceComputations/interface_states.jl b/src/EarthSystemModels/InterfaceComputations/interface_states.jl index 4f1ae6c2d..55c76245b 100644 --- a/src/EarthSystemModels/InterfaceComputations/interface_states.jl +++ b/src/EarthSystemModels/InterfaceComputations/interface_states.jl @@ -255,60 +255,87 @@ end Jᵀ = Qa * λ # Calculating the atmospheric temperature - # We use to compute the sensible heat flux Tᵃᵗ = surface_atmosphere_temperature(Ψₐ, ℙₐ) ΔT = Tᵃᵗ - Ψₛ.T - Ωc = ifelse(ΔT == 0, zero(ΔT), 𝒬ᵀ / ΔT * λ) # Sensible heat transfer coefficient (W/m²K) - # Computing the flux balance temperature - return (Ψᵢ.T * F.κ - (Jᵀ + Ωc * Tᵃᵗ) * F.δ) / (F.κ - Ωc * F.δ) + # Flux balance: T★ = (Tᵢ κ - (Jᵀ + Ωc Tᵃᵗ) δ) / (κ - Ωc δ) + # where Ωc = 𝒬ᵀ λ / ΔT. Multiply through by ΔT to avoid Inf when ΔT → 0. + Ωᵀ = 𝒬ᵀ * λ # unnormalized sensible heat coefficient (= Ωc * ΔT) + D = F.κ * ΔT - Ωᵀ * F.δ + T★ = (Ψᵢ.T * F.κ * ΔT - (Jᵀ * ΔT + Ωᵀ * Tᵃᵗ) * F.δ) / D + + return ifelse(D == 0, Ψₛ.T, T★) end -# 𝒬ᵛ + ℐꜛˡʷ + Qd + Ωc * (Tᵃᵗ - Tˢ) + k / h * (Tˢ - Tˢⁱ) = 0 -# where Ωc (the sensible heat transfer coefficient) is given by Ωc = 𝒬ᵀ / (Tᵃᵗ - Tˢ) -# ⟹ Tₛ = (Tˢⁱ * k - (𝒬ᵛ + ℐꜛˡʷ + Qd + Ωc * Tᵃᵗ) * h / (k - Ωc * h) -@inline function flux_balance_temperature(st::SkinTemperature{<:ClimaSeaIce.ConductiveFlux}, Ψₛ, ℙₛ, 𝒬ᵀ, 𝒬ᵛ, ℐꜛˡʷ, Qd, Ψᵢ, ℙᵢ, Ψₐ, ℙₐ) - F = st.internal_flux - k = F.conductivity - h = Ψᵢ.h - hc = Ψᵢ.hc # Critical thickness for ice consolidation - - # Bottom temperature at the melting temperature - Tˢⁱ = ClimaSeaIce.SeaIceThermodynamics.melting_temperature(ℙᵢ.liquidus, Ψᵢ.S) - Tˢⁱ = convert_to_kelvin(ℙᵢ.temperature_units, Tˢⁱ) +# Solve the surface flux balance equation: +# Qa(Tₛ) + Ωc (Tᵃᵗ - Tₛ) + (Tₛ - Tᵦ) / R = 0 +# where R is the total thermal resistance (h/k for bare ice, hₛ/kₛ + hᵢ/kᵢ with snow), +# Ωc = 𝒬ᵀ/(Tᵃᵗ-Tₛ) is the linearized sensible heat coefficient, and Qa = 𝒬ᵛ + ℐꜛˡʷ + Qd. +# The upward longwave ℐꜛˡʷ = σ ε Tₛ⁴ is strongly nonlinear in Tₛ; a pure Picard +# iteration (treating Qa constant) is unstable when 4σεTₛ³ ≳ 1/R (radiation +# dominated). We linearize: Qa(Tₛ) ≈ Qa(Tₛ⁻) + β (Tₛ − Tₛ⁻) with β = 4σεTₛ⁻³, +# yielding the Newton-like semi-implicit update: +# Tₛ = [Tᵦ + β R Tₛ⁻ - Ωc R Tᵃᵗ - Qa R] / [1 + β R - Ωc R] +@inline function conductive_flux_balance_temperature(st, R, hᵢ, Ψₛ, ℙₛ, 𝒬ᵀ, 𝒬ᵛ, ℐꜛˡʷ, Qd, Ψᵢ, ℙᵢ, Ψₐ, ℙₐ) + hc = Ψᵢ.hc + + # Bottom temperature at the melting point + Tᵦ = ClimaSeaIce.SeaIceThermodynamics.melting_temperature(ℙᵢ.liquidus, Ψᵢ.S) + Tᵦ = convert_to_kelvin(ℙᵢ.temperature_units, Tᵦ) Tₛ⁻ = Ψₛ.T - # Calculating the atmospheric temperature - # We use to compute the sensible heat flux Tᵃᵗ = surface_atmosphere_temperature(Ψₐ, ℙₐ) ΔT = Tᵃᵗ - Tₛ⁻ - Ωc = ifelse(ΔT == 0, zero(h), 𝒬ᵀ / ΔT) # Sensible heat transfer coefficient (W/m²K) - Qa = (𝒬ᵛ + ℐꜛˡʷ + Qd) # Net flux excluding sensible heat (positive out of the ocean) + Qa = 𝒬ᵛ + ℐꜛˡʷ + Qd + + # Sensible transfer coefficient Ωc = 𝒬ᵀ/ΔT, safely handling ΔT → 0. + Ωc = ifelse(ΔT == zero(ΔT), zero(Tₛ⁻), 𝒬ᵀ / ΔT) - # Computing the flux balance temperature - T★ = (Tˢⁱ * k - (Qa + Ωc * Tᵃᵗ) * h) / (k - Ωc * h) + # Newton linearization of upwelling longwave: ℐꜛˡʷ(Tₛ) ≈ ℐꜛˡʷ(Tₛ⁻) + β (Tₛ − Tₛ⁻). + σ = ℙₛ.radiation.σ + ϵ = ℙₛ.radiation.ϵ + β = 4 * σ * ϵ * Tₛ⁻^3 - # Fix a NaN - T★ = ifelse(isnan(T★), Tₛ⁻, T★) + # Flux balance solution with T⁴ linearization (stable even at ΔT = 0): + D = 1 + β * R - Ωc * R + T★ = (Tᵦ + β * R * Tₛ⁻ - Ωc * R * Tᵃᵗ - Qa * R) / D + T★ = ifelse(D == 0, Tₛ⁻, T★) - # To prevent instabilities in the fixed point iteration - # solver we cap the maximum temperature difference with `max_ΔT` + # Cap the temperature step for iteration stability ΔT★ = T★ - Tₛ⁻ max_ΔT = convert(typeof(T★), st.max_ΔT) - abs_ΔT = min(max_ΔT, abs(ΔT★)) - Tₛ⁺ = Tₛ⁻ + abs_ΔT * sign(ΔT★) + Tₛ⁺ = Tₛ⁻ + clamp(ΔT★, -max_ΔT, max_ΔT) - # Under heating fluxes, cap surface temperature by melting temperature + # Cap at melting temperature Tₘ = ℙᵢ.liquidus.freshwater_melting_temperature Tₘ = convert_to_kelvin(ℙᵢ.temperature_units, Tₘ) Tₛ⁺ = min(Tₛ⁺, Tₘ) - # If the ice is not consolidated, use the bottom temperature - Tₛ⁺ = ifelse(h ≥ hc, Tₛ⁺, Tˢⁱ) - + # If ice is not consolidated, use the bottom temperature + Tₛ⁺ = ifelse(hᵢ ≥ hc, Tₛ⁺, Tᵦ) + return Tₛ⁺ end +# Bare ice: R = hᵢ / kᵢ +@inline function flux_balance_temperature(st::SkinTemperature{<:ClimaSeaIce.ConductiveFlux}, + Ψₛ, ℙₛ, 𝒬ᵀ, 𝒬ᵛ, ℐꜛˡʷ, Qd, Ψᵢ, ℙᵢ, Ψₐ, ℙₐ) + k = st.internal_flux.conductivity + hᵢ = Ψᵢ.hi + R = hᵢ / k + return conductive_flux_balance_temperature(st, R, hᵢ, Ψₛ, ℙₛ, 𝒬ᵀ, 𝒬ᵛ, ℐꜛˡʷ, Qd, Ψᵢ, ℙᵢ, Ψₐ, ℙₐ) +end + +# Snow + ice: R = hₛ / kₛ + hᵢ / kᵢ +@inline function flux_balance_temperature(st::SkinTemperature{<:ClimaSeaIce.SeaIceThermodynamics.IceSnowConductiveFlux}, + Ψₛ, ℙₛ, 𝒬ᵀ, 𝒬ᵛ, ℐꜛˡʷ, Qd, Ψᵢ, ℙᵢ, Ψₐ, ℙₐ) + F = st.internal_flux + hᵢ = Ψᵢ.hi + hₛ = Ψᵢ.hs + R = hₛ / F.snow_conductivity + hᵢ / F.ice_conductivity + return conductive_flux_balance_temperature(st, R, hᵢ, Ψₛ, ℙₛ, 𝒬ᵀ, 𝒬ᵛ, ℐꜛˡʷ, Qd, Ψᵢ, ℙᵢ, Ψₐ, ℙₐ) +end + @inline function compute_interface_temperature(st::SkinTemperature, interface_state, atmosphere_state, diff --git a/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_fluxes.jl index a715493d3..f8d5e30a4 100644 --- a/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_fluxes.jl @@ -36,7 +36,7 @@ function compute_sea_ice_ocean_fluxes!(interface, ocean, sea_ice, ocean_properti hˢⁱ = sea_ice.model.ice_thickness hc = sea_ice.model.ice_consolidation_thickness - phase_transitions = sea_ice.model.ice_thermodynamics.phase_transitions + phase_transitions = sea_ice.model.phase_transitions liquidus = phase_transitions.liquidus L = phase_transitions.reference_latent_heat @@ -175,12 +175,12 @@ end qᶠ = δ𝒬ᶠʳᶻ / ℰ @inbounds begin - Tᴺ = Tᵒᶜ[i, j, Nz] - Sᴺ = Sᵒᶜ[i, j, Nz] + Tᴺ = Tᵒᶜ[i, j, Nz] + Sᴺ = Sᵒᶜ[i, j, Nz] Sˢⁱ = ice_salinity[i, j, 1] hˢⁱ = ice_thickness[i, j, 1] - ℵᵢ = ice_concentration[i, j, 1] - hc = ice_consolidation_thickness[i, j, 1] + ℵᵢ = ice_concentration[i, j, 1] + hc = ice_consolidation_thickness[i, j, 1] end # Extract internal temperature (for ConductiveFluxTEF, zero otherwise) @@ -198,8 +198,8 @@ end # ============================================= # Returns interfacial heat flux, melt rate qᵐ, and interface T, S 𝒬ⁱᵒ, qᵐ, Tᵦ, Sᵦ = compute_interface_heat_flux(flux_formulation, - ocean_surface_state, ice_state, - liquidus, ocean_properties, ℰ, u★) + ocean_surface_state, ice_state, + liquidus, ocean_properties, ℰ, u★) # Store interface values and heat flux @inbounds 𝒬ⁱⁿᵗ[i, j, 1] = 𝒬ⁱᵒ diff --git a/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_heat_flux_formulations.jl b/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_heat_flux_formulations.jl index ebf4f4ce8..915fa64b8 100644 --- a/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_heat_flux_formulations.jl +++ b/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_heat_flux_formulations.jl @@ -113,8 +113,8 @@ References - [holland1999modeling](@citet): Holland, D. M., & Jenkins, A. (1999). Modeling thermodynamic ice–ocean interactions at the base of an ice shelf. *Journal of Physical Oceanography*, 29(8), 1787-1800. -- [hieronymus2021comparison](@citet): Hieronymus, M., et al. (2021). A comparison of ocean-ice flux parametrizations. - *Geosci. Model Dev.*, 14, 4891-4908. +- [shi2021sensitivity](@citet): Shi, X., Notz, D., Liu, J., Yang, H., & Lohmann, G. (2021). Sensitivity of Northern + Hemisphere climate to ice-ocean interface heat flux parameterizations. *Geosci. Model Dev.*, 14, 4891-4908. """ struct ThreeEquationHeatFlux{F, T, FT, U} conductive_flux :: F @@ -139,7 +139,7 @@ Adapt.adapt_structure(to, f::ThreeEquationHeatFlux) = Construct a `ThreeEquationHeatFlux` with the specified parameters. -Default values follow [hieronymus2021comparison](@citet) with ``R = \\alpha_h / \\alpha_s = 35``. +Default values follow [shi2021sensitivity](@citet) with ``R = \\alpha_h / \\alpha_s = 35``. Keyword Arguments ================= diff --git a/src/EarthSystemModels/earth_system_model.jl b/src/EarthSystemModels/earth_system_model.jl index 113cce6b5..2823db249 100644 --- a/src/EarthSystemModels/earth_system_model.jl +++ b/src/EarthSystemModels/earth_system_model.jl @@ -317,7 +317,7 @@ function above_freezing_ocean_temperature!(ocean, grid, sea_ice) T = ocean_temperature(ocean) S = ocean_salinity(ocean) ℵ = sea_ice_concentration(sea_ice) - liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus + liquidus = sea_ice.model.phase_transitions.liquidus arch = architecture(grid) launch!(arch, grid, :xy, _above_freezing_ocean_temperature!, T, grid, S, ℵ, liquidus) diff --git a/src/SeaIces/SeaIces.jl b/src/SeaIces/SeaIces.jl index 57909bdc8..2f16a4f46 100644 --- a/src/SeaIces/SeaIces.jl +++ b/src/SeaIces/SeaIces.jl @@ -44,24 +44,31 @@ interpolate_state!(exchanger, grid, ::FreezingLimitedOceanTemperature, coupled_m # ComponentExchangers ComponentExchanger(sea_ice::FreezingLimitedOceanTemperature, grid) = nothing -function ComponentExchanger(sea_ice::Simulation{<:SeaIceModel}, grid) +function ComponentExchanger(sea_ice::Simulation{<:SeaIceModel}, grid) sea_ice_grid = sea_ice.model.grid - + if sea_ice_grid == grid u = sea_ice.model.velocities.u v = sea_ice.model.velocities.v - h = sea_ice.model.ice_thickness + hi = sea_ice.model.ice_thickness hc = sea_ice.model.ice_consolidation_thickness ℵ = sea_ice.model.ice_concentration + hs = sea_ice.model.snow_thickness else u = Field{Center, Center, Nothing}(grid) v = Field{Center, Center, Nothing}(grid) - h = Field{Center, Center, Nothing}(grid) + hi = Field{Center, Center, Nothing}(grid) hc = Field{Center, Center, Nothing}(grid) ℵ = Field{Center, Center, Nothing}(grid) + hs = Field{Center, Center, Nothing}(grid) + end + + # When there's no snow model, use ZeroField so kernels can read hs[i,j,1] = 0 + if isnothing(hs) + hs = ZeroField(eltype(grid)) end - return ComponentExchanger((; u, v, h, hc, ℵ), nothing) + return ComponentExchanger((; u, v, hi, hc, ℵ, hs), nothing) end end diff --git a/src/SeaIces/assemble_net_sea_ice_fluxes.jl b/src/SeaIces/assemble_net_sea_ice_fluxes.jl index 17d1e10c0..e7141bd41 100644 --- a/src/SeaIces/assemble_net_sea_ice_fluxes.jl +++ b/src/SeaIces/assemble_net_sea_ice_fluxes.jl @@ -26,13 +26,14 @@ function update_net_fluxes!(coupled_model, sea_ice::Simulation{<:SeaIceModel}) ℐꜜˡʷ = atmosphere_fields.ℐꜜˡʷ.data) freshwater_flux = atmosphere_fields.Jᶜ.data + snowfall_flux = atmosphere_fields.Jˢⁿ.data atmos_sea_ice_properties = coupled_model.interfaces.atmosphere_sea_ice_interface.properties sea_ice_properties = coupled_model.interfaces.sea_ice_properties sea_ice_surface_temperature = coupled_model.interfaces.atmosphere_sea_ice_interface.temperature ice_concentration = sea_ice_concentration(sea_ice) - + launch!(arch, grid, :xy, _assemble_net_sea_ice_fluxes!, top_fluxes, @@ -42,6 +43,7 @@ function update_net_fluxes!(coupled_model, sea_ice::Simulation{<:SeaIceModel}) atmosphere_sea_ice_fluxes, sea_ice_ocean_fluxes, freshwater_flux, + snowfall_flux, ice_concentration, sea_ice_surface_temperature, downwelling_radiation, @@ -57,7 +59,8 @@ end clock, atmosphere_sea_ice_fluxes, sea_ice_ocean_fluxes, - freshwater_flux, # Where do we add this one? + freshwater_flux, + snowfall_flux, ice_concentration, surface_temperature, downwelling_radiation, @@ -79,6 +82,7 @@ end 𝒬ᵛ = atmosphere_sea_ice_fluxes.latent_heat[i, j, 1] # latent heat flux 𝒬ᶠʳᶻ = sea_ice_ocean_fluxes.frazil_heat[i, j, 1] # frazil heat flux 𝒬ⁱⁿᵗ = sea_ice_ocean_fluxes.interface_heat[i, j, 1] # interfacial heat flux + Jˢⁿ = snowfall_flux[i, j, 1] end ρτˣ = atmosphere_sea_ice_fluxes.x_momentum # zonal momentum flux @@ -98,8 +102,9 @@ end # Mask fluxes over land for convenience inactive = inactive_node(i, j, kᴺ, grid, Center(), Center(), Center()) - @inbounds top_fluxes.heat[i, j, 1] = ifelse(inactive, zero(grid), ΣQt) - @inbounds top_fluxes.u[i, j, 1] = ifelse(inactive, zero(grid), ℑxᶠᵃᵃ(i, j, 1, grid, ρτˣ)) - @inbounds top_fluxes.v[i, j, 1] = ifelse(inactive, zero(grid), ℑyᵃᶠᵃ(i, j, 1, grid, ρτʸ)) - @inbounds bottom_heat_flux[i, j, 1] = ifelse(inactive, zero(grid), ΣQb) + @inbounds top_fluxes.heat[i, j, 1] = ifelse(inactive, zero(grid), ΣQt) + @inbounds top_fluxes.snowfall[i, j, 1] = ifelse(inactive, zero(grid), Jˢⁿ) + @inbounds top_fluxes.u[i, j, 1] = ifelse(inactive, zero(grid), ℑxᶠᵃᵃ(i, j, 1, grid, ρτˣ)) + @inbounds top_fluxes.v[i, j, 1] = ifelse(inactive, zero(grid), ℑyᵃᶠᵃ(i, j, 1, grid, ρτʸ)) + @inbounds bottom_heat_flux[i, j, 1] = ifelse(inactive, zero(grid), ΣQb) end diff --git a/src/SeaIces/sea_ice_simulation.jl b/src/SeaIces/sea_ice_simulation.jl index dc0681041..c718ab665 100644 --- a/src/SeaIces/sea_ice_simulation.jl +++ b/src/SeaIces/sea_ice_simulation.jl @@ -1,28 +1,47 @@ using ClimaSeaIce -using ClimaSeaIce: SeaIceModel, SlabThermodynamics, PhaseTransitions, ConductiveFlux -using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium +using ClimaSeaIce: SeaIceModel, PhaseTransitions, ConductiveFlux, + sea_ice_slab_thermodynamics, snow_slab_thermodynamics +using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium, IceSnowConductiveFlux using ClimaSeaIce.SeaIceDynamics: SplitExplicitSolver, SemiImplicitStress, SeaIceMomentumEquation, StressBalanceFreeDrift using ClimaSeaIce.Rheologies: IceStrength, ElastoViscoPlasticRheology +using Oceananigans.TimeSteppers: SplitRungeKuttaTimeStepper + using NumericalEarth.EarthSystemModels: ocean_surface_salinity, ocean_surface_velocities -using NumericalEarth.Oceans: Default +using NumericalEarth.Oceans: Default, reference_density default_rotation_rate = Oceananigans.defaults.planet_rotation_rate +ocean_reference_density(ocean::Simulation, FT) = convert(FT, reference_density(ocean)) +ocean_reference_density(::Nothing, FT) = convert(FT, 1026.0) + +function default_snow_thermodynamics(grid) + FT = eltype(grid) + snow_conductivity = FT(0.31) + # Use PrescribedTemperature so ClimaSeaIce does NOT run its own surface solve; + # the coupled flux solver in NumericalEarth handles the snow surface temperature. + snow_surface_temperature = Field{Center, Center, Nothing}(grid) + top_heat_boundary_condition = PrescribedTemperature(snow_surface_temperature.data) + return snow_slab_thermodynamics(grid; conductivity = snow_conductivity, top_heat_boundary_condition) +end + function sea_ice_simulation(grid, ocean=nothing; Δt = 5minutes, ice_salinity = 4, # psu - advection = nothing, # for the moment + advection = nothing, tracers = (), ice_heat_capacity = 2100, # J kg⁻¹ K⁻¹ ice_consolidation_thickness = 0.05, # m - ice_density = 900, # kg m⁻³ + sea_ice_density = 900, # kg m⁻³ + snow_density = 330, # kg m⁻³ dynamics = sea_ice_dynamics(grid, ocean), bottom_heat_boundary_condition = nothing, top_heat_boundary_condition = nothing, - phase_transitions = PhaseTransitions(; heat_capacity = ice_heat_capacity, density = ice_density), - conductivity = 2, # kg m s⁻³ K⁻¹ - internal_heat_flux = ConductiveFlux(; conductivity)) + timestepper = :SplitRungeKutta3, + phase_transitions = PhaseTransitions(eltype(grid); heat_capacity=ice_heat_capacity, density=sea_ice_density), + conductivity = 2, # W m⁻¹ K⁻¹ + internal_heat_flux = ConductiveFlux(; conductivity), + snow_thermodynamics = default_snow_thermodynamics(grid)) # Build consistent boundary conditions for the ice model: # - bottom -> flux boundary condition @@ -37,20 +56,19 @@ function sea_ice_simulation(grid, ocean=nothing; if isnothing(ocean) surface_ocean_salinity = 0 else - kᴺ = size(grid, 3) surface_ocean_salinity = ocean_surface_salinity(ocean) end bottom_heat_boundary_condition = IceWaterThermalEquilibrium(surface_ocean_salinity) end - ice_thermodynamics = SlabThermodynamics(grid; - internal_heat_flux, - phase_transitions, - top_heat_boundary_condition, - bottom_heat_boundary_condition) + ice_thermodynamics = sea_ice_slab_thermodynamics(grid; + internal_heat_flux, + top_heat_boundary_condition, + bottom_heat_boundary_condition) bottom_heat_flux = Field{Center, Center, Nothing}(grid) top_heat_flux = Field{Center, Center, Nothing}(grid) + snowfall = Field{Center, Center, Nothing}(grid) # Build the sea ice model sea_ice_model = SeaIceModel(grid; @@ -58,30 +76,51 @@ function sea_ice_simulation(grid, ocean=nothing; advection, tracers, ice_consolidation_thickness, + sea_ice_density, + snow_density, + phase_transitions, ice_thermodynamics, + snow_thermodynamics, + snowfall, dynamics, + timestepper, bottom_heat_flux, top_heat_flux) verbose = false - - # Build the simulation sea_ice = Simulation(sea_ice_model; Δt, verbose) return sea_ice end +default_coriolis(ocean::Simulation) = ocean.model.coriolis +default_coriolis(ocean::Nothing) = HydrostaticSphericalCoriolis(; rotation_rate=default_rotation_rate) + +default_solver(grid, ocean) = SplitExplicitSolver(grid; substeps=120) + +# We assume RK3 has a larger timestep +function default_solver(grid, ocean::Simulation) + substeps = if ocean.model.timestepper isa SplitRungeKuttaTimeStepper + 240 + else + 120 + end + return SplitExplicitSolver(grid; substeps) +end + function sea_ice_dynamics(grid, ocean=nothing; - sea_ice_ocean_drag_coefficient = 5.5e-3, + sea_ice_ocean_drag_coefficient = 3.24e-3, rheology = ElastoViscoPlasticRheology(), - coriolis = HydrostaticSphericalCoriolis(; rotation_rate=default_rotation_rate), + coriolis = default_coriolis(ocean), free_drift = nothing, - solver = SplitExplicitSolver(grid; substeps=120)) + solver = default_solver(grid, ocean)) SSU, SSV = ocean_surface_velocities(ocean) - sea_ice_ocean_drag_coefficient = convert(eltype(grid), sea_ice_ocean_drag_coefficient) + FT = eltype(grid) + sea_ice_ocean_drag_coefficient = convert(FT, sea_ice_ocean_drag_coefficient) + ρₑ = ocean_reference_density(ocean, FT) - τo = SemiImplicitStress(uₑ=SSU, vₑ=SSV, Cᴰ=sea_ice_ocean_drag_coefficient) + τo = SemiImplicitStress(uₑ=SSU, vₑ=SSV, Cᴰ=sea_ice_ocean_drag_coefficient, ρₑ=ρₑ) τua = Field{Face, Center, Nothing}(grid) τva = Field{Center, Face, Nothing}(grid) @@ -105,8 +144,10 @@ end sea_ice_thickness(sea_ice::Simulation{<:SeaIceModel}) = sea_ice.model.ice_thickness sea_ice_concentration(sea_ice::Simulation{<:SeaIceModel}) = sea_ice.model.ice_concentration -heat_capacity(sea_ice::Simulation{<:SeaIceModel}) = sea_ice.model.ice_thermodynamics.phase_transitions.heat_capacity -reference_density(sea_ice::Simulation{<:SeaIceModel}) = sea_ice.model.ice_thermodynamics.phase_transitions.density +heat_capacity(sea_ice::Simulation{<:SeaIceModel}) = sea_ice.model.phase_transitions.heat_capacity +# `sea_ice.model.sea_ice_density` is wrapped as a `ConstantField` by `SeaIceModel`; +# the scalar value lives on `phase_transitions.density`. +reference_density(sea_ice::Simulation{<:SeaIceModel}) = sea_ice.model.phase_transitions.density function net_fluxes(sea_ice::Simulation{<:SeaIceModel}) net_momentum_fluxes = if isnothing(sea_ice.model.dynamics) @@ -119,15 +160,21 @@ function net_fluxes(sea_ice::Simulation{<:SeaIceModel}) (; u, v) end - net_top_sea_ice_fluxes = merge((; heat=sea_ice.model.external_heat_fluxes.top), net_momentum_fluxes) + net_top_sea_ice_fluxes = merge((; heat=sea_ice.model.external_heat_fluxes.top, snowfall=sea_ice.model.snowfall), net_momentum_fluxes) net_bottom_sea_ice_fluxes = (; heat=sea_ice.model.external_heat_fluxes.bottom) return (; bottom = net_bottom_sea_ice_fluxes, top = net_top_sea_ice_fluxes) end function default_ai_temperature(sea_ice::Simulation{<:SeaIceModel}) - conductive_flux = sea_ice.model.ice_thermodynamics.internal_heat_flux.parameters.flux - return SkinTemperature(conductive_flux) + ice_flux = sea_ice.model.ice_thermodynamics.internal_heat_flux + snow_thermo = sea_ice.model.snow_thermodynamics + internal_flux = if isnothing(snow_thermo) + ice_flux + else + IceSnowConductiveFlux(snow_thermo.internal_heat_flux.conductivity, ice_flux.conductivity) + end + return SkinTemperature(internal_flux) end # Constructor that accepts the sea-ice model @@ -136,7 +183,7 @@ function ThreeEquationHeatFlux(sea_ice::Simulation{<:SeaIceModel}, FT::DataType salt_transfer_coefficient = heat_transfer_coefficient / 35, friction_velocity = convert(FT, 0.002)) - conductive_flux = sea_ice.model.ice_thermodynamics.internal_heat_flux.parameters.flux + conductive_flux = sea_ice.model.ice_thermodynamics.internal_heat_flux ice_temperature = sea_ice.model.ice_thermodynamics.top_surface_temperature return ThreeEquationHeatFlux(conductive_flux, diff --git a/test/Project.toml b/test/Project.toml index 54e377ccc..8824c296e 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -28,8 +28,8 @@ WorldOceanAtlasTools = "04f20302-f1b9-11e8-29d9-7d841cb0a64a" XESMF = "2e0b0046-e7a1-486f-88de-807ee8ffabe5" [sources] +ClimaSeaIce = {url = "https://github.com/CliMA/ClimaSeaIce.jl", rev = "ss/refactor-thermodynamics"} NumericalEarth = {path = ".."} -Breeze = {url = "https://github.com/NumericalEarth/Breeze.jl/", rev = "main"} [compat] Breeze = "0.4" diff --git a/test/runtests.jl b/test/runtests.jl index 08a4ce9b9..95e300072 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -72,6 +72,7 @@ function __init__() try atmosphere = JRA55PrescribedAtmosphere(backend=JRA55NetCDFBackend(2)) + land = JRA55PrescribedLand(backend=JRA55NetCDFBackend(2)) catch e @warn "Original JRA55 download failed, trying NumericalEarthArtifacts fallback..." exception=(e, catch_backtrace()) emit_ci_warning("Broken JRA55 download", "Original source failed during init") @@ -79,7 +80,6 @@ function __init__() datum = Metadatum(name; dataset=JRA55.RepeatYearJRA55()) download_from_artifacts(metadata_path(datum)) end - atmosphere = JRA55PrescribedAtmosphere(backend=JRA55NetCDFBackend(2)) end ##### diff --git a/test/test_ocean_sea_ice_model.jl b/test/test_ocean_sea_ice_model.jl index 4188435ec..acfc1246f 100644 --- a/test/test_ocean_sea_ice_model.jl +++ b/test/test_ocean_sea_ice_model.jl @@ -44,7 +44,7 @@ using ClimaSeaIce.Rheologies ocean = ocean_simulation(grid; free_surface) sea_ice = sea_ice_simulation(grid, ocean; advection=WENO(order=7)) - liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus + liquidus = sea_ice.model.phase_transitions.liquidus # Set the ocean temperature and salinity set!(ocean.model, T=temperature_metadata[1], S=salinity_metadata[1]) diff --git a/test/test_snow_model_integration.jl b/test/test_snow_model_integration.jl new file mode 100644 index 000000000..a5860e6d9 --- /dev/null +++ b/test/test_snow_model_integration.jl @@ -0,0 +1,218 @@ +include("runtests_setup.jl") + +using ClimaSeaIce: SeaIceModel, ConductiveFlux +using ClimaSeaIce.SeaIceThermodynamics: IceSnowConductiveFlux +using NumericalEarth.SeaIces: default_snow_thermodynamics +using NumericalEarth.EarthSystemModels.InterfaceComputations: + ComponentInterfaces, + SkinTemperature, + InterfaceProperties, + conductive_flux_balance_temperature + +using Oceananigans.Fields: ZeroField +using Oceananigans.Units: hours, days + +##### +##### Unit tests +##### + +@testset "Snow model unit tests" begin + for arch in test_architectures + A = typeof(arch) + + grid = RectilinearGrid(arch; + size = (4, 4, 1), + extent = (1, 1, 1), + topology = (Periodic, Periodic, Bounded)) + + @testset "sea_ice_simulation with_snow=false [$A]" begin + sea_ice = sea_ice_simulation(grid; dynamics=nothing, snow_thermodynamics=nothing) + @test sea_ice isa Simulation + @test sea_ice.model isa SeaIceModel + @test sea_ice.model.snow_thermodynamics === nothing + @test sea_ice.model.ice_thermodynamics.internal_heat_flux isa ConductiveFlux + end + + @testset "sea_ice_simulation with_snow=true [$A]" begin + sea_ice = sea_ice_simulation(grid; dynamics=nothing) + @test sea_ice isa Simulation + @test sea_ice.model.snow_thermodynamics !== nothing + @test sea_ice.model.snow_thickness isa Field + end + + @testset "PhaseTransitions API [$A]" begin + sea_ice = sea_ice_simulation(grid; dynamics=nothing) + pt = sea_ice.model.phase_transitions + @test pt.heat_capacity == 2100 + @test pt.density == 900 + end + + @testset "ComponentExchanger includes hs [$A]" begin + using NumericalEarth.EarthSystemModels.InterfaceComputations: ComponentExchanger + + # Without snow: hs should be ZeroField + sea_ice = sea_ice_simulation(grid; dynamics=nothing, snow_thermodynamics=nothing) + exchanger = ComponentExchanger(sea_ice, grid) + @test haskey(exchanger.state, :hs) + @test haskey(exchanger.state, :hi) + @test exchanger.state.hs isa ZeroField + + # With snow: hs should be a Field + sea_ice_snow = sea_ice_simulation(grid; dynamics=nothing) + exchanger_snow = ComponentExchanger(sea_ice_snow, grid) + @test exchanger_snow.state.hs isa Field + end + + @testset "default_ai_temperature dispatches on snow [$A]" begin + using NumericalEarth.SeaIces: default_ai_temperature + + sea_ice = sea_ice_simulation(grid; dynamics=nothing, snow_thermodynamics=nothing) + st = default_ai_temperature(sea_ice) + @test st isa SkinTemperature + @test st.internal_flux isa ConductiveFlux + + sea_ice_snow = sea_ice_simulation(grid; dynamics=nothing) + st_snow = default_ai_temperature(sea_ice_snow) + @test st_snow isa SkinTemperature + @test st_snow.internal_flux isa IceSnowConductiveFlux + end + + @testset "net_fluxes includes snowfall [$A]" begin + using NumericalEarth.SeaIces: net_fluxes + + sea_ice = sea_ice_simulation(grid; dynamics=nothing) + fluxes = net_fluxes(sea_ice) + @test haskey(fluxes.top, :snowfall) + end + + @testset "Snow insulates: warmer surface temperature [$A]" begin + # Build two coupled models — one without snow, one with snow and + # nonzero snow thickness — then compare surface temperatures after + # one coupled time step. Snow adds thermal resistance, so the + # surface should be warmer (closer to the warmer atmosphere). + # + # Radiation is disabled (ε=0) so the surface energy balance + # reduces to conductive + turbulent fluxes; otherwise the + # Stefan–Boltzmann loss swamps the small snow-insulation signal. + ocean_grid = RectilinearGrid(arch; + size = (1, 1, 2), + extent = (1, 1, 1), + topology = (Periodic, Periodic, Bounded)) + + function build_coupled(; with_snow) + ocean = ocean_simulation(ocean_grid; + momentum_advection = nothing, + tracer_advection = nothing, + closure = nothing, + coriolis = nothing) + set!(ocean.model, T = -1.8, S = 34) + + snow_thermodynamics = with_snow ? default_snow_thermodynamics(ocean_grid) : nothing + sea_ice = sea_ice_simulation(ocean_grid, ocean; + dynamics = nothing, + snow_thermodynamics) + set!(sea_ice.model, h = 1.0, ℵ = 1.0) + if with_snow + set!(sea_ice.model, hs = 0.2) + end + + atmosphere = PrescribedAtmosphere(ocean_grid, [0.0]) + parent(atmosphere.velocities.u) .= 2.0 + radiation = Radiation(ocean_emissivity = 0, sea_ice_emissivity = 0) + return OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) + end + + bare = build_coupled(with_snow = false) + snowy = build_coupled(with_snow = true) + + time_step!(bare, 1) + time_step!(snowy, 1) + + Ts_bare = bare.interfaces.atmosphere_sea_ice_interface.temperature + Ts_snowy = snowy.interfaces.atmosphere_sea_ice_interface.temperature + + @allowscalar begin + # Snow insulation → warmer (or equal) surface temperature + @test Ts_snowy[1, 1, 1] ≥ Ts_bare[1, 1, 1] + end + end + end +end + +##### +##### Integration tests +##### + +@testset "Snow model integration tests" begin + for arch in test_architectures + A = typeof(arch) + + grid = RectilinearGrid(arch; + size = (4, 4, 2), + extent = (1, 1, 1), + topology = (Periodic, Periodic, Bounded)) + + @testset "Coupled model with snow [$A]" begin + ocean = ocean_simulation(grid; + momentum_advection = nothing, + tracer_advection = nothing, + closure = nothing, + coriolis = nothing) + + sea_ice = sea_ice_simulation(grid, ocean; + dynamics = nothing) + + atmosphere = PrescribedAtmosphere(grid, [0.0]) + radiation = Radiation() + + @test begin + coupled = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) + time_step!(coupled, 1) + true + end + end + + @testset "Snowfall routing [$A]" begin + ocean = ocean_simulation(grid; + momentum_advection = nothing, + tracer_advection = nothing, + closure = nothing, + coriolis = nothing) + + sea_ice = sea_ice_simulation(grid, ocean; + dynamics = nothing, + snow_thermodynamics = nothing) + + atmosphere = PrescribedAtmosphere(grid, [0.0]) + radiation = Radiation() + + coupled = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) + + # The snowfall field should exist in the exchanger + exchanger = coupled.interfaces.exchanger + @test haskey(exchanger.atmosphere.state, :Jˢⁿ) + + # The net fluxes should include snowfall + top_fluxes = coupled.interfaces.net_fluxes.sea_ice.top + @test haskey(top_fluxes, :snowfall) + end + + @testset "Snow insulation effect [$A]" begin + # With snow, the IceSnowConductiveFlux is used for the surface + # temperature solve. Verify this dispatch is wired correctly. + ocean = ocean_simulation(grid; + momentum_advection = nothing, + tracer_advection = nothing, + closure = nothing, + coriolis = nothing) + + sea_ice = sea_ice_simulation(grid, ocean; + dynamics = nothing) + + ai_temp = NumericalEarth.SeaIces.default_ai_temperature(sea_ice) + @test ai_temp.internal_flux isa IceSnowConductiveFlux + @test ai_temp.internal_flux.ice_conductivity == 2.0 + @test ai_temp.internal_flux.snow_conductivity == 0.31 + end + end +end