ConstantForcing{unit,T} <: Forcing{unit,T}

Simple Forcing type that just returns a constant value. This type is primarily intended for writing tests but could still be used in applications if necessary.

CryoGridOutput{TSol}

Helper type that stores the raw output from a CryoGrid run along with DimArray views of all logged variables. CryoGridOutput overrides Base.getproperty to allow for direct dot-syntax access of state variables. For example, if your model has a grid variable named T, out.T returns a DimArray with indexed time and depth axes. For OrdinaryDiffEq.jl outputs, the ODESolution can be accessed via out.sol, or for convenience, the continuous solution at time t can be computed via out(t) which is equivalent to withaxes(out.sol(t)).

Evaluates the continuous solution at time t.

CryoGridOutput(sol::TSol, tspan::NTuple{2,Float64}=(-Inf,Inf)) where {TSol<:SciMLBase.AbstractODESolution}

Constructs a CryoGridOutput from the given ODESolution. Optional argument tspan restricts the time span of the output.

CryoGridParams{T,TM} <: DenseArray{T,1}

Wraps a ModelParameters.Model parameter handler for CryoGrid types. It is recommended not to use this type directly in math or linear algebra operations but rather to use Base.values to obtain a normal array of parameter values.

CryoGridParams(obj)

Constructs a modelParameters.Model wrapped with CryoGridParams from obj. If full_metadata is true, additonal fields for nested Parameterization types will be added.

DateNumber <: Real

Represents a MATLAB DateNumber type.

  Forcing{unit,T}

Abstract type representing a generic external forcing term.

Represents an externally specified format for forcing inputs. IO functions should dispatch on specific types T<:ForcingFormat that they implement.

JSON forcing input format (from CryoGridLite) with specified version indicator.

NetCDF forcing input format.

  Forcings{names,TF,TMeta}

Generic container for forcing types with optional metadata.

  InterpolatedForcing{unit,T,TI}

Forcing data provided by a discrete time series of data.

Get interpolated forcing value at t seconds from t0.

NCDFormat

Base type for NetCDF based forcing formats.

NCDGeneric <: NCDFormat

Generic NetCDF forcing format.

NCDTopoPyScale

TopoPyScale NetCDF forcing format.

JSON parameter input format (from CryoGridLite). Not yet implemented.

YAML parameter input format matching that of the CryoGrid community model. Not yet implemented.

Resource{T}

Simple representation of a local or remote resource with a name, type, and url. The name and type will be used to specify a file name and suffix respectively.

  TimeVaryingForcing{unit,T} <: Forcing{unit,T}

Simple Forcing type that evaluates a given function f(t) of where t is the current time.

  TransformedForcing(transformed_unit,orig_unit,TF,T,TO<:Forcing{orig_unit,T}) <: Forcing{transformed_unit,T}

Wraps another Forcing and applies an arbitrary transformation f when evaluated. The transformed unit can either be automatically determined on construction or specified directly.

fetch(resource::Resource, dir::String)

Downloads the remote resource from the location given by its URL to directory dir.

loadforcings(filename::String)::Forcings
loadforcings(resource::Resource; outdir=DEFAULT_FORCINGS_DIR)::Forcings
loadforcings([format::ForcingFormat], filename::String; outdir=DEFAULT_FORCINGS_DIR)::Forcings

Loads forcing data from the given file according to the format specified by format. By default, the forcing format is automatically detected via detectformat. Returns a Forcings struct containing all forcing data and metadata

time_derivative_forcing(
    f::InterpolatedForcing{unit},
    new_name::Symbol;
    interp=Numerics.Linear()
) where {unit}

Computes the finite difference time derivative of the given InterpolatedForcing time series and returns a new forcing with units [unit].s⁻¹

datenum(d::Dates.DateTime)

Converts a Julia DateTime to a MATLAB style DateNumber. MATLAB represents time as DateNumber, a double precision floating point number being the the number of days since January 0, 0000

todatetime(d::DateNumber)

Converts a MATLAB DateNumber to a Julia DateTime.

Utils

Common utility functions, constants, and macros used throughout the CryoGrid.jl codebase.

Optional{T}

Option type, alias for Union{Nothing,T}. Intended for cases where a property/field may be optionally provided and is nothing otherwise.

Named{name,T}

Wraps a value of type T with a name type parameter.

NamedTupleWrapper

Base type for container types that hold a NamedTuple of arbitrary field values. NamedTupleWrapper provides dispatches for getproperty and propertynames that forward property name queries to the NamedTuple container. Subtypes are by default assumed to have a field named values that corresponds to the NamedTuple container, but this can be overriden by providing a dispatch for Base.values.

applyunits(u::Unitful.Units, x::Number)

Conditionally applies unit u to x if and only if x is a unit-free quantity. If x is a unitful quantity, asserts that the unit matches u.

convert_t(t::DateTime)
convert_t(t::Float64)

Convenience method for converting between Dates.DateTime and solver time.

convert_tspan(tspan::Tuple{DateTime,DateTime})
convert_tspan(tspan::Tuple{Float64,Float64})

Convenience method for converting between Dates.DateTime and solver time.

fastiterate(f!::F, iters::NTuple{N,Any}...) where {F,N}

Same as fastmap but simply invokes f! on each argument set without constructing a tuple.

fastmap(f::F, iter::NTuple{N,Any}...) where {F,N}

Same as map for NTuples but with guaranteed type stability. fastmap is a @generated function which unrolls calls to f into a loop-free tuple construction expression.

ffill!(x::AbstractVector{T}) where {E,T<:Union{Missing,E}}

Forward fills missing values in vector x.

getscalar(x)
getscalar(x, i)

Helper method for generalizing between arrays and scalars. Without an index, retrieves the first element of x if x is an array, otherwise simply returning x. If an index i, is specified, returns the ith value of x if x is an array, or x otherwise. Note that this method is not strictly necessary since Julia allows for scalar quantities to be accessed at the first index like an array; however, the point is to make it expliclty clear in scalar-typed code that a state variable is treated as such and is not a vector valued quantity.

groupby(f, xs)

Simple implementation of group-by operations which groups together items in a collection by the return-value of f. Copied from Lazy.jl.

pstrip(obj; keep_units=false)

Strips Param types and units from obj. If keep_units=true, then Param types will be stripped but units preserved.

tuplejoin([x, y], z...)

Concatenates one or more tuples together; should generally be type stable.

ModelParameters.stripunits(obj)

Additional override for stripunits which reconstructs obj with all fields that have unitful quantity types converted to base SI units and then stripped to be unit free.

Similar to Unitful.@u_str (i.e. u"kg") but produces the type of the quantity rather than the instance. NOT conditional on debug mode.

Similar to Unitful.@u_str (i.e. u"kg") but produces the type of the unit rather than the instance. NOT conditional on debug mode.

properties(expr)

Defines a new container type for properties which subtypes NamedTupleWrapper. Usage:

@properties MyProperties(
    prop1 = 1.0,
    prop2 = 2.0,
)

outputs:

struct MyProperties{TV} <: NamedTupleWrapper
    values::TV
    function HydraulicProperties(;
        prop1 = 1.0,
        prop2 = 2.0,
        additional_kwargs...
    )
        props = (;prop1, prop2, additional_kwargs...)
        return new{typeof(props)}(props)
    end
end

@pstrip(expr, kwargs...)

Convenience macro for Utils.pstrip; equivalent to pstrip(expr, kwargs...).

Convenience macro for setting scalar (single-element) arrays/vectors. It turns an expression of the form: a.b = ... into a.b[1] = ...

This is primarily intended for code clarity, i.e to clearly discern scalar and non-scalar values.

sym_str(val)

Convenience macro, sym"val", for creating a Symbol from val. Equivalent to Symbol(val). Use in situations where normal Julia :val syntax is not possible, e.g. sym"1" instead of Symbol(1) or sym"my.var" instead of Symbol("my.var").

Prepends expr with Threads.@threads if and only if Threads.nthreads() > 1, thus avoiding the overhead of @threads when running in single-threaded mode.

