Soil heat with free water freeze curve

In this example, we use the preset SoilHeatTile to construct a Tile consisting of a soil column with heat conduction forced using air temperatures from Samoylov Island. The enthalpy-based HeatBalance process defaults to the so-called "free water" freezing characteristic which assumes that water only freezes and thaws at a melting temperature of 0°C.

First we load the built-in forcing file from Nitzbon et al. 2020 (CryoGrid 3). Note that this will download the forcing files from the AWI NextCloud if they are not already present in the input/ folder.

using CryoGrid
forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);

We use the provided default soil and temperature profiles for Samoylov.

soilprofile, tempprofile = CryoGrid.Presets.SamoylovDefault;

We choose the default grid discretization with 5 cm spacing at the surface.

grid = CryoGrid.Presets.DefaultGrid_5cm;

We construct a state variable initializer for temperature T from the temperature profile.

initT = initializer(:T, tempprofile)
tile = CryoGrid.Presets.SoilHeatTile(
    :H,
    TemperatureBC(forcings.Tair, NFactor(nf=0.6)),
    GeothermalHeatFlux(0.053u"W/m^2"),
    soilprofile,
    initT;
    grid=grid
);

Here we define the time span:

tspan = (DateTime(2010,10,30),DateTime(2011,10,30));

Evaluate the initial condition

u0, du0 = initialcondition!(tile, tspan);

Here we construct a CryoGridProblem with tile, initial condition, and timespan.

prob = CryoGridProblem(tile, u0, tspan, saveat=24*3600.0, savevars=(:T,));

Solve the configured problem with the built-in forward Euler method. note that, due to compile time, this may take 1-2 minutes when executed in a fresh Julia session. Subsequent solves will be much faster.

sol = @time solve(prob);
out = CryoGridOutput(sol)

CryoGridOutput with 2921 time steps (2010-10-30T00:00:00 to 2011-10-30T00:00:00) and 2 variables:
    H => DimArray of Quantity{Float64, 𝐌 𝐋^-1 𝐓^-2, Unitful.FreeUnits{(J, m^-3), 𝐌 𝐋^-1 𝐓^-2, nothing}} with dimensions (218, 2921)
    T => DimArray of Quantity{Float64, 𝚯, Unitful.FreeUnits{(K,), 𝚯, Unitful.Affine{-5463//20}}} with dimensions (218, 2921)

Now we plot the reuslts!

import Plots
zs = [1,10,20,30,50,100,200,500,1000]u"cm"
cg = Plots.cgrad(:copper,rev=true);
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)

This page was generated using Literate.jl.