Project.toml

 name = "CryoGrid"
 uuid = "a535b82e-5f3d-4d97-8b0b-d6483f5bebd5"
 authors = ["Brian Groenke <brian.groenke@awi.de>", "Moritz Langer <moritz.langer@awi.de>"]
-version = "0.10.3"
+version = "0.11.0"
 
 [deps]
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
@@ -14,6 +14,7 @@ DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 Flatten = "4c728ea3-d9ee-5c9a-9642-b6f7d7dc04fa"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+FreezeCurves = "71e4ad71-e4f2-45a3-aa0a-91ffaa9676be"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
@@ -23,7 +24,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 ModelParameters = "4744a3fa-6c31-4707-899e-a3298e4618ad"
-NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
@@ -31,6 +32,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SimulationLogs = "e3c6736c-93d9-479d-93f3-96ff5e8836d5"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StructTypes = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
@@ -38,8 +40,10 @@ TimeSeries = "9e3dc215-6440-5c97-bce1-76c03772f85e"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
-julia = ">= 1.6"
-PreallocationTools = ">= 0.3.1"
+FreezeCurves = "0.2"
+OrdinaryDiffEq = "6"
+PreallocationTools = "0.4"
+julia = "1.6,1.7"
 
 [extras]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"

src/IO/InputOutput.jl

 module InputOutput
 
 import CryoGrid
+
 using CryoGrid.Numerics
+using CryoGrid.Utils
 
 using Base: @propagate_inbounds
 using ComponentArrays
+using ConstructionBase
 using DataStructures: DefaultDict
 using Dates
 using DimensionalData
@@ -23,10 +26,10 @@ import DimensionalData: stack
 export loadforcings
 
 include("ioutils.jl")
+export CryoGridParams, parameterize
+include("params.jl")
 export Forcing, TimeSeriesForcing
 include("forcings.jl")
-export CryoGridParams
-include("params.jl")
 export CryoGridOutput
 include("output.jl")
 

src/IO/forcings.jl

@@ -8,6 +8,7 @@ abstract type Forcing{T,N} end
 @inline @propagate_inbounds (forcing::Forcing)(t::DateTime) = forcing(ustrip(u"s", float(Dates.datetime2epochms(t))u"ms"))
 # disable flattening for all fields of forcing types by default
 Flatten.flattenable(::Type{<:Forcing}, ::Type) = false
+InputOutput.parameterize(f::Forcing; fields...) = f
 
 """
       TimeSeriesForcing{T,A,I}

src/IO/params.jl

@@ -40,14 +40,6 @@ function Base.show(io::IO, ::MIME"text/plain", ps::CryoGridParams{T}) where T
 end
 Tables.columns(ps::CryoGridParams) = Tables.columns(ps.table)
 Tables.rows(ps::CryoGridParams) = Tables.rows(ps.table)
-function CryoGridParams(obj; full_metadata=false)
-    m = Model(obj)
-    if full_metadata
-        m[:idx] = 1:length(m)
-        m = _setparafields(m)
-    end
-    return CryoGridParams(m)
-end
 function _setparafields(m::Model)
     function _setparafield(name, type::Type, para::CryoGrid.Parameterization)
         if length(ModelParameters.params(para)) > 0
@@ -68,3 +60,37 @@ function _setparafields(m::Model)
     newparent = Flatten.reconstruct(parent(m), updated_parameterizations, CryoGrid.Parameterization)
     return Model(newparent)
 end
+function CryoGridParams(obj; full_metadata=false)
+    m = Model(obj)
+    if full_metadata
+        m[:idx] = 1:length(m)
+        m = _setparafields(m)
+    end
+    return CryoGridParams(m)
+end
+"""
+    parameterize(x::T) where {T}
+    parameterize(x::Unitful.AbstractQuantity; fields...)
+    parameterize(x::Number; fields...)
+    parameterize(p::Param; ignored...)
+
+If `x` is a numeric type, `x` will be wrapped in a `ModelParameters.Param`, including a `units` field if `x` has units.
+If `x` is a `Param` type, `x` will be returned as-is.
+If `x` is some other type, `x` will be recursively unpacked and `parameterize` called on each field. It may be necessary
+or desirable for some types to override `parameterize` to define custom behavior, e.g. if only some fields should be
+parameterized.
+"""
+function parameterize(x::Unitful.AbstractQuantity; fields...)
+    let x = normalize_units(x);
+        Param(ustrip(x); untis=unit(x), fields...)
+    end
+end
+function parameterize(x::T; fields...) where {T}
+    _parameterize(x) = parameterize(x; fields...)
+    new_fields = map(_parameterize, getfields(x))
+    ctor = ConstructionBase.constructorof(T)
+    return ctor(new_fields...)
+end
+parameterize(x::Number; fields...) = Param(x; fields...)
+parameterize(x::Param; ignored...) = x
+parameterize(x::AbstractArray; ignored...) = x

src/Numerics/init.jl

@@ -25,7 +25,7 @@ struct InterpInitializer{varname,P,I,E} <: VarInitializer{varname}
     extrap::E
     InterpInitializer(varname::Symbol, profile::P, interp::I=Linear(), extrap::E=Flat()) where {P<:Profile,I,E} = new{varname,P,I,E}(profile, interp, extrap)
 end
-function (init::InterpInitializer{var})(state) where var
+function (init::InterpInitializer{var})(layer, state) where var
     profile, interp, extrap = init.profile, init.interp, init.extrap
     depths = collect(map(knot -> dustrip(knot.depth), profile.knots))
     u = getproperty(state, var)
@@ -51,6 +51,6 @@ Base.getindex(init::InterpInitializer{var}, itrv::Interval) where var = InterpIn
 
 Convenience constructor for `VarInitializer` that selects the appropriate initializer type based on the arguments.
 """
-initializer(varname::Symbol, x::Number) = FunctionInitializer(varname, state -> getproperty(state, varname) .= x)
+initializer(varname::Symbol, x::Number) = FunctionInitializer(varname, (layer,state) -> getproperty(state, varname) .= x)
 initializer(varname::Symbol, f::Function) = FunctionInitializer(varname, f)
 initializer(varname::Symbol, profile::Profile, interp=Linear(), extrap=Flat()) = InterpInitializer(varname, profile, interp, extrap)

src/Numerics/math.jl

-∇(f, x) = ∇(typeof(x), f, x)
-function ∇(::Type{T}, f, x) where {T}
-    res = ForwardDiff.derivative!(ForwardDiff.DiffResult(zero(T), x), f, x)
+function ∇(f::F, x::Number) where {F}
+    res = ForwardDiff.derivative!(ForwardDiff.DiffResult(zero(x), zero(x)), f, x)
     return res.value, res.derivs[1]
 end
