src/Physics/HeatConduction/HeatConduction.jl

 module HeatConduction
 
-import CryoGrid: SubSurfaceProcess, BoundaryStyle, Dirichlet, Neumann, BoundaryProcess, Layer, Top, Bottom, SubSurface
-import CryoGrid: diagnosticstep!, prognosticstep!, interact!, initialcondition!, boundaryflux, boundaryvalue, variables
+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.InputOutput: Forcing
 using CryoGrid.Physics
 using CryoGrid.Physics.Boundaries
-using CryoGrid.Physics.Water: VanGenuchten
 using CryoGrid.Numerics
 using CryoGrid.Numerics: nonlineardiffusion!, harmonicmean!, harmonicmean, heaviside
-using CryoGrid.Layers: Soil, θp, θw, θm, θo
 using CryoGrid.Utils
 
-using DimensionalData
+using Base: @propagate_inbounds, @kwdef
 using IfElse
-using Interpolations: Linear, Flat
-using IntervalSets
 using ModelParameters
-using Parameters
 using Unitful
 
-import Flatten: @flattenable, flattenable
-
 export Heat, TemperatureProfile
 export FreeWater, FreezeCurve, freezecurve
 
 abstract type FreezeCurve end
 struct FreeWater <: FreezeCurve end
 
-TemperatureProfile(pairs::Pair{<:DistQuantity,<:TempQuantity}...) =
-    Profile(map(p -> uconvert(u"m", p[1]) => uconvert(u"°C", p[2]),pairs)...)
-    
-abstract type HeatVariable end
-struct Enthalpy <: HeatVariable end
-struct PartitionedEnthalpy <: HeatVariable end
-struct Temperature <: HeatVariable end
-
-@with_kw struct Heat{F<:FreezeCurve,S} <: SubSurfaceProcess
-    ρ::Float"kg/m^3" = 1000.0xu"kg/m^3" #[kg/m^3]
-    Lsl::Float"J/kg" = 334000.0xu"J/kg" #[J/kg] (latent heat of fusion)
-    L::Float"J/m^3" = ρ*Lsl             #[J/m^3] (specific latent heat of fusion)
-    freezecurve::F = FreeWater()        # freeze curve, defautls to free water fc
-    sp::S = Enthalpy()
+"""
+    TemperatureProfile(pairs::Pair{<:Union{DistQuantity,Param},<:Union{TempQuantity,Param}}...)
+
+Convenience constructor for `Numerics.Profile` which automatically converts temperature quantities.
+"""
+TemperatureProfile(pairs::Pair{<:Union{DistQuantity,Param},<:Union{TempQuantity,Param}}...) = Profile(map(p -> uconvert(u"m", p[1]) => uconvert(u"°C", p[2]), pairs))
+
+"""
+    HeatImpl
+
+Base type for different numerical formulations of two-phase heat diffusion.
+"""
+abstract type HeatImpl end
+struct Enthalpy <: HeatImpl end
+struct Temperature <: HeatImpl 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
+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
 end
 # convenience constructors for specifying prognostic variable as symbol
-Heat(var::Union{Symbol,Tuple{Vararg{Symbol}}}; kwargs...) = Heat(Val{var}(); kwargs...)
+Heat(var::Symbol; kwargs...) = Heat(Val{var}(); kwargs...)
 Heat(::Val{:H}; kwargs...) = Heat(;sp=Enthalpy(), kwargs...)
-Heat(::Val{(:Hₛ,:Hₗ)}; kwargs...) = Heat(;sp=PartitionedEnthalpy(), kwargs...)
 Heat(::Val{:T}; kwargs...) = Heat(;sp=Temperature(), kwargs...)
+freezecurve(heat::Heat) = heat.freezecurve
+# Default implementation of `variables` for freeze curve
+variables(::SubSurface, ::Heat, ::FreezeCurve) = ()
 
-export enthalpy, heatcapacity, heatcapacity!, thermalconductivity, thermalconductivity!
-export heatconduction!
-include("heat.jl")
 export ConstantTemp, GeothermalHeatFlux, TemperatureGradient, NFactor, Damping
 include("heat_bc.jl")
