From 29035b9ea273fd48996de7e0958c47f88a19bbed Mon Sep 17 00:00:00 2001
From: Brian Groenke <brian.groenke@awi.de>
Date: Thu, 16 Dec 2021 11:28:45 +0100
Subject: [PATCH] Rework soil thermal properties to allow both scalar and
 vector forms

---
 Project.toml                                 |  2 +-
 src/Layers/Layers.jl                         |  3 +-
 src/Layers/soil.jl                           | 17 ++---
 src/Physics/HeatConduction/HeatConduction.jl | 57 +++++++++++++++-
 src/Physics/HeatConduction/heat.jl           | 36 ++++------
 src/Physics/HeatConduction/soil/sfcc.jl      | 23 ++++---
 src/Physics/HeatConduction/soil/soilheat.jl  | 71 ++++++++++++++++----
 src/Strat/stratigraphy.jl                    |  2 +-
 src/Utils/Utils.jl                           |  2 +
 test/Physics/HeatConduction/sfcc_tests.jl    | 20 +++---
 10 files changed, 158 insertions(+), 75 deletions(-)

diff --git a/Project.toml b/Project.toml
index 544ab68a..517d8d6c 100755
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "CryoGrid"
 uuid = "a535b82e-5f3d-4d97-8b0b-d6483f5bebd5"
 authors = ["Brian Groenke <brian.groenke@awi.de>", "Moritz Langer <moritz.langer@awi.de>"]
-version = "0.5.7"
+version = "0.5.8"
 
 [deps]
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
diff --git a/src/Layers/Layers.jl b/src/Layers/Layers.jl
index 227bfb1b..c3defedb 100644
--- a/src/Layers/Layers.jl
+++ b/src/Layers/Layers.jl
@@ -12,7 +12,8 @@ using ModelParameters
 using Parameters
 using Unitful
 
-export Soil, SoilParameterization, SoilCharacteristicFractions, SoilProfile, SoilType, Sand, Silt, Clay, soilparameters, soilcomp, porosity
+export Soil, SoilParameterization, SoilCharacteristicFractions, SoilProfile, SoilType, Sand, Silt, Clay
+export soilparameters, soilcomp, porosity
 include("soil.jl")
 
 end
\ No newline at end of file
diff --git a/src/Layers/soil.jl b/src/Layers/soil.jl
index 57e597dd..a7e42799 100644
--- a/src/Layers/soil.jl
+++ b/src/Layers/soil.jl
@@ -64,17 +64,8 @@ Basic Soil layer.
     hc::SoilHCParams = SoilHCParams()
     sp::S = nothing # user-defined specialization
 end
-# Methods
+# Volumetric content methods
+totalwater(soil::Soil{T,<:SoilCharacteristicFractions}) where T = soilcomp(Val{:θw}(), soil.para)
 porosity(soil::Soil{T,<:SoilCharacteristicFractions}) where T = soilcomp(Val{:θp}(), soil.para)
-variables(::Soil) = (
-    Diagnostic(:θp, Float"1", OnGrid(Cells)), # porosity
-    Diagnostic(:θw, Float"1", OnGrid(Cells)), # total water content (xice + saturated pore space)
-    Diagnostic(:θm, Float"1", OnGrid(Cells)), # mineral content
-    Diagnostic(:θo, Float"1", OnGrid(Cells)), # organic content
-)
-function initialcondition!(soil::Soil{T,<:SoilCharacteristicFractions}, state) where T
-    state.θp .= soilcomp(Val{:θp}(), soil.para)
-    state.θw .= soilcomp(Val{:θw}(), soil.para)
-    state.θm .= soilcomp(Val{:θm}(), soil.para)
-    state.θo .= soilcomp(Val{:θo}(), soil.para)
-end
+mineral(soil::Soil{T,<:SoilCharacteristicFractions}) where T = soilcomp(Val{:θm}(), soil.para)
+organic(soil::Soil{T,<:SoilCharacteristicFractions}) where T = soilcomp(Val{:θo}(), soil.para)
diff --git a/src/Physics/HeatConduction/HeatConduction.jl b/src/Physics/HeatConduction/HeatConduction.jl
index 067fd43a..b688de7a 100644
--- a/src/Physics/HeatConduction/HeatConduction.jl
+++ b/src/Physics/HeatConduction/HeatConduction.jl
@@ -2,13 +2,13 @@ module HeatConduction
 
 import CryoGrid: SubSurfaceProcess, BoundaryStyle, Dirichlet, Neumann, BoundaryProcess, Layer, Top, Bottom, SubSurface, Callback
 import CryoGrid: diagnosticstep!, prognosticstep!, interact!, initialcondition!, boundaryflux, boundaryvalue, variables, callbacks, criterion, affect!