Credit to @ranocha (Hendrik Ranocha) https://discourse.julialang.org/t/overhead-of-threads-threads/53964/22

active_layer_thickness(T::AbstractDimArray{<:TempQuantity})

Computes active layer thickness annually. The active layer thickness is defined here as the maximum thaw depth throughout the calendar year. Assumes T to have dimensions Ti (time) and Z (depth) in any order.

build_dummy_state(grid::Grid, layer::Layer; t=0.0, with_units=true)

Collects variables defined on layer and initializes a dummy state named tuple with all state variables initialized from the given grid. This is intended to be used for unit tests and debugging in order to avoid the full-fledged construction of a Tile/Stratigraphy and associated state types.

computejac(tile::Tile, u, p, t)

Helper function that computes the Jacobian of the given tile at u with parameters p and time t.

integrate(X::AbstractDimArray, grid::Grid{Edges}; upper_limit=0u"m", lower_limit=10u"m")

Integrates the quantity X over the given grid, which is assumed to be spatially alligned, i.e. length(grid) == length(dims(X,Z)) + 1 and cells(grid) .≈ dims(X,Z) are necessary preconditions.

mean_annual_ground_temperature(T::AbstractDimArray; upper_limit=0u"m", lower_limit=10u"m")

Computes mean annual ground temperature between upper_limit and lower_limit. Assumes T to have dimensions Ti (time) and Z (depth) in any order.

permafrostbase(T::AbstractDimArray{<:TempQuantity})

Computes depth of permafrost base for all years, i.e. the closest depth to the "bottom" at which the maximum annual temperature is strictly < 0°C. Assumes T to have dimensions Ti (time) and Z (depth) in any order.

permafrosttable(T::AbstractDimArray{<:TempQuantity})

Computes depth of permafrost table for all years, i.e. the closest depth to the surface at which the maximum annual temperature is strictly less than Tmelt. Assumes T to have dimensions Ti (time) and Z (depth) in any order.

plot_alt(out::CryoGridOutput; kwargs...)

Plots the active layer thickness from the given CryoGridOutput.

plot_at_depths(var::Symbol, out::CryoGridOutput, depths::AbstractArray; kwargs...)

Plots a time series of var from out at the given depths.

plot_temperature_heatmap(out::CryoGridOutput; kwargs...)

Plots a heatmap of temperatures along with the zero degree isotherm.

spinup(setup::Tile, tspan::NTuple{2,DateTime}, p, tol, layername; kwargs...)

Implements a simple, iterative spin-up procedure. Runs the model specified by setup over tspan until the profile mean up to maxdepth over the whole time span changes only within the given tolerance tol. Returns the ODESolution generated by the final iteration.

thawdepth(T::AbstractDimArray{<:TempQuantity}; Tmelt=0.0u"°C")

Computes sub-grid thaw depth (a.k.a freezing front) from temperature at all time steps. The sub-grid depth of the zero degree isotherm is determined by linearly interpolating between grid cell temperatures in T. Note that T is assumed to have units °C and dimensions Ti (time) and Z (depth) in any order.

zero_annual_amplitude(T::AbstractDimArray{<:TempQuantity}; threshold=0.5u"K")

Computes annual depth of zero amplitude (where |max - min| < threshold) and returns the result for each year. Assumes T to have dimensions Ti (time) and Z (depth) in any order.

Analytical

Analytical formulation from Byun 1990.

Businger 1971

Høgstrøm 1988 (unstable conditions) SHEBA, Uttal et al., 2002, Grachev et al. 2007 (stable conditions)

Iterative

Westermann 2016, use state from last time step.

Numerical

Equations from Westermann 2016, but use nonlinear solver.

SEBInputs{TT1,TT2,TR,TP,TQ,TW,TZ,Tsurf}

Non-prognostic input variables to the SEB that are assumed to be known a priori.

SEBOutputs{TQs,TQg,TQnet,TSout,TLout,TL,TU}

Outputs of surface energy balance computation that include both the updated state as well as the net radiation Qnet and outgoing radiation components Sout and Lout.

SEBState{TQ,TL,TU,TIn}

Unknown state variables that constitute a candidate solution to the SEB.

SurfaceEnergyBalance{TSolution,TStabFun,TPara,F} <: BoundaryProcess{HeatBalance}

Surface energy balance upper boundary condition.

SurfaceWaterBalance{TR,TS} <: BoundaryProcess{Union{WaterBalance, SnowMassBalance}}

The SurfaceWaterBalance represents the closure of the water balance at the surface and acts as both a Neumann-type upper boundary condition for snow and water fluxes as well as an accountant for the water mass balance at the surface.

SurfaceWaterBalance(forcings::Forcings)

Automatically attempts to extract precipitation forcings from forcings.

SurfaceWaterBalance(precip::Forcing{u"m/s",T1}, Tair::Forcing{u"°C",T2}) where {T1,T2}

SWB constructor which derives rainfall and snowfall forcings from total precipitation and air temperature.

Type alias for WaterHeatBC{TSWB,TSEB} where {TSWB<:SurfaceWaterBalance,TSEB<:SurfaceEnergyBalance}.

Latent heat of evaporation/condensation of water in [J/kg] according to https://en.wikipedia.org/wiki/Latentheat#citenote-RYfit-11

Latent heat of sublimation/resublimation of water in [J/kg] accodring to https://en.wikipedia.org/wiki/Latentheat#citenote-RYfit-11

Obukhov length according to Monin-Obukhov theory, iterative determination as in CryoGrid3 / Westermann et al. 2016

• uses the turubulent fluxes Qe and Qh as well as the friction velocity of the previous time step

Obukhov length, analytical solution of Monin-Obukhov theory according to Byun 1990

• implicitly assumes the Businger 1971 stability functions
• only uses the current surface temperature as input

Latent heat flux, defined as positive if it is a flux towards the surface. Represents evapo(transpi)ration/condensation at positive surface temperatures and sublimation/resublimation at negative surface temperatures

Sensible heat flux, defined as positive if it is a flux towards the surface

Density of air at given tempeature and pressure

Saturation pressure of water/ice according to the empirical August-Roche-Magnus formula

surfaceproperties(::SurfaceEnergyBalance, ::SubSurface)

Retrieves the SurfaceProperties for the given SubSurface layer.

Friction velocity according to Monin-Obukhov theory

Integrated stability function for heat/water transport Businger 1971

Integrated stability function for heat/water transport Högström, 1988 SHEBA, Uttal et al., 2002, Grachev et al. 2007

Integrated stability function for momentum transport Busigner 1971

Integrated stability function for momentum transport Högström, 1988 SHEBA, Uttal et al., 2002, Grachev et al. 2007

Algebraic{name,S,T,units,domain} <: Var{name,S,T,units,domain}

Defines an algebraic (implicit) state variable.

BCKind

Trait that specifies the kind of boundary condition. This can be used to write generic implementations of interact! that are (relatively) agnostic to specific implementations of BoundaryProcess. A good example of this can be found in the boundaryflux method interface.

BCKind(::Type{T})

Can be overriden by BoundaryProcess types to indicate the type of boundary condition, e.g:

BCKind(::Type{BP}) = Dirichlet()

where BP is a BoundaryProcess that provides the boundary conditions.

Bottom{TProc} <: Layer

Generic "bottom" layer that marks the lower boundary of the subsurface grid.

BoundaryEffect

Base type for boundary "effects" which modify boundary conditions based on some given parameterization.

BoundaryProcess{T<:SubSurfaceProcess}

Abstract base type for boundary processes, i.e. processes that operate at the boundaries of the subsurface. A BoundaryProcess represents the boundary conditions of one or more SubSurfaceProcesses but may include its own diagnostic (or even prognostic) variables, if necessary.

CFL{Tmax<:MaxDelta}

Courant-Fredrichs-Lewy condition (where defined) with the given Courant number and embedded MaxDelta condition.

CGEuler <: SciMLBase.AbstractODEAlgorithm

Simple, lightweight implementation of the forward Euler integration algorithm. Does not include support for fancier features such as interpolation, adaptive timestepping, or event handling. In order to get these features, you must install the OrdinaryDiffEq package.

