Update Richards eq example and change reduction factor

Now using the original reduction factor formulation from CryoGrid CM again.

examples/08_heat_sfcc_richardseq_samoylov.jl

@@ -18,30 +18,30 @@ initT = initializer(:T, tempprofile);
 
 # Here we conigure the water retention curve and freeze curve. The van Genuchten parameters coorespond to that
 # which would be reasonable for a silty soil.
-swrc = VanGenuchten(α=0.1, n=1.8)
+swrc = VanGenuchten(α=1.0, n=1.5)
 sfcc = PainterKarra(ω=0.0; swrc)
-waterflow = RichardsEq(swrc=swrc);
+waterflow = RichardsEq(;swrc);
 
 # We use the enthalpy-based heat diffusion with high accuracy Newton-based solver for inverse enthalpy mapping
-heatop = Heat.MOLEnthalpy(SFCCNewtonSolver())
+heatop = Heat.MOLEnthalpy(SFCCPreSolver())
 upperbc = WaterHeatBC(SurfaceWaterBalance(forcings), TemperatureGradient(forcings.Tair, NFactor(nf=0.6, nt=0.9)));
 
 # We will use a simple stratigraphy with three subsurface soil layers.
 # Note that the @Stratigraphy macro lets us list multiple subsurface layers without wrapping them in a tuple.
-heat = HeatBalance(heatop, freezecurve=sfcc, advection=false)
-water = WaterBalance(RichardsEq(; swrc))
+heat = HeatBalance(heatop, freezecurve=sfcc)
+water = WaterBalance(waterflow)
 strat = @Stratigraphy(
     -2.0u"m" => Top(upperbc),
-    0.0u"m" => Ground(MineralOrganic(por=0.80,sat=0.7,org=0.75); heat, water),
-    0.2u"m" => Ground(MineralOrganic(por=0.40,sat=0.8,org=0.10); heat, water),
+    0.0u"m" => Ground(MineralOrganic(por=0.80,sat=0.8,org=0.75); heat, water),
+    0.2u"m" => Ground(MineralOrganic(por=0.40,sat=0.9,org=0.10); heat, water),
     2.0u"m" => Ground(MineralOrganic(por=0.10,sat=1.0,org=0.0); heat, water),
     1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"W/m^2"))
 );
 grid = CryoGrid.Presets.DefaultGrid_2cm
 tile = Tile(strat, grid, initT);
 u0, du0 = initialcondition!(tile, tspan)
 prob = CryoGridProblem(tile, u0, tspan, saveat=3*3600, savevars=(:T,:θw,:θwi,:kw));
-integrator = init(prob, Euler(), dt=60.0)
+integrator = init(prob, CGEuler())
 
 using BenchmarkTools
 @btime $tile($du0, $u0, $prob.p, $prob.tspan[1])
@@ -62,15 +62,17 @@ state = getstate(integrator);
 @time while integrator.t < prob.tspan[end]
     @assert all(isfinite.(integrator.u))
     @assert all(0 .<= integrator.u.sat .<= 1)
-    ## run the integrator forward in daily increments
+    ## run the integrator forward in 24 hour increments
     step!(integrator, 24*3600.0)
     t = convert_t(integrator.t)
     @info "t=$t, current dt=$(integrator.dt*u"s")"
+    if integrator.dt < 0.01
+        @warn "dt=$(integrator.dt) < 0.01 s; this may indicate a problem!"
+        break
+    end
 end;
 out = CryoGridOutput(integrator.sol)
 
-@run step!(integrator)
-
 # Check mass conservation...
 water_added = values(sum(upreferred.(forcings.rainfall.(tspan[1]:Hour(3):tspan[2]).*u"m/s".*3u"hr")))[1]
 water_mass = Diagnostics.integrate(out.θwi, tile.grid)