+import CryoGrid.Layers: Soil, totalwater, porosity, mineral, organic
 
 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, porosity
 using CryoGrid.Utils
 
 using Base: @propagate_inbounds
@@ -24,6 +24,7 @@ import Flatten: @flattenable, flattenable
 
 export Heat, TemperatureProfile
 export FreeWater, FreezeCurve, freezecurve
+export enthalpy, heatcapacity, heatcapacity!, thermalconductivity, thermalconductivity!
 
 abstract type FreezeCurve end
 struct FreeWater <: FreezeCurve end
@@ -46,8 +47,60 @@ end
 Heat(var::Union{Symbol,Tuple{Vararg{Symbol}}}; kwargs...) = Heat(Val{var}(); kwargs...)
 Heat(::Val{:H}; kwargs...) = Heat(;sp=Enthalpy(), 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) = ()
+# Fallback (error) implementation for freeze curve
+(fc::FreezeCurve)(sub::SubSurface, heat::Heat, state) = error("freeze curve $(typeof(fc)) not implemented for $(typeof(heat)) on layer $(typeof(sub))")
+# 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*θ
+"""
+    totalwater(sub::SubSurface, heat::Heat, state)
+    totalwater(sub::SubSurface, heat::Heat, 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, heat::Heat, state) = totalwater(sub)
+@inline totalwater(sub::SubSurface, heat::Heat, state, i) = Utils.getscalar(totalwater(sub, heat, state), i)
+"""
+    heatcapacity(sub::SubSurface, heat::Heat, state, i)
+
+Computes the heat capacity at grid cell `i` for the given layer from the current state.
+"""
+heatcapacity(sub::SubSurface, heat::Heat, state, i) = error("heatcapacity not defined for $(typeof(heat)) on $(typeof(sub))")
+"""
+    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(sub::SubSurface, heat::Heat, state, i)
+
+Computes the thrmal conductivity at grid cell `i` for the given layer from the current state.
+"""
+thermalconductivity(sub::SubSurface, heat::Heat, state, i) = error("thermalconductivity not defined for $(typeof(heat)) on $(typeof(sub))")
+"""
+    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
 
-export enthalpy, heatcapacity, heatcapacity!, thermalconductivity, thermalconductivity!
 export heatconduction!
 include("heat.jl")
 export ConstantTemp, GeothermalHeatFlux, TemperatureGradient, NFactor, Damping
diff --git a/src/Physics/HeatConduction/heat.jl b/src/Physics/HeatConduction/heat.jl
index ab1f8b46..9949f11c 100644
--- a/src/Physics/HeatConduction/heat.jl
+++ b/src/Physics/HeatConduction/heat.jl
@@ -31,29 +31,19 @@ function heatconduction!(∂H,T,ΔT,k,Δk)
     end
     return nothing
 end