CompositeBoundaryProcess{B1,B2,F,S} <: BoundaryProcess

Represents a composition of two boundary processes, B1 and B2, via an operator F. A typical use case is combining ConstantBC with a forcing-driven boundary process to scale or shift the forcing.

ConstantBC{P,S,T} <: BoundaryProcess{P}

Constant boundary condition (of any type/unit) specified by value.

ContinuousTrigger

Base type for continuous trigger flags, Increasing and Decreasing, which indicate an upcrossing of the function root (negative to positive) and a downcrossing (positive to negative) respectively. Both subtypes have a field idx which, for GridContinuousEvent is set to the grid cell index for which the event was triggered.

Coupled2{P1,P2} = CoupledProcesses{Tuple{T1,T2}} where {T1,T2}

Type alias for coupled processes, i.e. CoupledProcesses{Tuple{P1,P2}}. Coupled provides a simple mechanism for defining new behaviors on multi-processes systems.

CoupledProcesses{TProcs} <: Process

Represents an explicitly or implicitly coupled system of processes. TProcs is always a Tuple of other processes.

CryoGridProblem{iip,Tu,Tt,Tp,TT,Tsv,Tsf,Tcb,Tdf,Tkw} <: SciMLBase.AbstractODEProblem{Tu,Tt,iip}

Represents a CryoGrid discretized PDE forward model configuration using the SciMLBase/DiffEqBase problem interface.

CryoGridProblem(
    tile::Tile,
    u0::ComponentVector,
    tspan::NTuple{2,Float64},
    p=nothing;
    saveat=3600.0,
    savevars=(),
    save_everystep=false,
    save_start=true,
    save_end=true,
    step_limiter=timestep,
    safety_factor=1,
    max_step=true,
    callback=nothing,
    isoutofdomain=Tiles.domain(tile),
    specialization=SciMLBase.AutoSpecialize,
    function_kwargs=(),
    prob_kwargs...
)

Constructor for CryoGridProblem that automatically generates all necessary callbacks.

CryoGridProblem(tile::Tile, u0::ComponentVector, tspan::NTuple{2,DateTime}, args...;kwargs...)

DVar{dname,name,S,T,units,domain} <: Var{dname,S,T,units,domain}

Defines a delta/differential term du for variable u, which is the time-derivative/divergence for prognostic variables or the residual for algebraic variables.

Trigger for the criterion function crossing zero from positive to negative.

Diagnostic{name,S,T,units,domain} <: Var{name,S,T,units,domain}

Defines a diagnostic variable which is allocated and cached per timestep but not integrated over time.

BCKind instance for Dirichlet boundary conditions.

DynamicParameterization

Base type for dynamic parameterizations whose values may be time dependent. Implementations of DynamicParameterization should be callable structs that accept a single parameter t which corresponds to the current timestamp in seconds.

Event{name}

Base type for integration "events" (i.e. callbacks). name should be a Symbol or type which identifies the event for the purposes of dispatch on criterion and trigger!.

FunctionBC{P,S,F} <: BoundaryProcess{P}

Generic boundary condition type that invokes a user supplied function func(state) where state is the layer state associated with the boundary layer.

FunctionInitializer{varname,F} <: VarInitializer{varname}

Initializes a scalar or vector-valued state variable using an arbitrary method/function.

Base type for grid alignment, i.e. either "edges" or "cells".

Initializer

Base type for all state initializers.

InterpInitializer{varname,P,I,E} <: VarInitializer{varname}

Initializer for on-grid variables that takes a Profile as initial values and interpolates along the model grid. The interpolation mode is linear by default, but can also be any other Gridded interpolation mode supported by Interpolations.jl.

Layer

Abstract base type for all layers.

MaxDelta{T}

Allow a maximum change of Δmax in the integrated quantity.

BCKind instance for Neumann boundary conditions.

Parameterization

Base type for generic parameterizations of processes or components.

PeriodicBC{P,S,T1,T2,T3,T4} <: BoundaryProcess{P}

Periodic boundary condition (of any type/unit) specified by period, amplitude, and phaseshift.

Abstract base type for all dynamical processes.

Prognostic{name,S,T,units,domain} <: Var{name,S,T,units,domain}

Defines a prognostic (time-integrated) state variable.

SubSurface <: Layer

Abstract base type for layers in the stratigraphy, e.g. soil, snow, pond, etc.

SubSurfaceProcess <: Process

Abstract base type for subsurface processes, i.e. processes that operate at or below the surface, such as heat conduction, water infiltration, etc.

Top{TProc} <: Layer

Generic "top" layer that marks the upper boundary of the subsurface grid.

Var{name,S<:VarDim,T,units,domain}

Base type for symbolic state variables in the model.

VarInitializer{varname}

Base type for state variable initializers. Initializers of this type are dedicated to a specific state variable with name varname and are run before all other initializers.

(p::DiffEqCallbacks.StepsizeLimiterAffect{typeof(CryoGrid.timestep)})(integrator)

Custom implementation of StepsizeLimiterAffect function for CryoGrid.timestep that invokes the timestep function with tile,du,u,p,t as arguments.

Coupled(ps::SubSurfaceProcess...)

Constructs a composite/coupled process from one or more subsurface processes. Alias for CoupledProcesses(ps...).

Coupled(types::Type{<:Process}...)

Convenince method which constructs a CoupledProcesses type corresponding to each type in types, e.g:

Coupled(SnowMassBalance, HeatBalance) = CoupledProcesses{Tuple{T1,T2}} where {T1<:SnowMassBalance, T2<:HeatBalance}

also equivalent to Coupled2{<:SnowMassBalance,<:HeatBalance}.

boundaryflux(bc::BoundaryProcess, b::Union{Top,Bottom}, p::SubSurfaceProcess, sub::SubSurface, sbc, ssub)
boundaryflux(s::BCKind, bc::BoundaryProcess, b::Union{Top,Bottom}, p::SubSurfaceProcess, sub::SubSurface, sbc, ssub)

Computes the flux dH/dt at the boundary layer. Calls boundaryflux(BCKind(B),...) to allow for generic implementations by boundary condition type. Note that this method uses a different argument order convention than interact!. This is intended to faciliate stratigraphy independent implementations of certain boundary conditions (e.g. a simple Dirichlet boundary could be applied in the same manner to both the upper and lower boundary).

boundaryvalue(bc::BoundaryProcess, state)

Computes the value of the boundary condition specified by bc for the given layer/process combinations. Note that this method uses a different argument order convention than interact!. This is intended to faciliate stratigraphy independent implementations of certain boundary conditions (e.g. a simple Dirichlet boundary could be applied in the same manner to both the upper and lower boundary).

caninteract(layer1::Layer, layer2::Layer, state1, state2)
caninteract(l1::Layer, ::Process, l2::Layer, ::Process, state1, state2)

Returns true if and only if the given layer/process types are able to interact based on the current state. Defaults to checking whether both layers are currently active. This behavior should be overridden by subtypes where necessary.

computediagnostic!(l::Layer, state)
computediagnostic!(l::Layer, p::Process, state)

Updates all diagnostic/non-flux state variables for the given Layer based on the current prognostic state.

computefluxes!(l::Layer, p::Process, state)

Calculates all internal fluxes for a given layer. Note that an instance of computefluxes! must be provided for all non-boundary (subsurface) processes/layers.

criterion!(out::AbstractArray, ev::GridContinuousEvent, ::Layer, ::Process, state)

Event criterion for on-grid (i.e. multi-valued) continuous events. The condition for each grid cell should be stored in out.

criterion(::Event, ::Layer, ::Process, state)

Event criterion/condition. Should return a Bool for discrete events. For continuous events, this should be a real-valued function where the event is fired at the zeros/roots.

diagnosticstep!(layer::Layer, state)

Optionally performs discontinuous/discrete-time updates to the layer state. Should return true if the prognostic state was modified and false otherwise. Defaults to returning false.

events(::Layer, ::Process)

Defines "events" for a given Process on the given Layer. Implementations should return a Tuple of Events.