-
+function ∇(f::F, x::AbstractArray) where {F}
+    res = ForwardDiff.gradient!(ForwardDiff.DiffResult(eltype(x), zero(x)), f, x)
+    return res.value, res.derivs
+end
+# Flux calculations
+@propagate_inbounds @inline _flux_kernel(x₁, x₂, Δx, k) = -k*(x₂ - x₁)/ Δx
+@propagate_inbounds @inline function _flux!(j::AbstractVector, x₁::AbstractVector, x₂::AbstractVector, Δx::AbstractVector, k::AbstractVector, ::Val{use_turbo}) where {use_turbo}
+    @. j = _flux_kernel(x₁, x₂, Δx, k)
+end
+@propagate_inbounds @inline function _flux!(j::AbstractVector{Float64}, x₁::AbstractVector{Float64}, x₂::AbstractVector{Float64}, Δx::AbstractVector{Float64}, k::AbstractVector{Float64}, ::Val{true})
+    @turbo @. j = _flux_kernel(x₁, x₂, Δx, k)
+end
 """
-    finitediff!(∂x::AbstractVector, x::AbstractVector, Δ::AbstractVector)
+    flux!(j::AbstractVector, x::AbstractVector, Δx::AbstractVector, k::AbstractVector)
 
-First order forward finite difference operator.
+Calculates the first-order, non-linear spatial flux over a discretized variable `x` with conductivity `k`.
+`x` is assumed to have shape `(N,)`, `Δx` shape `(N-1,)`, and `j` and `k` shape `(N+1,)` such that `j[2:end-1]` represents
+the fluxes over the inner grid cell faces. Fluxes are added to existing values in `j`.
 """
-function finitediff!(∂x::AbstractVector, x::AbstractVector, Δ::AbstractVector)
-    @inbounds let x₁ = (@view x[1:end-1]),
-        x₂ = (@view x[2:end]);
-        @. ∂x = (x₂ - x₁) / Δ
+function flux!(j::AbstractVector, x::AbstractVector, Δx::AbstractVector, k::AbstractVector)
+    @inbounds let x₁ = @view(x[1:end-1]),
+        x₂ = @view(x[2:end]),
+        k = @view(k[2:end-1]),
+        j = @view(j[2:end-1]);
+        _flux!(j, x₁, x₂, Δx, k, Val{USE_TURBO}())
     end
 end
+# Divergence
+@propagate_inbounds @inline _div_kernel(j₁, j₂, Δj) = (j₁ - j₂) / Δj
+@propagate_inbounds @inline function _div!(dx::AbstractVector, j₁::AbstractVector, j₂::AbstractVector, Δj::AbstractVector, ::Val{use_turbo}) where {use_turbo}
+    @. dx = _div_kernel(j₁, j₂, Δj)
+end
+@propagate_inbounds @inline function _div!(dx::AbstractVector{Float64}, j₁::AbstractVector{Float64}, j₂::AbstractVector{Float64}, Δj::AbstractVector{Float64}, ::Val{true})
+    @turbo @. dx = _div_kernel(j₁, j₂, Δj)
+end
 """
-    lineardiffusion!(∂y::AbstractVector, x::AbstractVector, Δ::AbstractVector, k::Number)
+    divergence!(dx::AbstractVector, j::AbstractVector, Δj::AbstractVector)
 
-Second order Laplacian with constant diffusion k.
+Calculates the first-order divergence over a 1D flux vector field `j` and grid cell lengths `Δj`. Divergences are added to existing values in `dx`.
 """
-function lineardiffusion!(∂y::AbstractVector, x::AbstractVector, Δ::AbstractVector, k::Number)
-    @inbounds let x₁ = (@view x[1:end-2]),
-        x₂ = (@view x[2:end-1]),
-        x₃ = (@view x[3:end]),
-        Δ₁ = (@view Δ[1:end-1]),
-        Δ₂ = (@view Δ[2:end]);
-        @. ∂y = k*((x₃ - x₂)/Δ₂ - (x₂-x₁)/Δ₁)/Δ₁
+function divergence!(dx::AbstractVector, j::AbstractVector, Δj::AbstractVector)
+    @inbounds let j₁ = @view(j[1:end-1]),
+        j₂ = @view(j[2:end]);
+        _div!(dx, j₁, j₂, Δj, Val{USE_TURBO}())
     end
 end
-"""
-    nonlineardiffusion!(∂y, x, Δx, k, Δk)
-
-Second order Laplacian with non-linear diffusion operator, `k`. Accelerated using `LoopVectorization.@turbo` for `Float64` vectors.
-"""
-function nonlineardiffusion!(∂y::AbstractVector, x::AbstractVector, Δx::AbstractVector, k::AbstractVector, Δk::AbstractVector)
-    @inbounds let x₁ = (@view x[1:end-2]),
-        x₂ = (@view x[2:end-1]),
-        x₃ = (@view x[3:end]),
-        k₁ = (@view k[1:end-1]),
-        k₂ = (@view k[2:end]),
-        Δx₁ = (@view Δx[1:end-1]),
-        Δx₂ = (@view Δx[2:end]);
-        nonlineardiffusion!(∂y, x₁, x₂, x₃, k₁, k₂, Δx₁, Δx₂, Δk)
+# non-linear diffusion
+@propagate_inbounds @inline function _nonlineardiffusion_kernel(x₁, x₂, Δx, j₁, k, Δj)
+    j₂ = _flux_kernel(x₁, x₂, Δx, k)
+    div = _div_kernel(j₁, j₂, Δj)
+    return j₂, div
+end
+@propagate_inbounds @inline function _nonlineardiffusion!(dx::AbstractVector, j::AbstractVector, x::AbstractVector, Δx::AbstractVector, k::AbstractVector, Δk::AbstractVector, i::Integer)
+    let x₁ = x[i-1],
+        x₂ = x[i],
+        j₁ = j[i-1],
+        k = k[i],
+        δx = Δx[i-1],
+        δj = Δk[i-1];
+        j₂, divx₁ = _nonlineardiffusion_kernel(x₁, x₂, δx, j₁, k, δj)
+        j[i] += j₂
+        dx[i-1] += divx₁
     end
 end
 """
-    nonlineardiffusion(x₁, x₂, x₃, k₁, k₂, Δx₁, Δx₂, Δk)
+    nonlineardiffusion!(dx::AbstractVector, j::AbstractVector, x::AbstractVector, Δx::AbstractVector, k::AbstractVector, Δk::AbstractVector)
 
-Scalar second order Laplacian with non-linear diffusion operator, `k`.
+Fast alternative to `flux!` and `divergence!` which computes fluxes and divergences (via `_flux_kernel` and `_div_kernel`) in a single pass. Note, however, that
+loop vectorization with `@turbo` is not possible because of necessary loop-carried dependencies. Fluxes and divergences are added to the existing values stored
+in `j` and `dx`.
 """
-@inline nonlineardiffusion(x₁, x₂, x₃, k₁, k₂, Δx₁, Δx₂, Δk) = (k₂*(x₃ - x₂)/Δx₂ - k₁*(x₂ - x₁)/Δx₁)/Δk
-@propagate_inbounds function nonlineardiffusion!(∂y, x₁, x₂, x₃, k₁, k₂, Δx₁, Δx₂, Δk)
-    @. ∂y = nonlineardiffusion(x₁, x₂, x₃, k₁, k₂, Δx₁, Δx₂, Δk)
-end
-@propagate_inbounds function nonlineardiffusion!(
-    ∂y::AbstractVector{Float64},
-    x₁::AbstractVector{Float64},
-    x₂::AbstractVector{Float64},
-    x₃::AbstractVector{Float64},
-    k₁::AbstractVector{Float64},
-    k₂::AbstractVector{Float64},
-    Δx₁::AbstractVector{Float64},
-    Δx₂::AbstractVector{Float64},
-    Δk::AbstractVector{Float64},
-)
-    if USE_TURBO
-        @turbo @. ∂y = nonlineardiffusion(x₁, x₂, x₃, k₁, k₂, Δx₁, Δx₂, Δk)
-    else
-        @. ∂y = nonlineardiffusion(x₁, x₂, x₃, k₁, k₂, Δx₁, Δx₂, Δk)
+function nonlineardiffusion!(dx::AbstractVector, j::AbstractVector, x::AbstractVector, Δx::AbstractVector, k::AbstractVector, Δk::AbstractVector) 
+    @inbounds for i in 2:length(j)-1
+        _nonlineardiffusion!(dx, j, x, Δx, k, Δk, i)
     end
