Fix freeW init bug and make porosity a state var

src/Layers/soil.jl

@@ -67,11 +67,13 @@ end
 # Methods
 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)

src/Physics/HeatConduction/heat.jl

@@ -75,12 +75,6 @@ variables(heat::Heat{<:FreezeCurve,Temperature}) = (
     Diagnostic(:kc, Float"W/m/K", OnGrid(Cells)),
     Diagnostic(:θl, Float"1", OnGrid(Cells)),
 )
-""" Initial condition for heat conduction (all state configurations) on any subsurface layer. """
-function initialcondition!(layer::SubSurface, heat::Heat, state)
-    L = heat.L
-    heatcapacity!(layer, heat, state)
-    @. state.H = enthalpy(state.T, state.C, L, state.θl)
-end
 """ Diagonstic step for heat conduction (all state configurations) on any subsurface layer. """
 function diagnosticstep!(layer::SubSurface, heat::Heat, state)
     # Reset energy flux to zero; this is redundant when H is the prognostic variable
@@ -88,10 +82,10 @@ function diagnosticstep!(layer::SubSurface, heat::Heat, state)
     @. state.dH = zero(eltype(state.dH))
     # Evaluate the freeze curve (updates T, C, and θl)
     fc! = freezecurve(heat);
-    fc!(layer,heat,state)
+    fc!(layer, heat, state)
     # Update thermal conductivity
     thermalconductivity!(layer, heat, state)
-    # Harmonic mean of conductivities
+    # Harmonic mean of inner conductivities
     @inbounds let k = (@view state.k[2:end-1]),
         Δk = Δ(state.grids.k);
         harmonicmean!(k, state.kc, Δk)
@@ -144,7 +138,7 @@ Generic top interaction. Computes flux dH at top cell.
 function interact!(top::Top, bc::BoundaryProcess, sub::SubSurface, heat::Heat, stop, ssub)
     # thermal conductivity at boundary
     # assumes (1) k has already been computed, (2) surface conductivity = cell conductivity
-    @inbounds ssub.k[1] = ssub.k[2]
+    @inbounds ssub.k[1] = ssub.kc[1]
     # boundary flux
     @inbounds ssub.dH[1] += boundaryflux(bc, top, heat, sub, stop, ssub)
     return nothing # ensure no allocation
@@ -155,7 +149,7 @@ Generic bottom interaction. Computes flux dH at bottom cell.
 function interact!(sub::SubSurface, heat::Heat, bot::Bottom, bc::BoundaryProcess, ssub, sbot)
     # thermal conductivity at boundary
     # assumes (1) k has already been computed, (2) bottom conductivity = cell conductivity
-    @inbounds ssub.k[end] = ssub.k[end-1]
+    @inbounds ssub.k[end] = ssub.kc[end]
     # boundary flux
     @inbounds ssub.dH[end] += boundaryflux(bc, bot, heat, sub, sbot, ssub)
     return nothing # ensure no allocation
@@ -192,7 +186,7 @@ total water content (θw), and liquid water content (θl).
 """
 @inline function (fc::FreeWater)(layer::SubSurface, heat::Heat{FreeWater,Enthalpy}, state)
     @inline function enthalpyinv(H, C, L, θtot)
-        let θtot = max(1.0e-8,θtot),
+        let θtot = max(1.0e-8, θtot),
             Lθ = L*θtot,
             I_t = H > Lθ,
             I_f = H <= 0.0;
@@ -200,7 +194,7 @@ total water content (θw), and liquid water content (θl).
         end
     end
     @inline function freezethaw(H, L, θtot)
-        let θtot = max(1.0e-8,θtot),
+        let θtot = max(1.0e-8, θtot),
             Lθ = L*θtot,
             I_t = H > Lθ,
             I_c = (H > 0.0) && (H <= Lθ);

src/Physics/HeatConduction/soil/sfcc.jl

@@ -93,7 +93,7 @@ end
 sfccparams(f::DallAmico, soil::Soil, heat::Heat, state) = (
     f.Tₘ,
     f.θres,
-    porosity(soil), # θ saturated = porosity
+    state.θp, # θ saturated = porosity
     state.θw, # total water content
     heat.L, # specific latent heat of fusion, L
     f.α,
@@ -129,7 +129,7 @@ end
 sfccparams(f::McKenzie, soil::Soil, heat::Heat, state) = (
     f.Tₘ,
     f.θres,
-    porosity(soil), # θ saturated = porosity
+    state.θp, # θ saturated = porosity
     state.θw, # total water content
     f.γ,
 )
@@ -154,7 +154,7 @@ end
 sfccparams(f::Westermann, soil::Soil, heat::Heat, state) = (
     f.Tₘ,
     f.θres,
-    porosity(soil), # θ saturated = porosity
+    state.θp, # θ saturated = porosity
     state.θw, # total water content
     f.δ,
 )

src/Physics/HeatConduction/soil/soilheat.jl

@@ -26,11 +26,23 @@ end
 include("sfcc.jl")
 
 """ Initial condition for heat conduction (all state configurations) on soil layer. """
-function initialcondition!(soil::Soil, heat::Heat{<:SFCC}, state)
+function initialcondition!(soil::Soil, heat::Heat, state)
     initialcondition!(soil, state)
     L = heat.L
-    sfcc = freezecurve(heat)
-    state.θl .= sfcc.f.(state.T, sfccparams(sfcc.f, soil, heat, state)...)
+    fc! = freezecurve(heat)
+    fc!(soil, heat, state)
     heatcapacity!(soil, heat, state)
     @. state.H = enthalpy(state.T, state.C, L, state.θl)
 end
+""" Initial condition for heat conduction (all state configurations) with free water freeze curve on soil layer. """
+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, 0.0, state.θp)
+    heatcapacity!(soil, heat, state)
+    @. state.H = enthalpy(state.T, state.C, L, state.θl)
+    # evaluate freeze curve
+    fc! = freezecurve(heat)
+    fc!(soil, heat, state)
+end

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))
+                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))
                 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))
+                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))
                 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))
+                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))
                 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))
+                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))
                 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))
+                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))
                 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))
+                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))
                 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))
+                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))
                 @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))
+                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))
                 @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))
+                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))
                 @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))
+            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))
             SFCC(McKenzie(γ=Param(p.γ[1])))(soil, heat, state)
             state.T[1]
         end