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, :snowpack, :ground1, :ground2, :ground3, :ground4, :ground5, :bottom), Tuple{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{CryoGrid.Snow.Bulk{CryoGrid.Snow.ConstantDensity{Float64}, Float64, CryoGrid.Snow.SnowThermalProperties{CryoGrid.Snow.SturmQuadratic{ClosedInterval{Float64}, NamedTuple{(:a, :b, :c), Tuple{Float64, Float64, Float64}}, NamedTuple{(:a, :b, :c), Tuple{Float64, Float64, Float64}}}, ThermalProperties{NamedTuple{(:kh_w, :kh_i, :kh_a, :ch_w, :ch_i, :ch_a), NTuple{6, Float64}}}}, HydraulicProperties{NamedTuple{(:kw_sat,), Tuple{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, :Lsl, :L), Tuple{Float64, Float64, Float64}}}}, WaterBalance{NoFlow, Nothing, CryoGrid.MaxDelta{Float64}, Nothing, WaterBalanceProperties{Float64, Float64, Float64}}, Nothing}, Ground{MineralOrganic{Float64, Float64, Float64, ThermalProperties{NamedTuple{(:kh_w, :kh_i, :kh_a, :ch_w, :ch_i, :ch_a, :kh_m, :kh_o, :ch_m, :ch_o), NTuple{10, Float64}}}, HydraulicProperties{NamedTuple{(:kw_sat, :fieldcapacity), Tuple{Float64, 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, :Lsl, :L), Tuple{Float64, Float64, Float64}}}}, WaterBalance{NoFlow, Nothing, CryoGrid.MaxDelta{Float64}, Nothing, WaterBalanceProperties{Float64, Float64, Float64}}, Nothing, Nothing}, Ground{MineralOrganic{Float64, Float64, Float64, ThermalProperties{NamedTuple{(:kh_w, :kh_i, :kh_a, :ch_w, :ch_i, :ch_a, :kh_m, :kh_o, :ch_m, :ch_o), NTuple{10, Float64}}}, HydraulicProperties{NamedTuple{(:kw_sat, :fieldcapacity), Tuple{Float64, 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, :Lsl, :L), Tuple{Float64, Float64, Float64}}}}, WaterBalance{NoFlow, Nothing, CryoGrid.MaxDelta{Float64}, Nothing, WaterBalanceProperties{Float64, Float64, Float64}}, Nothing, Nothing}, Ground{MineralOrganic{Float64, Float64, Float64, ThermalProperties{NamedTuple{(:kh_w, :kh_i, :kh_a, :ch_w, :ch_i, :ch_a, :kh_m, :kh_o, :ch_m, :ch_o), NTuple{10, Float64}}}, HydraulicProperties{NamedTuple{(:kw_sat, :fieldcapacity), Tuple{Float64, 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, :Lsl, :L), Tuple{Float64, Float64, Float64}}}}, WaterBalance{NoFlow, Nothing, CryoGrid.MaxDelta{Float64}, Nothing, WaterBalanceProperties{Float64, Float64, Float64}}, Nothing, Nothing}, Ground{MineralOrganic{Float64, Float64, Float64, ThermalProperties{NamedTuple{(:kh_w, :kh_i, :kh_a, :ch_w, :ch_i, :ch_a, :kh_m, :kh_o, :ch_m, :ch_o), NTuple{10, Float64}}}, HydraulicProperties{NamedTuple{(:kw_sat, :fieldcapacity), Tuple{Float64, 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, :Lsl, :L), Tuple{Float64, Float64, Float64}}}}, WaterBalance{NoFlow, Nothing, CryoGrid.MaxDelta{Float64}, Nothing, WaterBalanceProperties{Float64, Float64, Float64}}, Nothing, Nothing}, Ground{MineralOrganic{Float64, Float64, Float64, ThermalProperties{NamedTuple{(:kh_w, :kh_i, :kh_a, :ch_w, :ch_i, :ch_a, :kh_m, :kh_o, :ch_m, :ch_o), NTuple{10, Float64}}}, HydraulicProperties{NamedTuple{(:kw_sat, :fieldcapacity), Tuple{Float64, 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, :Lsl, :L), Tuple{Float64, Float64, Float64}}}}, WaterBalance{NoFlow, Nothing, CryoGrid.MaxDelta{Float64}, Nothing, WaterBalanceProperties{Float64, Float64, Float64}}, Nothing, Nothing}, 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,)))

CryoGridProblem with uType ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(top = ViewAxis(1:1, Axis(runoff = 1:1,)), snowpack = ViewAxis(2:3, Axis(swe = 1:1, Δz = 2:2)), sat = 4:2:440, H = 5:2:441)}}} and tType Float64. In-place: true
timespan: (6.3453024e10, 6.35161824e10)
u0: ComponentVector{Float64}(top = (runoff = [0.0]), snowpack = (swe = [0.0], Δz = [0.0]), sat = [0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], H = [0.0, -1.995e6, -1.995e6, -1.945e6, -1.945e6, -1.945e6, -1.945e6, -1.945e6, -1.945e6, -2.00125e6  …  -1.773e7, -1.5628666666666668e7, -1.1426000000000002e7, -7.223333333333333e6, -3.0206666666666674e6, 1.01796e8, 1.0747066666666667e8, 1.1314533333333333e8, 1.1882e8, 1.2449466666666666e8])

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)

CryoGridOutput with 5849 time steps (2010-09-30T00:00:00 to 2012-09-30T00:00:00) and 8 variables:
    sat => DimArray of Float64 with dimensions (219, 5849)
    H => DimArray of Quantity{Float64, 𝐌 𝐋^-1 𝐓^-2, Unitful.FreeUnits{(J, m^-3), 𝐌 𝐋^-1 𝐓^-2, nothing}} with dimensions (219, 5849)
    T => DimArray of Quantity{Float64, 𝚯, Unitful.FreeUnits{(K,), 𝚯, Unitful.Affine{-5463//20}}} with dimensions (219, 5849)
    top => 
        runoff => DimArray of Quantity{Float64, 𝐋^3, Unitful.FreeUnits{(m^3,), 𝐋^3, nothing}} with dimensions (1, 5849)
        T_ub => DimArray of Quantity{Float64, 𝚯, Unitful.FreeUnits{(K,), 𝚯, nothing}} with dimensions (1, 5849)
    snowpack => 
        swe => DimArray of Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}} with dimensions (1, 5849)
        Δz => DimArray of Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}} with dimensions (1, 5849)
        dsn => DimArray of Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}} with dimensions (1, 5849)

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.