+    # compute divergence for last grid cell
+    @inbounds dx[end] += _div_kernel(j[end-1], j[end], Δk[end])
+    return nothing
 end
-
+# other helper functions
 """
     harmonicmean(x₁, x₂, w₁, w₂)
 
@@ -169,38 +179,3 @@ softplusinv(x) = let x = clamp(x, eps(), Inf); IfElse.ifelse(x > 34, x, log(exp(
 # sqrt ∘ plusone == x -> sqrt(x. + 1.0)
 minusone(x) = x .- one.(x)
 plusone(x) = x .+ one.(x)
-
-# Function tabulation
-"""
-    Tabulated(f, argknots...)
-
-Alias for `tabulate` intended for function types.
-"""
-Tabulated(f, argknots...) = tabulate(f, argknots...)
-"""
-    tabulate(f, argknots::Pair{Symbol,<:Union{Number,AbstractArray}}...)
-
-Tabulates the given function `f` using a linear, multi-dimensional interpolant.
-Knots should be given as pairs `:arg => A` where `A` is a `StepRange` or `Vector`
-of input values (knots) at which to evaluate the function. `A` may also be a
-`Number`, in which case a pseudo-point interpolant will be used (i.e valid on
-`[A,A+ϵ]`). No extrapolation is provided by default but can be configured via
-`Interpolations.extrapolate`.
-"""
-function tabulate(f, argknots::Pair{Symbol,<:Union{Number,AbstractArray}}...)
-    initknots(a::AbstractArray) = Interpolations.deduplicate_knots!(a)
-    initknots(x::Number) = initknots([x,x])
-    interp(::AbstractArray) = Gridded(Linear())
-    interp(::Number) = Gridded(Constant())
-    extrap(::AbstractArray) = Flat()
-    extrap(::Number) = Throw()
-    names = map(first, argknots)
-    # get knots for each argument, duplicating if only one value is provided
-    knots = map(initknots, map(last, argknots))
-    f_argnames = Utils.argnames(f)
-    @assert all(map(name -> name ∈ names, f_argnames)) "Missing one or more arguments $f_argnames in $f"
-    arggrid = Iterators.product(knots...)
-    # evaluate function construct interpolant
-    f = extrapolate(interpolate(Tuple(knots), map(Base.splat(f), arggrid), map(interp ∘ last, argknots)), map(extrap ∘ last, argknots))
-    return f
-end

src/Numerics/variables.jl

@@ -37,16 +37,25 @@ end
 """
     Delta{dname,name,S,T,units,domain} <: Var{dname,S,T,units,domain}
 
-Defines a "delta" term `du` for variable `u`, which is the time-derivative or flux for prognostic variables and
+Defines a "delta" term `du` for variable `u`, which is the time-derivative/divergence for prognostic variables and
 the residual for algebraic variables.
 """
 struct Delta{dname,name,S,T,units,domain} <: Var{dname,S,T,units,domain}
     dim::S
-    Delta(::Symbol, dims::OnGrid, args...) = error("Off-cell prognostic/algebraic spatial variables are not currently supported.")
     Delta(dname::Symbol, name::Symbol, dims::Union{<:Shape,OnGrid{Cells,typeof(identity)}}, units=NoUnits, ::Type{T}=Float64; domain=-Inf..Inf) where {T} = new{dname,name,typeof(dims),T,units,domain}(dims)
     Delta(var::Prognostic{name,S,T,units,domain}) where {name,S,T,units,domain} = let dims=vardims(var); new{Symbol(:d,name),name,typeof(dims),T,upreferred(units)/u"s",domain}(dims) end
     Delta(var::Algebraic{name,S,T,units,domain}) where {name,S,T,units,domain} = let dims=vardims(var); new{Symbol(:d,name),name,typeof(dims),T,units,domain}(dims) end
 end
+"""
+    Flux{jname,name,S,T,units,domain} <: Var{jname,S,T,units,domain}
+
+Defines a "flux" term `ju` for prognostic variable `u`, which is typically a function of the gradient w.r.t space.
+"""
+struct Flux{jname,name,S,T,units,domain} <: Var{jname,S,T,units,domain}
+    dim::S
+    Flux(jname::Symbol, name::Symbol, dims::Union{<:Shape,OnGrid{Cells,typeof(identity)}}, units=NoUnits, ::Type{T}=Float64; domain=-Inf..Inf) where {T} = new{jname,name,typeof(dims),T,units,domain}(dims)
+    Flux(var::Prognostic{name,OnGrid{Cells,typeof(identity)},T,units,domain}) where {name,T,units,domain} = let dims=OnGrid(Edges, vardims(var).f); new{Symbol(:j,name),name,typeof(dims),T,upreferred(units)/u"s",domain}(dims) end
+end
 """
     Diagnostic{name,S,T,units,domain} <: Var{name,S,T,units,domain}
 
@@ -61,6 +70,7 @@ ConstructionBase.constructorof(::Type{Prognostic{name,S,T,units,domain}}) where
 ConstructionBase.constructorof(::Type{Diagnostic{name,S,T,units,domain}}) where {name,S,T,units,domain} = s -> Diagnostic(name, s, units, T; domain)
 ConstructionBase.constructorof(::Type{Algebraic{name,S,T,units,domain}}) where {name,S,T,units,domain} = s -> Algebraic(name, s, units, T; domain)
 ConstructionBase.constructorof(::Type{Delta{dname,name,S,T,units,domain}}) where {dname,name,S,T,units,domain} = s -> Delta(dname, name, s, units, T; domain)
+ConstructionBase.constructorof(::Type{Flux{jname,name,S,T,units,domain}}) where {jname,name,S,T,units,domain} = s -> Flux(jname, name, s, units, T; domain)
 ==(::Var{N1,S1,T1,u1,d1},::Var{N2,S2,T2,u2,d2}) where {N1,N2,S1,S2,T1,T2,u1,u2,d1,d2} = (N1 == N2) && (S1 == S2) && (T1 == T2) && (u1 == u2) && (d1 == d2)
 varname(::Var{name}) where {name} = name
 varname(::Type{<:Var{name}}) where {name} = name

src/Numerics/varstates.jl

