Soil heat with SFCC and constant BCs

In this example, we use the preset SoilHeatTile to construct a Tile consisting of a soil column with heat conduction and zero-flux boundary conditions. This is a useful test case for checking energy conservation since we can guarantee that no energy is being added or removed at the boundaries.

using CryoGrid

Select default grid and initial temperature profile.

grid = CryoGrid.Presets.DefaultGrid_10cm
tempprofile = TemperatureProfile(
    0.0u"m" => -20.0u"°C",
    1000.0u"m" => 10.0u"°C",

Here we use a simple single layer model with default soil parameters (50% porosity, no organic).

soilprofile = SoilProfile(
    0.0u"m" => MineralOrganic()

Here we specify the soil freezing characteristic curve (SFCC) formulation of Painter and Karra (2014). The van Genuchten parameters α=0.1 and n=1.8 correspond to a silty soil.

sfcc = PainterKarra(swrc=VanGenuchten(α=1.0, n=2.0));

import Plots
Plots.plot(-2.0u"°C":0.01u"K":0.0u"°C", sfcc)
Example block output

Enthalpy form of the heat transfer operator (i.e. prognostic :H). In this case, this is equivalent to the shorthand SoilHeatTile(:H, ...). However, it's worth demonstrating how the operator can be explicitly specified.

heatop = Heat.Diffusion1D(:H)
initT = initializer(:T, tempprofile)
tile = CryoGrid.Presets.SoilHeatTile(
    # 10 W/m^2 in and out, i.e. net zero flux
    ConstantBC(HeatBalance, CryoGrid.Neumann, 10.0u"W/m^2"),
    ConstantBC(HeatBalance, CryoGrid.Neumann, 10.0u"W/m^2"),

Define the simulation time span.

tspan = (DateTime(2010,1,1),DateTime(2010,12,31))
u0, du0 = initialcondition!(tile, tspan);

Construct and solve the CryoGridProblem using the built-in forward Euler integrator.

prob = CryoGridProblem(tile, u0, tspan, saveat=3600.0, savevars=(:T,:H))
sol = @time solve(prob, dtmax=900.0);
out = CryoGridOutput(sol)
CryoGridOutput with 2913 time steps (2010-01-01T00:00:00 to 2010-12-31T00:00:00) and 2 variables:
    H => DimArray of Quantity{Float64, 𝐌 𝐋^-1 𝐓^-2, Unitful.FreeUnits{(J, m^-3), 𝐌 𝐋^-1 𝐓^-2, nothing}} with dimensions (178, 2913)
    T => DimArray of Quantity{Float64, 𝚯, Unitful.FreeUnits{(K,), 𝚯, Unitful.Affine{-5463//20}}} with dimensions (178, 2913)

Compute total energy balance error.

Htot = Diagnostics.integrate(out.H, grid)
mass_balance_error = Htot[end] - Htot[1]
0.0006561279296875 J

Plot it!

import Plots
zs = [1,45,95,195,495,795,995]u"cm"
Diagnostics.plot_at_depths(:T, out, zs, ylabel="Temperature", leg=false, size=(800,500), dpi=150)
Example block output

Plot the energy balance error over time.

Plots.plot(uconvert.(u"MJ", Htot .- Htot[1]), title="", ylabel="Energy balance error")
Example block output

This page was generated using Literate.jl.