Coupled soil heat and water transport

In this example, we construct a Tile from a Stratigraphy of three soil layers with coupled water and heat transport. We use the Richards-Richardson equation for unsaturated flow in porous media to account for flow due to capillary suction and pressure gradients.

using CryoGrid


forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
tspan = (DateTime(2011,1,1),DateTime(2011,12,31))
T0 = forcings.Tair(tspan[1])
tempprofile = TemperatureProfile(
    0.0u"m" => 0.0*u"°C",
    1.0u"m" => -8.0u"°C",
    20.0u"m" => -10.0u"°C",
    1000.0u"m" => 1.0u"°C",
)
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(α=1.0, n=1.5)
sfcc = PainterKarra(ω=0.0; swrc)
waterflow = RichardsEq(;swrc);

We use the enthalpy-based heat diffusion with high accuracy Newton-based solver for inverse enthalpy mapping

heatop = Heat.Diffusion1D(:H)
upperbc = WaterHeatBC(SurfaceWaterBalance(), TemperatureBC(Input(: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)
water = WaterBalance(waterflow)
strat = @Stratigraphy(
    -2.0u"m" => Top(upperbc),
    0.0u"m" => Ground(SimpleSoil(por=0.80,sat=0.9,org=0.75,freezecurve=sfcc); heat, water),
    0.5u"m" => Ground(SimpleSoil(por=0.40,sat=1.0,org=0.10,freezecurve=sfcc); heat, water),
    1.0u"m" => Ground(SimpleSoil(por=0.30,sat=1.0,org=0.0,freezecurve=sfcc); heat, water),
    2.0u"m" => Ground(SimpleSoil(por=0.10,sat=1.0,org=0.0,freezecurve=sfcc); heat, water),
    1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"W/m^2"))
);
grid = CryoGrid.DefaultGrid_2cm
tile = Tile(strat, grid, forcings, initT);
u0, du0 = @time initialcondition!(tile, tspan)
prob = CryoGridProblem(tile, u0, tspan, saveat=3*3600, savevars=(:T,:θw,:θwi,:kw));
integrator = init(prob, CGEuler())

t: 6.34610592e10
u: ComponentVector{Float64}(top = (runoff = [0.0], infil = [0.0]), sat = [0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], H = [8.15049472830431e7, 4.781032348767249e7, 3.676461736314711e7, 3.068581187047787e7, 2.6642404873069827e7, 2.3665070390244327e7, 2.132790997123952e7, 1.9410593495478917e7, 1.7786184147987116e7, 1.6375715790648941e7  …  -1.72986924095188e7, -1.6032937325758928e7, -1.3720467877056617e7, -1.1388499640093846e7, -9.027026621657975e6, -6.616916590730443e6, -4.1153292066077325e6, -1.3964128860094277e6, 2.23844582003107e6, 3.4374081632653065e7])

Run the integrator forward in time until the end of the tspan; Note that this is currently somewhat slow since the integrator must take very small time steps during the thawed season; expect it to take about 3-5 minutes per year on a typical workstation/laptop. Minor speed-ups might be possible by tweaking the dt limiters or by using the SFCCPreSolver for the freeze curve.

@time while integrator.t < prob.tspan[end]
    @assert all(isfinite.(integrator.u))
    @assert all(0 .<= integrator.u.sat .<= 1)
    # 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)

CryoGridOutput with 2913 time steps (2011-01-01T00:00:00 to 2011-12-31T00:00:00) and 8 variables:
    sat => DimArray of Float64 with dimensions (278, 2913)
    H => DimArray of Quantity{Float64, 𝐌 𝐋^-1 𝐓^-2, Unitful.FreeUnits{(J, m^-3), 𝐌 𝐋^-1 𝐓^-2, nothing}} with dimensions (278, 2913)
    θwi => DimArray of Float64 with dimensions (278, 2913)
    θw => DimArray of Float64 with dimensions (278, 2913)
    kw => DimArray of Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(m, s^-1), 𝐋 𝐓^-1, nothing}} with dimensions (279, 2913)
    T => DimArray of Quantity{Float64, 𝚯, Unitful.FreeUnits{(K,), 𝚯, Unitful.Affine{-5463//20}}} with dimensions (278, 2913)
    top => 
        runoff => DimArray of Quantity{Float64, 𝐋^3, Unitful.FreeUnits{(m^3,), 𝐋^3, nothing}} with dimensions (1, 2913)
        infil => DimArray of Quantity{Float64, 𝐋^3, Unitful.FreeUnits{(m^3,), 𝐋^3, nothing}} with dimensions (1, 2913)

Check mass conservation...

water_added = out.top.infil[end]
water_mass = Diagnostics.integrate(out.θwi, grid)
water_resid = water_mass[end] - water_mass[1] - water_added

-1.5532920193275896 m^3

Plot it!

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(Near([1,5,10,15,20,25,30]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))

This page was generated using Literate.jl.