src/Physics/HeatConduction/HeatConduction.jl

 module HeatConduction
 
-import CryoGrid: SubSurfaceProcess, BoundaryStyle, Dirichlet, Neumann, BoundaryProcess, Layer, Top, Bottom, SubSurface, Callback
-import CryoGrid: diagnosticstep!, prognosticstep!, interact!, initialcondition!, boundaryflux, boundaryvalue, variables, callbacks, criterion, affect!, thickness, midpoints
-
+using CryoGrid
 using CryoGrid.InputOutput: Forcing
 using CryoGrid.Physics
 using CryoGrid.Physics.Boundaries
@@ -15,6 +13,9 @@ using IfElse
 using ModelParameters
 using Unitful
 
+import CryoGrid
+import CryoGrid.Physics
+
 export Heat, TemperatureProfile
 export FreeWater, FreezeCurve, freezecurve
 
@@ -29,55 +30,57 @@ Convenience constructor for `Numerics.Profile` which automatically converts temp
 TemperatureProfile(pairs::Pair{<:Union{DistQuantity,Param},<:Union{TempQuantity,Param}}...) = Profile(map(p -> uconvert(u"m", p[1]) => uconvert(u"°C", p[2]), pairs))
 
 """
-    HeatImpl
+    HeatParameterization
 
 Base type for different numerical formulations of two-phase heat diffusion.
 """
-abstract type HeatImpl end
-struct Enthalpy <: HeatImpl end
-struct Temperature <: HeatImpl end
+abstract type HeatParameterization end
+struct Enthalpy <: HeatParameterization end
+struct Temperature <: HeatParameterization end
 
-"""
-    ThermalProperties
-
-Base type for defining thermal properties.
-"""
-abstract type ThermalProperties <: IterableStruct end
-"""
-    HydroThermalProperties{Tρ,TLf,Tkw,Tki,Tcw,Tci}
-
-Thermal properties of water used in two-phase heat conduction.
-"""
-@kwdef struct HydroThermalProperties{Tρw,TLf,Tkw,Tki,Tka,Tcw,Tci,Tca} <: ThermalProperties
-    ρw::Tρw = Param(1000.0, units=u"kg/m^3") # density of water
-    Lf::TLf = Param(334000.0, units=u"J/kg") # latent heat of fusion of water
-    kw::Tkw = Param(0.57, units=u"W/m/K") # thermal conductivity of water [Hillel(1982)]
-    ki::Tki = Param(2.2, units=u"W/m/K") # thermal conductivity of ice [Hillel(1982)]
-    ka::Tka = Param(0.025, units=u"W/m/K") # air [Hillel(1982)]
-    cw::Tcw = Param(4.2e6, units=u"J/K/m^3") # heat capacity of water
-    ci::Tci = Param(1.9e6, units=u"J/K/m^3") # heat capacity of ice
-    ca::Tca = Param(0.00125e6, units=u"J/K/m^3") # heat capacity of air
+abstract type StepLimiter end
+@kwdef struct CFL <: StepLimiter
+    fallback_dt::Float64 = 60.0 # fallback dt [s]
 end
 
-@kwdef struct Heat{Tfc<:FreezeCurve,Tsp,Tinit,TProp<:HydroThermalProperties,TL} <: SubSurfaceProcess
-    prop::TProp = HydroThermalProperties()
-    L::TL = prop.ρw*prop.Lf # [J/m^3] (specific latent heat of fusion of water)
-    freezecurve::Tfc = FreeWater() # freeze curve, defautls to free water fc
-    sp::Tsp = Enthalpy() # specialization
-    init::Tinit = nothing # optional initialization scheme
+ThermalProperties(
+    consts=Physics.Constants();
+    ρw = consts.ρw,
+    Lf = consts.Lf,
+    L = consts.ρw*consts.Lf,
+    kw = Param(0.57, units=u"W/m/K"), # thermal conductivity of water [Hillel(1982)]
+    ki = Param(2.2, units=u"W/m/K"), # thermal conductivity of ice [Hillel(1982)]
+    ka = Param(0.025, units=u"W/m/K"), # air [Hillel(1982)]
+    cw = Param(4.2e6, units=u"J/K/m^3"), # heat capacity of water
+    ci = Param(1.9e6, units=u"J/K/m^3"), # heat capacity of ice
+    ca = Param(0.00125e6, units=u"J/K/m^3"), # heat capacity of air
+) = (; ρw, Lf, L, kw, ki, ka, cw, ci, ca)
+
+struct Heat{Tfc<:FreezeCurve,TPara<:HeatParameterization,Tdt,Tinit,TProp} <: SubSurfaceProcess
+    para::TPara
+    prop::TProp
+    freezecurve::Tfc
+    dtlim::Tdt  # timestep limiter
+    init::Tinit # optional initialization scheme
 end
 # convenience constructors for specifying prognostic variable as symbol
-Heat(var::Symbol; kwargs...) = Heat(Val{var}(); kwargs...)
-Heat(::Val{:H}; kwargs...) = Heat(;sp=Enthalpy(), kwargs...)
-Heat(::Val{:T}; kwargs...) = Heat(;sp=Temperature(), kwargs...)
+Heat(var::Symbol=:H; kwargs...) = Heat(Val{var}(); kwargs...)
+Heat(::Val{:H}; kwargs...) = Heat(Enthalpy(); kwargs...)
+Heat(::Val{:T}; kwargs...) = Heat(Temperature(); kwargs...)
+Heat(para::Enthalpy; freezecurve=FreeWater(), prop=ThermalProperties(), dtlim=nothing, init=nothing) = Heat(para, prop, deepcopy(freezecurve), dtlim, init)
+Heat(para::Temperature; freezecurve, prop=ThermalProperties(), dtlim=CFL(), init=nothing) = Heat(para, prop, deepcopy(freezecurve), dtlim, init)
+
+# getter functions
+thermalproperties(heat::Heat) = heat.prop
 freezecurve(heat::Heat) = heat.freezecurve
+
 # Default implementation of `variables` for freeze curve
-variables(::SubSurface, ::Heat, ::FreezeCurve) = ()
+CryoGrid.variables(::SubSurface, ::Heat, ::FreezeCurve) = ()
 
-export ConstantTemp, GeothermalHeatFlux, TemperatureGradient, NFactor, Damping
+export HeatBC, ConstantTemp, GeothermalHeatFlux, TemperatureGradient, NFactor
 include("heat_bc.jl")
 
-export heatconduction!, enthalpy, totalwater, liquidwater, freezethaw!, heatcapacity, heatcapacity!, thermalconductivity, thermalconductivity!
+export heatconduction!, enthalpy, enthalpyinv, waterice, liquidwater, freezethaw!, heatcapacity, heatcapacity!, thermalconductivity, thermalconductivity!
 include("heat.jl")
 
 end

src/Physics/HeatConduction/heat.jl

@@ -18,12 +18,6 @@ Computes the apparent or "effective" heat capacity `dHdT` as a function of tempe
 latent heat of fusion, derivative of the freeze curve `dθdT`, and the constituent heat capacities of water and ice.
 """
 @inline C_eff(T, C, L, dθdT, cw, ci) = C + dθdT*(L + T*(cw - ci))
-"""
-    liquidwater(sub::SubSurface, heat::Heat, state)
-
-Retrieves the liquid water content for the given layer.
-"""
-@inline liquidwater(::SubSurface, ::Heat, state) = state.θl
 """
     thermalconductivities(::SubSurface, heat::Heat)
 
@@ -36,23 +30,23 @@ Get thermal conductivities for generic `SubSurface` layer.
 Get heat capacities for generic `SubSurface` layer.
 """
 @inline heatcapacities(::SubSurface, heat::Heat) = (heat.prop.cw, heat.prop.ci, heat.prop.ca)
+
+# Generic heat conduction implementation
+
 """
-    volumetricfractions(::SubSurface, heat::Heat)
+    volumetricfractions(::SubSurface, heat::Heat, state, i)
 
 Get constituent volumetric fractions for generic `SubSurface` layer. Default implementation assumes
-the only constitutents are liquid water, ice, and air: `(θl,θi,θa)`.
-"""
-@inline function volumetricfractions(sub::SubSurface, heat::Heat, state, i)
-    return let θw = totalwater(sub, state, i),
-        θl = state.θl[i],
-        θa = 1.0 - θw,
-        θi = θw - θl;
-        (θl, θi, θa)
+the only constitutents are liquid water, ice, and air: `(θw,θi,θa)`.
+"""
+@inline function Physics.volumetricfractions(sub::SubSurface, heat::Heat, state, i)
+    let θwi = waterice(sub, heat, state, i),
+        θw = liquidwater(sub, heat, state, i),
+        θa = 1.0 - θwi,
+        θi = θwi - θw;
+        return (θw, θi, θa)
     end
 end
-
-# Generic heat conduction implementation
-
 """
     freezethaw!(sub::SubSurface, heat::Heat, state)
 
@@ -63,15 +57,13 @@ should also be computed depending on the thermal scheme being implemented.
 """
 freezethaw!(sub::SubSurface, heat::Heat, state) = error("missing implementation of freezethaw!")
 """
-    heatcapacity(capacities::NTuple{N}, fracs::NTuple{N}) where N
+    heatcapacity(sub::SubSurface, heat::Heat, θfracs...)
 
-Computes the heat capacity as a weighted average over constituent `capacities` with volumetric fractions `fracs`.
+Computes the heat capacity as a weighted average over constituent capacities with volumetric fractions `θfracs`.
 """
-heatcapacity(capacities::NTuple{N,Any}, fracs::NTuple{N,Any}) where N = sum(map(*, capacities, fracs))
-@inline function heatcapacity(sub::SubSurface, heat::Heat, state, i)
-    θs = volumetricfractions(sub, heat, state, i)
+@inline function heatcapacity(sub::SubSurface, heat::Heat, θfracs...)
     cs = heatcapacities(sub, heat)
-    return heatcapacity(cs, θs)
+    return sum(map(*, cs, θfracs))
 end
 """
     heatcapacity!(sub::SubSurface, heat::Heat, state)
@@ -80,19 +72,18 @@ Computes the heat capacity for the given layer from the current state and stores
 """
 @inline function heatcapacity!(sub::SubSurface, heat::Heat, state)
     @inbounds for i in 1:length(state.T)
-        state.C[i] = heatcapacity(sub, heat, state, i)
+        θfracs = volumetricfractions(sub, heat, state, i)
+        state.C[i] = heatcapacity(sub, heat, θfracs...)
     end
 end
 """
-    thermalconductivity(conductivities::NTuple{N}, fracs::NTuple{N}) where N
+    thermalconductivity(sub::SubSurface, heat::Heat, θfracs...)
 
-Computes the thermal conductivity as a squared weighted sum over constituent `conductivities` with volumetric fractions `fracs`.
+Computes the thermal conductivity as a squared weighted sum over constituent conductivities with volumetric fractions `θfracs`.
 """
-thermalconductivity(conductivities::NTuple{N,Any}, fracs::NTuple{N,Any}) where N = sum(map(*, map(sqrt, conductivities), fracs))^2
-@inline function thermalconductivity(sub::SubSurface, heat::Heat, state, i)
-    θs = volumetricfractions(sub, heat, state, i)
+@inline function thermalconductivity(sub::SubSurface, heat::Heat, θfracs...)
     ks = thermalconductivities(sub, heat)
-    return thermalconductivity(ks, θs)
+    return sum(map(*, map(sqrt, ks), θfracs))^2
 end
 """
     thermalconductivity!(sub::SubSurface, heat::Heat, state)
@@ -101,14 +92,15 @@ Computes the thermal conductivity for the given layer from the current state and
 """
 @inline function thermalconductivity!(sub::SubSurface, heat::Heat, state)
     @inbounds for i in 1:length(state.T)
-        state.kc[i] = thermalconductivity(sub, heat, state, i)
+        θfracs = volumetricfractions(sub, heat, state, i)
+        state.kc[i] = thermalconductivity(sub, heat, θfracs...)
     end
 end
 """
 Variable definitions for heat conduction on any subsurface layer. Joins variables defined on
 `layer` and `heat` individually as well as variables defined by the freeze curve.
 """