-export SFCC, DallAmico, Westermann, McKenzie, SFCCNewtonSolver
-include("soil/soilheat.jl")
+
+export heatconduction!, enthalpy, totalwater, liquidwater, freezethaw!, heatcapacity, heatcapacity!, thermalconductivity, thermalconductivity!
+include("heat.jl")
 
 end

src/Physics/HeatConduction/heat.jl

-# Default implementation of `variables` for freeze curve
-variables(::SubSurface, ::Heat, ::FreezeCurve) = ()
-# Fallback (error) implementation for freeze curve
-(fc::FreezeCurve)(layer::SubSurface, heat::Heat, state) =
-    error("freeze curve $(typeof(fc)) not implemented for $(typeof(heat)) on layer $(typeof(layer))")
-
-freezecurve(heat::Heat) = heat.freezecurve
-enthalpy(T::Number"°C", C::Number"J/K/m^3", L::Number"J/m^3", θ::Real) = T*C + L*θ
-totalwater(::SubSurface, ::Heat, state) = state.θw
-heatcapacity!(layer::SubSurface, heat::Heat, state) = error("heatcapacity not defined for $(typeof(heat)) on $(typeof(layer))")
-thermalconductivity!(layer::SubSurface, heat::Heat, state) = error("thermalconductivity not defined for $(typeof(heat)) on $(typeof(layer))")
+# Thermal properties (generic)
+"""
+    enthalpy(T, C, L, θ) = T*C + L*θ
+
+Discrete enthalpy function on temperature, heat capacity, specific latent heat of fusion, and liquid water content.
+"""
+@inline enthalpy(T, C, L, θ) = T*C + L*θ
+"""
+    enthalpyinv(H, C, L, θ) = (H - L*θ) / C
+
+Discrete inverse enthalpy function given H, C, L, and θ.
+"""
+@inline enthalpyinv(H, C, L, θ) = (H - L*θ) / C
+"""
+    C_eff(T, C, L, dθdT, cw, ci) = C + dθdT*(L + T*(cw - ci))
+
+Computes the apparent or "effective" heat capacity `dHdT` as a function of temperature, volumetric heat capacity,
+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)
+
+Get thermal conductivities for generic `SubSurface` layer.
+"""
+@inline thermalconductivities(::SubSurface, heat::Heat) = (heat.prop.kw, heat.prop.ki, heat.prop.ka)
+"""
+    heatcapacities(::SubSurface, heat::Heat)
+
+Get heat capacities for generic `SubSurface` layer.
+"""
+@inline heatcapacities(::SubSurface, heat::Heat) = (heat.prop.cw, heat.prop.ci, heat.prop.ca)
+"""
+    volumetricfractions(::SubSurface, heat::Heat)
+
+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)
+    end
+end
+
+# Generic heat conduction implementation
+
+"""
+    freezethaw!(sub::SubSurface, heat::Heat, state)
+
+Calculates freezing and thawing effects, including evaluation of the freeze curve.
+In general, this function should compute at least the liquid/frozen water contents
+and the corresponding heat capacity. Other variables such as temperature or enthalpy
+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
+
+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)
+    cs = heatcapacities(sub, heat)
+    return heatcapacity(cs, θs)
+end
+"""
+    heatcapacity!(sub::SubSurface, heat::Heat, state)
+
+Computes the heat capacity for the given layer from the current state and stores the result in-place in the state variable `C`.
+"""
+@inline function heatcapacity!(sub::SubSurface, heat::Heat, state)
+    @inbounds for i in 1:length(state.T)
+        state.C[i] = heatcapacity(sub, heat, state, i)
+    end
+end
+"""
+    thermalconductivity(conductivities::NTuple{N}, fracs::NTuple{N}) where N
+
+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)
+    ks = thermalconductivities(sub, heat)
+    return thermalconductivity(ks, θs)
+end
+"""
+    thermalconductivity!(sub::SubSurface, heat::Heat, state)
+
+Computes the thermal conductivity for the given layer from the current state and stores the result in-place in the state variable `C`.
+"""
+@inline function thermalconductivity!(sub::SubSurface, heat::Heat, state)
+    @inbounds for i in 1:length(state.T)
+        state.kc[i] = thermalconductivity(sub, heat, state, i)
+    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) = (
+    variables(sub)..., # layer variables
+    variables(heat)...,  # heat variables
+    variables(sub, heat, freezecurve(heat))..., # freeze curve variables
+)
+"""
+Variable definitions for heat conduction (enthalpy).
+"""
+variables(heat::Heat{<:FreezeCurve,Enthalpy}) = (
+    Prognostic(:H, OnGrid(Cells), u"J/m^3"),
+    Diagnostic(:T, OnGrid(Cells), u"°C"),
+    basevariables(heat)...,
+)
+"""
+Variable definitions for heat conduction (temperature).
+"""
+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)...,
+)
+"""
+Common variable definitions for all heat implementations.
+"""
+basevariables(::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(: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)),
+)
 """
     heatconduction!(∂H,T,ΔT,k,Δk)
 