-
-# 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, C, L, θ) = T*C + L*θ
-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))")
 """
 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(layer::SubSurface, heat::Heat) = (
-    variables(layer)..., # layer variables
+variables(sub::SubSurface, heat::Heat) = (
+    variables(sub)..., # layer variables
     variables(heat)...,  # heat variables
-    variables(layer, heat, freezecurve(heat))..., # freeze curve variables
+    variables(sub, heat, freezecurve(heat))..., # freeze curve variables
 )
 """
 Variable definitions for heat conduction (enthalpy).
 """
-variables(heat::Heat{<:FreezeCurve,Enthalpy}) = (
+variables(::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)),
@@ -65,7 +55,7 @@ variables(heat::Heat{<:FreezeCurve,Enthalpy}) = (
 """
 Variable definitions for heat conduction (temperature).
 """
-variables(heat::Heat{<:FreezeCurve,Temperature}) = (
+variables(::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)),
@@ -76,15 +66,15 @@ variables(heat::Heat{<:FreezeCurve,Temperature}) = (
     Diagnostic(:θl, Float"1", OnGrid(Cells)),
 )
 """ Diagonstic step for heat conduction (all state configurations) on any subsurface layer. """
-function diagnosticstep!(layer::SubSurface, heat::Heat, state)
+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)
+    fc!(sub, heat, state)
     # Update thermal conductivity
-    thermalconductivity!(layer, heat, state)
+    thermalconductivity!(sub, heat, state)
     # Harmonic mean of inner conductivities
     @inbounds let k = (@view state.k[2:end-1]),
         Δk = Δ(state.grids.k);
@@ -97,14 +87,14 @@ 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)
+    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)
     Δk = Δ(state.grids.k) # cell sizes
     ΔT = Δ(state.grids.T)
     # Diffusion on non-boundary cells
-    heatconduction!(state.dH,state.T,ΔT,state.k,Δk)
+    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
@@ -184,7 +174,7 @@ Implementation of "free water" freeze curve for any subsurface layer. Assumes th
 '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 (fc::FreeWater)(sub::SubSurface, heat::Heat{FreeWater,Enthalpy}, state)
     @inline function enthalpyinv(H, C, L, θtot)
         let θtot = max(1.0e-8, θtot),
             Lθ = L*θtot,
@@ -202,11 +192,11 @@ total water content (θw), and liquid water content (θl).
         end
     end
     @inbounds for i in 1:length(state.H)
-        let θw = state.θw[i],
+        let θw = totalwater(sub, heat, state, i),
             H = state.H[i],
             L = heat.L;
             state.θl[i] = θw*freezethaw(H, L, θw)
-            state.C[i] = heatcapacity(layer, heat, state, i) # update heat capacity, C
+            state.C[i] = heatcapacity(sub, heat, state, i) # update heat capacity, C
             state.T[i] = enthalpyinv(H, state.C[i], L, θw)
         end
     end
diff --git a/src/Physics/HeatConduction/soil/sfcc.jl b/src/Physics/HeatConduction/soil/sfcc.jl
index 530905bf..75fdc97c 100644
--- a/src/Physics/HeatConduction/soil/sfcc.jl
+++ b/src/Physics/HeatConduction/soil/sfcc.jl
@@ -57,8 +57,11 @@ function (sfcc::SFCC)(soil::Soil, heat::Heat{<:SFCC,Temperature}, state)
         f_args = tuplejoin((state.T,),sfccparams(f,soil,heat,state));
         for i in 1:length(state.T)
             f_argsáµ¢ = Utils.selectat(i, identity, f_args)
+            θw = totalwater(soil, heat, state, i)
+            θm = mineral(soil, heat, i)
+            θo = organic(soil, heat, i)
             state.θl[i] = f(f_argsᵢ...)
-            state.C[i] = heatcapacity(soil, state.θw[i], state.θl[i], state.θm[i], state.θo[i])
+            state.C[i] = heatcapacity(soil, θw, state.θl[i], θm, θo)
             state.H[i] = enthalpy(state.T[i], state.C[i], L, state.θl[i])
             state.Ceff[i] = L*∇f(f_argsᵢ) + state.C[i]
         end
@@ -93,8 +96,8 @@ end
 sfccparams(f::DallAmico, soil::Soil, heat::Heat, state) = (
     f.Tₘ,
     f.θres,
-    state.θp, # θ saturated = porosity
-    state.θw, # total water content
+    porosity(soil, heat, state), # θ saturated = porosity
+    totalwater(soil, heat, state), # total water content
     heat.L, # specific latent heat of fusion, L
     f.α,
     f.n,
@@ -129,8 +132,8 @@ end
 sfccparams(f::McKenzie, soil::Soil, heat::Heat, state) = (
     f.Tₘ,
     f.θres,
-    state.θp, # θ saturated = porosity
-    state.θw, # total water content
+    porosity(soil, heat, state), # θ saturated = porosity
+    totalwater(soil, heat, state), # total water content
     f.γ,
 )
 function (f::McKenzie)(T,Tₘ,θres,θsat,θtot,γ)
@@ -154,8 +157,8 @@ end
 sfccparams(f::Westermann, soil::Soil, heat::Heat, state) = (
     f.Tₘ,
     f.θres,
-    state.θp, # θ saturated = porosity
-    state.θw, # total water content
+    porosity(soil, heat, state), # θ saturated = porosity
+    totalwater(soil, heat, state), # total water content
     f.δ,
 )
 function (f::Westermann)(T,Tₘ,θres,θsat,θtot,δ)
@@ -225,9 +228,9 @@ function (s::SFCCNewtonSolver)(soil::Soil, heat::Heat{<:SFCC,Enthalpy}, state, f
             H = state.H[i] |> Utils.adstrip, # enthalpy
             C = state.C[i] |> Utils.adstrip, # heat capacity
             θl = state.θl[i] |> Utils.adstrip, # liquid water content
-            θw = state.θw[i] |> Utils.adstrip, # total water content
-            θm = state.θm[i] |> Utils.adstrip, # mineral content
-            θo = state.θo[i] |> Utils.adstrip, # organic content
+            θw = totalwater(soil, heat, state, i) |> Utils.adstrip, # total water content
+            θm = mineral(soil, heat, state, i) |> Utils.adstrip, # mineral content
+            θo = organic(soil, heat, state, i) |> Utils.adstrip, # organic content
             L = heat.L, # specific latent heat of fusion
             cw = soil.hc.cw, # heat capacity of liquid water
             α₀ = s.α₀,
diff --git a/src/Physics/HeatConduction/soil/soilheat.jl b/src/Physics/HeatConduction/soil/soilheat.jl
index 8d200fab..ed11bac9 100644
--- a/src/Physics/HeatConduction/soil/soilheat.jl
+++ b/src/Physics/HeatConduction/soil/soilheat.jl
@@ -1,4 +1,43 @@
-# Thermal properties
+# === Thermal properties ===
+# We use methods with optioinal index arguments `i` to allow for implementations both
+# where these variables are treated as constants and as state variables.
+# In the latter case, specializations should override only the index-free form
+# and return a state vector instead of a scalar. The `getscalar` function will
+# handle both the scalar and vector case!
+"""
+    mineral(soil::Soil, ::Heat, state)
+    mineral(soil::Soil, ::Heat, state, i)
+
+Retrieves the mineral content for the given layer at grid cell `i`, if provided.
+Defaults to using the scalar mineral content defined on `soil`.
+"""
+@inline mineral(soil::Soil, ::Heat, state) = mineral(soil)
+@inline mineral(soil::Soil, heat::Heat, state, i) = Utils.getscalar(mineral(soil, heat, state), i)
+"""
+    organic(soil::Soil, ::Heat, state)
+    organic(soil::Soil, ::Heat, state, i)
+
+Retrieves the organic content for the given layer at grid cell `i`, if provided.
+Defaults to using the scalar organic content defined on `soil`.
+"""
+@inline organic(soil::Soil, ::Heat, state) = organic(soil)
+@inline organic(soil::Soil, heat::Heat, state, i) = Utils.getscalar(organic(soil, heat, state), i)
+"""
+    porosity(soil::Soil, ::Heat, state)
+    porosity(soil::Soil, ::Heat, state, i)
+
+Retrieves the porosity for the given layer at grid cell `i`, if provided.
+Defaults to using the scalar porosity defined on `soil`.
+"""
+@inline porosity(soil::Soil, ::Heat, state) = porosity(soil)
+@inline porosity(soil::Soil, heat::Heat, state, i) = Utils.getscalar(porosity(soil, heat, state), i)
+@inline heatcapacity(soil::Soil, heat::Heat, state, i) = heatcapacity(
+    soil,
+    totalwater(soil, heat, state, i),
+    state.θl[i],
+    mineral(soil, heat, state, i),
+    organic(soil, heat, state, i),
+)
 @inline function heatcapacity(soil::Soil, totalwater, liquidwater, mineral, organic)
     @unpack cw,co,cm,ca,ci = soil.hc
     let air = 1.0 - totalwater - mineral - organic,
@@ -7,6 +46,13 @@
         liq*cw + ice*ci + mineral*cm + organic*co + air*ca
     end
 end
+@inline thermalconductivity(soil::Soil, heat::Heat, state, i) = thermalconductivity(
+    soil,
+    totalwater(soil, heat, state, i),
+    state.θl[i],
+    mineral(soil, heat, state, i),
+    organic(soil, heat, state, i),
+)
 @inline function thermalconductivity(soil::Soil, totalwater, liquidwater, mineral, organic)
     @unpack kw,ko,km,ka,ki = soil.tc
     let air = 1.0 - totalwater - mineral - organic,
@@ -15,14 +61,6 @@ end
         (liq*kw^0.5 + ice*ki^0.5 + mineral*km^0.5 + organic*ko^0.5 + air*ka^0.5)^2
     end
 end
-@inline @propagate_inbounds heatcapacity(soil::Soil, ::Heat, state, i) = heatcapacity(soil, state.θw[i], state.θl[i], state.θm[i], state.θo[i])
-@inline @propagate_inbounds function heatcapacity!(soil::Soil, ::Heat, state)
-    @. state.C = heatcapacity(soil, state.θw, state.θl, state.θm, state.θo)
-end
-@inline @propagate_inbounds thermalconductivity(soil::Soil, ::Heat, state, i) = thermalconductivity(soil, state.θw[i], state.θl[i], state.θm[i], state.θo[i])
-@inline @propagate_inbounds function thermalconductivity!(soil::Soil, ::Heat, state)
-    @. state.kc = thermalconductivity(soil, state.θw, state.θl, state.θm, state.θo)
-end
 
 # SFCC
 include("sfcc.jl")
@@ -35,8 +73,10 @@ function initialcondition!(soil::Soil, heat::Heat, state)
     L = heat.L
     fc! = freezecurve(heat)
     fc!(soil, heat, state)
-    heatcapacity!(soil, heat, state)
-    @. state.H = enthalpy(state.T, state.C, L, state.θl)
+    @inbounds for i in 1:length(state.T)
+        state.C[i] = heatcapacity(soil, heat, state, i)
+        state.H[i] = enthalpy(state.T[i], state.C[i], L, state.θl[i])
+    end
 end
 """
 Initial condition for heat conduction (all state configurations) with free water freeze curve on soil layer.