initialcondition!(::Layer, state)
initialcondition!(::Layer, ::Process, state)
initialcondition!(::VarInitializer, ::Layer, state)

Defines the initial condition for a given Layer and possibly an initializer. initialcondition! should compute initial values into all relevant state variables in state.

initialcondition!(layer1::Layer, layer2::Layer, state1, state2)
initialcondition!(::Layer, ::Process, ::Layer, ::Process, state1, state2)

Defines the initial condition for two processes on adjacent layers. initialcondition! should write initial values into all relevant state variables in state.

initializer(varname::Symbol, args...) = initializer(Val{varname}(), args...)
initializer(::Val{varname}, x::Number) => FunctionInitializer w/ constant
initializer(::Val{varname}, f::Function) => FunctionInitializer
initializer(::Val{varname}, profile::Profile, interp=Linear(), extrap=Flat()) => InterpInitializer

Convenience constructor for VarInitializer that selects the appropriate initializer type based on the arguments.

initializers(::Layer)
initializers(::Layer, ::Process)

Optional method that can be used to provide default initializers for state variables that will be run before user provided ones.

interact!(::Layer, ::Process, ::Layer, ::Process, state1, state2)

Defines a boundary interaction between two processes on adjacent layers. For any interaction, the order of the arguments follows decreasing depth, i.e. the first layer/process is always on top of the second layer/process. This ordering matters and separate dispatches must be provided for interactions in reverse order.

interactmaybe!(layer1::Layer, layer2::Layer, state1, state2)
interactmaybe!(layer1::Layer, p1::Process, layer2::Layer, p2::Process, state1, state2)

Conditionally invokes interact! if and only if caninteract is true.

isactive(::Layer, state)

Returns a boolean whether or not this layer is currently active in the stratigraphy and should interact with other layers. Note that computediagnostic! and computefluxes! are always invoked regardless of the current state of isactive. The default implementation of isactive always returns true.

midpoint(::Layer, state)
midpoint(::Layer, state, i)
midpoint(::Layer, state, ::typeof(first))
midpoint(::Layer, state, ::typeof(last))

Get midpoint (m) of layer or grid cell i.

odefunction(setup::Tile, u0, p, tspan; kwargs...)

Constructs a SciML ODEFunction given the model setup, initial state u0, parameters p, and tspan. Can (and should) be overridden by users to provide customized ODEFunction configurations for specific problem setups, e.g:

tile = Tile(strat,grid)
function CryoGrid.Setup.odefunction(::DefaultJac, setup::typeof(tile), u0, p, tspan)
    ...
    # make sure to return an instance of ODEFunction
end
...
prob = CryoGridProblem(tile, tspan, p)

JacobianStyle can also be extended to create custom traits which can then be applied to compatible Tiles.

parameterize(x::T) where {T}
parameterize(x::Unitful.AbstractQuantity; props...)
parameterize(p::AbstractParam; ignored...)

Recursively wraps x or nested numeric quantities in x with Param to mark them as parameters. If x is already a Param type, x will be returned as-is. If x is a numeric type, x will be wrapped in Param with associated properties props. If x is a struct type, x will be recursively unpacked and parameterize called on each field.

processes(l::Layer)

Fetches the process(es) attached to this layer, if any. Returned value must be of type Process. If the layer has more than one process, they should be combined together with Coupled(procs...).

resetfluxes!(layer::Layer, state)
resetfluxes!(layer::Layer, ::Process, state)

Resets all flux variables for the given layer/process to zero.

thickness(::Layer, state)
thickness(::Layer, state, i)
thickness(l::Layer, state, ::typeof(first))
thickness(l::Layer, state, ::typeof(last))

Get thickness (m) of layer or grid cell i.

timestep(::Layer, ::Process, state)

Retrieves the recommended timestep for the given Process defined on the given Layer. The default implementation returns Inf which indicates no timestep restriction. The actual chosen timestep will depend on the integrator being used and other user configuration options.

timestep(l::Layer, ps::CoupledProcesses{P}, state) where {P}

Default implementation of timestep for coupled process types. Calls each process in sequence.

trigger!(::Event, ::Layer, ::Process, state)
trigger!(ev::ContinuousEvent, ::ContinuousTrigger, ::Layer, ::Process, state)
trigger!(ev::GridContinuousEvent, ::ContinuousTrigger, ::Layer, ::Process, state)

Event action executed when criterion is met.

variables(layer::Layer, process::Process)
variables(::Layer)
variables(::Any)

Defines variables for a given Layer, Process, or arbitrary user-defined type. Implementations should return a Tuple of Vars.

volumetricfractions(::SubSurface, state)
volumetricfractions(::SubSurface, state, i)

Get the volumetric fractions of each constituent in the volume (at grid cell i, if specificed). All implementations of volumetricfractions are expected to obey a semi-consistent order in the returned Tuple of fractions; the first three consituents should always be θw,θi,θa, i.e. water, ice, and air, followed by any number of additional constituents which may be defined by the specific layer. There is no feasible way to verify that client code actually obeys this ordering, so be sure to double check your implementation, otherwise this can cause very subtle bugs!

Driver module SciML diffeq solvers.

CryoGridParameterEnsemble(
    prob::CryoGridProblem,
    Θ::AbstractMatrix;
    transform=identity,
    output_func=(sol, i) -> CryoGridOutput(sol),
    reduction=(u,data,i) -> (append!(u,data),false),
    ensprob_kwargs...
)

Constructs an EnsembleProblem for a m x N parameter matrix Θ, where N is the size of the ensemble and m is the dimensionality of the ensmeble parameter space (e.g. the number of parameters). transform should be a transform function that accepts an m-dimensional vector and produces a parameter vector (or CryoGridParams instance) which matches the output of CryoGrid.parameters. By default, param_map is the identity function; however, it may be customized to permit the construction of reduced-rank or reparameterized ensembles for which the parameter space differs from the full CryoGrid model parameter space.

Keyword arguments:

output_func: a function (sol,i)::Any which processes the ODESolution for ensemble member i and returns the result. It is recommended to save output to disk for non-trivial time spans to avoid slowdowns from serialization time when running the ensemble using parallel workers.

reduction: a function (u,data,i) which accumulates the result of output_func in u. Defaults to just appending data to u.

output_dir: Only used to specify the output directory for the default implementation of output_func. If a custom output_func is provided, this arugment is ignored.

All additional keyword arguments will be passed to EnsembleProblem.

See also SciMLBase.EnsembleProblem

ArrayCache{T,TA} <: StateVarCache

Simple state variable cache that writes directly into the given array.

AutoGrid <: DiscretizationStrategy

DiffCache{T,TCache}

Wrapper around PreallocationTools.DiffCache that stores state variables in forward-diff compatible cache arrays.

Base type for discretization "strategies" that generate a spatial discretization (typically a grid).

struct Grid{S,G,Q,A} <: AbstractDiscretization{Q,1}

Represents the 1D spatial discretization on which time integration is performed. S is a GridOffset, either Edges or Cells (always edges upon initial construction). The grid representation can be converted (allocation free) between grid edges and cells via the cells and edges methods. G represents the geometry/volume on which the vertical 1D discretization is applied. A is the underlying array type, and Q is the numerical type (e.g. Float64 or a Unitful.Quantity).

Simple discretization strategy that just supplies a pre-specified grid.

Profile{N,idxTypes,valTypes}

Represents a "profile" of values indexed typically over depth, i.e:

x₁: value 1 x₂: value 2 ...

where xᵢ are the indices.

StateVars{names,griddvars,TU,TV,TI,DF,DG}

Generic container for holding discretized state arrays from symbolic variables (Var types). The uproto field represetns a "prototype" of the prognostic state array which should fully define the state of the system at any given point.

StateVars(vars::NamedTuple, D::Numerics.AbstractDiscretization, cachector::Function, arrayproto::Type{A}=Vector) where {A<:AbstractVector}

Constructs a StateVars container from vars and the given discretization D. The function cachector must accept argments f(name::Symbol, x::AbstractVector) where name is the name of the state variable and x is the array generated by the discretization.