@@ -21,9 +155,9 @@ function heatconduction!(∂H,T,ΔT,k,Δk)
     @inbounds ∂H[1] += let T₂=T[2],
         T₁=T[1],
         k=k[2],
-        δ=ΔT[1],
-        a=Δk[1];
-        k*(T₂-T₁)/δ/a
+        ϵ=ΔT[1],
+        δ=Δk[1];
+        k*(T₂-T₁)/ϵ/δ
     end
     # diffusion on non-boundary cells
     @inbounds let T = T,
@@ -36,160 +170,110 @@ function heatconduction!(∂H,T,ΔT,k,Δk)
     @inbounds ∂H[end] += let T₂=T[end],
         T₁=T[end-1],
         k=k[end-1],
-        δ=ΔT[end],
-        a=Δk[end];
-        -k*(T₂-T₁)/δ/a
+        ϵ=ΔT[end],
+        δ=Δk[end];
+        -k*(T₂-T₁)/ϵ/δ
     end
     return nothing
 end
-""" Variable definitions for heat conduction (enthalpy) on any subsurface layer. """
-variables(layer::SubSurface, heat::Heat{<:FreezeCurve,Enthalpy}) = (
-    Prognostic(:H, Float"J/m^3", OnGrid(Cells)),
-    Diagnostic(:T, Float"°C", OnGrid(Cells)),
-    Diagnostic(:C, Float"J//K/m^3", OnGrid(Cells)),
-    Diagnostic(:Ceff, Float"J/K/m^3", OnGrid(Cells)),
-    Diagnostic(:k, Float"W/m/K", OnGrid(Edges)),
-    Diagnostic(:kc, Float"W//m/K", OnGrid(Cells)),
-    Diagnostic(:θl, Float"1", OnGrid(Cells)),
-    # add freeze curve variables (if any are present)
-    variables(layer, heat, freezecurve(heat))...,
-)
-""" Variable definitions for heat conduction (partitioned enthalpy) on any subsurface layer. """
-variables(layer::SubSurface, heat::Heat{<:FreezeCurve,PartitionedEnthalpy}) = (
-    Prognostic(:Hₛ, Float"J/m^3", OnGrid(Cells)),
-    Prognostic(:Hₗ, Float"J/m^3", OnGrid(Cells)),
-    Diagnostic(:dH, Float"J/s/m^3", OnGrid(Cells)),
-    Diagnostic(:H, Float"J", OnGrid(Cells)),
-    Diagnostic(:T, Float"°C", OnGrid(Cells)),
-    Diagnostic(:C, Float"J/K/m^3", OnGrid(Cells)),
-    Diagnostic(:Ceff, Float"J/K/m^3", OnGrid(Cells)),
-    Diagnostic(:dθdT, Float"m/m", OnGrid(Cells)),
-    Diagnostic(:k, Float"W/m/K", OnGrid(Edges)),
-    Diagnostic(:kc, Float"W/m/K", OnGrid(Cells)),
-    Diagnostic(:θl, Float"1", OnGrid(Cells)),
-    # add freeze curve variables (if any are present)
-    variables(layer, heat, freezecurve(heat))...,
-)
-""" Variable definitions for heat conduction (temperature) on any subsurface layer. """
-variables(layer::SubSurface, heat::Heat{<:FreezeCurve,Temperature}) = (
-    Prognostic(:T, Float"°C", OnGrid(Cells)),
-    Diagnostic(:H, Float"J/m^3", OnGrid(Cells)),
-    Diagnostic(:dH, Float"J/s/m^3", OnGrid(Cells)),
-    Diagnostic(:C, Float"J/K/m^3", OnGrid(Cells)),
-    Diagnostic(:Ceff, Float"J/K/m^3", OnGrid(Cells)),
-    Diagnostic(:k, Float"W/m/K", OnGrid(Edges)),
-    Diagnostic(:kc, Float"W/m/K", OnGrid(Cells)),
-    Diagnostic(:θl, Float"1", OnGrid(Cells)),
-    # add freeze curve variables (if any are present)
-    variables(layer, heat, freezecurve(heat))...,
-)
-""" Initial condition for heat conduction (all state configurations) on any subsurface layer. """
-function initialcondition!(layer::SubSurface, heat::Heat, state)
-    L = heat.L
-    heatcapacity!(layer, heat, state)
-    @. state.H = enthalpy(state.T, state.C, L, state.θl)
-end
-""" Diagonstic step for heat conduction (all state configurations) on any subsurface layer. """
-function diagnosticstep!(layer::SubSurface, heat::Heat, state)
+"""
+Diagonstic step for heat conduction (all state configurations) on any subsurface layer.
+"""
+function 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))
-    # Evaluate the freeze curve (updates T, C, and θl)
-    fc! = freezecurve(heat);
-    fc!(layer,heat,state)
+    @. state.dH_upper = zero(eltype(state.dH_upper))
+    @. state.dH_lower = zero(eltype(state.dH_lower))
+    # Evaluate freeze/thaw processes
+    freezethaw!(sub, heat, state)
     # Update thermal conductivity