@@ -81,5 +83,5 @@ import Plots
 zs = [1,5,10,15,20,30,40,50,100,150,200]u"cm"
 cg = Plots.cgrad(:copper,rev=true);
 p1 = Plots.plot(out.T[Z(Near(zs))], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", leg=false, size=(800,500), dpi=150)
-p2 = Plots.plot(out.sat[Z(1:10)], color=cg[LinRange(0.0,1.0,10)]', ylabel="Saturation", leg=false, size=(800,500), dpi=150)
+p2 = Plots.plot(out.sat[Z(Near([1,5,10,15,30,50,100]u"cm"))], color=cg[LinRange(0.0,1.0,10)]', ylabel="Saturation", leg=false, size=(800,500), dpi=150)
 Plots.plot(p1, p2, size=(1200,400))

src/Physics/Hydrology/water_balance.jl

@@ -9,8 +9,8 @@ function limit_upper_flux(water::WaterBalance, jw, θw, θwi, θsat, sat, Δz)
     # case (ii): jw > 0 -> inflow -> limit based on available space
     # case (iii): jw = 0 -> no flow, limit has no effect
     jw = min(max(jw, -θw*Δz), (θsat - θwi)*Δz)
-    # influx reduction factor
-    r = reductionfactor(water, sat)*(jw > zero(jw)) + 1.0*(jw <= zero(jw))
+    # influx reduction factors
+    r = reduction_factor_in(water, sat)*(jw > zero(jw)) + reduction_factor_out(water, sat)*(jw <= zero(jw))
     return r*jw
 end
 
@@ -25,8 +25,8 @@ function limit_lower_flux(water::WaterBalance, jw, θw, θwi, θsat, sat, Δz)
     # case (ii): jw > 0 -> outflow -> limit based on available water
     # case (iii): jw = 0 -> no flow, limit has no effect
     jw = min(max(jw, (θwi - θsat)*Δz), θw*Δz)
-    # influx reduction factor
-    r = reductionfactor(water, sat)*(jw < zero(jw)) + 1.0*(jw >= zero(jw))
+    # influx reduction factors
+    r = reduction_factor_in(water, sat)*(jw < zero(jw)) + reduction_factor_out(water, sat)*(jw >= zero(jw))
     return r*jw
 end
 
@@ -128,16 +128,18 @@ function waterprognostic!(::SubSurface, ::WaterBalance{<:BucketScheme}, state)
 end
 
 # Helper methods
-"""
-    reductionfactor(water::WaterBalance, x)
+function reduction_factor_out(water::WaterBalance, sat)
+    smoothness = water.prop.rf_smoothness
+    rf = 1 - exp(-sat/smoothness)
+    return rf
+end
+function reduction_factor_in(water::WaterBalance, sat)
+    smoothness = water.prop.rf_smoothness
+    rf = 1 - exp((sat-1)/smoothness)
+    return rf
+end
 
-Hydraulic conductivity reduction factor for near-saturated conditions:
-```math
-r(x) = (1-exp(-βx^2))*(1-exp(-β*(1-x)^2))
-```
-where β is a smoothness parameter and c is the "center" or shift parameter.
-"""
-reductionfactor(water::WaterBalance, x) = 1 - exp(-water.prop.r_β*(1-x)^2)
+# reduction_factor(water::WaterBalance, x) = 1 - exp(-water.prop.r_β*(1-x)^2)
 
 """
     resetfluxes!(::SubSurface, water::WaterBalance, state)
@@ -220,7 +222,7 @@ function CryoGrid.timestep(
 ) where {TFlow,TET}
     dtmax = Inf
     @inbounds for i in 1:length(state.sat)
-        dt = water.dtlim(state.∂sat∂t[i], state.sat[i], state.t, zero(state.t), one(state.t))
+        dt = water.dtlim(state.∂θwi∂t[i], state.θwi[i], state.t, zero(state.t), state.θsat[i])
         dt = isfinite(dt) && dt > zero(dt) ? dt : Inf # make sure it's +Inf
         dtmax = min(dtmax, dt)
     end

src/Physics/Hydrology/water_types.jl

@@ -3,11 +3,12 @@
 
 Numerical constants shared across water balance implementations.
 """
-Base.@kwdef struct WaterBalanceProperties{Tρw,TLsg,Trb,Trc}
+Base.@kwdef struct WaterBalanceProperties{Tρw,TLsg,Trfs}
     ρw::Tρw = CryoGrid.Constants.ρw
     Lsg::TLsg = CryoGrid.Constants.Lsg
-    r_β::Trb = 1e3 # reduction factor scale parameter
-    r_c::Trc = 0.96325 # reduction factor shift parameter
+    rf_smoothness::Trfs = 0.3
+    # r_β::Trb = 1e3 # reduction factor scale parameter
+    # r_c::Trc = 0.96325 # reduction factor shift parameter
 end
 
 # do not parameterize water balance constants

src/Physics/Soils/soil_water.jl

@@ -78,6 +78,7 @@ function Hydrology.watercontent!(soil::Soil, water::WaterBalance{<:RichardsEq{Sa
         @inbounds for i in 1:length(state.ψ₀)
             state.θsat[i] = θsat = Hydrology.maxwater(soil, water, state, i)
             state.θwi[i] = θwi = state.sat[i]*θsat
+            state.θw[i] = state.θwi[i] # initially set liquid water content to total water content (coupling with HeatBalance will overwrite this)
             # this is a bit shady because we're allowing for incorrect/out-of-bounds values of θwi, but this is necessary
             # for solving schemes that might attempt to use illegal state values
             state.ψ₀[i] = f⁻¹(max(θres, min(θwi, θsat)); θsat)