@@ -45,7 +85,10 @@ function initialcondition!(soil::Soil, heat::Heat{FreeWater}, state)
     initialcondition!(soil, state)
     L = heat.L
     # initialize liquid water content based on temperature
-    @. state.θl = ifelse.(state.T > 0.0, state.θw, 0.0)
-    heatcapacity!(soil, heat, state)
-    @. state.H = enthalpy(state.T, state.C, L, state.θl)
+    @inbounds for i in 1:length(state.T)
+        θw = totalwater(soil, heat, state, i)
+        state.θl[i] = ifelse(state.T[i] > 0.0, θw, 0.0)
+        state.C[i] = heatcapacity(soil, heat, state, i)
+        state.H[i] = enthalpy(state.T[i], state.C[i], L, state.θl[i])
+    end
 end
diff --git a/src/Strat/stratigraphy.jl b/src/Strat/stratigraphy.jl
index d296a561..25471835 100644
--- a/src/Strat/stratigraphy.jl
+++ b/src/Strat/stratigraphy.jl
@@ -21,7 +21,7 @@ Base.show(io::IO, node::StratComponent{L,P,name}) where {L,P,name} = print(io, "
 # Constructors for stratigraphy nodes
 top(bcs::BoundaryProcess...) = StratComponent(:top, Top(), CoupledProcesses(bcs...))
 bottom(bcs::BoundaryProcess...) = StratComponent(:bottom, Bottom(), CoupledProcesses(bcs...))
-subsurface(name::Symbol, layer::SubSurface, processes::SubSurfaceProcess...) = StratComponent(name, layer, CoupledProcesses(processes...))
+subsurface(name::Symbol, sub::SubSurface, processes::SubSurfaceProcess...) = StratComponent(name, sub, CoupledProcesses(processes...))
 
 """
     Stratigraphy{N,TComponents,Q}
diff --git a/src/Utils/Utils.jl b/src/Utils/Utils.jl
index 49fcf409..05df88ac 100644
--- a/src/Utils/Utils.jl
+++ b/src/Utils/Utils.jl
@@ -83,7 +83,9 @@ Base.iterate(p::Params, state) = structiterate(p,state)
 @inline tuplejoin(x, y, z...) = (x..., tuplejoin(y, z...)...)
 
 getscalar(x::Number) = x
+getscalar(x::Number, i) = x
 getscalar(a::AbstractArray) = a[1]
+getscalar(a::AbstractArray, i) = a[i]
 
 """
     convert_tspan(tspan::Tuple{DateTime,DateTime})
diff --git a/test/Physics/HeatConduction/sfcc_tests.jl b/test/Physics/HeatConduction/sfcc_tests.jl
index 2916e4bc..d1c6748a 100644
--- a/test/Physics/HeatConduction/sfcc_tests.jl
+++ b/test/Physics/HeatConduction/sfcc_tests.jl
@@ -36,7 +36,7 @@ using ComponentArrays
                     C = heatcapacity.(soil,θw,θl,θm,θo);
                    enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
                 end
-                state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
+                state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,)
                 sfcc(soil, heat, state)
                 @test all(abs.(T.-(H .- L.*θl)./C) .<= tol)
             end
