Soil heat with bulk snow scheme

In this example, we construct a Tile consisting of a soil column with (i) heat conduction and (ii) a bulk (single-layer) snow scheme. The snowfall data comes from the ERA-Interim reanalysis product.

using CryoGrid
using OrdinaryDiffEq

First we set up the model:

forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
soilprofile, tempprofile = CryoGrid.Presets.SamoylovDefault
initT = initializer(:T, tempprofile)
initsat = initializer(:sat, 1.0)
z_top = -2.0u"m"
z_sub = keys(soilprofile)
z_bot = 1000.0u"m"
upperbc = WaterHeatBC(
    SurfaceWaterBalance(forcings),
    TemperatureBC(forcings.Tair)
)
snowmass = SnowMassBalance(
    ablation = Snow.DegreeDayMelt(factor=5.0u"mm/K/d")
)
snowpack = Snowpack(
    para=Snow.Bulk(),
    mass=snowmass,
    heat=HeatBalance(),
    water=WaterBalance(),
)
ground_layers = map(soilprofile) do para
    Ground(para, heat=HeatBalance(), water=WaterBalance())
end
strat = @Stratigraphy(
    z_top => Top(upperbc),
    z_top => snowpack,
    ground_layers...,
    z_bot => Bottom(GeothermalHeatFlux(0.053u"J/s/m^2"))
);
modelgrid = CryoGrid.Presets.DefaultGrid_5cm
tile = Tile(strat, modelgrid, initT, initsat)

Tile (iip=true) with layers (:top, :snowpack, :ground1, :ground2, :ground3, :ground4, :ground5, :bottom), Grid{Edges, CryoGrid.Numerics.UnitRectangle, Float64, Vector{Float64}}, Stratigraphy{8, @NamedTuple{top::Top{CoupledProcesses{Tuple{SurfaceWaterBalance{InterpolatedForcing{m s^-1, Float64, Interpolations.GriddedInterpolation{Float64, 1, Vector{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}}}}, InterpolatedForcing{m s^-1, Float64, Interpolations.GriddedInterpolation{Float64, 1, Vector{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}}}}}, TemperatureBC{Nothing, InterpolatedForcing{°C, Float64, Interpolations.GriddedInterpolation{Float64, 1, Vector{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}}}}}}}}, snowpack::Snowpack{CryoGrid.Snow.Bulk{CryoGrid.Snow.ConstantDensity{Float64}, Float64, CryoGrid.Snow.SnowThermalProperties{CryoGrid.Snow.SturmQuadratic{ClosedInterval{Float64}, @NamedTuple{a::Float64, b::Float64, c::Float64}, @NamedTuple{a::Float64, b::Float64, c::Float64}}, ThermalProperties{@NamedTuple{kh_w::Float64, kh_i::Float64, kh_a::Float64, ch_w::Float64, ch_i::Float64, ch_a::Float64}}}, HydraulicProperties{@NamedTuple{kw_sat::Float64}}}, SnowMassBalance{CryoGrid.Snow.LinearAccumulation{Float64}, CryoGrid.Snow.DegreeDayMelt{Float64, Float64}, Nothing, CryoGrid.MaxDelta{Float64}}, HeatBalance{FreeWater, CryoGrid.Heat.Diffusion1D{:H, typeof(CryoGrid.Heat.quadratic_parallel_conductivity), typeof(CryoGrid.Heat.weighted_average_heatcapacity)}, CryoGrid.MaxDelta{Int64}, CryoGrid.Heat.HeatBalanceProperties{@NamedTuple{ρw::Float64, Lsl::Float64, L::Float64}}}, WaterBalance{NoFlow, Nothing, CryoGrid.MaxDelta{Float64}, Nothing, WaterBalanceProperties{Float64, Float64, Float64}}, Nothing}, ground1::Ground{MineralOrganic{Float64, Float64, Float64, ThermalProperties{@NamedTuple{kh_w::Float64, kh_i::Float64, kh_a::Float64, ch_w::Float64, ch_i::Float64, ch_a::Float64, kh_m::Float64, kh_o::Float64, ch_m::Float64, ch_o::Float64}}, HydraulicProperties{@NamedTuple{kw_sat::Float64, fieldcapacity::Float64}}}, HeatBalance{FreeWater, CryoGrid.Heat.Diffusion1D{:H, typeof(CryoGrid.Heat.quadratic_parallel_conductivity), typeof(CryoGrid.Heat.weighted_average_heatcapacity)}, CryoGrid.MaxDelta{Int64}, CryoGrid.Heat.HeatBalanceProperties{@NamedTuple{ρw::Float64, Lsl::Float64, L::Float64}}}, WaterBalance{NoFlow, Nothing, CryoGrid.MaxDelta{Float64}, Nothing, WaterBalanceProperties{Float64, Float64, Float64}}, Nothing, Nothing}, ground2::Ground{MineralOrganic{Float64, Float64, Float64, ThermalProperties{@NamedTuple{kh_w::Float64, kh_i::Float64, kh_a::Float64, ch_w::Float64, ch_i::Float64, ch_a::Float64, kh_m::Float64, kh_o::Float64, ch_m::Float64, ch_o::Float64}}, HydraulicProperties{@NamedTuple{kw_sat::Float64, fieldcapacity::Float64}}}, HeatBalance{FreeWater, CryoGrid.Heat.Diffusion1D{:H, typeof(CryoGrid.Heat.quadratic_parallel_conductivity), typeof(CryoGrid.Heat.weighted_average_heatcapacity)}, CryoGrid.MaxDelta{Int64}, CryoGrid.Heat.HeatBalanceProperties{@NamedTuple{ρw::Float64, Lsl::Float64, L::Float64}}}, WaterBalance{NoFlow, Nothing, CryoGrid.MaxDelta{Float64}, Nothing, WaterBalanceProperties{Float64, Float64, Float64}}, Nothing, Nothing}, ground3::Ground{MineralOrganic{Float64, Float64, Float64, ThermalProperties{@NamedTuple{kh_w::Float64, kh_i::Float64, kh_a::Float64, ch_w::Float64, ch_i::Float64, ch_a::Float64, kh_m::Float64, kh_o::Float64, ch_m::Float64, ch_o::Float64}}, HydraulicProperties{@NamedTuple{kw_sat::Float64, fieldcapacity::Float64}}}, HeatBalance{FreeWater, CryoGrid.Heat.Diffusion1D{:H, typeof(CryoGrid.Heat.quadratic_parallel_conductivity), typeof(CryoGrid.Heat.weighted_average_heatcapacity)}, CryoGrid.MaxDelta{Int64}, CryoGrid.Heat.HeatBalanceProperties{@NamedTuple{ρw::Float64, Lsl::Float64, L::Float64}}}, WaterBalance{NoFlow, Nothing, CryoGrid.MaxDelta{Float64}, Nothing, WaterBalanceProperties{Float64, Float64, Float64}}, Nothing, Nothing}, ground4::Ground{MineralOrganic{Float64, Float64, Float64, ThermalProperties{@NamedTuple{kh_w::Float64, kh_i::Float64, kh_a::Float64, ch_w::Float64, ch_i::Float64, ch_a::Float64, kh_m::Float64, kh_o::Float64, ch_m::Float64, ch_o::Float64}}, HydraulicProperties{@NamedTuple{kw_sat::Float64, fieldcapacity::Float64}}}, HeatBalance{FreeWater, CryoGrid.Heat.Diffusion1D{:H, typeof(CryoGrid.Heat.quadratic_parallel_conductivity), typeof(CryoGrid.Heat.weighted_average_heatcapacity)}, CryoGrid.MaxDelta{Int64}, CryoGrid.Heat.HeatBalanceProperties{@NamedTuple{ρw::Float64, Lsl::Float64, L::Float64}}}, WaterBalance{NoFlow, Nothing, CryoGrid.MaxDelta{Float64}, Nothing, WaterBalanceProperties{Float64, Float64, Float64}}, Nothing, Nothing}, ground5::Ground{MineralOrganic{Float64, Float64, Float64, ThermalProperties{@NamedTuple{kh_w::Float64, kh_i::Float64, kh_a::Float64, ch_w::Float64, ch_i::Float64, ch_a::Float64, kh_m::Float64, kh_o::Float64, ch_m::Float64, ch_o::Float64}}, HydraulicProperties{@NamedTuple{kw_sat::Float64, fieldcapacity::Float64}}}, HeatBalance{FreeWater, CryoGrid.Heat.Diffusion1D{:H, typeof(CryoGrid.Heat.quadratic_parallel_conductivity), typeof(CryoGrid.Heat.weighted_average_heatcapacity)}, CryoGrid.MaxDelta{Int64}, CryoGrid.Heat.HeatBalanceProperties{@NamedTuple{ρw::Float64, Lsl::Float64, L::Float64}}}, WaterBalance{NoFlow, Nothing, CryoGrid.MaxDelta{Float64}, Nothing, WaterBalanceProperties{Float64, Float64, Float64}}, Nothing, Nothing}, bottom::Bottom{GeothermalHeatFlux{Float64}}}, NTuple{8, Float64}}