build_mass_matrix(u::ComponentVector, states::StateVars)

Constructs a mass matrix M⋅∂u∂t = f(u) suitable for the prognostic state vector u based on the defined variable types.

divergence!(dx::AbstractVector, j::AbstractVector, Δj::AbstractVector)

Calculates the first-order divergence over a 1D flux vector field j and grid cell lengths Δj. Divergences are added to existing values in dx.

dual(x::Number, ::Type{tag}) where {tag}
dual(x::A, ::Type{tag}) where {N,T,A<:SVector{N,T},tag}

Constructs a ForwardDiff.Dual number (or static array thereof) with tag from x.

flux!(j::AbstractVector, x::AbstractVector, Δx::AbstractVector, k::AbstractVector)

Calculates the first-order, non-linear spatial flux over a discretized variable x with conductivity k. x is assumed to have shape (N,), Δx shape (N-1,), and j and k shape (N+1,) such that j[2:end-1] represents the fluxes over the inner grid cell faces. Fluxes are added to existing values in j.

fpzero(x::Real)

If abs(x) is below the machine epsilon value eps(x), return zero. Otherwise return x.

harmonicmean!(h::AbstractVector, x::AbstractVector, w::AbstractVector)

Vectorized harmonic mean of elements in x with weights w. Output is stored in h, which should have size length(x)-1.

harmonicmean(x₁, x₂, w₁, w₂)

Simple weighted harmonic mean of two values, x₁ and x₂.

heaviside(x)

Differentiable implementation of heaviside step function, i.e:

$h(x) = \begin{cases} 1 & x ≥ 0 \\ 0 & x < 0 \end{cases}$

instantiate(::Var, ::T, ::Type{A}) where {T,N,D<:AbstractDiscretization{T,N},A<:AbstractArray{T,N}}

Produces a discretization of the given variable based on T and array type A.

linearmean(x, w)

Calculates a weighted average of x with weights w.

logistic(x)

Numerically stable logistic function.

$σ(x) = \begin{cases} \frac{1}{1+\exp(-x)} & x ≥ 0 \\ \frac{\exp(x)}{1+\exp(x)} & x < 0 \end{cases}$

logit(x)

Numerically stable logit function. True domain is (0,1) but inputs are clamped to (ϵ,1-ϵ) for numerical convenience, making the effective domain (-∞,∞).

makegrid(strategy::DiscretizationStrategy, bounds::NTuple{2,<:DistQuantity})
makegrid(layer::Layer, strategy::DiscretizationStrategy, bounds::NTuple{2,<:DistQuantity})

Constructs a Grid spanning bounds using the given strategy. makegrid can also be specified for specific Layer types when layers have specific discretization requirements.

nonlineardiffusion!(dx::AbstractVector, j::AbstractVector, x::AbstractVector, Δx::AbstractVector, k::AbstractVector, Δk::AbstractVector)

Fast alternative to flux! and divergence! which computes fluxes and divergences (via _flux_kernel and _div_kernel) in a single pass. Note, however, that loop vectorization with @turbo is not possible because of necessary loop-carried dependencies. Fluxes and divergences are added to the existing values stored in j and dx.

softplus(x)

Numerically stable softplus function.

$s(x) = \log(1+\exp(-|x|)) + \max(x,ϵ)$

softplusinv(x)

Numerically stable softplus inverse function. True domain is (0,∞) but inputs are clamped to (ϵ,∞) for numerical convenience, making the effective domain (-∞,∞).

tdma_solve!(x, a, b, c, d)

Tridiagonal matrix solver via the Thomas algorithm; adapted from original CryoGridLite implementation. Modifies all input vectors in-place.

updategrid!(grid::Grid{Edges}, edges)
updategrid!(grid::Grid{Edges}, z0, thick::AbstractVector)

Updates all grid values based on new grid edges or an initial bottom z0 + cell thick.

∇(f::F, x::AbstractArray) where {F}

Takes a function y = f(x) and vector-valued argument x and returns a tuple: (y, ∇ₓy). The gradient is calculated using forward-mode automatic differentiation.

∇(f::F, x::Number) where {F}

Takes a function y = f(x) and argument x and returns a tuple: (y, ∂y∂x). The derivative is calculated using forward-mode automatic differentiation.

Pre-built CryoGrid configurations for rapid prototyping.

High resolution grid spacing used by the GIPLv2 model: https://github.com/Elchin/GIPL

SoilHeatTile([heatop=:H], upperbc::BoundaryProcess, soilprofile::Profile, init::CryoGrid.Initializer...; grid::Grid=DefaultGrid_10cm, tile_kwargs...) where {F<:FreezeCurve}

Builds a simple one-layer soil/heat-conduction model with the given grid and configuration.

AbstractGround{Tpara<:GroundParameterization,Theat<:HeatBalance,Twater<:WaterBalance} <: SubSurface

Base type for all ground layers defining heat and water balance schemes.

Ground{Tpara,Theat<:HeatBalance,Twater<:WaterBalance,Taux} <: Soil{Tpara,Theat,Twater}

Generic representation of a Ground layer with material parameterization para.

Ground(soilprofile::SoilProfile; kwargs...)

Heterogeneous{V,N,D,Taux} <: SoilParameterization

Special SoilParameterization which wraps a Profile of another soil parameterization type to indicate that it should be heterogeneous with over depth.

Pressure <: RREqForm

Pressure head based form of Richards equation:

∂θ∂ψ⋅∇ₜΨ = ∇ₓ[k(θ)[1 + ∇ₓΨ]]

RREqForm

Base type for different formulations of Richardson-Richard's equation.

RichardsEq{Tform<:RichardsEqFormulation,Tswrc<:SWRC,Taux,TΩ} <: Hydrology.WaterFlow

The Richardson-Richards equation describes the flow of water in porous media under unsaturated condition.

Saturation <: RREqForm

Sautration based form of Richards equation:

∇ₜθ = ∇ₓ[k(θ)[1 + ∇ₓΨ(θ)]]

SimpleSoil{Tfc,Tpor,Tsat,Torg,Thp,Twp} <: SoilParameterization

Represents a simple organic/mineral soil mixutre in terms of its characteristic fractions: i.e. natural porosity, saturation, and organic solid fraction. This is the standard CryoGrid representation of a discrete soil volume.

Soil{Tpara,Theat,Twater}

Type alias for any AbstractGround layer with parameterization of type SoilParameterization.

SoilParameterization

Base type for parameterizations of soil consituents.

SoilProfile{N,IT,VT} = Profile{N,IT,VT} where {N,IT<:NTuple{N},VT<:NTuple{N,SoilParameterization}}

Alias for depthwise Profile where the values are SoilParameterization types.

SoilProfile(pairs::Pair{<:DistQuantity,<:SoilParameterization}...)

Alias for Profile(pairs...) specific for SoilProfiles.

SoilTexture{Tsand,Tclay,Tsilt}

Represents soil "texture" as a simple mixture of sand, silt, and clay.

freezethaw!(sfcc::SFCC, soil::Soil, heat::HeatBalance{<:TemperatureBased}, state)
freezethaw!(sfcc::SFCC, soil::Soil, heat::HeatBalance{<:EnthalpyBased}, state)

Updates state variables according to the specified SFCC function and solver. For heat conduction with enthalpy, evaluation of the inverse enthalpy function is performed using the given solver. For heat conduction with temperature, we can simply evaluate the freeze curve to get dHdT, θw, and H.

Gets the ThermalProperties for the given soil layer.

Gets the ThermalProperties for the given soil layer.

Gets the HydraulicProperties for the given soil layer.

Gets the HydraulicProperties for the given soil layer.

SoilHydraulicProperties{::Type{SoilParameterization}; kw_sat, fieldcapacity}

Material hydraulic properties for the given SoilParameterization type.

impedencefactor(water::WaterBalance{<:RichardsEq}, θw, θwi)

Impedence factor which represents the blockage of water-filled pores by ice (see Hansson et al. 2004 and Westermann et al. 2022).

initialize_sfccsolver!([::FreezeCurve], ::Soil, state)