@@ -51,7 +51,7 @@ using ComponentArrays
                     C = heatcapacity.(soil,θw,θl,θm,θo);
                    enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
                 end
-                state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
+                state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,)
                 sfcc(soil, heat, state)
                 @test all(abs.(T.-(H .- L.*θl)./C) .<= tol)
             end
@@ -66,7 +66,7 @@ using ComponentArrays
                     C = heatcapacity.(soil,θw,θl,θm,θo);
                    enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
                 end
-                state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
+                state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,)
                 sfcc(soil, heat, state)
                 @test all(abs.(T.-(H .- L.*θl)./C) .<= tol)
             end
@@ -100,7 +100,7 @@ using ComponentArrays
                     C = heatcapacity.(soil,θw,θl,θm,θo);
                    enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
                 end
-                state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
+                state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,)
                 sfcc(soil, heat, state)
                 @test all(abs.(T.-(H .- L.*θl)./C) .<= tol)
             end
@@ -113,7 +113,7 @@ using ComponentArrays
                     C = heatcapacity.(soil,θw,θl,θm,θo);
                    enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
                 end
-                state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
+                state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,)
                 sfcc(soil, heat, state)
                 @test all(abs.(T.-(H .- L.*θl)./C) .<= tol)
             end
@@ -126,7 +126,7 @@ using ComponentArrays
                     C = heatcapacity.(soil,θw,θl,θm,θo);
                    enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
                 end