-    thermalconductivity!(layer, heat, state)
-    # Harmonic mean of conductivities
+    thermalconductivity!(sub, heat, state)
+    # thermal conductivity at boundaries
+    # assumes boundary conductivities = cell conductivities
+    @inbounds state.k[1] = state.kc[1]
+    @inbounds state.k[end] = state.kc[end]
+    # Harmonic mean of inner conductivities
     @inbounds let k = (@view state.k[2:end-1]),
         Δk = Δ(state.grids.k);
         harmonicmean!(k, state.kc, Δk)
     end
     return nothing # ensure no allocation
 end
-""" Prognostic step for heat conduction (enthalpy) on subsurface layer. """
+"""
+Prognostic step for heat conduction (enthalpy) on subsurface layer.
+"""
 function prognosticstep!(::SubSurface, ::Heat{<:FreezeCurve,Enthalpy}, state)
     Δk = Δ(state.grids.k) # cell sizes
     ΔT = Δ(state.grids.T)
     # Diffusion on non-boundary cells
-    heatconduction!(state.dH,state.T,ΔT,state.k,Δk)
-end
-""" Prognostic step for heat conduction (partitioned enthalpy) on subsurface layer."""
-function prognosticstep!(::SubSurface, heat::Heat{<:FreezeCurve,PartitionedEnthalpy}, state)
-    Δk = Δ(state.grids.k) # cell sizes
-    ΔT = Δ(state.grids.T)
-    # Diffusion on non-boundary cells
-    heatconduction!(state.dH,state.T,ΔT,state.k,Δk)
-    let L = heat.L;
-        @. state.dHₛ = state.dH / (L/state.C*state.dθdT + 1)
-        # This could also be expressed via a mass matrix with 1
-        # in the upper right block diagonal. But this is easier.
-        @. state.dHₗ = state.dH - state.dHₛ
-    end
+    heatconduction!(state.dH, state.T, ΔT, state.k, Δk)
 end
