Fix heterogeneous soil water/ice initialization

Author: bgroenks96

examples/heat_freeW_lite_implicit.jl

@@ -9,6 +9,13 @@ using CryoGrid.LiteImplicit
 
 # Load forcings and build stratigraphy like before.
 forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
+soilprofile = SoilProfile(
+    0.0u"m" => MineralOrganic(por=0.80,sat=1.0,org=0.75),
+    0.1u"m" => MineralOrganic(por=0.80,sat=1.0,org=0.25),
+    0.4u"m" => MineralOrganic(por=0.55,sat=1.0,org=0.25),
+    3.0u"m" => MineralOrganic(por=0.50,sat=1.0,org=0.0),
+    10.0u"m" => MineralOrganic(por=0.30,sat=1.0,org=0.0),
+)
 tempprofile_linear = TemperatureProfile(
     0.0u"m" => -30.0u"°C",
     10.0u"m" => -10.0u"°C", 
@@ -20,22 +27,19 @@ upperbc = TemperatureBC(forcings.Tair, NFactor())
 initT = initializer(:T, tempprofile_linear)
 heatop = Heat.EnthalpyImplicit()
 freezecurve = FreeWater()
-strat = @Stratigraphy(
+heat = HeatBalance(heatop; freezecurve)
+soil_layers = map(para -> Ground(para; heat), soilprofile)
+strat = Stratigraphy(
     z_top => Top(upperbc),
-    0.0u"m" => Ground(MineralOrganic(por=0.80,sat=1.0,org=0.75), heat=HeatBalance(heatop; freezecurve)),
-    0.1u"m" => Ground(MineralOrganic(por=0.80,sat=1.0,org=0.25), heat=HeatBalance(heatop; freezecurve)),
-    0.4u"m" => Ground(MineralOrganic(por=0.55,sat=1.0,org=0.25), heat=HeatBalance(heatop; freezecurve)),
-    3.0u"m" => Ground(MineralOrganic(por=0.50,sat=1.0,org=0.0), heat=HeatBalance(heatop; freezecurve)),
-    10.0u"m" => Ground(MineralOrganic(por=0.30,sat=1.0,org=0.0), heat=HeatBalance(heatop; freezecurve)),
+    soil_layers,
     z_bot => Bottom(GeothermalHeatFlux(0.053u"W/m^2"))
 );
 modelgrid = CryoGrid.Presets.DefaultGrid_2cm
 tile = Tile(strat, modelgrid, initT);
 
 # Since the solver can take daily timesteps, we can easily specify longer simulation time spans at minimal cost.
-# Here we specify a time span of 5 years.
-tspan = (DateTime(2010,1,1), DateTime(2015,1,1))
-tspan_sol = convert_tspan(tspan)
+# Here we specify a time span of 20 years.
+tspan = (DateTime(2000,1,1), DateTime(2020,1,1))
 u0, du0 = initialcondition!(tile, tspan);
 prob = CryoGridProblem(tile, u0, tspan, saveat=24*3600.0, savevars=(:T,))
 sol = @time solve(prob, LiteImplicitEuler(), dt=24*3600.0)

src/Physics/Soils/ground.jl

@@ -16,14 +16,14 @@ Base.@kwdef struct Ground{Tpara,Theat<:Optional{HeatBalance},Twater<:Optional{Wa
     para::Tpara = MineralOrganic() # ground parameterization
     heat::Theat = HeatBalance() # heat conduction
     water::Twater = nothing # water balance
-    fcsolver::Tsolver = default_fcsolver(heat, water)
+    fcsolver::Tsolver = default_fcsolver(para, heat, water)
     aux::Taux = nothing # user-defined specialization
 end
 # Convenience constructors
 Ground(para::GroundParameterization; kwargs...) = Ground(; para, kwargs...)
 
-default_fcsolver(::Any, ::Any) = nothing
-default_fcsolver(::HeatBalance{<:SFCC}, ::Nothing) = SFCCPreSolver(FreezeCurves.SFCCPreSolverCache1D())
-default_fcsolver(::HeatBalance{<:SFCC}, ::WaterBalance) = SFCCPreSolver(FreezeCurves.SFCCPreSolverCacheND())
+default_fcsolver(::Any, ::Any, ::Any) = nothing
+default_fcsolver(::GroundParameterization, ::HeatBalance{<:SFCC}, ::Nothing) = SFCCPreSolver(FreezeCurves.SFCCPreSolverCache1D())
+default_fcsolver(::GroundParameterization, ::HeatBalance{<:SFCC}, ::WaterBalance) = SFCCPreSolver(FreezeCurves.SFCCPreSolverCacheND())
 
 fcsolver(ground::Ground) = ground.fcsolver

src/Physics/Soils/soil_para.jl

@@ -103,6 +103,9 @@ end
 """
 Ground(soilprofile::SoilProfile; kwargs...) = Ground(Heterogeneous(soilprofile); kwargs...)
 
+# add dispatch for default_fcsolver that selects the ND presolver
+default_fcsolver(::Heterogeneous, ::HeatBalance{<:SFCC}, ::Nothing) = SFCCPreSolver(FreezeCurves.SFCCPreSolverCacheND())
+
 saturation(::Soil{<:Heterogeneous{<:MineralOrganic}}, state) = state.sat
 porosity(::Soil{<:Heterogeneous{<:MineralOrganic}}, state) = state.por
 mineral(::Soil{<:Heterogeneous{<:MineralOrganic}}, state) = state.θm
@@ -138,6 +141,7 @@ CryoGrid.initializers(soil::Soil{<:Heterogeneous{<:MineralOrganic}}) = (
     initializer(:por, map(para -> para.por, soil.para.profile)),
     initializer(:sat, map(para -> para.sat, soil.para.profile)),
     initializer(:org, map(para -> para.org, soil.para.profile)),
+    initializer(:θwi, map(para -> para.por*para.sat, soil.para.profile)),
     initializer(:kh_m, map(para -> para.heat.kh_m, soil.para.profile)),
     initializer(:kh_o, map(para -> para.heat.kh_o, soil.para.profile)),
     initializer(:kh_w, map(para -> para.heat.kh_w, soil.para.profile)),

src/Tiles/stratigraphy.jl

@@ -22,6 +22,11 @@ struct Stratigraphy{N,TLayers<:NamedTuple,TBoundaries}
         sub::AbstractVector{<:Pair{<:DistQuantity}},
         bot::Pair{<:DistQuantity,<:Bottom}
     ) = Stratigraphy(top, Tuple(sub), bot)
+    Stratigraphy(
+        top::Pair{<:DistQuantity,<:Top},
+        sub::Numerics.Profile{N,<:Tuple{Vararg{DistQuantity}},<:Tuple{Vararg{SubSurface}}},
+        bot::Pair{<:DistQuantity,<:Bottom}
+    ) where N = Stratigraphy(top, Tuple(sub), bot)
     function Stratigraphy(
         # use @nospecialize to (hopefully) reduce compilation overhead
         @nospecialize(top::Pair{<:DistQuantity,<:Top}),