-                state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
+                state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,)
                 sfcc(soil, heat, state)
                 @test all(abs.(T.-(H .- L.*θl)./C) .<= tol)
             end
@@ -163,7 +163,7 @@ using ComponentArrays
                     C = heatcapacity.(soil,θw,θl,θm,θo);
                    enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
                 end
-                state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
+                state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,)
                 @inferred sfcc(soil, heat, state)
                 @test all(abs.(T.-(H .- L.*θl)./C) .<= tol)
             end
@@ -176,7 +176,7 @@ using ComponentArrays
                     C = heatcapacity.(soil,θw,θl,θm,θo);
                    enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
                 end
-                state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
+                state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,)
                 @inferred sfcc(soil, heat, state)
                 @test all(abs.(T.-(H .- L.*θl)./C) .<= tol)
             end
@@ -189,7 +189,7 @@ using ComponentArrays
                     C = heatcapacity.(soil,θw,θl,θm,θo);
                    enthalpy.(T,C,L,θl) # compute enthalpy at "true" temperature
                 end
-                state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
+                state = (T=T,C=C,Ceff=similar(C),H=H,θl=θl,)
                 @inferred sfcc(soil, heat, state)
                 @test all(abs.(T.-(H .- L.*θl)./C) .<= tol)
             end
@@ -219,7 +219,7 @@ using ComponentArrays
             T_ .= T
             C = similar(C,eltype(p))
             θl = similar(θl,eltype(p))
-            state = (T=T_,C=C,Ceff=similar(C),H=H,θl=θl,θw=θw.+zero(T),θm=θm.+zero(T),θo=θo.+zero(T),θp=θp.+zero(T))
+            state = (T=T_,C=C,Ceff=similar(C),H=H,θl=θl,)
             SFCC(McKenzie(γ=Param(p.γ[1])))(soil, heat, state)
             state.T[1]
         end
-- 
GitLab