-""" Prognostic step for heat conduction (temperature) on subsurface layer. """
-function prognosticstep!(::SubSurface, ::Heat{<:FreezeCurve,Temperature}, state)
+"""
+Prognostic step for heat conduction (temperature) on subsurface layer.
+"""
+function prognosticstep!(sub::SubSurface, ::Heat{<:FreezeCurve,Temperature}, state)
     Δk = Δ(state.grids.k) # cell sizes
     ΔT = Δ(state.grids.T)
     # Diffusion on non-boundary cells
-    heatconduction!(state.dH,state.T,ΔT,state.k,Δk)
-    # Compute temperature flux by dividing by C_eff;
-    # C_eff should be computed by the freeze curve.
-    @inbounds @. state.dT = state.dH / state.Ceff
+    heatconduction!(state.dH, state.T, ΔT, state.k, Δk)
+    # Compute temperature flux by dividing by dHdT;
+    # dHdT should be computed by the freeze curve.
+    @inbounds @. state.dT = state.dH / state.dHdT
     return nothing
 end
-boundaryflux(::Neumann, bc::BoundaryProcess, top::Top, heat::Heat, sub::SubSurface, stop, ssub) = boundaryvalue(bc,top,heat,sub,stop,ssub)
-boundaryflux(::Neumann, bc::BoundaryProcess, bot::Bottom, heat::Heat, sub::SubSurface, sbot, ssub) = boundaryvalue(bc,bot,heat,sub,sbot,ssub)
-function boundaryflux(::Dirichlet, bc::BoundaryProcess, top::Top, heat::Heat, sub::SubSurface, stop, ssub)
-    Δk = Δ(ssub.grids.k)
+# 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
     @inbounds let Tupper=boundaryvalue(bc,top,heat,sub,stop,ssub),
         Tsub=ssub.T[1],
         k=ssub.k[1],
-        d=Δk[1],
-        δ₀=(d/2); # distance to boundary
-        -k*(Tsub-Tupper)/δ₀/d
+        δ=Δk[1]/2; # distance to boundary
+        -k*(Tsub-Tupper)/δ
     end
 end
-function boundaryflux(::Dirichlet, bc::BoundaryProcess, bot::Bottom, heat::Heat, sub::SubSurface, sbot, ssub)
-    Δk = Δ(ssub.grids.k)
+@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
     @inbounds let Tlower=boundaryvalue(bc,bot,heat,sub,sbot,ssub),
         Tsub=ssub.T[end],
         k=ssub.k[end],
-        d=Δk[1],
-        δ₀=(d/2); # distance to boundary
-        -k*(Tsub-Tlower)/δ₀/d
+        δ=Δk[end]/2; # distance to boundary
+        -k*(Tsub-Tlower)/δ
     end
 end
 """
 Generic top interaction. Computes flux dH at top cell.
 """
-function interact!(top::Top, bc::BoundaryProcess, sub::SubSurface, heat::Heat, stop, ssub)
-    # thermal conductivity at boundary
-    # assumes (1) k has already been computed, (2) surface conductivity = cell conductivity
-    @inbounds ssub.k[1] = ssub.k[2]
+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
-    @inbounds ssub.dH[1] += boundaryflux(bc, top, heat, sub, stop, ssub)
+    @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::BoundaryProcess, ssub, sbot)
-    # thermal conductivity at boundary
-    # assumes (1) k has already been computed, (2) bottom conductivity = cell conductivity
-    @inbounds ssub.k[end] = ssub.k[end-1]
+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
-    @inbounds ssub.dH[end] += boundaryflux(bc, bot, heat, sub, sbot, ssub)
+    @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!(::SubSurface, ::Heat, ::SubSurface, ::Heat, s1, s2)
+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],
-        Δ₁ = Δ(s1.grids.k)[end],
-        Δ₂ = Δ(s2.grids.k)[1];
+        Δ₁ = Δk₁[end],
+        Δ₂ = Δk₂[1];
         k = harmonicmean(k₁, k₂, Δ₁, Δ₂);
         s1.k[end] = s2.k[1] = k
     end