define time span, 2 years + 3 months

tspan = (DateTime(2010,9,30), DateTime(2012,9,30))
u0, du0 = @time initialcondition!(tile, tspan)
prob = CryoGridProblem(tile, u0, tspan, saveat=3*3600.0, savevars=(:T, :top => (:T_ub), :snowpack => (:dsn,)))

solve full tspan with forward Euler and initial timestep of 5 minutes

@info "Running model ..."
sol = @time solve(prob, Euler(), dt=300.0);
out = CryoGridOutput(sol)

Plot it!

using Plots: plot, plot!, heatmap, cgrad, Measures
zs = [1,10,20,30,50,100,200,500]u"cm"
cg = cgrad(:copper,rev=true);
plot(ustrip(out.T[Z(Near(zs))]), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature (°C)", leg=false, dpi=150)
plot!(ustrip(out.T[1,:]), color=:darkgray, ylabel="Temperature (°C)", leg=false, dpi=150)
plt1 = plot!(ustrip.(out.top.T_ub), color=:skyblue, linestyle=:dash, alpha=0.5, leg=false, dpi=150)

Plot snow water equivalent and depth:

plot(ustrip(out.snowpack.swe), ylabel="Depth (m)", label="Snow water equivalent", dpi=150)
plt2 = plot!(ustrip.(out.snowpack.dsn), label="Snow depth", ylabel="Depth (m)", legendtitle=nothing, dpi=150)
plot(plt1, plt2, size=(1600,700), margins=5*Measures.mm)

Temperature heatmap:

T_sub = out.T[Z(Between(0.0u"m",10.0u"m"))]
heatmap(T_sub, yflip=true, size=(1200,600), dpi=150)

Thaw depth:

td = Diagnostics.thawdepth(out.T[Z(Where(>=(0.0u"m")))])
plot(td, yflip=true, ylabel="Thaw depth (m)", size=(1200,600))

This page was generated using Literate.jl.