Retrieves and initializes the SFCC solver for the given layer.

mineral(soil::Soil, state, i)
mineral(soil::Soil, state)
mineral(soil::Soil)

Retrieves the volumetric mineral content for the given layer at grid cell i, if provided.

organic(soil::Soil, state, i)
organic(soil::Soil, state)
organic(soil::Soil)

Retrieves the volumetric organic content for the given layer at grid cell i, if provided.

porosity(soil::Soil, state, i)
porosity(soil::Soil, state)
porosity(soil::Soil)

Retrieves the porosity for the given layer at grid cell i, if provided.

saturation(soil, state, i)
saturation(soil::Soil, state)
saturation(soil::Soil)

Retrieves the saturation level for the given layer at grid cell i, if provided.

sfccheatcap(soil::Soil, state, i)

Constructs a function heatcap(θw,θwi,θsat) that comptues the heat capacity given the current state for grid cell i.

sfcckwargs(f::SFCC, soil::Soil, state, i)

Builds a named tuple of values corresponding to each keyword arguments of the SFCC f which should be set according to the layer/process properties or state. The default implementation sets only the saturated water content, θsat = porosity.

(tile::AbstractTile{true})(du,u,p,t,dt=1.0)
(tile::AbstractTile{false})(u,p,t,dt=1.0)

Invokes computefluxes! on tile to compute the time derivative du/dt.

AbstractTile{iip}

Base type for 1D tiles. iip is a boolean value that indicates, if true, whether the model operates on state variables in-place (overwriting arrays) or if false, out-of-place (copying arrays). Current only in-place is supported.

JacobianStyle

Trait for indicating Jacobian sparsity of a CryoGrid ODEProblem.

JacobianStyle(::AbstractTile)

Can be overriden/extended to specify Jacobian structure for specific Tiles.

LayerState{TStates<:NamedTuple,Tt,Tdt}

State for a single layer, typically constructed from a parent TileState.

Stratigraphy{N,TLayers<:NamedTuple,TBoundaries}

Defines a 1-dimensional stratigraphy by connecting a top and bottom layer to one or more subsurface layers.

Tile{TStrat,TGrid,TStates,TInits,TEvents,iip} <: AbstractTile{iip}

Defines the full specification of a single CryoGrid tile; i.e. stratigraphy, grid, and state variables.

Tile(integrator::SciMLBase.DEIntegrator)

Constructs a Tile from a SciMLBase DEIntegrator.

Tile(
    @nospecialize(strat::Stratigraphy),
    @nospecialize(discretization_strategy::DiscretizationStrategy),
    @nospecialize(inits::CryoGrid.Initializer...);
    metadata::Dict=Dict(),
    arraytype::Type{A}=Vector,
    iip::Bool=true,
    chunk_size=nothing,
)

Constructs a Tile from the given stratigraphy and discretization strategy. arraytype keyword arg should be an array instance (of any arbitrary length, including zero, contents are ignored) that will determine the array type used for all state vectors.

TileState{TGrid,TStates,Tu,Tt,Tdt}

Represents the state of a CryoGrid Tile.

getvar(var::Symbol, integrator::SciMLBase.DEIntegrator)

getvar(name::Symbol, tile::Tile, u; interp=true)
getvar(::Val{name}, tile::Tile, u; interp=true)

Retrieves the (diagnostic or prognostic) grid variable from tile given prognostic state u. If name is not a variable in the tile, or if it is not a grid variable, nothing is returned.

Rebuilds the obj adding name to the layer field to all Params, if any are defined.

Initialize StateVars which holds the caches for all defined state variables.

computefluxes!(_tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,true}, _du, _u, p, t) where {TStrat,TGrid,TStates,TInits,TEvents}

Time derivative step function (i.e. du/dt) for any arbitrary Tile. Specialized code is generated and compiled on the fly via the @generated macro to ensure type stability. The generated code updates each layer in the stratigraphy in sequence, i.e for each layer 1 <= i <= N:

computediagnostic!(layer[i], ...)
interact!(layer[i], ..., layer[i+1], ...)
computefluxes!(layer[i], ...)

computefluxes!(::T,du,u,p,t) where {T<:AbstractTile}

In-place update function for tile T. Computes du/dt and stores the result in du.

computefluxes(::T,u,p,t) where {T<:AbstractTile}

Out-of-place update function for tile T. Computes and returns du/dt as vector with same size as u.

domain(tile::Tile)

Returns a function isoutofdomain(u,p,t) which checks whether any prognostic variable values in u are outside of their respective domains. Only variables which have at least one finite domain endpoint are checked; variables with unbounded domains are ignored.

getstate(tile::Tile, u, du, t, dt=1.0)

Constructs a LayerState representing the full state of layername given tile, state vectors u and du, and the time step t.

getstate(integrator::SciMLBase.DEIntegrator)

Builds the TileState given an initialized integrator.

isinplace(tile::AbstractTile{iip}) where {iip}

Returns true if tile uses in-place mode, false if out-of-place.

parameters(tile::Tile; kwargs...)

Extracts all parameters from tile.

resolve(tile::Tile, p, t)

Resolves/updates the given tile by: (1) Replacing all ModelParameters.AbstractParam values in tile with their (possibly updated) value from p. (2) Resolving the boundary depths of the Stratigraphy layers by invoking resolveboundaries. (3) Replacing all instances of DynamicParameterization with their resolved values given the current state.

Returns the reconstructed Tile instance.

withaxes(u::AbstractArray, ::Tile)

Constructs a ComponentArray with labeled axes from the given state vector u. Assumes u to be of the same type/shape as setup.uproto.

initialcondition!(tile::Tile, tspan::NTuple{2,Float64}[, p])
initialcondition!(tile::Tile, tspan::NTuple{2,DateTime}[, p])

Calls initialcondition! on all layers/processes and returns the fully constructed u0 and du0 states.

interact!(strat::Stratigraphy, state)

Special implementation of interact! that iterates over each pair of layers in the stratigraphy which are adjacent and "active" based on the current state. state must have properties defined corresponding to the name of each layer such that getproperty would return the appropriate state object for the i'th layer in the stratigraphy.

parameterize(tile::Tile)

Adds parameter information to all nested types in tile by recursively calling parameterize.

timestep(_tile::Tile, _du, _u, p, t)

Computes the maximum permissible forward timestep for this Tile given the current u, p, and t.

variables(tile::Tile)

Returns a tuple of all variables defined in the tile.

Convenience macro for defining stratigraphies with multiple subsurface layers.

BucketScheme <: WaterFlow

"Bucket" water scheme for downward advective flow due to gravity.

DampedET{Tftr,Tfev,Tdtr,Tdev}

Corresponds to evapotranspiration scheme 2 described in section 2.2.4 of Westermann et al. (2022).

EvapTop <: Evapotranspiration

Represents a simple evaporation-only scheme where water is drawn only from the top-most grid cell. Corresponds to evapotranspiration scheme 3 described in section 2.2.4 of Westermann et al. (2022).

Evapotranspiration

Base type for parameterizations of evapotranspiration (ET).

HydraulicProperties

Default material hydraulic properties.

NoFlow <: WaterFlow

Represents a zero flow scheme

WaterBalance{TFlow<:WaterFlow,TET<:Union{Nothing,Evapotranspiration},Tdt,Taux,TProp} <: CryoGrid.SubSurfaceProcess

Represents subsurface water transport processes.

WaterBalanceProperties

Numerical constants shared across water balance implementations.

WaterFlow

Base type for different formulations of water flow in WaterBalance.

WaterTableInitializer <: VarInitializer{:sat}

Simple, piecewise constant initializer for saturation state that takes a surface-level saturation sat0 and water table depth z_tab and produces a two-segment, piecewise constant profile with the saturation level set to (sat0 + 1.0) / 2 from the halfway point down to the water table.

ETflux(::SubSurface, water::WaterBalance{<:WaterFlow,<:Evapotranspiration}, state)

Computes the ET base flux as Qe / (Lsg*ρw) where state.Qe is typically provided as a boundary condition.