@@ -197,42 +281,58 @@ function interact!(::SubSurface, ::Heat, ::SubSurface, ::Heat, s1, s2)
     Qᵢ = @inbounds let k₁ = s1.k[end],
         k₂ = s2.k[1],
         k = k₁ = k₂,
-        δ = s2.grids.T[1] - s1.grids.T[end];
+        z₁ = midpoints(sub1, s1),
+        z₂ = midpoints(sub2, s2),
+        δ = z₂[1] - z₁[end];
         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ᵢ / Δ(s1.grids.k)[end]
-    @inbounds s2.dH[1] += -Qᵢ / Δ(s2.grids.k)[1]
+    @inbounds s1.dH[end] += Qᵢ / Δk₁[end]
+    @inbounds s2.dH[1] += -Qᵢ / Δk₂[1]
     return nothing # ensure no allocation
 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
+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,
+        Lθ = L*θtot,
+        I_t = H > Lθ,
+        I_c = (H > 0.0) && (H <= Lθ);
+        (I_c*(H/Lθ) + I_t)θtot
+    end
+end
 """
-Implementation of "free water" freeze curve for any subsurface layer. Assumes that
-'state' contains at least temperature (T), enthalpy (H), heat capacity (C),
+    freezethaw!(sub::SubSurface, heat::Heat{FreeWater,Enthalpy}, state)
+
+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).
 """
-@inline function (fc::FreeWater)(layer::SubSurface, heat::Heat{FreeWater,Enthalpy}, state)
-    @inline function enthalpyinv(H, C, L, θtot)
-        let θtot = max(1.0e-8,θtot),
-            Lθ = L*θtot,
-            I_t = H > Lθ,
-            I_f = H <= 0.0;
-            (I_t*(H-Lθ) + I_f*H)/C
-        end
-    end
-    @inline function freezethaw(H, L, θtot)
-        let θtot = max(1.0e-8,θtot),
-            Lθ = L*θtot,
-            I_t = H > Lθ,
-            I_c = (H > 0.0) && (H <= Lθ);
-            I_c*(H/Lθ) + I_t
-        end
-    end
-    let L = heat.L,
-        θw = totalwater(layer, heat, state);
-        @. state.θl = freezethaw(state.H, L, θw)*θw
-        heatcapacity!(layer, heat, state) # update heat capacity, C
-        @. state.T = enthalpyinv(state.H, state.C, L, θ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)
+        # set dHdT (a.k.a dHdT)
+        state.dHdT[i] = state.T[i] ≈ 0.0 ? 1e8 : 1/C
     end
     return nothing
 end

src/Physics/HeatConduction/heat_bc.jl

 # Boundary condition type aliases
-const ConstantTemp = Constant{Dirichlet,Float"°C"}
-ConstantTemp(value::UFloat"K") = Constant(Dirichlet, dustrip(u"°C", value))
-ConstantTemp(value::UFloat"°C") = Constant(Dirichlet, dustrip(value))
-const GeothermalHeatFlux = Constant{Neumann,Float"J/s/m^2"}
-GeothermalHeatFlux(value::UFloat"J/s/m^2"=0.053xu"J/s/m^2") = Constant(Neumann, dustrip(value))
+const HeatBC = BoundaryProcess{T} where {Heat<:T<:SubSurfaceProcess}
+const ConstantTemp = ConstantBC{Heat, Dirichlet,Float"°C"}
+ConstantTemp(value::UFloat"K") = ConstantBC(Heat, Dirichlet, dustrip(u"°C", value))
+ConstantTemp(value::UFloat"°C") = ConstantBC(Heat, Dirichlet, dustrip(value))
+const GeothermalHeatFlux = ConstantBC{Heat, Neumann, Float"J/s/m^2"}
+GeothermalHeatFlux(value::UFloat"J/s/m^2"=0.053xu"J/s/m^2") = ConstantBC(Heat, Neumann, dustrip(value))
 