@@ -97,8 +97,10 @@ function VarStates(vars::GroupedVars, D::Numerics.AbstractDiscretization, chunks
     diagvars = map(group -> filter(isdiagnostic, group), vars)
     progvars = map(group -> filter(isprognostic, group), vars)
     algvars = map(group -> filter(isalgebraic, group), vars)
-    # create variables for gradients/fluxes
+    # create variables for time delta variables (divergence/residual)
     dpvars = map(group -> map(Delta, filter(var -> isalgebraic(var) || isprognostic(var), group)), vars)
+    # create variables for fluxes (spatial gradients)
+    jpvars = map(group -> map(Flux, filter(var -> isprognostic(var) && isongrid(var), group)), vars)
     allprogvars = tuplejoin(_flatten(progvars), _flatten(algvars)) # flattened prognostic/algebraic variable group
     # TODO: not a currently necessary use case, but could be supported by partitioning the state array further by layer/group name;
     # this would require passing the name to discretize(...)
@@ -111,9 +113,9 @@ function VarStates(vars::GroupedVars, D::Numerics.AbstractDiscretization, chunks
     freediagvars = map(group -> filter(!isongrid, group), diagvars)
     freediagstate = map(group -> (;map(v -> varname(v) => DiffCache(varname(v), discretize(A, D, v), chunksize), group)...), freediagvars)
     # build gridded diagnostic state vectors
-    griddiagvars = filter(isongrid, _flatten(diagvars))
+    griddiagvars = Tuple(unique(filter(isongrid, _flatten(map(tuplejoin, diagvars,  jpvars)))))
     griddiagstate = map(v -> varname(v) => DiffCache(varname(v), discretize(A, D, v), chunksize), griddiagvars)
-    # join prognostic variables with flux variables, then build nested named tuples in each group with varnames as keys
-    allvars = map(vars -> NamedTuple{map(varname, vars)}(vars), map(tuplejoin, vars, dpvars))
+    # join prognostic variables with delta and flux variables, then build nested named tuples in each group with varnames as keys
+    allvars = map(vars -> NamedTuple{map(varname, vars)}(vars), map(tuplejoin, vars, dpvars, jpvars))
     VarStates(uproto, D, allvars, (;freediagstate...), (;griddiagstate...))
 end

src/Physics/HeatConduction/HeatConduction.jl

@@ -10,6 +10,7 @@ using CryoGrid.Utils
 
 using Base: @propagate_inbounds, @kwdef
 using IfElse
+using FreezeCurves: FreezeCurves, FreezeCurve, FreeWater
 using ModelParameters
 using Unitful
 
@@ -19,9 +20,6 @@ import CryoGrid.Physics
 export Heat, TemperatureProfile
 export FreeWater, FreezeCurve, freezecurve
 
-abstract type FreezeCurve end
-struct FreeWater <: FreezeCurve end
-
 """
     TemperatureProfile(pairs::Pair{<:Union{DistQuantity,Param},<:Union{TempQuantity,Param}}...)
 
@@ -38,24 +36,17 @@ abstract type HeatParameterization end
 struct Enthalpy <: HeatParameterization end
 struct Temperature <: HeatParameterization end
 
-abstract type StepLimiter end
-@kwdef struct CFL <: StepLimiter
-    fallback_dt::Float64 = 60.0 # fallback dt [s]
+@Base.kwdef struct ThermalProperties{Tconsts,TL,Tkw,Tki,Tka,Tcw,Tci,Tca}
+    consts::Tconsts = Physics.Constants()
+    L::TL = consts.ρw*consts.Lsl
+    kw::Tkw = 0.57u"W/m/K" # thermal conductivity of water [Hillel(1982)]
+    ki::Tki = 2.2u"W/m/K" # thermal conductivity of ice [Hillel(1982)]
+    ka::Tka = 0.025u"W/m/K" # air [Hillel(1982)]
+    cw::Tcw = 4.2e6u"J/K/m^3" # heat capacity of water
+    ci::Tci = 1.9e6u"J/K/m^3" # heat capacity of ice
+    ca::Tca = 0.00125e6u"J/K/m^3" # heat capacity of air
 end
 
-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
@@ -68,7 +59,7 @@ 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)
+Heat(para::Temperature; freezecurve, prop=ThermalProperties(), dtlim=Physics.CFL(), init=nothing) = Heat(para, prop, deepcopy(freezecurve), dtlim, init)
 
 # getter functions
 thermalproperties(heat::Heat) = heat.prop

src/Physics/HeatConduction/heat.jl

@@ -127,56 +127,23 @@ CryoGrid.variables(heat::Heat{<:FreezeCurve,Temperature}) = (
 Common variable definitions for all heat implementations.
 """
 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", 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(:θw, OnGrid(Cells), domain=0..1),
 )
-"""
-    heatconduction!(∂H,T,ΔT,k,Δk)
-
-1-D heat conduction/diffusion given T, k, and their deltas. Resulting enthalpy gradient is stored in ∂H.
-Note that this function does not perform bounds checking. It is up to the user to ensure that all variables are
-arrays of the correct length.
-"""
-function heatconduction!(∂H,T,ΔT,k,Δk)
-    # upper boundary
-    @inbounds ∂H[1] += let T₂=T[2],
-        T₁=T[1],
-        k=k[2],
-        ϵ=ΔT[1],
-        δ=Δk[1];
-        k*(T₂-T₁)/ϵ/δ
-    end
-    # diffusion on non-boundary cells
-    @inbounds let T = T,
-        k = (@view k[2:end-1]),
-        Δk = (@view Δk[2:end-1]),
-        ∂H = (@view ∂H[2:end-1]);
-        nonlineardiffusion!(∂H, T, ΔT, k, Δk)
-    end
-    # lower boundary
-    @inbounds ∂H[end] += let T₂=T[end],
-        T₁=T[end-1],
-        k=k[end-1],
-        ϵ=ΔT[end],
-        δ=Δk[end];
-        -k*(T₂-T₁)/ϵ/δ
-    end
-    return nothing
+function resetfluxes!(sub::SubSurface, heat::Heat, state)
+    # Reset energy fluxes to zero; this is redundant when H is the prognostic variable
+    # but necessary when it is not.
+    @. state.dH = zero(eltype(state.dH))
+    @. state.jH = zero(eltype(state.jH))
 end
 """
 Diagonstic step for heat conduction (all state configurations) on any subsurface layer.
 """
 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))
-    @. state.dH_upper = zero(eltype(state.dH_upper))
-    @. state.dH_lower = zero(eltype(state.dH_lower))
+    resetfluxes!(sub, heat, state)
     # Evaluate freeze/thaw processes
     freezethaw!(sub, heat, state)
     # Update thermal conductivity
@@ -196,20 +163,16 @@ 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
+    ssub.jH[1] += boundaryflux(bc, top, heat, sub, stop, ssub)
     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
+    # boundary flux; here we flip the sign since a positive flux is by convention downward
+    ssub.jH[end] -= boundaryflux(bc, bot, heat, sub, sbot, ssub)
     return nothing # ensure no allocation
 end
 """
@@ -226,18 +189,17 @@ function CryoGrid.interact!(sub1::SubSurface, ::Heat, sub2::SubSurface, ::Heat,
             Δ₂ = Δk₂[1];
             harmonicmean(k₁, k₂, Δ₁, Δ₂)
         end
-    # calculate heat flux between cells
+    # calculate heat flux between cells (positive downward)
     Qᵢ = @inbounds let z₁ = CryoGrid.midpoint(sub1, s1, last),
         z₂ = CryoGrid.midpoint(sub2, s2, first),
         δ = z₂ - z₁;
-        k*(s2.T[1] - s1.T[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ᵢ / Δk₁[end]
-    @inbounds s2.dH[1] += -Qᵢ / Δk₂[1]
+    # set inner boundary flux
+    s1.jH[end] += Qᵢ
+    # these should already be equal since jH[1] on s2 and jH[end] on s1 should point to the same array index;
+    # but we set them explicitly just to be safe.
+    s2.jH[1] = s1.jH[end]
     return nothing # ensure no allocation
 end
 """