advectiveflux(θw_up, θwi_lo, θsat_lo, θmin, kw)

Computes the advective downward water flux given liquid water content in the upper grid cell (θw_up), the minimum water content (θmin), and hydraulic conductivity at the boundary (kw).

balanceflux([water_up::WaterBalance], water_lo::WaterBalance, jw, θw_up, θw_lo, θwi_up, θwi_lo, θsat_up, θsat_lo, sat_up, sat_lo, Δz_up, Δz_lo)

Rescales the water flux jw to satisfy mass conservation based on the given hydrological states for two adjacent grid cells ('up' refers to upper and 'lo' refers to lower). Also applies nonlinear reduction factors for near-saturated conditions.

balancefluxes!(sub::SubSurface, water::WaterBalance, state)

Sums vertical and evapotranspirative flux components jw_v and jw_ET and applies balanceflux to all grid cell faces.

evapotranspiration!(::SubSurface, ::WaterBalance, state)

Computes diagnostic evapotranspiration quantities for the given layer and water balance configuration, storing the results in state. This method should generally be called before interact! for WaterBalance, e.g. in computediagnostic!.

evapotranspirative_fluxes!(::SubSurface, ::WaterBalance, state)

Computes diagnostic evapotranspiration quantities for the given layer and water balance configuration, storing the results in state. This method should generally be called in or after the surface interaction for WaterBalance.

hydraulicconductivity!(sub::SubSurface, water::WaterBalance, state)

Computes hydraulic conductivities for the given subsurface layer and water balance scheme.

hydraulicconductivity(water::WaterBalance, kw_sat, θw, θwi, θsat)

Computes the hydraulic conductivity for the given water balance configuration, current unfrozen water content θw, total water/ice content θwi, and saturated (maximum) water content θsat.

hydraulicproperties(::SubSurface)

Retrieves the hydraulic properties from the given subsurface layer.

interact_ET!(::Top, ::WaterBC, ::SubSurface, ::WaterBalance, state1, state2)
interact_ET!(::SubSurface, ::WaterBalance, ::Bottom, ::WaterBC, state1, state2)
interact_ET!(::SubSurface, ::WaterBalance, ::SubSurface, ::WaterBalance, state1, state2)

Specialized layer interaction for evapotranspirative processes. Default implementation does nothing. Can be overridden by ET schemes and invoked in the relevant interact! implementation.

kwsat(::SubSurface, state, i)

Hydraulic conductivity at saturation.

limit_lower_flux(water::WaterBalance, jw, θw, θwi, θsat, sat, Δz)

Flux limiter for fluxes coming from "below"; limits jw based on the available water in the current cell as well as free pore space.

limit_upper_flux(water::WaterBalance, jw, θw, θwi, θsat, sat, Δz)

Flux limiter for fluxes coming from "above"; limits jw based on the available water in the current cell as well as free pore space.

maxinfiltration(::SubSurface, ::WaterBalance, state)

Returns the current maximum infiltration rate for the given SubSurface layer and water balance scheme. Defaults to the hydraulic conductivity at the upper boundary.

maxwater(sub::SubSurface, ::WaterBalance) 
maxwater(sub::SubSurface, water::WaterBalance, state)
maxwater(::SubSurface, ::WaterBalance, state, i)

Returns the maximum volumetric water content (saturation point) for grid cell i.

minwater(::SubSurface, water::WaterBalance)
minwater(::SubSurface, water::WaterBalance, state, i)

Returns the minimum volumetric water content (typically field capacity for simplified schemes) for grid cell i. Defaults to zero.

wateradvection!(sub::SubSurface, water::WaterBalance, state)

Computes the advective component of water fluxes due to gravity and stores the result in state.jw.

watercontent!(::SubSurface, ::WaterBalance, state)

Computes the volumetric water content from current saturation or pressure state.

watercontent(::SubSurface, state)
watercontent(::SubSurface, state, i)

Returns the total water content θwi from the given subsurface layer and/or current state.

waterdensity(sub::SubSurface)

Retrieves the density of water ρw from the given SubSurface layer. Default implementation assumes that WaterBalance is provided as a field water on sub; this can of course, however, be overridden.

waterdiffusion!(::SubSurface, ::WaterBalance, state)

Computes diffusive fluxes for water balance, if defined. Default implementation does nothing.

waterprognostic!(::SubSurface, ::WaterBalance, state)

Computes the prognostic time derivative for the water balance, usually based on dθwi. Implementation depends on which water flow scheme is being used.

resetfluxes!(::SubSurface, water::WaterBalance, state)

Resets flux terms (jw and dθwi) for WaterBalance.

Constant <: SourceTerm

Parametric source term that is constant through time and space.

Periodic <: SourceTerm

Parametric source term that is periodic through time and constant through space.

Source{P,T,S} <: SubSurfaceProcess

Generic "source" process to provide additive fluxes for one or more prognostic variables in process P. The behavior is governed by the SourceTerm, T, and further user specialization can be implemented via aux::S.

Abstract base type for source terms which define a Source process.

Diffusion1D{progvar} <: HeatOperator{progvar}

Represents a standard method-of-lines (MOL) forward diffusion operator in 1 dimension.

Type alias for HeatOperator{:H}, i.e. enthalpy-based heat transfer operators.

EnthalpyImplicit <: HeatOperator{:H}

Implicit enthalpy formulation of Swaminathan and Voller (1992) and Langer et al. (2022). Note that this heat operator formulation does not compute a divergence dH but only computes the necessary diffusion coefficients for use by an appropriate solver. See the LiteImplicit module for the appropriate solver algorithms.

GeothermalHeatFlux{TQ} <: BoundaryProcess{HeatBalance}

Represents a simple, forced Neumann heat flux boundary condition for HeatBalance processes.

GroundHeatFlux{TE,TQ} <: BoundaryProcess{HeatBalance}

Represents a simple, forced Neumann heat flux boundary condition for HeatBalance processes.

HeatBalance{THeatOp<:HeatOperator,Tdt,Tprop} <: SubSurfaceProcess

Represents subsurface heat transfer processes. The formulation of heat transfer is governed by the HeatOperator, op.

Type alias for the implicit enthalpy formulation of HeatBalance.

Numerical constants for pararameterizing heat processes.

HeatOperator{progvar}

Base type for different numerical formulations of heat transfer.

LinearTwoPhaseTempProfile{TT1,TT2,TTb,TTm,Tz1,Tz2,Tz3,Tz4} <: VarInitializer{:T}

Simple, piecewise linear temprature initializer that uses three temperature values and four characteristic depths to initialize the temperature profile. T0 is the initial surface temperature, T1 is the permafrost temperature at z_deep, z_thaw and z_base are the top and bottom freezing fronts which are both assumed to be have temperature equal to Tm.

(init::LinearTwoPhaseTempProfile)(T::AbstractVector, grid::Grid)

Evaluate the initializer on the given temperature vector T which should match the length of cells(grid).

StefanProblem{Tp<:StefanParameters,Tx,Tt}

Represents the simple two-phase Stefan problem defined on a semi-infinite slab. The one-phase Stefan problem can be computed by setting the parameters T_s = T_m.

(sol::StefanSolution)(x,t)

Temperature solution, T(x,t).

(sol::StefanSolution)(t)

Moving boundary solution, x_m(t).

TemperatureBC{E,F} <: BoundaryProcess{HeatBalance}

Represents a simple, forced Dirichlet temperature boundary condition for HeatBalance processes.

Type alias for HeatOperator{:T}, i.e. temperature-based heat transfer operators.

ThermalProperties

Material thermal properties, e.g. conductivity and heat capacity. By default, this includes the thermal properties of water, ice, and air. This can be extended by passing additional properties into the constructor.

TemperatureProfile(pairs::Pair{<:Union{DistQuantity,Param},<:Union{TempQuantity,Param}}...)

Convenience constructor for Numerics.Profile which automatically converts temperature quantities.

advectiveflux(jw, T₁, T₂, cw, L)

Computes the advective energy flux between grid cells with temperatures T₁ and T₂ given the heat capacity of water cw and latent heat of fusion L.