-struct TemperatureGradient{E,F} <: BoundaryProcess
+struct TemperatureGradient{E,F} <: BoundaryProcess{Heat}
     T::F # temperature forcing
     effect::E # effect
     TemperatureGradient(T::F, effect::E=nothing) where {F<:Forcing,E} = new{E,F}(T, effect)
 end
 BoundaryStyle(::Type{<:TemperatureGradient}) = Dirichlet()
-BoundaryStyle(::Type{<:TemperatureGradient{<:Damping}}) = Neumann()
-@inline boundaryvalue(bc::TemperatureGradient,l2,p2,s1,s2) where {F} = bc.T(s1.t)
+@inline boundaryvalue(bc::TemperatureGradient, l1, ::Heat, l2, s1, s2) = getscalar(s1.T_ub)
 
-@with_kw struct NFactor{W,S} <: BoundaryEffect
-    winterfactor::W = Param(1.0, bounds=(0.0,1.0)) # applied when Tair <= 0
-    summerfactor::S = Param(1.0, bounds=(0.0,1.0)) # applied when Tair > 0
+variables(::Top, bc::TemperatureGradient) = (
+    Diagnostic(:T_ub, Scalar, u"K"),
+)
+function diagnosticstep!(::Top, bc::TemperatureGradient, state)
+    @setscalar state.T_ub = bc.T(state.t)
 end
-@inline function boundaryvalue(bc::TemperatureGradient{<:NFactor}, l1, ::Heat, l2, s1, s2)
-    let nfw = bc.effect.winterfactor,
-        nfs = bc.effect.summerfactor,
-        Tair = bc.T(s1.t),
-        nf = (Tair <= zero(Tair))*nfw + (Tair > zero(Tair))*nfs;
-        nf*Tair # apply damping factor to air temperature
-    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
 end
-# damped temperature gradient
-@inline function boundaryvalue(bc::TemperatureGradient{<:Damping}, l1, ::Heat, l2, s1,s2)
-    let Ttop = bc.T(s1.t),
-        Tsub = s2.T[1],
-        δ = Δ(s2.grids.k)[1], # grid cell size
-        k_sub = s2.kc[1];
-        bc.effect(Ttop, Tsub, k_sub, δ, s1.t)
-    end
+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)
+    nfw = bc.effect.nf
+    nfs = bc.effect.nt
+    Tair = bc.T(state.t)
+    @setscalar state.nfactor = nfactor(Tair, nfw, nfs)
+    @setscalar state.T_ub = getscalar(state.nfactor)*Tair
 end

src/Physics/HeatConduction/soil/soilheat.jl

-totalwater(soil::Soil, heat::Heat, state) = θw(soil.para)
-function heatcapacity(soil::Soil, totalwater, liquidwater, mineral, organic)
-    @unpack cw,co,cm,ca,ci = soil.hc
-    let air = 1.0 - totalwater - mineral - organic,
-        ice = totalwater - liquidwater,
-        liq = liquidwater;
-        liq*cw + ice*ci + mineral*cm + organic*co + air*ca
-    end
-end
-function thermalconductivity(soil::Soil, totalwater, liquidwater, mineral, organic)
-    @unpack kw,ko,km,ka,ki = soil.tc
-    let air = 1.0 - totalwater - mineral - organic,
-        ice = totalwater - liquidwater,
-        liq = liquidwater;
-        (liq*kw^0.5 + ice*ki^0.5 + mineral*km^0.5 + organic*ko^0.5 + air*ka^0.5)^2
-    end
-end
-""" Heat capacity for soil layer """
-function heatcapacity!(soil::Soil, ::Heat, state)
-    let w = θw(soil.para),
-        m = θm(soil.para),
-        o = θo(soil.para);
-        @. state.C = heatcapacity(soil, w, state.θl, m, o)
-    end
-end
-""" Thermal conductivity for soil layer """
-function thermalconductivity!(soil::Soil, ::Heat, state)
-    let w = θw(soil.para),
-        m = θm(soil.para),
-        o = θo(soil.para);
-        @. state.kc = thermalconductivity(soil, w, state.θl, m, o)
-    end
-end
-
-include("sfcc.jl")
-
-""" Initial condition for heat conduction (all state configurations) on soil layer. """
-function initialcondition!(soil::Soil, heat::Heat{<:SFCC}, state)
-    L = heat.L
-    sfcc = freezecurve(heat)
-    state.θl .= sfcc.f.(state.T, sfccparams(sfcc.f, soil, heat, state)...)
-    heatcapacity!(soil, heat, state)
-    @. state.H = enthalpy(state.T, state.C, L, state.θl)
-end
-""" Diagonstic step for partitioned heat conduction on soil layer. """
-function initialcondition!(soil::Soil, heat::Heat{<:SFCC,PartitionedEnthalpy}, state)
-    L = heat.L
-    sfcc = freezecurve(heat)
-    state.θl .= sfcc.f.(state.T, sfccparams(sfcc.f, soil, heat, state)...)
-    heatcapacity!(soil, heat, state)
-    @. state.Hₛ = state.T*state.C
-    @. state.Hₗ = state.θl*L
-end