@@ -245,18 +207,19 @@ Prognostic step for heat conduction (enthalpy) on subsurface layer.
 """
 function CryoGrid.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)
+    ΔT = Δ(state.grids.T) # midpoint distances
+    # compute internal fluxes and non-linear diffusion assuming boundary fluxes have been set
+    nonlineardiffusion!(state.dH, state.jH, state.T, ΔT, state.k, Δk)
+    return nothing
 end
 """
 Prognostic step for heat conduction (temperature) on subsurface layer.
 """
 function CryoGrid.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)
+    ΔT = Δ(state.grids.T) # midpoint distances
+    # compute internal fluxes and non-linear diffusion assuming boundary fluxes have been set
+    nonlineardiffusion!(state.dH, state.jH, 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
@@ -280,11 +243,12 @@ end
         Tsub=ssub.T[end],
         k=ssub.k[end],
         δ=Δk/2; # distance to boundary
-        -k*(Tsub-Tlower)/δ
+        # note again the inverted sign; positive here means *upward from* the bottom boundary
+        k*(Tlower-Tsub)/δ
     end
 end
 # CFL not defined for free-water freeze curve
-CryoGrid.timestep(::SubSurface, heat::Heat{FreeWater,Enthalpy,CFL}, state) = Inf
+CryoGrid.timestep(::SubSurface, heat::Heat{FreeWater,Enthalpy,Physics.CFL}, state) = Inf
 """
     timestep(::SubSurface, ::Heat{Tfc,TPara,CFL}, state) where {TPara}
 
@@ -292,7 +256,7 @@ Implementation of `timestep` for `Heat` using the Courant-Fredrichs-Lewy conditi
 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}
+@inline function CryoGrid.timestep(::SubSurface, heat::Heat{Tfc,TPara,Physics.CFL}, state) where {Tfc,TPara}
     Δx = Δ(state.grid)
     dtmax = Inf
     @inbounds for i in 1:length(Δx)
@@ -307,21 +271,14 @@ is taken to be the diffusivity: `dHdT / kc`.
 end
 # Free water freeze curve
 @inline function enthalpyinv(sub::SubSurface, heat::Heat{FreeWater,Enthalpy}, state, i)
-    f_hc = partial(heatcapacity, liquidwater, sub, heat, state, i)
-    return enthalpyinv(heat.freezecurve, f_hc, state.H[i], heat.prop.L, f_hc.θwi)
+    hc = partial(heatcapacity, liquidwater, sub, heat, state, i)
+    return enthalpyinv(heat.freezecurve, hc, state.H[i], hc.θwi, heat.prop.L)
 end
-@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θ);
-        # 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
+@inline function enthalpyinv(::FreeWater, hc::F, H, θwi, L) where {F}
+    θw, I_t, I_f, I_c, Lθ = FreezeCurves.freewater(H, θwi, L)
+    C = hc(θw)
+    T = (I_t*(H-Lθ) + I_f*H)/C
+    return T, θw, C
 end
 """
     freezethaw!(sub::SubSurface, heat::Heat{FreeWater,Enthalpy}, state)

src/Physics/HeatConduction/heat_bc.jl

@@ -24,9 +24,10 @@ function CryoGrid.diagnosticstep!(::Top, bc::TemperatureGradient, state)
 end
 
 Base.@kwdef struct NFactor{W,S} <: BoundaryEffect
-    nf::W = Param(1.0, domain=0..1) # applied when Tair <= 0
-    nt::S = Param(1.0, domain=0..1) # applied when Tair > 0
+    nf::W = 1.0 # applied when Tair <= 0
+    nt::S = 1.0 # applied when Tair > 0
 end
+InputOutput.parameterize(nf::NFactor; fields...) = NFactor(Param(nf.nf; domain=0..1, fields...), Param(nf.nt; domain=0..1, fields...))
 CryoGrid.variables(::Top, bc::TemperatureGradient{<:NFactor}) = (
     Diagnostic(:T_ub, Scalar, u"K"),
     Diagnostic(:nfactor, Scalar),

src/Physics/Hydrology/Hydrology.jl

@@ -5,7 +5,7 @@ import CryoGrid
 using CryoGrid
 using CryoGrid.Physics
 using CryoGrid.Numerics
-using CryoGrid.Numerics: nonlineardiffusion!
+using CryoGrid.Numerics: flux!, divergence!
 
 using IfElse
 using ModelParameters

src/Physics/Physics.jl

 module Physics
 
+import CryoGrid
+
 using CryoGrid
 using CryoGrid.Numerics
 using CryoGrid.Utils
 
-import CryoGrid
-
 using Reexport
 using Unitful
 

src/Physics/Snow/Snow.jl

@@ -8,6 +8,7 @@ using CryoGrid.Numerics
 using CryoGrid.Utils
 
 import CryoGrid
+import CryoGrid.InputOutput
 import CryoGrid.Physics
 import CryoGrid.Physics.HeatConduction
 
@@ -20,8 +21,8 @@ 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"),
+    ρsn_new = 250.0u"kg/m^3",
+    ρsn_old = 500.0u"kg/m^3",
 ) = (; ρw, ρsn_new, ρsn_old)
 
 """
@@ -50,8 +51,9 @@ 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
+    rate_scale::S = 1.0 # scaling factor for snowfall rate
 end
+InputOutput.parameterize(acc::LinearAccumulation{<:Real}; fields...) = LinearAccumulation(Param(acc.rate_sacle; domain=0..Inf, fields...))
 
 abstract type SnowDensityScheme end
 # constant density (using Snowpack properties)

src/Physics/Snow/snow_bulk.jl

@@ -75,14 +75,12 @@ 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]
+    ssnow.jH[1] += CryoGrid.boundaryflux(bc, top, heat, snow, stop, ssnow)
     return nothing # ensure no allocation
 end
 
@@ -101,6 +99,7 @@ function CryoGrid.diagnosticstep!(
 ) where {TAcc,TAbl<:DegreeDayMelt,TDen<:ConstantDensity}
     smb, heat = procs
     ρsn = snow.prop.ρsn_new
+    HeatConduction.resetfluxes!(snow, heat, state)
     @setscalar state.θwi = θwi = ρsn / snow.prop.ρw
     @setscalar state.ρsn = ρsn
     dsn = getscalar(state.swe) / θwi
