Overview

Setting up a model

At the highest level, a model in CryoGrid.jl is defined by one or more Tiles each consisting of a Grid and a Stratigraphy, constructed top-down from individual Layers, each of which has one or more Processes. Each layer in the Stratigraphy is assigned a depth, which then aligns it with the Grid. All models must consist of at least three layers/nodes: Top and Bottom layers with corresponding boundary conditions, as well as one or more SubSurface layers. Here we define a simple three-layer model (or one-layer, exlcuding the boundaries) with a single sub-surface process, i.e. HeatBalance (heat conduction):

# ... load forcings, set up profiles, etc.
# see examples/heat_vgfc_seb_saoylov_custom.jl for more details
strat = Stratigraphy(
    -2.0u"m" => Top(SurfaceEnergyBalance(Tair,pr,q,wind,Lin,Sin,z)),
    0.0u"m" => Ground(soilprofile, HeatBalance(:H; freezecurve=DallAmico())),
    1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"J/s/m^2"))
);
grid = CryoGrid.Presets.DefaultGrid_5cm
# define initial conditions for temperature using a given profile;
# The default initializer linearly interpolates between profile points.
initT = initializer(:T, tempprofile)
tile = Tile(strat, grid, initT);

This model can then be used to construct a CryoGridProblem:

tspan = (DateTime(2010,10,30),DateTime(2011,10,30))
p = parameters(tile)
u0 = initialcondition!(tile, tspan, p, initT)
prob = CryoGridProblem(tile, u0, tspan, p, saveat=24*3600.0, savevars=(:T,)) # produces an ODEProblem with problem type CryoGridODEProblem

It can then be solved/integrated using the solve function (from DiffEqBase and OrdinaryDiffEq):

# solve and construct CryoGridOutput from solution
sol = @time solve(prob, saveat=24*3600.0, progress=true);
out = CryoGridOutput(sol)

The resulting CryoGridOutput type provides DimArrays containing the model outputs over time and space:

julia> out.T
278×366 DimArray{Float64,2} with dimensions: 
  Z: Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}[0.01 m, 0.03 m, …, 850.0 m, 950.0 m] Sampled: Ordered Irregular Points,
  Ti (Time): DateTime[2010-10-30T00:00:00, …, 2011-10-30T00:00:00] Sampled: Ordered Irregular Points

Defining model behavior

Notice that, in the example above, it is types such as Ground, HeatBalance, DallAmico, etc. that specify which components the model should use. These components are defined by adding method dispatches to the CryoGrid interface methods. State variables are declared via the variables method, e.g:

variables(soil::Soil, heat::HeatBalance{<:EnthalpyBased}) = (
    Prognostic(:H, OnGrid(Cells), u"J/m^3"),
    Diagnostic(:T, OnGrid(Cells), u"°C"),
    Diagnostic(:C, OnGrid(Cells), u"J//K*/m^3"),
    Diagnostic(:∂H∂T, OnGrid(Cells), u"J/K/m^3"),
    Diagnostic(:k, OnGrid(Edges), u"W/m/K"),
    Diagnostic(:kc, OnGrid(Cells), u"W//m/K"),
)

When the HeatBalance process is assigned to a Soil layer, Tile will invoke this method and create state variables corresponding to each Var. Prognostic variables are assigned derivatives (in this case, dH, since H is the prognostic state variable) and integrated over time. Diagnostic variables provide in-place caches for derived/intermediary state variables.

Each variable definition consists of a name (a Julia Symbol), a type, and a shape. For variables discretized on the grid, the shape is specified by OnGrid, which will generate an array of the appropriate size when the model is compiled. The arguments Cells and Edges specify whether the variable should be defined on the grid cells or edges respecitvely.

The real work finally happens in computediagnostic! and computefluxes!, the latter of which should be used to compute the time derivatives (here dH). interact! defines the behavior at the boundaries and should be used to compute the derivatives (and any other necessary values) at the interface between layers.

We can take as an example the implementation of computefluxes! for enthalpy-based heat conduction (note that jH is a diagnostic variable representing the energy flux over each cell edge):

function CryoGrid.computefluxes!(::SubSurface, ::HeatBalance{<:FreezeCurve,<:EnthalpyBased}, state)
    Δk = Δ(state.grid) # cell sizes
    ΔT = Δ(cells(state.grid)) # midpoint distances
    # compute internal fluxes and non-linear diffusion assuming boundary fluxes have been set
    nonlineardiffusion!(state.dH, state.jH, state.T, ΔT, state.k, Δk)
    return nothing
end

Note that state is (typically) of type LayerState with properties corresponding to the state variables declared by the variables function for Soil and HeatBalance. Additionally, output arrays for the time derivatives are provided (here dH), as well as the current timestep, layer boundary depths, and variable grids (accessible via state.t, and state.grid respectively). Note that state will also contain other variables declared on this Soil layer by other SubSurfaceProcesses, allowing for implicit coupling between processes where appropriate.