src/Physics/Physics.jl

 module Physics
 
-import CryoGrid: Process, CompoundProcess, Coupled, Layer, Top, Bottom, SubSurface
-import CryoGrid: diagnosticstep!, initialcondition!, interact!, prognosticstep!, variables, observe
+import CryoGrid: Process, CoupledProcesses, BoundaryProcess, Coupled, Layer, Top, Bottom, SubSurface
+import CryoGrid: diagnosticstep!, initialcondition!, interact!, prognosticstep!, variables, callbacks, observe
 
-using CryoGrid.Layers
 using CryoGrid.Numerics
 using CryoGrid.Utils
 
@@ -11,15 +10,29 @@ using Lazy: @>>
 using Reexport
 using Unitful
 
-include("compound.jl")
+"""
+    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)
+
+include("coupled.jl")
 include("Boundaries/Boundaries.jl")
-include("Water/Water.jl")
 include("HeatConduction/HeatConduction.jl")
+include("WaterBalance/WaterBalance.jl")
+include("Soils/Soils.jl")
 include("SEB/SEB.jl")
 include("Sources/Sources.jl")
 
 @reexport using .Boundaries
 @reexport using .HeatConduction
+@reexport using .WaterBalance
+@reexport using .Soils
 @reexport using .SEB
 @reexport using .Sources
 

src/Physics/SEB/SEB.jl

@@ -3,16 +3,17 @@ module SEB
 using ..HeatConduction: Heat
 using ..Physics
 using ..Boundaries
-using CryoGrid.Layers: Soil
+using CryoGrid.InputOutput: Forcing
+using CryoGrid.Physics.Soils
 using CryoGrid.Numerics
 using CryoGrid.Utils
 
-using Parameters
-
 import CryoGrid: BoundaryProcess, BoundaryStyle, Neumann, Top
 import CryoGrid: initialcondition!, variables
 
-@with_kw struct SEBParams <: Params
+using Unitful
+
+Base.@kwdef struct SEBParams <: IterableStruct
     # surface properties --> should be associated with the Stratigraphy and maybe made state variables
     α::Float"1" = 0.2xu"1"                          # surface albedo [-]
     ϵ::Float"1" = 0.97xu"1"                         # surface emissivity [-]
@@ -31,7 +32,7 @@ import CryoGrid: initialcondition!, variables
     cₐ::Float"J/(m^3*K)" = 1005.7xu"J/(kg*K)" * ρₐ  # volumetric heat capacity of dry air at standard pressure and 0°C [J/(m^3*K)]
 end
 
-struct SurfaceEnergyBalance{F} <: BoundaryProcess
+struct SurfaceEnergyBalance{F} <: BoundaryProcess{Heat}
     forcing::F
     sebparams::SEBParams
     function SurfaceEnergyBalance(