@@ -112,7 +111,7 @@ function CryoGrid.diagnosticstep!(
     # 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)
+    T, θw, C = HeatConduction.enthalpyinv(heat.freezecurve, f_hc, getscalar(state.H), θwi_cap, heat.prop.L)
     # do not allow temperature to exceed 0°C
     @. state.T = min(T, zero(T))
     @. state.θw = θw
@@ -120,6 +119,7 @@ function CryoGrid.diagnosticstep!(
     # compute thermal conductivity
     HeatConduction.thermalconductivity!(snow, heat, state)
     @. state.k = state.kc
+    return nothing
 end
 # snowfall upper boundary
 function CryoGrid.interact!(
@@ -147,15 +147,15 @@ function CryoGrid.prognosticstep!(
     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]
+        jH_upper = state.jH[1] # [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]
+        dmelt = max(ddf*(T_ub-Tref), zero(eltype(state.dswe))) # [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))
+        # set upper heat flux to zero if dmelt > 0;
+        # this is due to the energy being (theoretically) "consumed" to melt the snow
+        state.jH[1] *= 1 - (dmelt > zero(dmelt))
     end
     return nothing
 end
@@ -212,10 +212,11 @@ function CryoGrid.diagnosticstep!(
 )
     smb, heat = procs
     ρw = snow.prop.ρw
+    HeatConduction.resetfluxes!(snow, heat, state)
     new_swe = swe(snow, smb, state)
     new_ρsn = snowdensity(snow, smb, state)
     new_dsn = new_swe*ρw/new_ρsn
-    ρw = heat.prop.ρw
+    ρw = heat.prop.consts.ρw
     if new_dsn > threshold(snow)
         # if new snow depth is above threshold, set state variables
         @setscalar state.swe = new_swe
@@ -238,14 +239,17 @@ function CryoGrid.diagnosticstep!(
         @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)
+function CryoGrid.prognosticstep!(snow::BulkSnowpack, ps::Coupled2{<:SnowMassBalance,<:Heat{FreeWater,Enthalpy}}, state)
+    smb, heat = ps
+    prognosticstep!(snow, smb, state)
     dsn = getscalar(state.dsn)
     if dsn < snow.para.thresh
-        # set energy flux to zero if there is no snow
+        # set divergence to zero if there is no snow
         @. state.dH = zero(eltype(state.H))
+    else
+        # otherwise call prognosticstep! for heat
+        prognosticstep!(snow, heat, state)
     end
 end
+``
\ No newline at end of file

src/Physics/Soils/Soils.jl

@@ -13,10 +13,13 @@ using IfElse
 using Interpolations: Interpolations
 using IntervalSets
 using ForwardDiff
+using FreezeCurves
 using ModelParameters
+using Setfield
 using Unitful
 
 import CryoGrid
+import CryoGrid.InputOutput
 import CryoGrid.Physics
 import CryoGrid.Physics.HeatConduction
 
@@ -25,6 +28,9 @@ import Flatten: @flattenable, flattenable
 export Soil, SoilParameterization, CharacteristicFractions, SoilProfile
 export soilparameters, soilcomponent, porosity, mineral, organic
 
+# from FreezeCurves
+export SFCC, DallAmico, DallAmicoSalt, Westermann, McKenzie, VanGenuchten
+
 const Enthalpy = HeatConduction.Enthalpy
 const Temperature = HeatConduction.Temperature
 
@@ -48,12 +54,19 @@ abstract type SoilParameterization end
 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 = 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
+    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
+end
+function 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)
+    )
 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...)
@@ -66,12 +79,12 @@ soilcomponent(::Val{:θo}, χ, ϕ, θ, ω) = (1-χ)*(1-ϕ)*ω
 """
 Soil thermal properties.
 """
-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)
+Base.@kwdef struct SoilThermalProperties{Tko,Tkm,Tco,Tcm}
+    ko::Tko = 0.25u"W/m/K" # organic [Hillel(1982)]
+    km::Tkm = 3.8u"W/m/K" # mineral [Hillel(1982)]
+    co::Tco = 2.5e6u"J/K/m^3" # heat capacity organic
+    cm::Tcm = 2.0e6u"J/K/m^3" # heat capacity mineral
+end
 """
 Basic Soil layer.
 """
@@ -80,21 +93,22 @@ Basic Soil layer.
     prop::TProp = SoilThermalProperties()
     sp::TSp = nothing # user-defined specialization
 end
+InputOutput.parameterize(soil::Soil; fields...) = Soil(para=InputOutput.parameterize(soil.para; fields...), prop=soil.prop, sp=soil.sp)
 HeatConduction.thermalproperties(soil::Soil) = soil.prop
 # SoilComposition trait impl
 SoilComposition(soil::Soil) = SoilComposition(typeof(soil))
 SoilComposition(::Type{<:Soil}) = Heterogeneous()
 SoilComposition(::Type{<:Soil{<:HomogeneousCharacteristicFractions}}) = Homogeneous()
 # 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)
+Physics.waterice(soil::Soil, state) = waterice(SoilComposition(soil), soil, state)
+porosity(soil::Soil, state) = porosity(SoilComposition(soil), soil, state)
+mineral(soil::Soil, state) = mineral(SoilComposition(soil), soil, state)
+organic(soil::Soil, state) = organic(SoilComposition(soil), soil, state)
 ## Homogeneous soils
-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)
+Physics.waterice(::Homogeneous, soil::Soil{<:CharacteristicFractions}, state=nothing) = soilcomponent(Val{:θwi}(), soil.para)
+porosity(::Homogeneous, soil::Soil{<:CharacteristicFractions}, state=nothing) = soilcomponent(Val{:θp}(), soil.para)
+mineral(::Homogeneous, soil::Soil{<:CharacteristicFractions}, state=nothing) = soilcomponent(Val{:θm}(), soil.para)
+organic(::Homogeneous, soil::Soil{<:CharacteristicFractions}, state=nothing) = soilcomponent(Val{:θo}(), soil.para)
 ## Heterogeneous soils
 Physics.waterice(::Heterogeneous, soil::Soil, state) = state.θwi   
 porosity(::Heterogeneous, soil::Soil, state) = state.θp
@@ -153,10 +167,6 @@ 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/sfcc.jl

