Integrate ClimaSeaIce snow model

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"

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"

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},

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

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,16 +148,19 @@ 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
 
 #####
 ##### 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...)

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

src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl

@@ -84,18 +84,17 @@ 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]
         Tₛ = convert_to_kelvin(sea_ice_properties.temperature_units, Tₛ)
     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)

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

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,

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] = 𝒬ⁱᵒ

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
 =================

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)