diff --git a/Project.toml b/Project.toml
index 544ab68a8539b99efe518c98abd007d24db50342..517d8d6c882dbb1d60d3fa7435181b0477e944f1 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 227bfb1b263a7ca38e74baedf51cb47eb3e37b21..c3defedb4f14ed2e04bea96afaf90919ac2de33c 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 57e597dd93eb9ff457c490e4f315e00efe952d0b..a7e4279901d5d2448f24209c2ae190126e15edf2 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 067fd43a01ec505446f187bf82aa5d37eab68f8e..b688de7a7dfba1b9ffed2a1fa01da611c7ad2b31 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 ab1f8b46769ea26c287d0b3fe9ff91b562f60eba..9949f11c583b2d45858e1054c791a9cb2c523bfd 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 530905bff865c112eb3db50e49abf877fe8d0229..75fdc97c0341010f117e13e69c970642ff3005b1 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 8d200fab51ee446c92b69da267b6c7d288b36412..ed11bac9a31cbe84a8a7357269e27a009be3674b 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 d296a561ee43ce7c13f4e213d604684caba006f1..25471835a76b625662951472a687a9cf32ac93b2 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 49fcf40974a262def55cf9925fa8fd69f0c4cfaf..05df88acf0199478232606d62cba089c45533adc 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 2916e4bc81741ef15db6a83cbfd619b2a50a17ff..d1c6748a95cd5713e601335c700be10de6f2bcd0 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