-"""
-Abstract representation of a soil freeze characteristic curve (SFCC) function.
-Subtypes should be callable structs that implement the freeze curve and contain
-any necessary additional constants or configuration options. User-specified parameters
-can either be supplied in the struct or declared as model parameters via the `variables`
-method.
-"""
-abstract type SFCCFunction end
-"""
-Abstract type for SFCC H <--> T solvers.
-"""
-abstract type SFCCSolver end
-"""
-    SFCC{F,S} <: FreezeCurve
-
-Generic representation of the soil freeze characteristic curve. The shape and parameters
-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.
-"""
-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
-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 CryoGrid.initialcondition!(soil::Soil, heat::Heat, sfcc::SFCC, state)
-    HeatConduction.freezethaw!(soil, heat, state)
-end
-
-"""
-    sfccargs(f::SFCCFunction, soil::Soil, heat::Heat, state)
-
-Retrieves a tuple of values corresponding to each parameter declared by SFCCFunction `f` given the
-Soil layer, Heat process, and model state. The order of parameters *must match* the argument order
-of the freeze curve function `f`.
-"""
-sfccargs(::SFCCFunction, ::Soil, ::Heat, state) = ()
-# Fallback implementation of variables for SFCCFunction
-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,Θ,G,Tvg<:VanGenuchten} <: SFCCFunction
-    Tₘ::T = Param(0.0, units=u"°C")
-    θres::Θ = Param(0.0, domain=0..1)
-    g::G = 9.80665u"m/s^2" # acceleration due to gravity
-    swrc::Tvg = VanGenuchten()
-end
-sfccargs(f::DallAmico, soil::Soil, heat::Heat, state) = (
-    porosity(soil, state), # θ saturated = porosity
-    waterice(soil, heat, state), # total water content
-    heat.prop.Lf, # specific latent heat of fusion, L
-    f.θres,
-    f.Tₘ,
-    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.swrc.α, n=f.swrc.n)
-    let θsat = max(θtot, θsat),
-        g = f.g,
-        Tₘ = normalize_temperature(Tₘ),
-        ψ₀ = f.swrc(inv, θtot, θres, θsat, α, n),
-        Tstar = Tₘ + g*Tₘ/Lf*ψ₀,
-        T = normalize_temperature(T),
-        ψ = ψ(T, Tstar, ψ₀, Lf, g);
-        f.swrc(ψ, θres, θsat, α, n)
-    end
-end
-
-"""
-    McKenzie <: SFCCFunction
-
-McKenzie JM, Voss CI, Siegel DI, 2007. Groundwater flow with energy transport and water-ice phase change:
-    numerical simulations, benchmarks, and application to freezing in peat bogs. Advances in Water Resources,
-    30(4): 966–983. DOI: 10.1016/j.advwatres.2006.08.008.
-"""
-Base.@kwdef struct McKenzie{T,Θ,Γ} <: SFCCFunction
-    Tₘ::T = Param(0.0, units=u"°C")
-    θ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
-    waterice(soil, heat, state), # total water content
-    f.θres,
-    f.Tₘ,
-    f.γ,
-)
-function (f::McKenzie)(T, θsat, θtot, θres=f.θres, Tₘ=f.Tₘ, γ=f.γ)
-    let T = normalize_temperature(T),
-        Tₘ = normalize_temperature(Tₘ),
-        θsat = max(θtot, θsat);
-        IfElse.ifelse(T<=Tₘ, θres + (θsat-θres)*exp(-((T-Tₘ)/γ)^2), θtot)
-    end
-end
-
-"""
-    Westermann <: SFCCFunction
-
-Westermann, S., Boike, J., Langer, M., Schuler, T. V., and Etzelmüller, B.: Modeling the impact of
-    wintertime rain events on the thermal regime of permafrost, The Cryosphere, 5, 945–959,
-    https://doi.org/10.5194/tc-5-945-2011, 2011. 
-"""
-Base.@kwdef struct Westermann{T,Θ,Δ} <: SFCCFunction
-    Tₘ::T = Param(0.0, units=u"°C")
-    θ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
-    waterice(soil, heat, state), # total water content
-    f.θres,
-    f.Tₘ,
-    f.δ,
-)
-function (f::Westermann)(T,θsat,θtot,θres=f.θres,Tₘ=f.Tₘ,δ=f.δ)
-    let T = normalize_temperature(T),
-        Tₘ = normalize_temperature(Tₘ),
-        θsat = max(θtot, θsat);
-        IfElse.ifelse(T<=Tₘ, θres - (θsat-θres)*(δ/(T-Tₘ-δ)), θtot)
-    end
-end
-struct SFCCTable{F,I} <: SFCCFunction
-    f::F
-    f_tab::I
-end
-(f::SFCCTable)(args...) = f.f_tab(args...)
-"""
-    Tabulated(f::SFCCFunction, args...)
-
-Produces an `SFCCTable` function which is a tabulation of `f`.
-"""
-Numerics.Tabulated(f::SFCCFunction, args...; kwargs...) = SFCCTable(f, Numerics.tabulate(f, args...; kwargs...))
-
-include("sfcc_solvers.jl")

src/Physics/Soils/sfcc/sfcc_solvers.jl

-using NLsolve
-
-"""
-    SFCCNewtonSolver <: SFCCSolver
-
-Specialized implementation of Newton's method with backtracking line search for resolving
-the energy conservation law, H = TC + Lθ. Attempts to find the root of the corresponding
-temperature residual: ϵ = T - (H - Lθ(T)) / C(θ(T)) and uses backtracking to avoid
-jumping over the solution. This prevents convergence issues that arise due to
-discontinuities and strong non-linearity in most common soil freeze curves.
-"""
-Base.@kwdef struct SFCCNewtonSolver <: SFCCSolver
-    maxiter::Int = 100 # maximum number of iterations
-    abstol::Float64 = 1e-2 # absolute tolerance for convergence
-    reltol::Float64 = 1e-4 # relative tolerance for convergence
-    α₀::Float64 = 1.0 # initial step size multiplier
-    τ::Float64 = 0.7 # step size decay for backtracking
-end
-# 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, 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, 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, θ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, θw, itercount, T_converged)
-        end
-        # derivative of freeze curve
-        ∂θ∂T = ForwardDiff.derivative(T -> f(T, f_args...), T)
-        # derivative of residual by quotient rule;
-        # 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 - θ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, θ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, θw, itercount, T_converged)
-            end
-            α = α*τ # decrease step size by τ
-            T̂ = T - α*Tres # new guess for T
-            T̂res, θw, C = sfccresidual(soil, heat, f, f_args, f_hc, T̂, H)
-            inneritercount += 1
-        end
-        T = T̂ # update T
-        Tres = T̂res # update residual
-        itercount += 1
-    end
-    T_converged = true
-    return (;T, Tres, θw, itercount, T_converged)
-end
-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
-    # eventually passed to the `residual` function; this is less than ideal but
-    # probably shouldn't incur too much of a performance hit, just a few extra stack pointers!
-    f_args = sfccargs(f, soil, heat, state)
-    solver = sfcc.solver
-    @inbounds let T₀ = i > 1 ? state.T[i-1] : nothing,
-        H = state.H[i] |> Utils.adstrip, # enthalpy
-        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ᵢ, f_hc, H, T₀)
-        return T
-    end
-end
-"""
-    SFCCPreSolver{TCache} <: SFCCSolver
-
-A fast SFCC "solver" which pre-builds an interpolant for the freeze curve in terms of enthalpy, θ(H).
-Note that this solver is **only valid when all freeze curve parameters are held constant** and will
-produce incorrect results otherwise.
-"""
-@flattenable struct SFCCPreSolver{TCache} <: SFCCSolver
-    cache::TCache | false
-    Tmin::Float64 | false
-    errtol::Float64 | false
-    SFCCPreSolver(cache, Tmin, errtol) = new{typeof(cache)}(cache, Tmin, errtol)
-    """
-        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. `errtol` determines the permitted local error in the interpolant.
-    """
-    function SFCCPreSolver(;Tmin=-60.0, errtol=1e-4)
-        cache = SFCCPreSolverCache()
-        new{typeof(cache)}(cache, Tmin, errtol)
-    end
-end
-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;
-        # 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),typeof(dummy_f)}(dummy_f, dummy_∇f, dummy_f)
-    end
-end
-_interpgrad(f) = (args...) -> Interpolations.gradient(f, args...)
-function _build_interpolant(Hs, θs)
-    return Interpolations.extrapolate(
-        Interpolations.interpolate((Vector(Hs),), θs, Interpolations.Gridded(Interpolations.Linear())),
-        Interpolations.Flat()
-    )
-end
-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.θw .= sfcc.f.(state.T, args...)
-    heatcapacity!(soil, heat, state)
-    @. 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
-    # 2) soil properties are uniform in `soil`
-    let θsat = args[1],
-        θtot = args[2],
-        tail_args = args[3:end],
-        θm = mineral(soil, state),
-        θo = organic(soil, state),
-        θa = 1-θm-θo-θtot,
-        L = heat.prop.L,
-        Tmin = sfcc.solver.Tmin,
-        Tmax = 0.0,
-        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_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/Soils/sfcc/swrc.jl