dHdT(T, C, L, ∂θw∂T, ch_w, ch_i) = C + ∂θw∂T*(L + T*(ch_w - ch_i))

Computes the apparent or "effective" heat capacity ∂H∂T as a function of temperature, volumetric heat capacity, latent heat of fusion, derivative of the freeze curve ∂θw∂T, and the constituent heat capacities of water and ice.

enthalpy(T, C, L, θ) = T*C + L*θ

Discrete enthalpy function on temperature, heat capacity, specific latent heat of fusion, and liquid water content.

enthalpyinv(H, C, L, θ) = (H - L*θ) / C

Discrete inverse enthalpy function given H, C, L, and θ.

enthalpyinv([::FreezeCurve], sub::SubSurface, heat::HeatBalance, state)

Evaluates the inverse enthalpy function (H -> T) on the current state.

freezecurve(sub::SubSurface)

Returns the soil freezing characteristic FreezeCurve for the given subsurface layer. Defautls to FreeWater.

Implementation of "free water" freezing characteristic for any subsurface layer. Assumes that state contains at least temperature (T), enthalpy (H), heat capacity (C), total water content (θwi), and liquid water content (θw).

freezethaw!(::FreezeCurve, ::SubSurface, ::Process, state)

Calculates freezing and thawing effects, including evaluation of the freeze curve. In general, this function should compute at least the liquid/frozen water contents and the corresponding heat capacity. Other variables such as temperature or enthalpy may also need to be computed depending on the thermal scheme being implemented.

geometric_conductivity(ks::NTuple{N}, θs::NTuple{N}) where {N}

Geometric mean of constituent thermal conductivities according to Woodside and Messmer (1961).

Woodside, W. & Messmer, J.H. 1961. Thermal conductivity of porous media. I. Unconsolidated sands. Journal of Applied Physics, 32, 1688–1699.

Analytical solution to linear form of the heat equation on a semi-infinite rod with periodic fluctuations at the upper boundary. Solution taken from Riseborough et al. 2008.

heatcapacities(::SubSurface)

Get heat capacities for generic SubSurface layer.

heatcapacity!(sub::SubSurface, state)

Computes the heat capacity for the given layer from the current state and stores the result in-place in the state variable C.

heatcapacity(sub::SubSurface, state, i)

Computes the heat capacity as a weighted average over constituent capacities with volumetric fractions θfracs.

Common variable definitions for all heat implementations.

quadratic_parallel_conductivity(ks, θs)

The "quadratic parallel" thermal conductivity formula as defined by Cosenza et al. 2003:

\[k = [\sum_{i=1}^N θᵢ\sqrt{kᵢ}]^2\]

Cosenza, P., Guérin, R., and Tabbagh, A.: Relationship between thermal conductivity and water content of soils using numerical modelling, European Journal of Soil Science, 54, 581–588, https://doi.org/10.1046/j.1365-2389.2003.00539.x, 2003.

stefan_boundary(x0, t0, λ, α, t)

Evaluates the position of the Stefan moving boundary at time t given the Stefan constant λ, diffusivity α, and initial time and position t0, x0.

stefan_number(c, ρ, Lf, ΔT)

Calculates the Stefan number for heat capacity c, desnity ρ, specific latent heat of fusion Lf, and temperature delta ΔT.

stefan_residual(λ, T_m, T_l, T_s, k_l, c_l, k_s, c_s, ρ, Lf)

Evaluates the residual function for the transcendental equation in the two-phase Stefan problem.

stefan_temperature_liquid(T_l, T_m, λ, α_l, x0, t0, x, t)

Evaluates the generic Stefan analytical solution for temperature in the liquid region.

stefan_temperature_solid(T_s, T_m, λ, α_s, α_l, x0, t0, x, t)

Evaluates the generic Stefan analytical solution for temperature in the solid region.

thermalconductivities(sub::SubSurface)

Get thermal conductivities for generic SubSurface layer.

thermalconductivity!(sub::SubSurface, state)

Computes the thermal conductivity for the given layer from the current state and stores the result in-place in the state variable k.

thermalconductivity(sub::SubSurface, heat::HeatBalance, state, i)

Computes the thermal conductivity as a squared weighted sum over constituent conductivities with volumetric fractions θfracs.

thermalproperties(::SubSurface)

Returns the thermal properties for the given subsurface layer.

water_energy_advection!(::SubSurface, ::Coupled(WaterBalance, HeatBalance), state)

Adds advective energy fluxes for all internal grid cell faces.

Diagnostic state update for heat conduction (all state configurations) on any subsurface layer.

timestep(::SubSurface, ::HeatBalance{THeatOp,CFL}, state) where {THeatOp}

Implementation of timestep for HeatBalance using the Courant-Fredrichs-Lewy condition defined as: Δt_max = u*Δx^2, whereu is the "characteristic velocity" which here is taken to be the diffusivity: ∂H∂T / kc.

Variable definitions for heat conduction (enthalpy) on any SubSurface layer.

Variable definitions for heat conduction (temperature).

Bulk{Tden,Tthresh,Theat,Twater} <: SnowpackParameterization

Simple, bulk ("single layer") snow scheme where snowpack is represented as a single grid cell with homogenous state.

BulkSnowpack = Snowpack{<:Bulk}

Type alias for Snowpack with Bulk parameterization.

LiteGridded{T1,T2,T3,T4} <: SnowpackParameterization

Simple bulk density snow scheme minimally adapted from CryoGridLite.

SnowAblationScheme

Base type for different snow ablation (i.e. melting or redistribution) schemes.

SnowAccumulationScheme

Base type for different snow accumulation schemes.

SnowBC

Type alias for any BoundaryProcess compatible with SnowMassBalance.

SnowDensityScheme

Base type for different snow density schemes.

SnowMassBalance{TAcc,TAbl} <: CryoGrid.SubSurfaceProcess

Subsurface process for snow layers governing how snow is accumulated and ablated.

SnowThermalConductivity

Base type for snow thermal conductivity parameterizations.

SnowThermalProperties{Tcond<:SnowThermalConductivity,Tprop}

Specifies the thermal properties of the snowpack.

Snowpack{Tpara<:SnowpackParameterization,Tmass<:SnowMassBalance,Twater<:WaterBalance,Theat<:HeatBalance,Taux} <: CryoGrid.SubSurface

Generic representation of a snowpack "subsurface" layer.

Snowpack(para::SnowpackParameterization; kwargs...)

Convenience constructor that accepts the parameterization as a positional argument.

SnowpackParameterization

Base type for snowpack paramterization schemes.

ablation!(::Top, ::SnowBC, ::Snowpack, ::SnowMassBalance, stop, ssnow)

Computes snow mass balance fluxes due to ablation (e.g. snow melt).

accumulation!(::Top, ::SnowBC, ::Snowpack, ::SnowMassBalance, stop, ssnow)

Computes snow mass balance fluxes due to accumulation (e.g. snowfall).

calculate_degree_day_snow_melt(ddm::DegreeDayMelt, T_ub::Number)

Implementation of degree day melting scheme.

compaction!(::Top, ::SnowBC, ::Snowpack, ::SnowMassBalance, stop, ssnow)

Computes snow density changes due to compaction, if defined.

snowdensity!(::Snowpack, ::SnowMassBalance, state)

Computes the current snow density (if necessary) and stores the result in state.ρsn.

snowdensity(::Snowpack, state)

Retrieve the current snow density.

snowdepth(::Snowpack, ::SnowMassBalance, state)

Retrieve the current snow depth.

snowfall(::SnowBC, state)

Retrieves the current snowfall flux for the given SnowBC. Defaults to state.jw_snow.

snowvariables(::Snowpack)

Default variables common to all snow schemes.

snowwater(::Snowpack, ::SnowMassBalance, state)

Retrieve the current snow water equivalent of the snowpack.

threshold(snow::BulkSnowpack)

Retrieves the minimum snow threshold for the bulk snow scheme.

threshold(::Snowpack)

Retrieves the snow cover threshold for this Snowpack layer to become active.