diff --git a/src/Physics/HeatConduction/HeatConduction.jl b/src/Physics/HeatConduction/HeatConduction.jl
index 1db0bf5e24d796ebc35de3b881c01411a2fe747b..067fd43a01ec505446f187bf82aa5d37eab68f8e 100644
--- a/src/Physics/HeatConduction/HeatConduction.jl
+++ b/src/Physics/HeatConduction/HeatConduction.jl
@@ -11,6 +11,7 @@ using CryoGrid.Numerics: nonlineardiffusion!, harmonicmean!, harmonicmean, heavi
 using CryoGrid.Layers: Soil, porosity
 using CryoGrid.Utils
 
+using Base: @propagate_inbounds
 using DimensionalData
 using IfElse
 using Interpolations: Linear, Flat
diff --git a/src/Physics/HeatConduction/heat.jl b/src/Physics/HeatConduction/heat.jl
index 79c07c58be79f74bb9999189c8ae5dbca8c3e4e5..ab1f8b46769ea26c287d0b3fe9ff91b562f60eba 100644
--- a/src/Physics/HeatConduction/heat.jl
+++ b/src/Physics/HeatConduction/heat.jl
@@ -201,8 +201,14 @@ total water content (θw), and liquid water content (θl).
             I_c*(H/Lθ) + I_t
         end
     end
-    @inbounds @. state.θl = freezethaw(state.H, heat.L, state.θw)*state.θw
-    heatcapacity!(layer, heat, state) # update heat capacity, C
-    @inbounds @. state.T = enthalpyinv(state.H, state.C, heat.L, state.θw)
+    @inbounds for i in 1:length(state.H)
+        let θw = state.θw[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.T[i] = enthalpyinv(H, state.C[i], L, θw)
+        end
+    end
     return nothing
 end
diff --git a/src/Physics/HeatConduction/heat_bc.jl b/src/Physics/HeatConduction/heat_bc.jl
index 9983809f3e155975b6c8ac919819cd4498b19279..e7ea2fdf7135e2c204e09bd0f34872c4b2863eaa 100644
--- a/src/Physics/HeatConduction/heat_bc.jl
+++ b/src/Physics/HeatConduction/heat_bc.jl
@@ -12,7 +12,7 @@ struct TemperatureGradient{E,F} <: BoundaryProcess
 end
 BoundaryStyle(::Type{<:TemperatureGradient}) = Dirichlet()
 BoundaryStyle(::Type{<:TemperatureGradient{<:Damping}}) = Neumann()
-@inline boundaryvalue(bc::TemperatureGradient, l1, p2, l2, s1, s2) where {F} = bc.T(s1.t)
+@inline boundaryvalue(bc::TemperatureGradient, l1, p2, l2, s1, s2) = bc.T(s1.t)
 
 @with_kw struct NFactor{W,S} <: BoundaryEffect
     winterfactor::W = Param(1.0, bounds=(0.0,1.0)) # applied when Tair <= 0
diff --git a/src/Physics/HeatConduction/soil/soilheat.jl b/src/Physics/HeatConduction/soil/soilheat.jl
index e1e08992f7d7b4dc4f0efc2ac8feae73a39128d6..8d200fab51ee446c92b69da267b6c7d288b36412 100644
--- a/src/Physics/HeatConduction/soil/soilheat.jl
+++ b/src/Physics/HeatConduction/soil/soilheat.jl
@@ -1,3 +1,4 @@
+# Thermal properties
 @inline function heatcapacity(soil::Soil, totalwater, liquidwater, mineral, organic)
     @unpack cw,co,cm,ca,ci = soil.hc
     let air = 1.0 - totalwater - mineral - organic,
@@ -14,13 +15,16 @@ 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 function heatcapacity!(soil::Soil, ::Heat, state)
-    @inbounds @. state.C = heatcapacity(soil, state.θw, state.θl, state.θm, state.θo)
+@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 function thermalconductivity!(soil::Soil, ::Heat, state)
-    @inbounds @. state.kc = thermalconductivity(soil, state.θw, state.θl, state.θm, state.θo)
+@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")
 
 """