-"""
-    SWRCFunction
-
-Base type for soil water retention curve SWRC function implementations.
-"""
-abstract type SWRCFunction end
-Base.inv(f::SWRCFunction) = (args...) -> f(inv, args...)
-"""
-    SWRC{F}
-
-Soil water retention curve with function type `F`.
-"""
-struct SWRC{F}
-    f::F # soil water retention curve function f: (ψ,...) -> θ
-    SWRC(f::F) where {F<:SWRCFunction} = new{F}(f)
-end
-"""
-    VanGenuchten <: SWRCFunction
-
-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.
-"""
-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

@@ -17,11 +17,44 @@
         (θw, θi, θa, θm, θo)
     end
 end
+# Default implementations of variables for SFCCFunction
+CryoGrid.variables(::Soil, ::Heat, f::SFCCFunction) = ()
+CryoGrid.variables(::Soil, ::Heat, s::SFCCSolver) = ()
+# Join the declared state variables of the SFCC function and the solver
+CryoGrid.variables(soil::Soil, heat::Heat, sfcc::SFCC) = tuplejoin(CryoGrid.variables(soil, heat, sfcc.f), CryoGrid.variables(soil, heat, sfcc.solver))
+"""
+    sfcckwargs(f::SFCCFunction, soil::Soil, heat::Heat, state, i)
 
+Builds a named tuple of values corresponding to each keyword arguments of the SFCCFunction `f`
+which should be set according to the layer/process properties or state. The default implementation
+sets only the total water content, θtot = θwi, and the saturated water content, θsat = θp.
+"""
+sfcckwargs(::SFCCFunction, soil::Soil, heat::Heat, state, i) = (
+    θtot = waterice(soil, heat, state, i), # total water content    
+    θsat = porosity(soil, state, i), # θ saturated = porosity
+)
 """
 Initial condition for heat conduction (all state configurations) on soil layer w/ SFCC.
 """
-CryoGrid.initialcondition!(soil::Soil, heat::Heat{<:SFCC}, state) = HeatConduction.initialcondition!(soil, heat, freezecurve(heat), state)
+function CryoGrid.initialcondition!(soil::Soil, heat::Heat{<:SFCC}, state)
+    fc = freezecurve(heat)
+    L = heat.prop.L
+    @inbounds for i in 1:length(state.T)
+        fc_kwargsᵢ = sfcckwargs(fc.f, soil, heat, state, i)
+        hc = partial(heatcapacity, liquidwater, soil, heat, state, i)
+        if i == 1
+            # TODO: this is currently only relevant for the pre-solver scheme and assumes that
+            # the total water content is uniform throughout the layer and does not change over time.
+            FreezeCurves.Solvers.initialize!(fc.solver, fc.f, hc; fc_kwargsᵢ...)
+        end
+        T = state.T[i]
+        θw, dθdT = ∇(T -> fc(T; fc_kwargsᵢ...), T)
+        state.θw[i] = θw
+        state.C[i] = hc(θw)
+        state.H[i] = enthalpy(state.T[i], state.C[i], L, state.θw[i])
+        state.dHdT[i] = HeatConduction.C_eff(T, state.C[i], L, dθdT, heat.prop.cw, heat.prop.ci)
+    end
+end
 """
 Initial condition for heat conduction (all state configurations) on soil layer w/ free water freeze curve.
 """
@@ -46,12 +79,11 @@ For heat conduction with temperature, we can simply evaluate the freeze curve to
 function HeatConduction.freezethaw!(soil::Soil, heat::Heat{<:SFCC,Temperature}, state)
     sfcc = freezecurve(heat)
     let L = heat.prop.L,
-        f = sfcc.f,
-        f_args = sfccargs(f,soil,heat,state);
+        f = sfcc.f;
         @inbounds @fastmath for i in 1:length(state.T)
             T = state.T[i]
-            f_argsᵢ = Utils.selectat(i, identity, f_args)
-            θw, dθdT = ∇(T -> f(T, f_argsᵢ...), T)
+            f_argsᵢ = sfcckwargs(f, soil, heat, state, i)
+            θ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, volumetricfractions(soil, heat, state, i)...)
@@ -60,42 +92,40 @@ function HeatConduction.freezethaw!(soil::Soil, heat::Heat{<:SFCC,Temperature},
         end
     end
 end
-function HeatConduction.freezethaw!(soil::Soil, heat::Heat{<:SFCC{F,SFCCNewtonSolver},Enthalpy}, state) where {F}
+function HeatConduction.freezethaw!(soil::Soil, heat::Heat{<:SFCC,Enthalpy}, state)
     sfcc = freezecurve(heat)
-    f_args = sfccargs(sfcc.f, soil, heat, state)
     @inbounds for i in 1:length(state.H)
-        let f_argsᵢ = Utils.selectat(i, identity, f_args),
-            f = sfcc.f;
-            # Evaluate inverse enthalpy function
-            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 θw, C, and T here with the tracked variables.
-            # Note that this results in one additional freeze curve function evaluation.
-            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 = 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
+        @inbounds let H = state.H[i], # enthalpy
+            L = heat.prop.L,
+            cw = heat.prop.cw,
+            ci = heat.prop.ci,
+            θwi = waterice(soil, heat, state, i), # total water content
+            T₀ = i > 1 ? state.T[i-1] : FreezeCurves.freewater(H, θwi, L), # initial guess for T
+            hc = partial(heatcapacity, liquidwater, soil, heat, state, i),
+            f = sfcc.f,
+            f_kwargsᵢ = sfcckwargs(f, soil, heat, state, i),
+            obj = FreezeCurves.SFCCInverseEnthalpyObjective(f, f_kwargsᵢ, hc, L, H);
+            res = FreezeCurves.sfccsolve(obj, sfcc.solver, T₀, Val{true}())
+            state.T[i] = res.T
+            state.θw[i] = res.θw
+            state.C[i] = res.C
+            state.dHdT[i] = HeatConduction.C_eff(state.T[i], state.C[i], L, res.dθwdT, cw, ci)
         end
     end
 end
-function HeatConduction.freezethaw!(soil::Soil{<:HomogeneousCharacteristicFractions}, heat::Heat{<:SFCC{F,<:SFCCPreSolver},Enthalpy}, state) where {F}
-    solver = freezecurve(heat).solver
-    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.θ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)
+function HeatConduction.enthalpyinv(soil::Soil, heat::Heat{<:SFCC,Enthalpy}, state, i)
+    sfcc = freezecurve(heat)
+    @inbounds let H = state.H[i], # enthalpy
+        L = heat.prop.L, # latent heat of fusion of water
+        θwi = waterice(soil, heat, state, i), # total water content
+        T₀ = i > 1 ? state.T[i-1] : freewater(H, θwi, L),
+        hc = partial(heatcapacity, liquidwater, soil, heat, state, i),
+        f = sfcc.f,
+        f_kwargsᵢ = sfcckwargs(f, soil, heat, state, i),
+        obj = FreezeCurves.SFCCInverseEnthalpyObjective(f, f_kwargsᵢ, hc, L, H);
+        T_sol = FreezeCurves.sfccsolve(obj, sfcc.solver, T₀, Val{false}())
+        return T_sol
     end
 end
+# Parameterization
+InputOutput.parameterize(solver::FreezeCurves.SFCCSolver; fields...) = solver