diff --git a/examples/08_heat_sfcc_richardseq_samoylov.jl b/examples/08_heat_sfcc_richardseq_samoylov.jl index 4d82bf00..c48afb16 100644 --- a/examples/08_heat_sfcc_richardseq_samoylov.jl +++ b/examples/08_heat_sfcc_richardseq_samoylov.jl @@ -18,22 +18,22 @@ 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")) ); @@ -41,7 +41,7 @@ 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)) diff --git a/src/Physics/Hydrology/water_balance.jl b/src/Physics/Hydrology/water_balance.jl index 959d287c..d5fe1e35 100644 --- a/src/Physics/Hydrology/water_balance.jl +++ b/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 diff --git a/src/Physics/Hydrology/water_types.jl b/src/Physics/Hydrology/water_types.jl index 02626dcb..5c1440d5 100644 --- a/src/Physics/Hydrology/water_types.jl +++ b/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 diff --git a/src/Physics/Soils/soil_water.jl b/src/Physics/Soils/soil_water.jl index 13ee5086..dc5f356b 100644 --- a/src/Physics/Soils/soil_water.jl +++ b/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)