-variables(sub::SubSurface, heat::Heat) = (
+CryoGrid.variables(sub::SubSurface, heat::Heat) = (
     variables(sub)..., # layer variables
     variables(heat)...,  # heat variables
     variables(sub, heat, freezecurve(heat))..., # freeze curve variables
@@ -116,32 +108,32 @@ variables(sub::SubSurface, heat::Heat) = (
 """
 Variable definitions for heat conduction (enthalpy).
 """
-variables(heat::Heat{<:FreezeCurve,Enthalpy}) = (
+CryoGrid.variables(heat::Heat{<:FreezeCurve,Enthalpy}) = (
     Prognostic(:H, OnGrid(Cells), u"J/m^3"),
     Diagnostic(:T, OnGrid(Cells), u"°C"),
-    basevariables(heat)...,
+    heatvariables(heat)...,
 )
 """
 Variable definitions for heat conduction (temperature).
 """
-variables(heat::Heat{<:FreezeCurve,Temperature}) = (
+CryoGrid.variables(heat::Heat{<:FreezeCurve,Temperature}) = (
     Prognostic(:T, OnGrid(Cells), u"°C"),
     Diagnostic(:H, OnGrid(Cells), u"J/m^3"),
     Diagnostic(:dH, OnGrid(Cells), u"W/m^3"),
-    Diagnostic(:dθdT, OnGrid(Cells)),
-    basevariables(heat)...,
+    Diagnostic(:dθdT, OnGrid(Cells), domain=0..Inf),
+    heatvariables(heat)...,
 )
 """
 Common variable definitions for all heat implementations.
 """
-basevariables(::Heat) = (
+heatvariables(::Heat) = (
     Diagnostic(:dH_upper, Scalar, u"J/K/m^2"),
     Diagnostic(:dH_lower, Scalar, u"J/K/m^2"),
-    Diagnostic(:dHdT, OnGrid(Cells), u"J/K/m^3"),
+    Diagnostic(:dHdT, OnGrid(Cells), u"J/K/m^3", domain=0..Inf),
     Diagnostic(:C, OnGrid(Cells), u"J/K/m^3"),
     Diagnostic(:k, OnGrid(Edges), u"W/m/K"),
     Diagnostic(:kc, OnGrid(Cells), u"W/m/K"),
-    Diagnostic(:θl, OnGrid(Cells)),
+    Diagnostic(:θw, OnGrid(Cells), domain=0..1),
 )
 """
     heatconduction!(∂H,T,ΔT,k,Δk)
@@ -179,7 +171,7 @@ end
 """
 Diagonstic step for heat conduction (all state configurations) on any subsurface layer.
 """
-function diagnosticstep!(sub::SubSurface, heat::Heat, state)
+function CryoGrid.diagnosticstep!(sub::SubSurface, heat::Heat, state)
     # Reset energy flux to zero; this is redundant when H is the prognostic variable
     # but necessary when it is not.
     @. state.dH = zero(eltype(state.dH))
@@ -201,9 +193,57 @@ function diagnosticstep!(sub::SubSurface, heat::Heat, state)
     return nothing # ensure no allocation
 end
 """
+Generic top interaction. Computes flux dH at top cell.
+"""
+function CryoGrid.interact!(top::Top, bc::HeatBC, sub::SubSurface, heat::Heat, stop, ssub)
+    Δk = CryoGrid.thickness(sub, ssub, first) # using `thickness` allows for generic layer implementations
+    # boundary flux
+    @setscalar ssub.dH_upper = boundaryflux(bc, top, heat, sub, stop, ssub)
+    @inbounds ssub.dH[1] += getscalar(ssub.dH_upper) / Δk
+    return nothing # ensure no allocation
+end
+"""
+Generic bottom interaction. Computes flux dH at bottom cell.
+"""
+function CryoGrid.interact!(sub::SubSurface, heat::Heat, bot::Bottom, bc::HeatBC, ssub, sbot)
+    Δk = CryoGrid.thickness(sub, ssub, last) # using `thickness` allows for generic layer implementations
+    # boundary flux
+    @setscalar ssub.dH_lower = boundaryflux(bc, bot, heat, sub, sbot, ssub)
+    @inbounds ssub.dH[end] += getscalar(ssub.dH_lower) / Δk
+    return nothing # ensure no allocation
+end
+"""
+Generic subsurface interaction. Computes flux dH at boundary between subsurface layers.
+"""
+function CryoGrid.interact!(sub1::SubSurface, ::Heat, sub2::SubSurface, ::Heat, s1, s2)
+    Δk₁ = CryoGrid.thickness(sub1, s1, last)
+    Δk₂ = CryoGrid.thickness(sub2, s2, first)
+    # thermal conductivity between cells
+    k = s1.k[end] = s2.k[1] =
+        @inbounds let k₁ = s1.kc[end],
+            k₂ = s2.kc[1],
+            Δ₁ = Δk₁[end],
+            Δ₂ = Δk₂[1];
+            harmonicmean(k₁, k₂, Δ₁, Δ₂)
+        end
+    # calculate heat flux between cells
+    Qᵢ = @inbounds let z₁ = CryoGrid.midpoint(sub1, s1, last),
+        z₂ = CryoGrid.midpoint(sub2, s2, first),
+        δ = z₂ - z₁;
+        k*(s2.T[1] - s1.T[end]) / δ
+    end
+    # diagnostics
+    @setscalar s1.dH_lower = Qᵢ
+    @setscalar s2.dH_upper = -Qᵢ
+    # add fluxes scaled by grid cell size
+    @inbounds s1.dH[end] += Qᵢ / Δk₁[end]
+    @inbounds s2.dH[1] += -Qᵢ / Δk₂[1]
+    return nothing # ensure no allocation
+end
+"""
 Prognostic step for heat conduction (enthalpy) on subsurface layer.
 """
-function prognosticstep!(::SubSurface, ::Heat{<:FreezeCurve,Enthalpy}, state)
+function CryoGrid.prognosticstep!(::SubSurface, ::Heat{<:FreezeCurve,Enthalpy}, state)
     Δk = Δ(state.grids.k) # cell sizes
     ΔT = Δ(state.grids.T)
     # Diffusion on non-boundary cells
@@ -212,7 +252,7 @@ end
 """
 Prognostic step for heat conduction (temperature) on subsurface layer.
 """
-function prognosticstep!(sub::SubSurface, ::Heat{<:FreezeCurve,Temperature}, state)
+function CryoGrid.prognosticstep!(sub::SubSurface, ::Heat{<:FreezeCurve,Temperature}, state)
     Δk = Δ(state.grids.k) # cell sizes
     ΔT = Δ(state.grids.T)
     # Diffusion on non-boundary cells
@@ -223,97 +263,64 @@ function prognosticstep!(sub::SubSurface, ::Heat{<:FreezeCurve,Temperature}, sta
     return nothing
 end
 # Boundary fluxes
-@inline boundaryflux(::Neumann, bc::HeatBC, top::Top, heat::Heat, sub::SubSurface, stop, ssub) = boundaryvalue(bc,top,heat,sub,stop,ssub)
-@inline boundaryflux(::Neumann, bc::HeatBC, bot::Bottom, heat::Heat, sub::SubSurface, sbot, ssub) = boundaryvalue(bc,bot,heat,sub,sbot,ssub)
-@inline function boundaryflux(::Dirichlet, bc::HeatBC, top::Top, heat::Heat, sub::SubSurface, stop, ssub)
-    Δk = thickness(sub, ssub) # using `thickness` allows for generic layer implementations
+@inline CryoGrid.boundaryflux(::Neumann, bc::HeatBC, top::Top, heat::Heat, sub::SubSurface, stop, ssub) = boundaryvalue(bc,top,heat,sub,stop,ssub)
+@inline CryoGrid.boundaryflux(::Neumann, bc::HeatBC, bot::Bottom, heat::Heat, sub::SubSurface, sbot, ssub) = boundaryvalue(bc,bot,heat,sub,sbot,ssub)
+@inline function CryoGrid.boundaryflux(::Dirichlet, bc::HeatBC, top::Top, heat::Heat, sub::SubSurface, stop, ssub)
+    Δk = CryoGrid.thickness(sub, ssub, first) # using `thickness` allows for generic layer implementations
     @inbounds let Tupper=boundaryvalue(bc,top,heat,sub,stop,ssub),
         Tsub=ssub.T[1],
         k=ssub.k[1],
-        δ=Δk[1]/2; # distance to boundary
+        δ=Δk/2; # distance to boundary
         -k*(Tsub-Tupper)/δ
     end
 end
-@inline function boundaryflux(::Dirichlet, bc::HeatBC, bot::Bottom, heat::Heat, sub::SubSurface, sbot, ssub)
-    Δk = thickness(sub, ssub) # using `thickness` allows for generic layer implementations
+@inline function CryoGrid.boundaryflux(::Dirichlet, bc::HeatBC, bot::Bottom, heat::Heat, sub::SubSurface, sbot, ssub)
+    Δk = CryoGrid.thickness(sub, ssub, last) # using `thickness` allows for generic layer implementations
     @inbounds let Tlower=boundaryvalue(bc,bot,heat,sub,sbot,ssub),
         Tsub=ssub.T[end],
         k=ssub.k[end],
-        δ=Δk[end]/2; # distance to boundary
+        δ=Δk/2; # distance to boundary
         -k*(Tsub-Tlower)/δ
     end
 end
+# CFL not defined for free-water freeze curve
+CryoGrid.timestep(::SubSurface, heat::Heat{FreeWater,Enthalpy,CFL}, state) = Inf
 """
-Generic top interaction. Computes flux dH at top cell.
-"""
-function interact!(top::Top, bc::HeatBC, sub::SubSurface, heat::Heat, stop, ssub)
-    Δk = thickness(sub, ssub) # using `thickness` allows for generic layer implementations
-    # boundary flux
-    @setscalar ssub.dH_upper = boundaryflux(bc, top, heat, sub, stop, ssub)
-    @inbounds ssub.dH[1] += getscalar(ssub.dH_upper) / Δk[1]
-    return nothing # ensure no allocation
-end
-"""
-Generic bottom interaction. Computes flux dH at bottom cell.
-"""
-function interact!(sub::SubSurface, heat::Heat, bot::Bottom, bc::HeatBC, ssub, sbot)
-    Δk = thickness(sub, ssub) # using `thickness` allows for generic layer implementations
-    # boundary flux
-    @setscalar ssub.dH_lower = boundaryflux(bc, bot, heat, sub, sbot, ssub)
-    @inbounds ssub.dH[end] += getscalar(ssub.dH_lower) / Δk[end]
-    return nothing # ensure no allocation
-end
-"""
-Generic subsurface interaction. Computes flux dH at boundary between subsurface layers.
-"""
-function interact!(sub1::SubSurface, ::Heat, sub2::SubSurface, ::Heat, s1, s2)
-    Δk₁ = thickness(sub1, s1)
-    Δk₂ = thickness(sub2, s2)
-    # thermal conductivity between cells
-    @inbounds let k₁ = s1.kc[end],
-        k₂ = s2.kc[1],
-        Δ₁ = Δk₁[end],
-        Δ₂ = Δk₂[1];
-        k = harmonicmean(k₁, k₂, Δ₁, Δ₂);
-        s1.k[end] = s2.k[1] = k
-    end
-    # calculate heat flux between cells
-    Qᵢ = @inbounds let k₁ = s1.k[end],
-        k₂ = s2.k[1],
-        k = k₁ = k₂,
-        z₁ = midpoints(sub1, s1),
-        z₂ = midpoints(sub2, s2),
-        δ = z₂[1] - z₁[end];
-        k*(s2.T[1] - s1.T[end]) / δ
+    timestep(::SubSurface, ::Heat{Tfc,TPara,CFL}, state) where {TPara}
+
+Implementation of `timestep` for `Heat` using the Courant-Fredrichs-Lewy condition
+defined as: Δt_max = u*Δx^2, where`u` is the "characteristic velocity" which here
+is taken to be the diffusivity: `dHdT / kc`.
+"""
+@inline function CryoGrid.timestep(::SubSurface, heat::Heat{Tfc,TPara,CFL}, state) where {Tfc,TPara}
+    Δx = Δ(state.grid)
+    dtmax = Inf
+    @inbounds for i in 1:length(Δx)
+        dtmax = let u = state.dHdT[i] / state.kc[i],
+            Δx = Δx[i],
+            Δt = 0.5*u*Δx^2;
+            min(dtmax, Δt)
+        end
     end
-    # diagnostics
-    @setscalar s1.dH_lower = Qᵢ
-    @setscalar s2.dH_upper = -Qᵢ
-    # add fluxes scaled by grid cell size
-    @inbounds s1.dH[end] += Qᵢ / Δk₁[end]
-    @inbounds s2.dH[1] += -Qᵢ / Δk₂[1]
-    return nothing # ensure no allocation
+    dtmax = isfinite(dtmax) && dtmax > 0 ? dtmax : heat.dtlim.fallback_dt
+    return dtmax
 end
 # Free water freeze curve
 @inline function enthalpyinv(sub::SubSurface, heat::Heat{FreeWater,Enthalpy}, state, i)
-    let θtot = max(1e-8, totalwater(sub, state, i)),
-        H = state.H[i],
-        C = state.C[i],
-        L = heat.L,
-        Lθ = L*θtot,
-        I_t = H > Lθ,
-        I_f = H <= 0.0;
-        (I_t*(H-Lθ) + I_f*H)/C
-    end
+    f_hc = partial(heatcapacity, liquidwater, sub, heat, state, i)
+    return enthalpyinv(heat.freezecurve, f_hc, state.H[i], heat.prop.L, f_hc.θwi)
 end
-@inline function liquidwater(sub::SubSurface, heat::Heat{FreeWater,Enthalpy}, state, i)
-    let θtot = max(1e-8, totalwater(sub, state, i)),
-        H = state.H[i],
-        L = heat.L,
+@inline function enthalpyinv(::FreeWater, f_hc::Function, H, L, θtot)
+    let θtot = max(1e-8, θtot),
         Lθ = L*θtot,
         I_t = H > Lθ,
+        I_f = H <= 0.0,
         I_c = (H > 0.0) && (H <= Lθ);
-        (I_c*(H/Lθ) + I_t)θtot
+        # compute liquid water content -> heat capacity -> temperature
+        θw = (I_c*(H/Lθ) + I_t)θtot
+        C = f_hc(θw)
+        T = (I_t*(H-Lθ) + I_f*H)/C
+        return T, θw, C
     end
 end
 """
@@ -321,18 +328,14 @@ end
 
 Implementation of "free water" freezing characteristic for any subsurface layer.
 Assumes that `state` contains at least temperature (T), enthalpy (H), heat capacity (C),
-total water content (θw), and liquid water content (θl).
+total water content (θwi), and liquid water content (θw).
 """
 @inline function freezethaw!(sub::SubSurface, heat::Heat{FreeWater,Enthalpy}, state)
     @inbounds for i in 1:length(state.H)
-        # liquid water content = (total water content) * (liquid fraction)
-        state.θl[i] = liquidwater(sub, heat, state, i)
-        # update heat capacity
-        state.C[i] = C = heatcapacity(sub, heat, state, i)
-        # enthalpy inverse function
-        state.T[i] = enthalpyinv(sub, heat, state, i)
+        # update T, θw, C
+        state.T[i], state.θw[i], state.C[i] = enthalpyinv(sub, heat, state, i)
         # set dHdT (a.k.a dHdT)
-        state.dHdT[i] = state.T[i] ≈ 0.0 ? 1e8 : 1/C
+        state.dHdT[i] = state.T[i] ≈ 0.0 ? 1e6 : 1/state.C[i]
     end
     return nothing
 end

src/Physics/HeatConduction/heat_bc.jl

+const Dirichlet = CryoGrid.Dirichlet
+const Neumann = CryoGrid.Neumann
 # Boundary condition type aliases
 const HeatBC = BoundaryProcess{T} where {Heat<:T<:SubSurfaceProcess}
 const ConstantTemp = ConstantBC{Heat, Dirichlet,Float"°C"}
@@ -11,26 +13,26 @@ struct TemperatureGradient{E,F} <: BoundaryProcess{Heat}
     effect::E # effect
     TemperatureGradient(T::F, effect::E=nothing) where {F<:Forcing,E} = new{E,F}(T, effect)
 end
-BoundaryStyle(::Type{<:TemperatureGradient}) = Dirichlet()
-@inline boundaryvalue(bc::TemperatureGradient, l1, ::Heat, l2, s1, s2) = getscalar(s1.T_ub)
+CryoGrid.BoundaryStyle(::Type{<:TemperatureGradient}) = Dirichlet()
+@inline CryoGrid.boundaryvalue(bc::TemperatureGradient, l1, ::Heat, l2, s1, s2) = getscalar(s1.T_ub)
 
-variables(::Top, bc::TemperatureGradient) = (
+CryoGrid.variables(::Top, bc::TemperatureGradient) = (
     Diagnostic(:T_ub, Scalar, u"K"),
 )
-function diagnosticstep!(::Top, bc::TemperatureGradient, state)
+function CryoGrid.diagnosticstep!(::Top, bc::TemperatureGradient, state)
     @setscalar state.T_ub = bc.T(state.t)
 end
 
 Base.@kwdef struct NFactor{W,S} <: BoundaryEffect
-    nf::W = Param(1.0, bounds=(0.0,1.0)) # applied when Tair <= 0
-    nt::S = Param(1.0, bounds=(0.0,1.0)) # applied when Tair > 0
+    nf::W = Param(1.0, domain=0..1) # applied when Tair <= 0
+    nt::S = Param(1.0, domain=0..1) # applied when Tair > 0
 end
-variables(::Top, bc::TemperatureGradient{<:NFactor}) = (
+CryoGrid.variables(::Top, bc::TemperatureGradient{<:NFactor}) = (
     Diagnostic(:T_ub, Scalar, u"K"),
     Diagnostic(:nfactor, Scalar),
 )
 nfactor(Tair, nfw, nfs) = (Tair <= zero(Tair))*nfw + (Tair > zero(Tair))*nfs
-function diagnosticstep!(::Top, bc::TemperatureGradient{<:NFactor}, state)
+function CryoGrid.diagnosticstep!(::Top, bc::TemperatureGradient{<:NFactor}, state)
     nfw = bc.effect.nf
     nfs = bc.effect.nt
     Tair = bc.T(state.t)

src/Physics/Hydrology/Hydrology.jl

+module Hydrology
+
+import CryoGrid
+
+using CryoGrid
+using CryoGrid.Physics
+using CryoGrid.Numerics
+using CryoGrid.Numerics: nonlineardiffusion!
+
+using IfElse
+using ModelParameters
+using Unitful
+
+abstract type WaterFlowParameterization end
+
+Base.@kwdef struct WaterFlow{TPara<:WaterFlowParameterization,Tinit,TProp} <: CryoGrid.SubSurfaceProcess
+    para::TPara = BucketScheme()
+    prop::TProp = HydrologicalProperties()
+    init::Tinit = nothing # optional initialization scheme
+end
+
+include("water_bucket.jl")
+
+end

src/Physics/Hydrology/water_bucket.jl

+Base.@kwdef struct BucketScheme{Tfc,Tkwsat} <: WaterFlowParameterization
+    θfc::Tfc = Param(0.5, domain=0..1)
+    kw_sat::Tkwsat = Param(1e-5, domain=0..Inf, units=u"m/s")
+end
+fieldcapacity(::SubSurface, water::WaterFlow{<:BucketScheme}) = water.para.θfc
+hydraulicconductivity(::SubSurface, water::WaterFlow{<:BucketScheme}) = water.para.kw_sat
+
+CryoGrid.variables(::WaterFlow{<:BucketScheme}) = (
+    Prognostic(:θwi, OnGrid(Cells)),
+    Diagnostic(:θw, OnGrid(Cells)),
+    Diagnostic(:θw_sat, OnGrid(Cells)),
+    Diagnostic(:jw, OnGrid(Edges), u"1/s"),
+    Diagnostic(:kw, OnGrid(Edges), u"m/s"),
+    Diagnostic(:kwc, OnGrid(Cells), u"m/s"),
+)
+function CryoGrid.diagnosticstep!(sub::SubSurface, water::WaterFlow{<:BucketScheme}, state)
+    kw_sat = hydraulicconductivity(sub, water)
+    por = porosity(sub, water, state)
+    @. state.θw_sat = state.θw / por
+    @. state.kwc = state.θw_sat*kw_sat
+    harmonicmean!(state.kw[2:end-1], state.kwc, Δ(state.grids.kw))
+    state.kw[1] = state.kwc[1] 
+    state.kw[end] = state.kwc[1]
+end
+function CryoGrid.prognosticstep!(sub::SubSurface, water::WaterFlow{<:BucketScheme}, state)
+    # set downward fluxes according to kw and field capacity
+    θfc = fieldcapacity(sub, water)
+    @inbounds for i in 2:length(state.kw)-1 # note that the index is over grid *edges*
+        let θwi_up = state.θwi[i-1], # cell above edge i
+            θw_up = state.θw[i-1],
+            θi_up = θwi_up - θw_up,
+            θwi_lo = state.θwi[i], # cell below edge i
+            por_lo = porosity(sub, water, state, i);
+            # downward advective flux due to gravity
+            jw = -state.kw[i]*(θw_up >= θfc)*10^(-7*θi_up)
+            # limit flux based on i) available water in cell above and ii) free pore space in cell below
+            state.jw[i] = min(min(jw, θw_up), min(por_lo - θwi_lo, zero(θwi_lo)))
+        end
+    end
+    @. state.dθwi += (state.jw[2:end] - state.jw[1:end-1]) / Δ(state.grids.jw)
+end

src/Physics/Physics.jl

 module Physics
 
-import CryoGrid: Process, CoupledProcesses, BoundaryProcess, Coupled, Layer, Top, Bottom, SubSurface
-import CryoGrid: diagnosticstep!, initialcondition!, interact!, prognosticstep!, variables, callbacks, observe
-
+using CryoGrid
 using CryoGrid.Numerics
 using CryoGrid.Utils
 
-using Lazy: @>>
+import CryoGrid
+
 using Reexport
 using Unitful
 
-"""
-    totalwater(sub::SubSurface, state)
-    totalwater(sub::SubSurface, state)
-    totalwater(sub::SubSurface, state, i)
-
-Retrieves the total water content for the given layer at grid cell `i`, if specified.
-Defaults to using the scalar total water content defined on layer `sub`.
-"""
-@inline totalwater(sub::SubSurface, state) = state.θw
-@inline totalwater(sub::SubSurface, state, i) = Utils.getscalar(totalwater(sub, state), i)
+export volumetricfractions, waterice, liquidwater, partial
 
-include("coupled.jl")
+include("common.jl")
 include("Boundaries/Boundaries.jl")
-include("HeatConduction/HeatConduction.jl")
-include("WaterBalance/WaterBalance.jl")
-include("Soils/Soils.jl")
-include("SEB/SEB.jl")
-include("Sources/Sources.jl")
-
 @reexport using .Boundaries
+include("HeatConduction/HeatConduction.jl")
 @reexport using .HeatConduction
-@reexport using .WaterBalance
+include("Hydrology/Hydrology.jl")
+@reexport using .Hydrology
+include("Snow/Snow.jl")
+@reexport using .Snow
+include("Soils/Soils.jl")
 @reexport using .Soils
+include("SEB/SEB.jl")
 @reexport using .SEB
+include("Sources/Sources.jl")
 @reexport using .Sources
 
 end
\ No newline at end of file

src/Physics/Snow/Snow.jl

+module Snow
+
+using CryoGrid
+using CryoGrid: ContinuousEvent, Increasing, Decreasing # for events/callbacks
+using CryoGrid.InputOutput: Forcing
+using CryoGrid.Physics.HeatConduction
+using CryoGrid.Numerics
+using CryoGrid.Utils
+
+import CryoGrid
+import CryoGrid.Physics
+import CryoGrid.Physics.HeatConduction
+
+using IfElse
+using ModelParameters
+using Unitful
+
+export Snowpack, SnowProperties, SnowMassBalance, Snowfall
+
+SnowProperties(
+    consts=Physics.Constants();
+    ρw = consts.ρw,
+    ρsn_new = Param(250.0, units=u"kg/m^3"),
+    ρsn_old = Param(500.0, units=u"kg/m^3"),
+) = (; ρw, ρsn_new, ρsn_old)
+
+"""
+    SnowpackParameterization
+
+Base type for snowpack paramterization schemes.
+"""
+abstract type SnowpackParameterization <: CryoGrid.Parameterization end
+
+"""
+    Snowpack{Tpara<:SnowpackParameterization,Tprop,Tsp} <: CryoGrid.SubSurface
+
+Generic representation of a ground surface snow pack.
+"""
+Base.@kwdef struct Snowpack{Tpara<:SnowpackParameterization,Tprop,Tsp} <: CryoGrid.SubSurface
+    para::Tpara = Bulk()
+    prop::Tprop = SnowProperties()
+    sp::Tsp = nothing
+end
+
+abstract type SnowAblationScheme end
+Base.@kwdef struct DegreeDayMelt{Tfactor,Tmax} <: SnowAblationScheme
+    factor::Tfactor = 5.0u"mm/K/d"
+    max_unfrozen::Tmax = 0.5
+end
+
+abstract type SnowAccumulationScheme end
+Base.@kwdef struct LinearAccumulation{S} <: SnowAccumulationScheme
+    rate_scale::S = Param(1.0, domain=0..Inf) # scaling factor for snowfall rate
+end
+
+abstract type SnowDensityScheme end
+# constant density (using Snowpack properties)
+struct ConstantDensity end
+
+abstract type SnowMassParameterization end
+Base.@kwdef struct Prescribed{Tswe,Tρsn} <: SnowMassParameterization
+    swe::Tswe = 0.0u"m" # depth snow water equivalent [m]
+    ρsn::Tρsn = 250.0u"kg/m^3" # snow density [kg m^-3]
+end
+Base.@kwdef struct Dynamic{TAcc,TAbl,TDen} <: SnowMassParameterization
+    accumulation::TAcc = LinearAccumulation()
+    ablation::TAbl = DegreeDayMelt()
+    density::TDen = ConstantDensity()
+end
+
+Base.@kwdef struct SnowMassBalance{Tpara} <: CryoGrid.SubSurfaceProcess
+    para::Tpara = Dynamic()
+end
+accumulation(snow::SnowMassBalance{<:Dynamic}) = snow.para.accumulation
+ablation(snow::SnowMassBalance{<:Dynamic}) = snow.para.ablation
+density(snow::SnowMassBalance{<:Dynamic}) = snow.para.density
+
+# convenience type aliases
+"""
+    PrescribedSnowMassBalance{Tswe,Tρsn} = SnowMassBalance{Prescribed{Tswe,Tρsn}} where {Tswe,Tρsn}
+"""
+const PrescribedSnowMassBalance{Tswe,Tρsn} = SnowMassBalance{Prescribed{Tswe,Tρsn}} where {Tswe,Tρsn}
+"""
+    DynamicSnowMassBalance{TAcc,TAbl,TDen} = SnowMassBalance{Dynamic{TAcc,TAbl,TDen}} where {TAcc,TAbl,TDen}
+"""
+const DynamicSnowMassBalance{TAcc,TAbl,TDen} = SnowMassBalance{Dynamic{TAcc,TAbl,TDen}} where {TAcc,TAbl,TDen}
+
+snowvariables(::Snowpack, ::SnowMassBalance) = (
+    Diagnostic(:dsn, Scalar, u"m", domain=0..Inf),
+    Diagnostic(:T_ub, Scalar, u"°C"),
+)
+
+swe(::Snowpack, ::SnowMassBalance, state) = state.swe
+swe(::Snowpack, smb::SnowMassBalance{<:Prescribed}, state) = smb.para.swe
+swe(::Snowpack, smb::SnowMassBalance{<:Prescribed{<:Forcing}}, state) = smb.para.swe(state.t)
+snowdensity(::Snowpack, ::SnowMassBalance, state) = state.ρsn
+snowdensity(::Snowpack, smb::SnowMassBalance{<:Prescribed}, state) = smb.para.ρsn
+snowdensity(::Snowpack, smb::SnowMassBalance{<:Prescribed{Tswe,<:Forcing}}, state) where {Tswe} = smb.para.ρsn(state.t)
+
+# Boundary conditions
+
+struct Snowfall{Tsn<:Forcing} <: BoundaryProcess{SnowMassBalance}
+    snowfall::Tsn
+end
+CryoGrid.BoundaryStyle(::Snowfall) = CryoGrid.Neumann()
+@inline boundaryvalue(bc::Snowfall, ::Top, ::SnowMassBalance, ::Snowpack, s1, s2) = bc.snowfall(s1.t)
+
+# Implementations
+
+# for prescribed snow depth/density, the mass balance is given so we do not need to do anything here
+CryoGrid.prognosticstep!(::Snowpack, ::SnowMassBalance{<:Prescribed}, ssnow) = nothing
+
+include("snow_bulk.jl")
+
+end
\ No newline at end of file

src/Physics/Snow/snow_bulk.jl

+"""
+    Bulk{Tthresh} <: SnowpackParameterization
+
+Simple, bulk ("single layer") snow scheme where snowpack is represented as a single grid cell with homogenous state.
+"""
+Base.@kwdef struct Bulk{Tthresh} <: SnowpackParameterization
+    thresh::Tthresh = 0.02u"m" # snow threshold
+end
+
+"""
+    BulkSnowpack = Snowpack{<:Bulk}
+
+Type alias for Snowpack with `Bulk` parameterization.
+"""
+const BulkSnowpack = Snowpack{<:Bulk}
+# Local alias for HeatConduction Enthalpy type
+const Enthalpy = HeatConduction.Enthalpy
+
+threshold(snow::BulkSnowpack) = snow.para.thresh
+
+CryoGrid.thickness(::BulkSnowpack, state, i::Integer=1) = getscalar(state.dsn)
+CryoGrid.midpoint(::BulkSnowpack, state, i::Integer=1) = -getscalar(state.dsn) / 2
+
+# Events
+CryoGrid.events(::BulkSnowpack, ::Coupled2{<:SnowMassBalance,<:Heat}) = (
+    ContinuousEvent(:snow_min),
+)
+# critierion for minimum snow threshold event
+function CryoGrid.criterion(
+    ::ContinuousEvent{:snow_min},
+    snow::BulkSnowpack,
+    smb::SnowMassBalance,
+    state,
+)
+    # get current snow water equivalent (from forcing or state variable)
+    new_swe = getscalar(swe(snow, smb, state))
+    # get current snow density
+    new_ρsn = getscalar(snowdensity(snow, smb, state))
+    # compute actual snow depth/height
+    new_dsn = new_swe*snow.prop.ρw/new_ρsn
+    # use threshold adjusted depth as residual;
+    # i.e. the event will fire when snow depth crosses this threshold.
+    return new_dsn - threshold(snow)
+end
+# triggers for minimum snow threshold event
+function CryoGrid.trigger!(
+    ::ContinuousEvent{:snow_min},
+    ::Decreasing,
+    snow::BulkSnowpack,
+    ::Coupled2{<:SnowMassBalance,<:Heat},
+    state
+)
+    # Case 1: Decreasing snow depth; set everything to zero to remove snowpack
+    state.H .= 0.0
+    state.T .= 0.0
+    state.θwi .= 0.0
+    state.swe .= 0.0
+    state.dsn .= 0.0
+end
+function CryoGrid.trigger!(
+    ::ContinuousEvent{:snow_min},
+    ::Increasing,
+    snow::BulkSnowpack,
+    procs::Coupled2{<:SnowMassBalance,<:Heat},
+    state
+)
+    # Case 2: Increasing snow depth; initialize temperature and enthalpy state
+    # using current upper boundary temperature.
+    _, heat = procs
+    θfracs = volumetricfractions(snow, heat, state, 1)
+    state.C .= C = HeatConduction.heatcapacity(snow, heat, θfracs...)
+    state.T .= state.T_ub
+    state.H .= state.T.*C
+end
+# heat upper boundary (for all bulk implementations)
+CryoGrid.interact!(top::Top, bc::HeatBC, snow::BulkSnowpack, heat::Heat, stop, ssnow) = CryoGrid.interact!(CryoGrid.BoundaryStyle(bc), top, bc, snow, heat, stop, ssnow)
+function CryoGrid.interact!(::CryoGrid.Dirichlet, top::Top, bc::HeatBC, snow::BulkSnowpack, heat::Heat, stop, ssnow)
+    Δk = CryoGrid.thickness(snow, ssnow) # using `thickness` allows for generic layer implementations
+    @setscalar ssnow.T_ub = CryoGrid.boundaryvalue(bc, top, heat, snow, stop, ssnow)
+    if getscalar(ssnow.dsn) < threshold(snow)
+        @setscalar ssnow.T = getscalar(ssnow.T_ub)
+    end
+    # boundary flux
+    @setscalar ssnow.dH_upper = CryoGrid.boundaryflux(bc, top, heat, snow, stop, ssnow)
+    @inbounds ssnow.dH[1] += getscalar(ssnow.dH_upper) / Δk[1]
+    return nothing # ensure no allocation
+end
+
+# ==== Dynamic bulk snow scheme ==== #
+
+CryoGrid.variables(snow::BulkSnowpack, smb::DynamicSnowMassBalance) = (
+    Prognostic(:swe, Scalar, u"m", domain=0..Inf),
+    Diagnostic(:ρsn, Scalar, u"kg/m^3", domain=0..Inf),
+    Diagnostic(:θwi, OnGrid(Cells), u"kg/m^3", domain=0..1), 
+    snowvariables(snow, smb)...,
+)
+function CryoGrid.diagnosticstep!(
+    snow::BulkSnowpack,
+    procs::Coupled2{<:DynamicSnowMassBalance{TAcc,TAbl,TDen},<:Heat{FreeWater,Enthalpy}},
+    state
+) where {TAcc,TAbl<:DegreeDayMelt,TDen<:ConstantDensity}
+    smb, heat = procs
+    ρsn = snow.prop.ρsn_new
+    @setscalar state.θwi = θwi = ρsn / snow.prop.ρw
+    @setscalar state.ρsn = ρsn
+    dsn = getscalar(state.swe) / θwi
+    # only update snowdepth if swe greater than threshold, otherwise, set to zero.
+    @setscalar state.dsn = IfElse.ifelse(getscalar(state.swe) >= threshold(snow)*θwi, dsn, zero(dsn))
+    # get heat capacity as a function of liquid water content
+    f_hc = partial(heatcapacity, liquidwater, snow, heat, state, 1)
+    # set temperature and liquid water content according to free water freeze curve,
+    # but capping the liquid fraction according to the 'max_unfrozen' parameter.
+    max_unfrozen = ablation(smb).max_unfrozen
+    θwi_cap = θwi*max_unfrozen
+    T, θw, C = HeatConduction.enthalpyinv(heat.freezecurve, f_hc, getscalar(state.H), heat.prop.L, θwi_cap)
+    # do not allow temperature to exceed 0°C
+    @. state.T = min(T, zero(T))
+    @. state.θw = θw
+    @. state.C = C
+    # compute thermal conductivity
+    HeatConduction.thermalconductivity!(snow, heat, state)
+    @. state.k = state.kc
+end
+# snowfall upper boundary
+function CryoGrid.interact!(
+    top::Top,
+    bc::Snowfall,
+    snow::BulkSnowpack,
+    smb::DynamicSnowMassBalance{<:LinearAccumulation},
+    stop,
+    ssnow
+)
+    # upper boundary condition for snow mass balance;
+    # apply snowfall to swe
+    rate_scale = accumulation(smb).rate_scale
+    snowfall_rate = boundaryvalue(bc, top, smb, snow, stop, ssnow)
+    @. ssnow.dswe += rate_scale*snowfall_rate
+end
+function CryoGrid.prognosticstep!(
+    snow::BulkSnowpack,
+    smb::DynamicSnowMassBalance{TAcc,TAbl},
+    state
+) where {TAcc,TAbl<:DegreeDayMelt}
+    if getscalar(state.dsn) < threshold(snow)
+        # set energy flux to zero if there is no snow
+        @. state.dH = zero(eltype(state.H))
+    end
+    if getscalar(state.swe) > 0.0 && getscalar(state.T_ub) > 0.0
+        ddf = ablation(smb).factor # [m/K/s]
+        dH_upper = getscalar(state.dH_upper) # [J/m^3]
+        T_ub = getscalar(state.T_ub) # upper boundary temperature
+        Tref = 0.0*unit(T_ub) # just in case T_ub has units
+        # calculate the melt rate per second via the degree day model
+        dmelt = max(ddf*(T_ub-Tref), zero(dH_upper)) # [m/s]
+        @. state.dswe += -dmelt
+        # remove upper heat flux from dH if dmelt > 0;
+        # this is due to the energy being "consumed" to melt the snow
+        @. state.dH += -dH_upper*(dmelt > zero(dmelt))
+    end
+    return nothing
+end
+
+# ==== Prescribed/forced snow scheme ==== #
+
+# Snow mass balance
+CryoGrid.variables(snow::BulkSnowpack, smb::PrescribedSnowMassBalance) = (
+    Diagnostic(:swe, Scalar, u"m", domain=0..Inf),
+    Diagnostic(:ρsn, Scalar, u"kg/m^3", domain=0..Inf),
+    Diagnostic(:θwi, OnGrid(Cells), u"kg/m^3", domain=0..1),
+    snowvariables(snow, smb)...,
+)
+CryoGrid.events(::BulkSnowpack, ::Coupled2{<:PrescribedSnowMassBalance,<:Heat}) = (
+    ContinuousEvent(:snow_min),
+)
+function CryoGrid.criterion(
+    ::ContinuousEvent{:snow_min},
+    snow::BulkSnowpack,
+    smb::PrescribedSnowMassBalance,
+    state,
+)
+    new_swe = swe(snow, smb, state)
+    new_ρsn = snowdensity(snow, smb, state)
+    new_dsn = new_swe*snow.prop.ρw/new_ρsn
+    return new_dsn - threshold(snow)
+end
+function CryoGrid.trigger!(
+    ::ContinuousEvent{:snow_min},
+    ::Decreasing,
+    snow::BulkSnowpack,
+    ::Coupled2{<:PrescribedSnowMassBalance,<:Heat},
+    state
+)
+    state.H .= 0.0
+end
+function CryoGrid.trigger!(
+    ::ContinuousEvent{:snow_min},
+    ::Increasing,
+    snow::BulkSnowpack,
+    procs::Coupled2{<:PrescribedSnowMassBalance,<:Heat},
+    state
+)
+    _, heat = procs
+    θfracs = volumetricfractions(snow, heat, state, 1)
+    C = HeatConduction.heatcapacity(snow, heat, θfracs...)
+    state.T .= state.T_ub
+    state.H .= state.T.*C
+end
+function CryoGrid.diagnosticstep!(
+    snow::BulkSnowpack,
+    procs::Coupled2{<:PrescribedSnowMassBalance,<:Heat{FreeWater,Enthalpy}},
+    state
+)
+    smb, heat = procs
+    ρw = snow.prop.ρw
+    new_swe = swe(snow, smb, state)
+    new_ρsn = snowdensity(snow, smb, state)
+    new_dsn = new_swe*ρw/new_ρsn
+    ρw = heat.prop.ρw
+    if new_dsn > threshold(snow)
+        # if new snow depth is above threshold, set state variables
+        @setscalar state.swe = new_swe
+        @setscalar state.ρsn = new_ρsn
+        @setscalar state.dsn = new_dsn
+        @. state.θwi = new_ρsn / ρw
+        HeatConduction.freezethaw!(snow, heat, state)
+        # cap temperature at 0°C
+        @. state.T = min(state.T, zero(eltype(state.T)))
+        HeatConduction.thermalconductivity!(snow, heat, state)
+        @. state.k = state.kc
+    else
+        # otherwise, set to zero
+        @setscalar state.swe = 0.0
+        @setscalar state.ρsn = 0.0
+        @setscalar state.dsn = 0.0
+        @setscalar state.θwi = 0.0
+        @setscalar state.θw = 0.0
+        @setscalar state.C = heat.prop.ca
+        @setscalar state.kc = heat.prop.ka
+    end
+end
+# override prognosticstep! for incompatible heat types to prevent incorrect usage;
+# this will actually result in an ambiguous dispatch error before this error is thrown.
+CryoGrid.prognosticstep!(snow::BulkSnowpack, heat::Heat, state) = error("prognosticstep! not implemented for $(typeof(heat)) on $(typeof(snow))")
+# prognosticstep! for free water, enthalpy based Heat on snow layer
+function CryoGrid.prognosticstep!(snow::BulkSnowpack, ::Heat{FreeWater,Enthalpy}, state)
+    dsn = getscalar(state.dsn)
+    if dsn < snow.para.thresh
+        # set energy flux to zero if there is no snow
+        @. state.dH = zero(eltype(state.H))
+    end
+end

src/Physics/Soils/Soils.jl

 module Soils
 
-import CryoGrid: SubSurface, Parameterization
-import CryoGrid: initialcondition!, variables
-import CryoGrid.Physics: totalwater
-import CryoGrid.Physics.HeatConduction: Enthalpy, Temperature, liquidwater, thermalconductivity, heatcapacity, freezethaw!, enthalpyinv
-
+using CryoGrid: SubSurface, Parameterization
 using CryoGrid.Numerics
 using CryoGrid.Numerics: heaviside
+using CryoGrid.Physics
 using CryoGrid.Physics.HeatConduction
-using CryoGrid.Physics.WaterBalance
+using CryoGrid.Physics.Hydrology
 using CryoGrid.Utils
 
 using Base: @propagate_inbounds, @kwdef
 using IfElse
+using Interpolations: Interpolations
+using IntervalSets
 using ForwardDiff
 using ModelParameters
 using Unitful
 
-import Interpolations
+import CryoGrid
+import CryoGrid.Physics
+import CryoGrid.Physics.HeatConduction
+
 import Flatten: @flattenable, flattenable
 
 export Soil, SoilParameterization, CharacteristicFractions, SoilProfile
-export soilcomponent, porosity, mineral, organic
+export soilparameters, soilcomponent, porosity, mineral, organic
+
+const Enthalpy = HeatConduction.Enthalpy
+const Temperature = HeatConduction.Temperature
 
 """
     SoilComposition
@@ -36,95 +41,122 @@ struct Heterogeneous <: SoilComposition end
 
 Abstract base type for parameterizations of soil properties.
 """
-abstract type SoilParameterization <: Parameterization end
+abstract type SoilParameterization end
 """
     CharacteristicFractions{P1,P2,P3,P4} <: SoilParameterization
 
 Represents uniform composition of a soil volume in terms of fractions: excess ice, natural porosity, saturation, and organic solid fraction.
 """
 Base.@kwdef struct CharacteristicFractions{P1,P2,P3,P4} <: SoilParameterization
-    xic::P1 = 0.0 # excess ice fraction
-    por::P2 = 0.5 # natural porosity
-    sat::P3 = 1.0 # saturation
-    org::P4 = 0.0 # organic fraction of solid; mineral fraction is 1-org
+    xic::P1 = Param(0.0, domain=0..1) # excess ice fraction
+    por::P2 = Param(0.5, domain=0..1) # natural porosity
+    sat::P3 = Param(1.0, domain=0..1) # saturation
+    org::P4 = Param(0.0, domain=0..1) # organic fraction of solid; mineral fraction is 1-org
 end
+soilparameters(::Type{CharacteristicFractions}=CharacteristicFractions; xic, por, sat, org) = CharacteristicFractions(xic=Param(xic, domain=0..1), por=Param(por, domain=0..1), sat=Param(sat, domain=0..1), org=Param(org, domain=0..1))
 # Type alias for CharacteristicFractions with all scalar/numeric constituents
 const HomogeneousCharacteristicFractions = CharacteristicFractions{<:Number,<:Number,<:Number,<:Number}
 SoilProfile(pairs::Pair{<:DistQuantity,<:SoilParameterization}...) = Profile(pairs...)
 # Helper functions for obtaining soil compositions from characteristic fractions.
 soilcomponent(::Val{var}, para::CharacteristicFractions) where var = soilcomponent(Val{var}(), para.xic, para.por, para.sat, para.org)
 soilcomponent(::Val{:θp}, χ, ϕ, θ, ω) = (1-χ)*ϕ
-soilcomponent(::Val{:θw}, χ, ϕ, θ, ω) = χ + (1-χ)*ϕ*θ
+soilcomponent(::Val{:θwi}, χ, ϕ, θ, ω) = χ + (1-χ)*ϕ*θ
 soilcomponent(::Val{:θm}, χ, ϕ, θ, ω) = (1-χ)*(1-ϕ)*(1-ω)
 soilcomponent(::Val{:θo}, χ, ϕ, θ, ω) = (1-χ)*(1-ϕ)*ω
 """
 Soil thermal properties.
 """
-@kwdef struct SoilThermalProperties{Tko,Tkm,Tco,Tcm} <: HeatConduction.ThermalProperties
-    ko::Tko = Param(0.25, units=u"W/m/K") # organic [Hillel(1982)]
-    km::Tkm = Param(3.8, units=u"W/m/K") # mineral [Hillel(1982)]
-    co::Tco = Param(2.5e6, units=u"J/K/m^3") # heat capacity organic
-    cm::Tcm = Param(2.0e6, units=u"J/K/m^3") # heat capacity mineral
-end
+SoilThermalProperties(;
+    ko = Param(0.25, units=u"W/m/K"), # organic [Hillel(1982)]
+    km = Param(3.8, units=u"W/m/K"), # mineral [Hillel(1982)]
+    co = Param(2.5e6, units=u"J/K/m^3"), # heat capacity organic
+    cm = Param(2.0e6, units=u"J/K/m^3"), # heat capacity mineral
+) = (; ko, km, co, cm)
 """
 Basic Soil layer.
 """
-@kwdef struct Soil{TPara<:SoilParameterization,TProp<:SoilThermalProperties,TSp} <: SubSurface
+@kwdef struct Soil{TPara<:SoilParameterization,TProp,TSp} <: SubSurface
     para::TPara = CharacteristicFractions()
     prop::TProp = SoilThermalProperties()
     sp::TSp = nothing # user-defined specialization
 end
+HeatConduction.thermalproperties(soil::Soil) = soil.prop
 # SoilComposition trait impl
 SoilComposition(soil::Soil) = SoilComposition(typeof(soil))
 SoilComposition(::Type{<:Soil}) = Heterogeneous()
 SoilComposition(::Type{<:Soil{<:HomogeneousCharacteristicFractions}}) = Homogeneous()
-# Fixed volumetric content methods
-totalwater(soil::Soil, state) = totalwater(SoilComposition(soil), soil)
+# Volumetric fraction methods
+Physics.waterice(soil::Soil, state) = waterice(SoilComposition(soil), soil)
 porosity(soil::Soil, state) = porosity(SoilComposition(soil), soil)
 mineral(soil::Soil, state) = mineral(SoilComposition(soil), soil)
 organic(soil::Soil, state) = organic(SoilComposition(soil), soil)
 ## Homogeneous soils
-totalwater(::Homogeneous, soil::Soil{<:CharacteristicFractions}) = soilcomponent(Val{:θw}(), soil.para)
+Physics.waterice(::Homogeneous, soil::Soil{<:CharacteristicFractions}) = soilcomponent(Val{:θwi}(), soil.para)
 porosity(::Homogeneous, soil::Soil{<:CharacteristicFractions}) = soilcomponent(Val{:θp}(), soil.para)
 mineral(::Homogeneous, soil::Soil{<:CharacteristicFractions}) = soilcomponent(Val{:θm}(), soil.para)
 organic(::Homogeneous, soil::Soil{<:CharacteristicFractions}) = soilcomponent(Val{:θo}(), soil.para)
 ## Heterogeneous soils
-totalwater(::Heterogeneous, soil::Soil, state) = state.θw   
+Physics.waterice(::Heterogeneous, soil::Soil, state) = state.θwi   
 porosity(::Heterogeneous, soil::Soil, state) = state.θp
 mineral(::Heterogeneous, soil::Soil, state) = state.θm
 organic(::Heterogeneous, soil::Soil, state) = state.θo
 
-variables(soil::Soil) = variables(SoilComposition(soil), soil)
-variables(::Homogeneous, ::Soil) = ()
-variables(::Heterogeneous, ::Soil) = (
-    Diagnostic(:θw, Float64, OnGrid(Cells)),
-    Diagnostic(:θp, Float64, OnGrid(Cells)),
-    Diagnostic(:θm, Float64, OnGrid(Cells)),
-    Diagnostic(:θo, Float64, OnGrid(Cells)),
+CryoGrid.variables(soil::Soil) = CryoGrid.variables(SoilComposition(soil), soil)
+CryoGrid.variables(::Homogeneous, ::Soil) = ()
+CryoGrid.variables(::Heterogeneous, ::Soil) = (
+    Diagnostic(:θwi, OnGrid(Cells), domain=0..1),
+    Diagnostic(:θp, OnGrid(Cells), domain=0..1),
+    Diagnostic(:θm, OnGrid(Cells), domain=0..1),
+    Diagnostic(:θo, OnGrid(Cells), domain=0..1),
 )
 
-initialcondition!(soil::Soil, state) = initialcondition!(SoilComposition(soil), soil, state)
-initialcondition!(::Homogeneous, ::Soil, state) = nothing
+CryoGrid.initialcondition!(soil::Soil, state) = CryoGrid.initialcondition!(SoilComposition(soil), soil, state)
+CryoGrid.initialcondition!(::Homogeneous, ::Soil, state) = nothing
 """
     initialcondition!(::Heterogeneous, soil::Soil{<:CharacteristicFractions}, state)
 
 Default implementation of initialcondition! for heterogeneous soils parameterized by characteristic fractions.
 Fields of `soil.para` may be either numbers (including `Param`s) or `Function`s of the form `f(::Soil, state)::AbstractVector`.
 """
-function initialcondition!(::Heterogeneous, soil::Soil{<:CharacteristicFractions}, state)
+function CryoGrid.initialcondition!(::Heterogeneous, soil::Soil{<:CharacteristicFractions}, state)
     evaluate(x::Number) = x
     evaluate(f::Function) = f(soil, state)
     χ = evaluate(soil.para.xic)
     ϕ = evaluate(soil.para.por)
     θ = evaluate(soil.para.sat)
     ω = evaluate(soil.para.org)
-    @. state.θw = soilcomponent(Val{:θw}(), χ, ϕ, θ, ω)
+    @. state.θwi = soilcomponent(Val{:θwi}(), χ, ϕ, θ, ω)
     @. state.θp = soilcomponent(Val{:θp}(), χ, ϕ, θ, ω)
     @. state.θm = soilcomponent(Val{:θm}(), χ, ϕ, θ, ω)
     @. state.θo = soilcomponent(Val{:θo}(), χ, ϕ, θ, ω)
 end
 
+"""
+    mineral(soil::Soil, state, i)
+
+Retrieves the mineral content for the given layer at grid cell `i`, if provided.
+Defaults to using the scalar mineral content defined on `soil`.
+"""
+@inline mineral(soil::Soil, state, i) = Utils.getscalar(mineral(soil, state), i)
+"""
+    organic(soil::Soil, state, i)
+
+Retrieves the organic content for the given layer at grid cell `i`, if provided.
+Defaults to using the scalar organic content defined on `soil`.
+"""
+@inline organic(soil::Soil, state, i) = Utils.getscalar(organic(soil, state), i)
+"""
+    porosity(soil::Soil, state, i)
+
+Retrieves the porosity for the given layer at grid cell `i`, if provided.
+Defaults to using the scalar porosity defined on `soil`.
+"""
+@inline porosity(soil::Soil, state, i) = Utils.getscalar(porosity(soil, state), i)
+
+export SWRC, VanGenuchten
+include("sfcc/swrc.jl")
 export SFCC, DallAmico, Westermann, McKenzie, SFCCNewtonSolver, SFCCPreSolver
+include("sfcc/sfcc.jl")
 include("soilheat.jl")
 
 end

src/Physics/Soils/sfcc.jl → src/Physics/Soils/sfcc/sfcc.jl

@@ -17,20 +17,17 @@ Generic representation of the soil freeze characteristic curve. The shape and pa
 of the curve are determined by the implementation of SFCCFunction `f`. Also requires
 an implementation of SFCCSolver which provides the solution to the non-linear mapping H <--> T.
 """
-@flattenable struct SFCC{F,S} <: FreezeCurve
-    f::F | true # freeze curve function f: (T,...) -> θ
-    solver::S | true # solver for H -> T or T -> H
+struct SFCC{F,S} <: FreezeCurve
+    f::F # freeze curve function f: (T,...) -> θ
+    solver::S # solver for H -> T or T -> H
     SFCC(f::F, s::S=SFCCPreSolver()) where {F<:SFCCFunction,S<:SFCCSolver} = new{F,S}(f,s)
 end
 
 # Join the declared state variables of the SFCC function and the solver
-variables(soil::Soil, heat::Heat, sfcc::SFCC) = tuplejoin(variables(soil, heat, sfcc.f), variables(soil, heat, sfcc.solver))
+CryoGrid.variables(soil::Soil, heat::Heat, sfcc::SFCC) = tuplejoin(CryoGrid.variables(soil, heat, sfcc.f), CryoGrid.variables(soil, heat, sfcc.solver))
 # Default SFCC initialization
-function initialcondition!(soil::Soil, heat::Heat, sfcc::SFCC, state)
-    L = heat.L
-    state.θl .= sfcc.f.(state.T, sfccargs(sfcc.f, soil, heat, state)...)
-    heatcapacity!(soil, heat, state)
-    @. state.H = enthalpy(state.T, state.C, L, state.θl)
+function CryoGrid.initialcondition!(soil::Soil, heat::Heat, sfcc::SFCC, state)
+    HeatConduction.freezethaw!(soil, heat, state)
 end
 
 """
@@ -42,38 +39,35 @@ of the freeze curve function `f`.
 """
 sfccargs(::SFCCFunction, ::Soil, ::Heat, state) = ()
 # Fallback implementation of variables for SFCCFunction
-variables(::Soil, ::Heat, f::SFCCFunction) = ()
-variables(::Soil, ::Heat, s::SFCCSolver) = ()
+CryoGrid.variables(::Soil, ::Heat, f::SFCCFunction) = ()
+CryoGrid.variables(::Soil, ::Heat, s::SFCCSolver) = ()
 """
     DallAmico <: SFCCFunction
 
 Dall'Amico M, 2010. Coupled water and heat transfer in permafrost modeling. Ph.D. Thesis, University of Trento, pp. 43.
 """
-Base.@kwdef struct DallAmico{T,Θ,A,N,G} <: SFCCFunction
+Base.@kwdef struct DallAmico{T,Θ,G,Tvg<:VanGenuchten} <: SFCCFunction
     Tₘ::T = Param(0.0, units=u"°C")
-    θres::Θ = Param(0.0, bounds=(0,1))
-    α::A = Param(4.0, bounds=(eps(),Inf), units=u"1/m")
-    n::N = Param(2.0, bounds=(1,Inf))
+    θres::Θ = Param(0.0, domain=0..1)
     g::G = 9.80665u"m/s^2" # acceleration due to gravity
-    swrc::VanGenuchten = VanGenuchten()
+    swrc::Tvg = VanGenuchten()
 end
 sfccargs(f::DallAmico, soil::Soil, heat::Heat, state) = (
     porosity(soil, state), # θ saturated = porosity
-    totalwater(soil, state), # total water content
+    waterice(soil, heat, state), # total water content
     heat.prop.Lf, # specific latent heat of fusion, L
     f.θres,
     f.Tₘ,
-    f.α,
-    f.n,
+    f.swrc.α,
+    f.swrc.n,
 )
 # pressure head at T
 @inline ψ(T,Tstar,ψ₀,Lf,g) = ψ₀ + Lf/(g*Tstar)*(T-Tstar)*heaviside(Tstar-T)
-@inline function (f::DallAmico)(T, θsat, θtot, Lf, θres=f.θres, Tₘ=f.Tₘ, α=f.α, n=f.n)
+@inline function (f::DallAmico)(T, θsat, θtot, Lf, θres=f.θres, Tₘ=f.Tₘ, α=f.swrc.α, n=f.swrc.n)
     let θsat = max(θtot, θsat),
         g = f.g,
-        m = 1-1/n,
         Tₘ = normalize_temperature(Tₘ),
-        ψ₀ = IfElse.ifelse(θtot < θsat, -1/α*(((θtot-θres)/(θsat-θres))^(-1/m)-1.0)^(1/n), zero(1/α)),
+        ψ₀ = f.swrc(inv, θtot, θres, θsat, α, n),
         Tstar = Tₘ + g*Tₘ/Lf*ψ₀,
         T = normalize_temperature(T),
         ψ = ψ(T, Tstar, ψ₀, Lf, g);
@@ -90,12 +84,12 @@ McKenzie JM, Voss CI, Siegel DI, 2007. Groundwater flow with energy transport an
 """
 Base.@kwdef struct McKenzie{T,Θ,Γ} <: SFCCFunction
     Tₘ::T = Param(0.0, units=u"°C")
-    θres::Θ = Param(0.0, bounds=(0,1))
-    γ::Γ = Param(0.1, bounds=(eps(),Inf), units=u"K")
+    θres::Θ = Param(0.0, domain=0..1)
+    γ::Γ = Param(0.1, domain=0..Inf, units=u"K")
 end
 sfccargs(f::McKenzie, soil::Soil, heat::Heat, state) = (
     porosity(soil, state), # θ saturated = porosity
-    totalwater(soil, state), # total water content
+    waterice(soil, heat, state), # total water content
     f.θres,
     f.Tₘ,
     f.γ,
@@ -117,12 +111,12 @@ Westermann, S., Boike, J., Langer, M., Schuler, T. V., and Etzelmüller, B.: Mod
 """
 Base.@kwdef struct Westermann{T,Θ,Δ} <: SFCCFunction
     Tₘ::T = Param(0.0, units=u"°C")
-    θres::Θ = Param(0.0, bounds=(0,1))
-    δ::Δ = Param(0.1, bounds=(eps(),Inf), units=u"K")
+    θres::Θ = Param(0.0, domain=0..1)
+    δ::Δ = Param(0.1, domain=0..Inf, units=u"K")
 end
 sfccargs(f::Westermann, soil::Soil, heat::Heat, state) = (
     porosity(soil, state), # θ saturated = porosity
-    totalwater(soil, state), # total water content
+    waterice(soil, heat, state), # total water content
     f.θres,
     f.Tₘ,
     f.δ,

src/Physics/Soils/sfcc_solvers.jl → src/Physics/Soils/sfcc/sfcc_solvers.jl

+using NLsolve
+
 """
     SFCCNewtonSolver <: SFCCSolver
 
@@ -14,32 +16,34 @@ Base.@kwdef struct SFCCNewtonSolver <: SFCCSolver
     α₀::Float64 = 1.0 # initial step size multiplier
     τ::Float64 = 0.7 # step size decay for backtracking
 end
-# Helper function for updating θl, C, and the residual.
-@inline function residual(soil::Soil, heat::Heat, T, H, L, f::F, f_args::Fargs, θw, θm, θo) where {F,Fargs}
-    θl = f(T, f_args...)
-    C = heatcapacity(soil, heat, θw, θl, θm, θo)
-    Tres = T - (H - θl*L) / C
-    return Tres, θl, C
+# Helper function for updating θw, C, and the residual.
+@inline function sfccresidual(soil::Soil, heat::Heat, f::F, f_args::Fargs, f_hc, T, H) where {F,Fargs}
+    L = heat.prop.L
+    θw = f(T, f_args...)
+    C = f_hc(θw)
+    Tres = T - (H - θw*L) / C
+    return Tres, θw, C
 end
 # Newton solver implementation
-function sfccsolve(solver::SFCCNewtonSolver, soil::Soil, heat::Heat, f, f_args, H, L, θw, θm, θo, T₀::Nothing=nothing)
-    T₀ = H / heatcapacity(soil, heat, θw, 0.0, θm, θo)
-    return sfccsolve(solver, soil, heat, f, f_args, H, L, θw, θm, θo, T₀)
+function sfccsolve(solver::SFCCNewtonSolver, soil::Soil, heat::Heat, f, f_args, f_hc, H, T₀::Nothing=nothing)
+    T₀ = IfElse.ifelse(H < zero(H), H / f_hc(0.0), zero(H))
+    return sfccsolve(solver, soil, heat, f, f_args, f_hc, H, T₀)
 end
 using ForwardDiff
-function sfccsolve(solver::SFCCNewtonSolver, soil::Soil, heat::Heat, f::F, f_args, H, L, θw, θm, θo, T₀) where {F}
+function sfccsolve(solver::SFCCNewtonSolver, soil::Soil, heat::Heat, f::F, f_args, f_hc, H, T₀) where {F}
     T = T₀
+    L = heat.prop.L
     cw = heat.prop.cw # heat capacity of liquid water
     ci = heat.prop.ci # heat capacity of ice
     α₀ = solver.α₀
     τ = solver.τ
     # compute initial residual
-    Tres, θl, C = residual(soil, heat, T, H, L, f, f_args, θw, θm, θo)
+    Tres, θw, C = sfccresidual(soil, heat, f, f_args, f_hc, T, H)
     itercount = 0
     T_converged = false
     while abs(Tres) > solver.abstol && abs(Tres) / abs(T) > solver.reltol
         if itercount >= solver.maxiter
-            return (;T, Tres, θl, itercount, T_converged)
+            return (;T, Tres, θw, itercount, T_converged)
         end
         # derivative of freeze curve
         ∂θ∂T = ForwardDiff.derivative(T -> f(T, f_args...), T)
@@ -47,22 +51,22 @@ function sfccsolve(solver::SFCCNewtonSolver, soil::Soil, heat::Heat, f::F, f_arg
         # note that this assumes heatcapacity to be a simple weighted average!
         # in the future, it might be a good idea to compute an automatic derivative
         # of heatcapacity in addition to the freeze curve function.
-        ∂Tres∂T = 1.0 - ∂θ∂T*(-L*C - (H - θl*L)*(cw-ci))/C^2
+        ∂Tres∂T = 1.0 - ∂θ∂T*(-L*C - (H - θw*L)*(cw-ci))/C^2
         α = α₀ / ∂Tres∂T
         T̂ = T - α*Tres
         # do first residual check outside of loop;
         # this way, we don't decrease α unless we have to.
-        T̂res, θl, C = residual(soil, heat, T̂, H, L, f, f_args, θw, θm, θo)
+        T̂res, θw, C = sfccresidual(soil, heat, f, f_args, f_hc, T̂, H)
         inneritercount = 0
         # simple backtracking line search to avoid jumping over the solution
         while sign(T̂res) != sign(Tres)
             if inneritercount > 100
                 @warn "Backtracking failed; this should not happen. Current state: α=$α, T=$T, T̂=$T̂, residual $(T̂res), initial residual: $(Tres)"
-                return (;T, Tres, θl, itercount, T_converged)
+                return (;T, Tres, θw, itercount, T_converged)
             end
             α = α*τ # decrease step size by τ
             T̂ = T - α*Tres # new guess for T
-            T̂res, θl, C = residual(soil, heat, T̂, H, L, f, f_args, θw, θm, θo)
+            T̂res, θw, C = sfccresidual(soil, heat, f, f_args, f_hc, T̂, H)
             inneritercount += 1
         end
         T = T̂ # update T
@@ -70,9 +74,9 @@ function sfccsolve(solver::SFCCNewtonSolver, soil::Soil, heat::Heat, f::F, f_arg
         itercount += 1
     end
     T_converged = true
-    return (;T, Tres, θl, itercount, T_converged)
+    return (;T, Tres, θw, itercount, T_converged)
 end
-function enthalpyinv(soil::Soil, heat::Heat{<:SFCC{F,SFCCNewtonSolver},Enthalpy}, state, i) where {F}
+function HeatConduction.enthalpyinv(soil::Soil, heat::Heat{<:SFCC{F,SFCCNewtonSolver},Enthalpy}, state, i) where {F}
     sfcc = freezecurve(heat)
     f = sfcc.f
     # get f arguments; note that this does create some redundancy in the arguments
@@ -82,12 +86,9 @@ function enthalpyinv(soil::Soil, heat::Heat{<:SFCC{F,SFCCNewtonSolver},Enthalpy}
     solver = sfcc.solver
     @inbounds let T₀ = i > 1 ? state.T[i-1] : nothing,
         H = state.H[i] |> Utils.adstrip, # enthalpy
-        L = heat.L, # specific latent heat of fusion
-        θw = totalwater(soil, state, i) |> Utils.adstrip, # total water content
-        θm = mineral(soil, state, i) |> Utils.adstrip, # mineral content
-        θo = organic(soil, state, i) |> Utils.adstrip, # organic content
+        f_hc = partial(heatcapacity, liquidwater, soil, heat, state, i),
         f_argsᵢ = Utils.selectat(i, Utils.adstrip, f_args);
-        T, _, _, _, _ = sfccsolve(solver, soil, heat, f, f_argsᵢ, H, L, θw, θm, θo, T₀)
+        T, _, _, _, _ = sfccsolve(solver, soil, heat, f, f_argsᵢ, f_hc, H, T₀)
         return T
     end
 end
@@ -101,31 +102,31 @@ produce incorrect results otherwise.
 @flattenable struct SFCCPreSolver{TCache} <: SFCCSolver
     cache::TCache | false
     Tmin::Float64 | false
-    dH::Float64 | false
-    SFCCPreSolver(cache, Tmin, dH) = new{typeof(cache)}(cache, Tmin, dH)
+    errtol::Float64 | false
+    SFCCPreSolver(cache, Tmin, errtol) = new{typeof(cache)}(cache, Tmin, errtol)
     """
-        SFCCPreSolver(Tmin=-50.0, dH=2e5)
+        SFCCPreSolver(Tmin=-60.0, errtol=1e-4)
 
     Constructs a new `SFCCPreSolver` with minimum temperature `Tmin` and integration step `dH`.
     Enthalpy values below `H(Tmin)` under the given freeze curve will be extrapolated with a
-    constant/flat function. `dH` determines the step size used when integrating `dθdH`; smaller
-    values will produce a more accurate interpolant at the cost of storing more knots and slower
-    initialization. The default value of `dH=2e5` should be sufficient for most use-cases.
+    constant/flat function. `errtol` determines the permitted local error in the interpolant.
     """
-    function SFCCPreSolver(;Tmin=-50.0, dH=2e5)
+    function SFCCPreSolver(;Tmin=-60.0, errtol=1e-4)
         cache = SFCCPreSolverCache()
-        new{typeof(cache)}(cache, Tmin, dH)
+        new{typeof(cache)}(cache, Tmin, errtol)
     end
 end
-mutable struct SFCCPreSolverCache{F,∇F}
-    f::F # H⁻¹ interpolant
-    ∇f::∇F # derivative of f
+mutable struct SFCCPreSolverCache{TFH,T∇FH,T∇FT}
+    f_H::TFH # θw(H) interpolant
+    ∇f_H::T∇FH # derivative of θw(H) interpolant
+    ∇f_T::T∇FT # ∂θ∂T interpolant
     function SFCCPreSolverCache()
-        # initialize with dummy functions to get type information
+        # initialize with dummy functions to get type information;
+        # this is just to make sure that the compiler can still infer all of the types.
         x  = -3e8:1e6:3e8
         dummy_f = _build_interpolant(x, zeros(length(x)))
         dummy_∇f = first ∘ _interpgrad(dummy_f)
-        return new{typeof(dummy_f),typeof(dummy_∇f)}(dummy_f, dummy_∇f)
+        return new{typeof(dummy_f),typeof(dummy_∇f),typeof(dummy_f)}(dummy_f, dummy_∇f, dummy_f)
     end
 end
 _interpgrad(f) = (args...) -> Interpolations.gradient(f, args...)
@@ -135,12 +136,12 @@ function _build_interpolant(Hs, θs)
         Interpolations.Flat()
     )
 end
-function initialcondition!(soil::Soil{<:HomogeneousCharacteristicFractions}, heat::Heat, sfcc::SFCC{F,<:SFCCPreSolver}, state) where {F}
-    L = heat.L
+function CryoGrid.initialcondition!(soil::Soil{<:HomogeneousCharacteristicFractions}, heat::Heat, sfcc::SFCC{F,<:SFCCPreSolver}, state) where {F}
+    L = heat.prop.L
     args = sfccargs(sfcc.f, soil, heat, state)
-    state.θl .= sfcc.f.(state.T, args...)
+    state.θw .= sfcc.f.(state.T, args...)
     heatcapacity!(soil, heat, state)
-    @. state.H = enthalpy(state.T, state.C, L, state.θl)
+    @. state.H = enthalpy(state.T, state.C, L, state.θw)
     # pre-solve freeze curve;
     # note that this is only valid given that the following assumptions hold:
     # 1) none of the freeze curve parameters (e.g. soil properties) change
@@ -150,29 +151,65 @@ function initialcondition!(soil::Soil{<:HomogeneousCharacteristicFractions}, hea
         tail_args = args[3:end],
         θm = mineral(soil, state),
         θo = organic(soil, state),
-        L = heat.L,
+        θa = 1-θm-θo-θtot,
+        L = heat.prop.L,
         Tmin = sfcc.solver.Tmin,
         Tmax = 0.0,
-        θ(T) = sfcc.f(T, θsat, θtot, tail_args...),
-        C(θl) = heatcapacity(soil, heat, θtot, θl, θm, θo),
-        Hmin = enthalpy(Tmin, C(θ(Tmin)), L, θ(Tmin)),
-        Hmax = enthalpy(Tmax, C(θ(Tmax)), L, θ(Tmax)),
-        dH = sfcc.solver.dH,
-        Hs = Hmin:dH:Hmax;
-        θs = Vector{eltype(state.θl)}(undef, length(Hs))
-        θs[1] = θ(Tmin)
-        Ts = Vector{eltype(state.T)}(undef, length(Hs))
-        Ts[1] = Tmin
-        solver = SFCCNewtonSolver(abstol=1e-3, reltol=1e-3, maxiter=100)
-        for i in 2:length(Hs)
-            Hᵢ = Hs[i]
-            T₀ = Ts[i-1] # use previous temperature value as initial guess
-            res = sfccsolve(solver, soil, heat, sfcc.f, args, Hᵢ, L, θtot, θm, θo, T₀)
-            @assert res.T_converged "solver failed to converge after $(res.itercount) iterations: H=$Hᵢ, T₀=$T₀, T=$(res.T), Tres=$(res.Tres), θl=$(res.θl)"
-            θs[i] = res.θl
-            Ts[i] = res.T
+        f(T) = sfcc.f(T, θsat, θtot, tail_args...),
+        C(θw) = heatcapacity(soil, heat, θw, θtot-θw, θa, θm, θo),
+        Hmin = enthalpy(Tmin, C(f(Tmin)), L, f(Tmin)),
+        Hmax = enthalpy(Tmax, C(f(Tmax)), L, f(Tmax));
+        # residual as a function of T and H
+        resid(T,H) = sfccresidual(soil, heat, sfcc.f, args, C, T, H)
+        function solve(H,T₀)
+            local T = nlsolve(T -> resid(T[1],H)[1], [T₀]).zero[1]
+            θw = f(T)
+            return (; T, θw)
+        end
+        function deriv(T) # implicit partial derivative w.r.t H as a function of T
+            θw, ∂θ∂T = ∇(f, T)
+            # get C_eff, i.e. dHdT
+            ∂H∂T = HeatConduction.C_eff(T, C(θw), heat.prop.L, ∂θ∂T, heat.prop.cw, heat.prop.ci)
+            # valid by chain rule and inverse function theorem
+            return ∂θ∂T / ∂H∂T, ∂θ∂T
+        end
+        function step(ΔH, H, θ, ∂θ∂H, T₀)
+            # get first order estimate
+            θest = θ + ΔH*∂θ∂H
+            # get true θ at H + ΔH
+            θsol = solve(H + ΔH, T₀).θw
+            err = abs(θsol - θest)
+            # return residual of error with target error
+            return err
+        end
+        T = [Tmin]
+        H = [Hmin]
+        θ = [f(T[1])]
+        dθdT₀, dθdH₀ = deriv(T[1])
+        dθdT = [dθdT₀]
+        dθdH = [dθdH₀]
+        @assert isfinite(H[1]) && isfinite(θ[1]) "H=$H, θ=$θ"
+        while H[end] < Hmax
+            # find the optimal step size
+            ϵ = Inf
+            ΔH = heat.prop.L*θtot/10 # initially set to large value
+            while abs(ϵ) > sfcc.solver.errtol
+                ϵ = step(ΔH, H[end], θ[end], dθdH[end], T[end])
+                # iteratively halve the step size until error tolerance is satisfied
+                ΔH *= 0.5
+            end
+            Hnew = H[end] + ΔH
+            @assert isfinite(Hnew) "isfinite(ΔH) failed; H=$(H[end]), T=$(T[end]), ΔH=$ΔH"
+            opt = solve(Hnew, T[end])
+            dθdHᵢ, dθdTᵢ = deriv(opt.T)
+            push!(H, Hnew)
+            push!(θ, opt.θw)
+            push!(T, opt.T)
+            push!(dθdT, dθdTᵢ)
+            push!(dθdH, dθdHᵢ)
         end
-        sfcc.solver.cache.f = _build_interpolant(Hs, θs)
-        sfcc.solver.cache.∇f = first ∘ _interpgrad(sfcc.solver.cache.f)
+        sfcc.solver.cache.f_H = _build_interpolant(H, θ)
+        sfcc.solver.cache.∇f_H = first ∘ _interpgrad(sfcc.solver.cache.f_H)
+        sfcc.solver.cache.∇f_T = _build_interpolant(T, dθdT)
     end
 end

src/Physics/WaterBalance/swrc.jl → src/Physics/Soils/sfcc/swrc.jl

@@ -4,6 +4,7 @@
 Base type for soil water retention curve SWRC function implementations.
 """
 abstract type SWRCFunction end
+Base.inv(f::SWRCFunction) = (args...) -> f(inv, args...)
 """
     SWRC{F}
 
@@ -19,9 +20,17 @@ end
 van Genuchten MT, 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils.
     Soil Science Society of America Journal, 44(5): 892–898. DOI: 10.2136/sssaj 1980.03615995004400050002x.
 """
-struct VanGenuchten <: SWRCFunction end
-function (f::VanGenuchten)(ψ,θres,θsat,α,n)
+Base.@kwdef struct VanGenuchten{Tα,Tn} <: SWRCFunction
+    α::Tα = Param(1.0, units=u"1/m")
+    n::Tn = Param(2.0)
+end
+function (f::VanGenuchten)(ψ, θres, θsat, α=f.α, n=f.n)
     let m = 1-1/n;
         IfElse.ifelse(ψ <= zero(ψ), θres + (θsat - θres)*(1 + abs(-α*ψ)^n)^(-m), θsat)
     end
 end
+function (f::VanGenuchten)(::typeof(inv), θ, θres, θsat, α=f.α, n=f.n)
+    let m = 1-1/n;
+        IfElse.ifelse(θ < θsat, -1/α*(((θ-θres)/(θsat-θres))^(-1/m)-1.0)^(1/n), zero(1/α))
+    end
+end

src/Physics/Soils/soilheat.jl

@@ -4,115 +4,63 @@
 # In the latter case, specializations should override only the index-free form
 # and return a state vector instead of a scalar. The `getscalar` function will
 # handle both the scalar and vector case!
-"""
-    mineral(soil::Soil, state, i)
-
-Retrieves the mineral content for the given layer at grid cell `i`, if provided.
-Defaults to using the scalar mineral content defined on `soil`.
-"""
-@inline mineral(soil::Soil, state, i) = Utils.getscalar(mineral(soil, state), i)
-"""
-    organic(soil::Soil, state, i)
-
-Retrieves the organic content for the given layer at grid cell `i`, if provided.
-Defaults to using the scalar organic content defined on `soil`.
-"""
-@inline organic(soil::Soil, state, i) = Utils.getscalar(organic(soil, state), i)
-"""
-    porosity(soil::Soil, state, i)
-
-Retrieves the porosity for the given layer at grid cell `i`, if provided.
-Defaults to using the scalar porosity defined on `soil`.
-"""
-@inline porosity(soil::Soil, state, i) = Utils.getscalar(porosity(soil, state), i)
 # Functions for retrieving constituents and volumetric fractions
-@inline thermalconductivities(soil::Soil, heat::Heat) = (heat.prop.kw, heat.prop.ki, heat.prop.ka, soil.prop.km, soil.prop.ko)
-@inline heatcapacities(soil::Soil, heat::Heat) = (heat.prop.cw, heat.prop.ci, heat.prop.ca, soil.prop.cm, soil.prop.co)
-@inline function volumetricfractions(soil::Soil, heat::Heat, state, i)
-    return let θw = totalwater(soil, state, i),
-        θl = liquidwater(soil, heat, state, i),
+@inline HeatConduction.thermalconductivities(soil::Soil, heat::Heat) = (heat.prop.kw, heat.prop.ki, heat.prop.ka, soil.prop.km, soil.prop.ko)
+@inline HeatConduction.heatcapacities(soil::Soil, heat::Heat) = (heat.prop.cw, heat.prop.ci, heat.prop.ca, soil.prop.cm, soil.prop.co)
+@inline function Physics.volumetricfractions(soil::Soil, heat::Heat, state, i)
+    return let θwi = waterice(soil, heat, state, i),
+        θw = liquidwater(soil, heat, state, i),
         θm = mineral(soil, state, i),
         θo = organic(soil, state, i),
-        θa = 1.0 - θw - θm - θo,
-        θi = θw - θl;
-        (θl, θi, θa, θm, θo)
-    end
-end
-@inline function heatcapacity(soil::Soil, heat::Heat, θw, θl, θm, θo)
-    let θa = 1.0 - θw - θm - θo,
-        θi = θw - θl;
-        return heatcapacity(heatcapacities(soil, heat), (θl, θi, θa, θm, θo))
-    end
-end
-@inline function thermalconductivity(soil::Soil, heat::Heat, θw, θl, θm, θo)
-    let θa = 1.0 - θw - θm - θo,
-        θi = θw - θl;
-        return thermalconductivity(thermalconductivities(soil, heat), (θl, θi, θa, θm, θo))
+        θa = 1.0 - θwi - θm - θo,
+        θi = θwi - θw;
+        (θw, θi, θa, θm, θo)
     end
 end
 
-# SFCC
-include("sfcc.jl")
-
 """
 Initial condition for heat conduction (all state configurations) on soil layer w/ SFCC.
 """
-function initialcondition!(soil::Soil, heat::Heat{<:SFCC}, state)
-    initialcondition!(soil, heat, freezecurve(heat), state)
-end
+CryoGrid.initialcondition!(soil::Soil, heat::Heat{<:SFCC}, state) = HeatConduction.initialcondition!(soil, heat, freezecurve(heat), state)
 """
 Initial condition for heat conduction (all state configurations) on soil layer w/ free water freeze curve.
 """
-function initialcondition!(soil::Soil, heat::Heat{FreeWater}, state)
-    L = heat.L
+function CryoGrid.initialcondition!(soil::Soil, heat::Heat{FreeWater,Enthalpy}, state)
+    L = heat.prop.L
     # initialize liquid water content based on temperature
     @inbounds for i in 1:length(state.T)
-        θw = totalwater(soil, state, i)
-        state.θl[i] = ifelse(state.T[i] > 0.0, θw, 0.0)
-        state.C[i] = heatcapacity(soil, heat, state, i)
-        state.H[i] = enthalpy(state.T[i], state.C[i], L, state.θl[i])
+        θwi = waterice(soil, heat, state, i)
+        state.θw[i] = ifelse(state.T[i] > 0.0, θwi, 0.0)
+        state.C[i] = heatcapacity(soil, heat, volumetricfractions(soil, heat, state, i)...)
+        state.H[i] = enthalpy(state.T[i], state.C[i], L, state.θw[i])
     end
 end
-function liquidwater(::SubSurface, heat::Heat{<:SFCC,Temperature}, state, i)
-    sfcc = freezecurve(heat)
-    f_args = tuplejoin((state.T,), sfccargs(sfcc.f, soil, heat, state))
-    f_argsᵢ = Utils.selectat(i, identity, f_args)
-    return sfcc.f(f_argsᵢ...)
-end
-function liquidwater(sub::SubSurface, heat::Heat{<:SFCC,Enthalpy}, state, i)
-    T = enthalpyinv(sub, heat, state, i)
-    sfcc = freezecurve(heat)
-    f_args = sfccargs(sfcc.f, soil, heat, state)
-    f_argsᵢ = Utils.selectat(i, identity, f_args)
-    return sfcc.f(T, f_argsᵢ...)
-end
 """
     freezethaw!(soil::Soil, heat::Heat{<:SFCC,Temperature}, state)
     freezethaw!(soil::Soil, heat::Heat{<:SFCC,Enthalpy}, state)
 
 Updates state variables according to the specified SFCC function and solver.
 For heat conduction with enthalpy, evaluation of the inverse enthalpy function is performed using the given solver.
-For heat conduction with temperature, we can simply evaluate the freeze curve to get C_eff, θl, and H.
+For heat conduction with temperature, we can simply evaluate the freeze curve to get C_eff, θw, and H.
 """
-function freezethaw!(soil::Soil, heat::Heat{<:SFCC,Temperature}, state)
+function HeatConduction.freezethaw!(soil::Soil, heat::Heat{<:SFCC,Temperature}, state)
     sfcc = freezecurve(heat)
-    @inbounds @fastmath let L = heat.L,
+    let L = heat.prop.L,
         f = sfcc.f,
-        f_args = sfccargs(f,soil,heat,state),
-        f_res = ForwardDiff.DiffResult(zero(eltype(state.T)), zero(eltype(state.T)));
-        for i in 1:length(state.T)
+        f_args = sfccargs(f,soil,heat,state);
+        @inbounds @fastmath for i in 1:length(state.T)
             T = state.T[i]
             f_argsᵢ = Utils.selectat(i, identity, f_args)
-            θl, dθdT = ∇(T -> f(T, f_args...), T)
-            state.θl[i] = θl
+            θw, dθdT = ∇(T -> f(T, f_argsᵢ...), T)
+            state.θw[i] = θw
             state.dθdT[i] = dθdT
-            state.C[i] = C = heatcapacity(soil, heat, state, i)
-            state.dHdT[i] = C + (L + T*(heat.prop.cw - heat.prop.ci))*dθdT
-            state.H[i] = enthalpy(T, C, L, θl)
+            state.C[i] = C = heatcapacity(soil, heat, volumetricfractions(soil, heat, state, i)...)
+            state.dHdT[i] = HeatConduction.C_eff(T, C, L, dθdT, heat.prop.cw, heat.prop.ci)
+            state.H[i] = enthalpy(T, C, L, θw)
         end
     end
 end
-function freezethaw!(soil::Soil, heat::Heat{<:SFCC{F,SFCCNewtonSolver},Enthalpy}, state) where {F}
+function HeatConduction.freezethaw!(soil::Soil, heat::Heat{<:SFCC{F,SFCCNewtonSolver},Enthalpy}, state) where {F}
     sfcc = freezecurve(heat)
     f_args = sfccargs(sfcc.f, soil, heat, state)
     @inbounds for i in 1:length(state.H)
@@ -122,36 +70,32 @@ function freezethaw!(soil::Soil, heat::Heat{<:SFCC{F,SFCCNewtonSolver},Enthalpy}
             T = enthalpyinv(soil, heat, state, i)
             # Here we apply the recovered temperature to the state variables;
             # Since we perform iteration on untracked variables, we need to
-            # recompute θl, C, and T here with the tracked variables.
+            # recompute θw, C, and T here with the tracked variables.
             # Note that this results in one additional freeze curve function evaluation.
-            state.θl[i] = f(T, f_argsᵢ...)
-            dθdT = ForwardDiff.derivative(T -> f(T, f_argsᵢ...), T)
-            let θl = state.θl[i],
-                H = state.H[i],
-                L = heat.L,
+            state.θw[i], dθdT = ∇(T -> f(T, f_argsᵢ...), T)
+            let H = state.H[i],
+                L = heat.prop.L,
                 cw = heat.prop.cw,
                 ci = heat.prop.ci,
-                θw = totalwater(soil, state, i), # total water content
-                θm = mineral(soil, state, i), # mineral content
-                θo = organic(soil, state, i), # organic content
-                θw = totalwater(soil, state, i);
-                state.C[i] = heatcapacity(soil, heat, θw, θl, θm, θo)
-                state.T[i] = (H - L*θl) / state.C[i]
+                θw = state.θw[i],
+                (_, θi, θa, θm, θo) = volumetricfractions(soil, heat, state, i);
+                state.C[i] = heatcapacity(soil, heat, θw, θi, θa, θm, θo)
+                state.T[i] = (H - L*θw) / state.C[i]
                 state.dHdT[i] = HeatConduction.C_eff(state.T[i], state.C[i], L, dθdT, cw, ci)
             end
         end
     end
 end
-function freezethaw!(soil::Soil{<:HomogeneousCharacteristicFractions}, heat::Heat{<:SFCC{F,<:SFCCPreSolver},Enthalpy}, state) where {F}
+function HeatConduction.freezethaw!(soil::Soil{<:HomogeneousCharacteristicFractions}, heat::Heat{<:SFCC{F,<:SFCCPreSolver},Enthalpy}, state) where {F}
     solver = freezecurve(heat).solver
-    f = solver.cache.f
-    ∇f = solver.cache.∇f
+    f_H = solver.cache.f_H
+    ∇f_T = solver.cache.∇f_T
+    L, cw, ci = heat.prop.L, heat.prop.cw, heat.prop.ci
     @inbounds for i in 1:length(state.H)
         H = state.H[i]
-        state.θl[i] = θl = f(H)
-        state.C[i] = C = heatcapacity(soil, heat, state, i)
-        state.T[i] = T = enthalpyinv(H, C, heat.L, θl)
-        dθdTᵢ = ∇f(state.H[i])
-        state.dHdT[i] = C + dθdTᵢ*(heat.L + T*(heat.prop.cw - heat.prop.ci))
+        state.θw[i] = θw = f_H(H)
+        state.C[i] = C = heatcapacity(soil, heat, volumetricfractions(soil, heat, state, i)...)
+        state.T[i] = T = enthalpyinv(H, C, heat.prop.L, θw)
+        state.dHdT[i] = HeatConduction.C_eff(T, C, L, ∇f_T(T), cw, ci)
     end
 end

src/Physics/WaterBalance/WaterBalance.jl

-module WaterBalance
-
-import CryoGrid: BoundaryStyle, diagnosticstep!, prognosticstep!, interact!, initialcondition!, variables
-
-using ..Physics
-using CryoGrid.Numerics
-using CryoGrid.Numerics: nonlineardiffusion!
-
-using IfElse
-
-export SWRC, VanGenuchten
-
-include("swrc.jl")
-
-# TODO: implement water fluxes
-
-end

src/Physics/common.jl

+Constants() = (
+    ρw = 1000.0u"kg/m^3", # density of water at standard conditions
+    Lf = 3.34e5u"J/kg", # specific latent heat of fusion of water [J/kg]
+    g = 9.80665u"m/s^2", # gravitational constant
+)
+# Composition
+"""
+    volumetricfractions(::SubSurface, ::SubSurfaceProcess, state)
+    volumetricfractions(::SubSurface, ::SubSurfaceProcess, state, i)
+
+Get the volumetric fractions of each constituent in the volume (at grid cell `i`, if specificed).
+All implementations of `volumetricfractions` are expected to obey a semi-consistent order
+in the returned `Tuple` of fractions; the first three consituents should always be `θw,θi,θa`,
+i.e. water, ice, and air, followed by any number of additional constituents which may be defined
+by the specific layer. There is no feasible way to verify that client code actually obeys this
+ordering, so be sure to double check your implementation, otherwise this can cause very subtle bugs!
+"""
+@inline volumetricfractions(::SubSurface, ::SubSurfaceProcess, state) = ()
+@inline volumetricfractions(sub::SubSurface, proc::SubSurfaceProcess, state, i) = volumetricfractions(sub, proc, state)
+"""
+    waterice(::SubSurface, state)
+    waterice(sub::SubSurface, ::SubSurfaceProcess, state)
+    waterice(sub::SubSurface, ::SubSurfaceProcess, state, i)
+
+Retrieves the total water content (water + ice) for the given layer at grid cell `i`, if specified.
+Defaults to retrieving the state variable `θwi` (assuming it exists).
+"""
+@inline waterice(::SubSurface, state) = state.θwi
+@inline waterice(sub::SubSurface, proc::SubSurfaceProcess, state) = waterice(sub, state)
+@inline waterice(sub::SubSurface, proc::SubSurfaceProcess, state, i) = Utils.getscalar(waterice(sub, proc, state), i)
+"""
+    liquidwater(sub::SubSurface, proc::SubSurfaceProcess, state)
+    liquidwater(sub::SubSurface, proc::SubSurfaceProcess, state, i)
+
+Retrieves the liquid water content for the given layer.
+"""
+@inline liquidwater(::SubSurface, ::SubSurfaceProcess, state) = state.θw
+@inline liquidwater(sub::SubSurface, proc::SubSurfaceProcess, state, i) = Utils.getscalar(liquidwater(sub, proc, state), i)
+"""
+    partial(f, ::typeof(liquidwater), sub::SubSurface, proc::SubSurfaceProcess, state, i)
+
+Returns a partially applied function `f` which takes liquid water `θw` as an argument and holds all other
+volumetric fractions constant. `f` must be a function of the form `f(::Layer, ::Process, θfracs...)` where
+`θfracs` are an arbitrary number of constituent volumetric fractions. Note that this method assumes that
+`volumetricfractions` and `f` both obey the implicit ordering convention: `(θw, θi, θa, θfracs...)` where
+`θfracs` are zero or more additional constituent fractions. The returned method is a closure which has the
+following properties available: `θw, θi, θa, θfracs, θwi` where `θwi` refers to the sum of `θw` and `θi`.
+"""
+function partial(f::F, ::typeof(liquidwater), sub::SubSurface, proc::SubSurfaceProcess, state, i) where F
+    (θw, θi, θa, θfracs...) = volumetricfractions(sub, proc, state, i)
+    θwi = θw + θi
+    return function apply(θw)
+        return f(sub, proc, θw, θwi - θw, θa, θfracs...)
+    end
+end

src/Presets/Presets.jl

@@ -31,9 +31,9 @@ end
 SoilHeatColumn(upperbc::BoundaryProcess, soilprofile::Profile, init::Numerics.VarInitializer; grid::Grid=DefaultGrid_2cm, freezecurve::F=FreeWater()) where {F<:FreezeCurve} = SoilHeatColumn(:H, upperbc, soilprofile, init; grid=grid, freezecurve=freezecurve)
 
 Forcings = (
-    Samoylov_ERA5_fitted_daily_1979_2020 = Resource("samoylov_era5_fitted_daily_1979-2020", "json", "https://nextcloud.awi.de/s/DSx3BWsfjCdy9MF"),
-    Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044 = Resource("Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044", "json", "https://nextcloud.awi.de/s/dp74KC2ceQKaG43"),
-    Samoylov_ERA_MkL3_CCSM4_long_term = Resource("Samoylov_ERA_MkL3_CCSM4_long_term", "json", "https://nextcloud.awi.de/s/gyoMTy9jpk2pMxL"),
+    Samoylov_ERA5_fitted_daily_1979_2020 = Resource("samoylov_era5_fitted_daily_1979-2020", "json", "https://nextcloud.awi.de/s/WJtT7CS7HtcoRDz/download/samoylov_era5_fitted_daily_1979-2020.json"),
+    Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044 = Resource("Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044", "json", "https://nextcloud.awi.de/s/dp74KC2ceQKaG43/download/samoylov_ERA_obs_fitted_1979_2014_spinup_extended2044.json"),
+    Samoylov_ERA_MkL3_CCSM4_long_term = Resource("Samoylov_ERA_MkL3_CCSM4_long_term", "json", "https://nextcloud.awi.de/s/gyoMTy9jpk2pMxL/download/FORCING_ULC_126_72.json"),
 )
 Parameters = (
     # Faroux et al. doi:10.1109/IGARSS.2007.4422971
@@ -42,11 +42,11 @@ Parameters = (
 
 const SamoylovDefault = (
     soilprofile = SoilProfile(
-        0.0u"m" => CharacteristicFractions(xic=0.0,por=0.80,sat=1.0,org=0.75), #(θw=0.80,θm=0.05,θo=0.15,ϕ=0.80),
-        0.1u"m" => CharacteristicFractions(xic=0.0,por=0.80,sat=1.0,org=0.25), #(θw=0.80,θm=0.15,θo=0.05,ϕ=0.80),
-        0.4u"m" => CharacteristicFractions(xic=0.30,por=0.55,sat=1.0,org=0.25), #(θw=0.80,θm=0.15,θo=0.05,ϕ=0.55),
-        3.0u"m" => CharacteristicFractions(xic=0.0,por=0.50,sat=1.0,org=0.0), #(θw=0.50,θm=0.50,θo=0.0,ϕ=0.50),
-        10.0u"m" => CharacteristicFractions(xic=0.0,por=0.30,sat=1.0,org=0.0), #(θw=0.30,θm=0.70,θo=0.0,ϕ=0.30),
+        0.0u"m" => soilparameters(xic=0.0,por=0.80,sat=1.0,org=0.75), #(θwi=0.80,θm=0.05,θo=0.15,ϕ=0.80),
+        0.1u"m" => soilparameters(xic=0.0,por=0.80,sat=1.0,org=0.25), #(θwi=0.80,θm=0.15,θo=0.05,ϕ=0.80),
+        0.4u"m" => soilparameters(xic=0.30,por=0.55,sat=1.0,org=0.25), #(θwi=0.80,θm=0.15,θo=0.05,ϕ=0.55),
+        3.0u"m" => soilparameters(xic=0.0,por=0.50,sat=1.0,org=0.0), #(θwi=0.50,θm=0.50,θo=0.0,ϕ=0.50),
+        10.0u"m" => soilparameters(xic=0.0,por=0.30,sat=1.0,org=0.0), #(θwi=0.30,θm=0.70,θo=0.0,ϕ=0.30),
     ),
     tempprofile = TemperatureProfile(
         0.0u"m" => -1.0u"°C",

src/Presets/presetgrids.jl

+const DefaultGrid_1cm = Grid(vcat([
+            0:0.01:0.5...,
+            0.55:0.02:2...,
+            2.05:0.05:4.0...,
+            4.1:0.1:10...,
+            10.2:0.2:20...,
+            21:1:30...,
+            35:5:50...,
+            60:10:100...,
+            200:100:1000...
+        ]...
+    )u"m"
+)
 const DefaultGrid_2cm = Grid(vcat([
             0:0.02:2...,
             2.05:0.05:4.0...,

src/Strat/Strat.jl

 module Strat
 
 import CryoGrid
-import CryoGrid: variables, initialcondition!, prognosticstep!, diagnosticstep!, interact!, observe
 
-using CryoGrid: Layer, Top, Bottom, SubSurface, Process, SubSurfaceProcess, BoundaryProcess, CoupledProcesses
+using CryoGrid
 using CryoGrid: Parameterization, DynamicParameterization
 using CryoGrid.InputOutput: Forcing, CryoGridParams
 using CryoGrid.Numerics
@@ -11,7 +10,6 @@ using CryoGrid.Physics
 using CryoGrid.Utils
 
 using ComponentArrays
-using ConstructionBase
 using DataStructures: OrderedDict
 using Dates
 using DimensionalData
@@ -26,6 +24,7 @@ using Setfield
 @reexport using ModelParameters
 @reexport using SimulationLogs
 
+import ConstructionBase
 import Flatten
 import ModelParameters: update
 

src/Strat/state.jl

@@ -4,75 +4,79 @@ Enumeration for in-place vs out-of-place mode.
 """
 @enum InPlaceMode inplace ooplace
 """
-    LayerState{iip,TGrid,TStates,TGrids,Tt,Tz,varnames}
+    LayerState{iip,TGrid,TStates,TGrids,Tt,Tdt,Tz,varnames}
 
 Represents the state of a single component (layer + processes) in the stratigraphy.
 """
-struct LayerState{iip,TGrid,TStates,TGrids,Tt,Tz,varnames}
+struct LayerState{iip,TGrid,TStates,TGrids,Tt,Tdt,Tz,varnames}
     grid::TGrid
     grids::NamedTuple{varnames,TGrids}
     states::NamedTuple{varnames,TStates}
     bounds::NTuple{2,Tz}
     t::Tt
-    LayerState(grid::TG, grids::NamedTuple{varnames,Tvg}, states::NamedTuple{varnames,Tvs}, t::Tt, bounds::NTuple{2,Tz}, ::Val{iip}=Val{inplace}()) where
-        {TG,Tvg,Tvs,Tt,Tz,varnames,iip} = new{iip,TG,Tvs,Tvg,Tt,Tz,varnames}(grid, grids, states, bounds, t)
+    dt::Tdt
+    LayerState(grid::TG, grids::NamedTuple{varnames,Tvg}, states::NamedTuple{varnames,Tvs}, bounds::NTuple{2,Tz}, t::Tt, dt::Tdt=nothing, ::Val{iip}=Val{inplace}()) where
+        {TG,Tvg,Tvs,Tt,Tdt,Tz,varnames,iip} = new{iip,TG,Tvs,Tvg,Tt,Tdt,Tz,varnames}(grid, grids, states, bounds, t, dt)
 end
 Base.getindex(state::LayerState, sym::Symbol) = getproperty(state, sym)
 function Base.getproperty(state::LayerState, sym::Symbol)
-    return if sym ∈ (:grid, :grids, :states, :t, :z)
+    return if sym ∈ (:grid, :grids, :states, :bounds, :t, :dt, :z)
         getfield(state, sym)
     else
         getproperty(getfield(state, :states), sym)
     end
 end
-@inline function LayerState(vs::VarStates, zs::NTuple{2}, u, du, t, ::Val{layername}, ::Val{iip}=Val{inplace}()) where {layername,iip}
+@inline function LayerState(vs::VarStates, zs::NTuple{2}, u, du, t, dt, ::Val{layername}, ::Val{iip}=Val{inplace}()) where {layername,iip}
     z_inds = subgridinds(edges(vs.grid), zs[1]..zs[2])
     return LayerState(
         vs.grid[z_inds],
         _makegrids(Val{layername}(), getproperty(vs.vars, layername), vs, z_inds),
         _makestates(Val{layername}(), getproperty(vs.vars, layername), vs, z_inds, u, du, t),
-        t,
         zs,
+        t,
+        dt,
         Val{iip}(),
     )
 end
 
 """
-    TileState{iip,TGrid,TStates,Tt,names}
+    TileState{iip,TGrid,TStates,Tt,Tdt,names}
 
 Represents the instantaneous state of a CryoGrid `Tile`.
 """
-struct TileState{iip,TGrid,TStates,Tt,names}
+struct TileState{iip,TGrid,TStates,Tt,Tdt,names}
     grid::TGrid
     states::NamedTuple{names,TStates}
     t::Tt
-    TileState(grid::TGrid, states::NamedTuple{names,TS}, t::Tt, ::Val{iip}=Val{inplace}()) where
-        {TGrid<:Numerics.AbstractDiscretization,TS<:Tuple{Vararg{<:LayerState}},Tt,names,iip} =
-            new{iip,TGrid,TS,Tt,names}(grid, states, t)
+    dt::Tdt
+    TileState(grid::TGrid, states::NamedTuple{names,TS}, t::Tt, dt::Tdt=nothing, ::Val{iip}=Val{inplace}()) where
+        {TGrid<:Numerics.AbstractDiscretization,TS<:Tuple{Vararg{<:LayerState}},Tt,Tdt,names,iip} =
+            new{iip,TGrid,TS,Tt,Tdt,names}(grid, states, t, dt)
 end
 Base.getindex(state::TileState, sym::Symbol) = Base.getproperty(state, sym)
 Base.getindex(state::TileState, i::Int) = state.states[i]
 function Base.getproperty(state::TileState, sym::Symbol)
-    return if sym ∈ (:grid,:states,:t)
+    return if sym ∈ (:grid,:states,:t,:dt)
         getfield(state, sym)
     else
         getproperty(getfield(state, :states), sym)
     end
 end
-@inline @generated function TileState(vs::VarStates{names}, zs::NTuple, u=copy(vs.uproto), du=similar(vs.uproto), t=0.0, ::Val{iip}=Val{inplace}()) where {names,iip}
+@inline @generated function TileState(vs::VarStates{names}, zs::NTuple, u=copy(vs.uproto), du=similar(vs.uproto), t=0.0, dt=nothing, ::Val{iip}=Val{inplace}()) where {names,iip}
     layerstates = (
         quote
-            bounds_i = (ustrip(bounds[$i][1]), ustrip(bounds[$i][2]))
-            LayerState(vs, bounds_i, u, du, t, Val{$(QuoteNode(names[i]))}(), Val{iip}())
+            bounds_i = (bounds[$i][1], bounds[$i][2])
+            LayerState(vs, bounds_i, u, du, t, dt, Val{$(QuoteNode(names[i]))}(), Val{iip}())
         end
         for i in 1:length(names)
     )
     quote
-        bounds = boundarypairs(zs, convert(eltype(zs), vs.grid[end]))
+        bounds = boundarypairs(map(ustrip, zs), vs.grid[end])
         return TileState(
             vs.grid,
             NamedTuple{tuple($(map(QuoteNode,names)...))}(tuple($(layerstates...))),
             t,
+            dt,
             Val{iip}(),
         )
     end

src/Strat/stratigraphy.jl

 const RESERVED_COMPONENT_NAMES = (:top, :bottom, :strat, :init)
 
 """
-    StratComponent{TLayer,TProcess,name}
+    StratComponent{TLayer,TProc,name}
 
-Represents a single component (layer + processes) in the stratigraphy.
+Represents a single component (layer + process(es)) in the stratigraphy.
 """
-struct StratComponent{TLayer,TProcesses,name}
+struct StratComponent{TLayer,TProc,name}
     layer::TLayer
-    processes::TProcesses
-    StratComponent(name::Symbol, layer::TLayer, processes::TProcesses) where {TLayer<:Layer,TProcesses<:CoupledProcesses} = new{TLayer,TProcesses,name}(layer,processes)
+    process::TProc
+    StratComponent(name::Symbol, layer::TLayer, proc::TProc) where {TLayer<:Layer,TProc<:Process} = new{TLayer,TProc,name}(layer, proc)
+    function StratComponent(name::Symbol, layer::TLayer, procs::CoupledProcesses; ignore_order=false) where {TLayer<:Layer}
+        # check coupling order
+        if !issorted(procs.processes) && !ignore_order
+            procs = CoupledProcesses(sort(procs.processes))
+            @warn "The ordering of the given coupled processes is inconsistent with the defined rules and has been automatically corrected: $(map(p -> typeof(p).name.wrapper, procs.processes)).
+            If this was on purpose, you can override the defined ordering and suppress this warning with `ignore_order=true` or define `isless` on your process types to explicitly declare the intended ordering."
+        end
+        new{TLayer,typeof(procs),name}(layer,procs)
+    end
 end
-ConstructionBase.constructorof(::Type{StratComponent{TLayer,TProcesses,name}}) where {TLayer,TProcesses,name} = (layer,processes) -> StratComponent(name, layer, processes)
+ConstructionBase.constructorof(::Type{StratComponent{TLayer,TProc,name}}) where {TLayer,TProc,name} = (layer,process) -> StratComponent(name, layer, process)
 """
 Get the name of the given stratigraphy node.
 """
@@ -21,11 +30,17 @@ componentnameval(::StratComponent{L,P,name}) where {L,P,name} = Val{name}
 Base.show(io::IO, node::StratComponent{L,P,name}) where {L,P,name} = print(io, "$name($L,$P)")
 
 # Constructors for stratigraphy nodes
-top(bcs::BoundaryProcess...) = StratComponent(:top, Top(), CoupledProcesses(bcs...))
-bottom(bcs::BoundaryProcess...) = StratComponent(:bottom, Bottom(), CoupledProcesses(bcs...))
-function subsurface(name::Symbol, sub::SubSurface, processes::SubSurfaceProcess...)
+top(bc::BoundaryProcess) = StratComponent(:top, Top(), bc)
+top(bcs::BoundaryProcess...; ignore_order=false) = StratComponent(:top, Top(), CoupledProcesses(bcs...); ignore_order)
+bottom(bc::BoundaryProcess) = StratComponent(:bottom, Bottom(), bc)
+bottom(bcs::BoundaryProcess...; ignore_order=false) = StratComponent(:bottom, Bottom(), CoupledProcesses(bcs...); ignore_order)
+function subsurface(name::Symbol, sub::SubSurface, proc::SubSurfaceProcess)
+    @assert name ∉ RESERVED_COMPONENT_NAMES "layer identifier $name is reserved"
+    return StratComponent(name, sub, proc)
+end
+function subsurface(name::Symbol, sub::SubSurface, procs::SubSurfaceProcess...; ignore_order=false)
     @assert name ∉ RESERVED_COMPONENT_NAMES "layer identifier $name is reserved"
-    return StratComponent(name, sub, CoupledProcesses(processes...))
+    return StratComponent(name, sub, CoupledProcesses(procs...); ignore_order)
 end
 
 """
@@ -57,10 +72,8 @@ struct Stratigraphy{N,TComponents,TBoundaries}
         @assert length(sub) > 0 "At least one subsurface layer must be specified"
         names = map(componentname, map(last, sub))
         @assert length(unique(names)) == length(names) "All layer names in Stratigraphy must be unique"
-        boundary(x) = x
-        boundary(x::DistQuantity) = Param(ustrip(u"m", x), units=u"m", bounds=(0.0,Inf), layer=:strat)
-        boundaries = Tuple(map(boundary ∘ first, (top, sub..., bot)))
-        @assert issorted(filter(Base.Fix2(isa, Param), boundaries), by=p -> p.val) "Stratigraphy boundary locations must be in strictly increasing order."
+        boundaries = Tuple(map(first, (top, sub..., bot)))
+        @assert issorted(boundaries) "Stratigraphy boundary locations must be in strictly increasing order."
         # get components
         components = Tuple(map(last, (top, sub..., bot)))
         # construct type