Tutorial

Terminology

ConceptualClimateModels.jl follows a process-based modelling approach to make differential equation systems from processes. A process is simply a particular equation defining the dynamics of a climate variable, while also explicitly defining which variable the equation defines. A vector of processes is composed by the user, and given to the main function processes_to_coupledodes which bundles them into a system of equations that creates a dynamical system. The dynamical system can then be further analyzed in terms of stability properties, multistability, tipping, periodic (or not) behavior, chaos, and many more aspects, via the DynamicalSystems.jl library (see the examples).

Note the distinction: a process is not the climate variable (such as "clouds" or "insolation"); rather it is the exact equation that defines the behavior of the climate variable, and could itself utilize many other already existing climate variables, or introduce new ones. In terminology of climate science a process is a generalization of the term "parameterization". Many different processes may describe the behavior of a particular variable and typically one wants to analyze what the model does when using one versus the other process.

Introductory example

Let's say that we want to create the most basic energy balance model,

\[c_T \frac{dT}{dt} = ASR - OLR + f\]

where $ASR$ is the absorbed solar radiation given by

\[ASR = S (1-\alpha)\]

with $\alpha$ the planetary albedo, $OLR$ is the outgoing longwave radiation given by the linearized expression

\[OLR = A + BT\]

and $f$ some radiative forcing at the top of the atmosphere, that is based on CO2 concentrations and given by

\[f = 3.7\log_2\left(\frac{CO_2}{400}\right)\]

with CO2 concentrations in ppm.

To create this model with ConceptualClimateModels.jl while providing the least information possible we can do:

using ConceptualClimateModels
using ConceptualClimateModels.GlobalMeanEBM # submodule for global mean models

processes = [
    BasicRadiationBalance(),
    LinearOLR(),
    ParameterProcess(α),
    CO2Forcing(), # for default CO2 valu this is zero forcing
]

ds = processes_to_coupledodes(processes, GlobalMeanEBM)
println(dynamical_system_summary(ds))

┌ Warning: Variable CO2(t) was introduced in process of variable f(t).
│ However, a process for CO2(t) was not provided,
│ and there is no default process for it either.
│ Since it has a default value, we make it a parameter by adding a process:
│ `ParameterProcess(CO2)`.
└ @ ProcessBasedModelling ~/.julia/packages/ProcessBasedModelling/itwz9/src/make.jl:98
┌ Warning: Variable S(t) was introduced in process of variable ASR(t).
│ However, a process for S(t) was not provided,
│ and there is no default process for it either.
│ Since it has a default value, we make it a parameter by adding a process:
│ `ParameterProcess(S)`.
└ @ ProcessBasedModelling ~/.julia/packages/ProcessBasedModelling/itwz9/src/make.jl:98
1-dimensional CoupledODEs
 deterministic: true
 discrete time: false
 in-place:      false

with equations:
 Differential(t)(T(t)) ~ (-0.0029390154298310064(ASR(t) + f(t) - OLR(t))) / (-τ_T)
 OLR(t) ~ -277.0 + 1.8T(t)
 S(t) ~ S_0
 α(t) ~ α_0
 CO2(t) ~ CO2_0
 ASR(t) ~ 340.25S(t)*(1 - α(t))
 f(t) ~ 3.7log2((1//400)*CO2(t))

and parameter values:
 τ_T ~ 1.4695077149155033e6
 α_0 ~ 0.3
 CO2_0 ~ 400.0
 S_0 ~ 1.0

The output is a dynamical system from DynamicalSystems.jl that is generated via symbolic expressions based on ModelingToolkit.jl, utilizing the process-based approach of ProcessBasedModelling.jl.

As the dynamical system is made by symbolic expressions, these can be obtained back at any time:

using DynamicalSystems
# access the MTK model that stores the symbolic bindings
mtk = referrenced_sciml_model(ds)
# show the equations of the dynamic state variables of the dynamical system
equations(mtk)

\[ \begin{align} \frac{\mathrm{d} T\left( t \right)}{\mathrm{d}t} =& \frac{ - 0.002939 \left( \mathrm{ASR}\left( t \right) + f\left( t \right) - \mathrm{OLR}\left( t \right) \right)}{ - \tau_{T}} \end{align} \]

# show state functions that are observable,
# i.e., they do not have a time derivative, they are not dynamic state variables
observed(mtk)

\[ \begin{align} \mathrm{OLR}\left( t \right) =& -277 + 1.8 T\left( t \right) \\ S\left( t \right) =& S_{0} \\ \alpha\left( t \right) =& \alpha_{0} \\ \mathrm{CO2}\left( t \right) =& CO2_{0} \\ \mathrm{ASR}\left( t \right) =& 340.25 S\left( t \right) \left( 1 - \alpha\left( t \right) \right) \\ f\left( t \right) =& 3.7 \log_{2}\left( \frac{1}{400} \mathrm{CO2}\left( t \right) \right) \end{align} \]

# show parameters
parameters(mtk)

4-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 τ_T
 α_0
 CO2_0
 S_0

The symbolic variables and parameters can also be used to query or alter the dynamical system. In the equations above we say that the symbolic variable CO2 was equated to the parameter CO2_0. There are multiple ways to obtain a symbolic index provided we know its name. First, we can re-create the symbolic parameter:

index = only(@parameters CO2_0)

\[ \begin{equation} CO2_{0} \end{equation} \]

Second, we can use the retrieved MTK model and access its CO2_0 field, which will return the symbolic variable:

index = mtk.CO2_0

\[ \begin{equation} CO2_{0} \end{equation} \]

Third, we can use a Symbol corresponding to the variable name. This is typically the simplest way.

index = :CO2_0

:CO2_0

we can query the value of this named parameter in the system,

current_parameter(ds, index)

400.0

or alter it:

# access symbolic parameter CO2_0 from the tracked symbolic list of the model
set_parameter!(ds, index, 800)

Similarly, we can obtain or alter values corresponding to the dynamic variables, or observed functions of the state of the system, using their symbolic indices. For example we can obtain the value corresponding to symbolic variable $T$ at the current state of the dynamical system with:

observe_state(ds, T) # binding `T` already exists in scope

290.0

We didn't have to create T because when we did using ConceptualClimateModels.GlobalMeanEBM, T was brought into scope. Similarly, we can obtain the $OLR$ (outgoing longwave radiation), using the Symbol representation of the symbolic variable:

observe_state(ds, :OLR)

245.0

The advantage of using the symbolic variable OLR itself instead of its Symbol representation, is that symbolic variables can do arbitrary symbolic computations. For example let's say that when we created the equations of the system we didn't have defined a variable representing the expression $T^2/OLR$. That's not a problem, we can query the expression nevertheless:

observe_state(ds, T^2/OLR) # symbolic expression

343.265306122449

Let's unpack the steps that led to this level of automation.

Processes

A process is conceptually an equation the defines a climate variable or observable. All processes that define the system are then composed into a set of differential equations via processes_to_coupledodes (or processes_to_mtkmodel) that represent the climate model.

For example, the process

T_process = BasicRadiationBalance()

BasicRadiationBalance(5.0e8, f(t), T(t), ASR(t), OLR(t))

is the process defining the variable $T$, representing temperature. We can learn this by either reading the documentation string of BasicRadiationBalance, or querying it directly:

using ProcessBasedModelling: lhs, rhs
# This is the equation created by the process
lhs(T_process) ~ rhs(T_process)

\[ \begin{equation} \frac{\mathrm{d} T\left( t \right)}{\mathrm{d}t} \tau_{T} = 0.002939 \left( \mathrm{ASR}\left( t \right) + f\left( t \right) - \mathrm{OLR}\left( t \right) \right) \end{equation} \]

Notice that this process does not further define e.g. outgoing longwave radiation OLR(t). That is why in the original example we also provided LinearOLR, which defines it:

OLR_process = LinearOLR()
lhs(OLR_process) ~ rhs(OLR_process)

\[ \begin{equation} \mathrm{OLR}\left( t \right) = -277 + 1.8 T\left( t \right) \end{equation} \]

Each physical "observable" or variable that can be configured in the system has its own process. This allows very easily exchanging the way processes are represented by equations without having to alter many equations. For example, if instead of LinearOLR we have provided

OLR_process = EmissivityStefanBoltzmanOLR()
lhs(OLR_process) ~ rhs(OLR_process)

\[ \begin{equation} \mathrm{OLR}\left( t \right) = 5.6704 \cdot 10^{-8} \left( T\left( t \right) \right)^{4} \varepsilon\left( t \right) \end{equation} \]

then we would have used a Stefan-Boltzmann grey-atmosphere representation for the outgoing longwave radiation. Notice how this introduced an additional variable $\varepsilon$.

Default processes

Hold on a minute though, because in the original processes we provided,

processes = [
    BasicRadiationBalance(),
    LinearOLR(),
    ParameterProcess(α),
    CO2Forcing(),
]

there was no process that defined for example the absorbed solar radiation $ASR$!

Well, ConceptualClimateModels.jl allows the concept of default processes. The package exports some submodules, and each submodule targets a particular application of conceptual climate models. In this Tutorial we are using the GlobalMeanEBM submodule, which provides functionality to model global mean energy balance models.

Each submodule defines and exports its own list of predefined symbolic variables. When we wrote

using ConceptualClimateModels.GlobalMeanEBM

we brought into scope all the variables that this (sub)module defines and exports, such as T, α, OLR.

Each (sub)module and provides a list of predefined processes for its predefined symbolic variables. These predefined default processes are loaded automatically when we provide the (sub)module as a second argument to processes_to_coupledodes, which we did above. In this (sub)module, the default process for the $ASR$ is $ASR = S(1-\alpha)$ with $S$ the solar constant.

The function processes_to_coupledodes goes through all processes the user provided and identifies variables that themselves do not have a process. It then checks the list of default processes and attempts to assign one to these variables. If there are no default processes for some variables, it makes the variables themselves parameters with the same name but with a subscript 0 if the variables have a default value.

For example, let's assume that we completely remove default processes and we don't specify a process for the albedo $α$:

processes = [
    BasicRadiationBalance(),
    LinearOLR(),
    CO2Forcing(), # for default CO2 value this is zero forcing
    ASR ~ S*(1-α), # add the processes for ASR, but not for S or α
]

4-element Vector{Any}:
 BasicRadiationBalance(5.0e8, f(t), T(t), ASR(t), OLR(t))
 OLR(t) ~ -277.0 + 1.8T(t)
 CO2Forcing(f(t), CO2(t), 3.7)
 ASR(t) ~ S(t)*(1 - α(t))

which we'll give to processes_to_mtkmodel. Notice how there is no second argument given, which would normally be the (sub)module to obtain default processes from.

side note: we use `processtomtkmodelhere, as we don't care yet for making aDynamicalSystem`; just the symbolic model representation

mtk = processes_to_mtkmodel(processes)
# we access the equations directly from the model
equations(mtk)

\[ \begin{align} \frac{\mathrm{d} T\left( t \right)}{\mathrm{d}t} \tau_{T} =& 0.002939 \left( \mathrm{ASR}\left( t \right) + f\left( t \right) - \mathrm{OLR}\left( t \right) \right) \\ \mathrm{OLR}\left( t \right) =& -277 + 1.8 T\left( t \right) \\ f\left( t \right) =& 3.7 \log_{2}\left( \frac{1}{400} \mathrm{CO2}\left( t \right) \right) \\ \mathrm{ASR}\left( t \right) =& S\left( t \right) \left( 1 - \alpha\left( t \right) \right) \\ \mathrm{CO2}\left( t \right) =& CO2_{0} \\ S\left( t \right) =& S_{0} \\ \alpha\left( t \right) =& \alpha_{0} \end{align} \]

You will notice the equation $α = α_0$ where $\alpha_0$ is now a parameter of the system (i.e., it can be altered after creating the system). The value of $\alpha_0$ is the default value of $\alpha$:

default_value(α)

0.3

# current value of the _parameter_ α_0 (not the variable!)
default_value(mtk.α_0)

0.3

When this automation occurs a warning is thrown:

┌ Warning: Variable α(t) was introduced in process of variable ASR(t).
│ However, a process for α(t) was not provided,
│ and there is no default process for it either.
│ Since it has a default value, we make it a parameter by adding a process:
└ `ParameterProcess(α)`.

ParameterProcess is the most trivial process: it simply means that the corresponding variable does not have any actual process describing it and rather it is a system parameter.

This automation does not occur if there is no default value. For example, the $ASR$ variable does not have a default value. If we have not assigned a process for $ASR$, the system construction would error instead:

processes = [
    BasicRadiationBalance(),
    LinearOLR(),
    CO2Forcing(),
]

mtk = processes_to_mtkmodel(processes)

ERROR: ArgumentError: Variable ASR(t) was introduced in process of
variable T(t). However, a process for ASR(t) was not provided,
there is no default process for ASR(t), and ASR(t) doesn't have a default value.
Please provide a process for variable ASR(t).

These warnings and errors are always "perfectly" informative. They tell us exactly which variable does not have a process, and exactly which other process introduced the process-less variable first. This drastically improves the modelling experience, especially when large and complex models are being created.

Adding your own processes

Each of the submodules of ConceptualClimateModels.jl provides an increasing list of predefined processes that you can use out of the box to compose climate models. The predefined processes all come from existing literature and cite their source via BiBTeX.

It is also very easy to make new processes on your own. The simplest way to make a process is to just provide an equation for it with the l.h.s. of the equation being the variable the process defines.

For example:

@variables x(t) = 0.5 # all variables must be functions of (t)
x_process = x ~ 0.5*T^2 # x is just a function of temperature

\[ \begin{equation} x\left( t \right) = 0.5 \left( T\left( t \right) \right)^{2} \end{equation} \]

A more re-usable approach however is to create a function that generates a process or create a new process type as we describe in making new processes.

A note on symbolic variable instances

Recall that when we wrote

using ConceptualClimateModels.GlobalMeanEBM

at the start of this tutorial, we brought into scope variables that this (sub)module defines and exports, such as T, α, OLR. They are listed on the (sub)module's documentation page, and are used in that module's default processes.

In some predefined processes documentation you will notice call signatures like

BasicRadiationBalance(; T, f, kwargs...)

There are keywords that do not have an assignment like T, f above. This means that use the (sub)module's predefined variables.

Crucially, these default variables are symbolic variables. They are defined as

@variables begin
    T(t) = 0.5 # ...
    # ...
end

which means that expressions that involve them result in symbolic expressions, for example

A2 = 0.5
B2 = 0.5
OLR2 = A2 + B2*T

\[ \begin{equation} 0.5 + 0.5 T\left( t \right) \end{equation} \]

In contrast, if we did instead

T2 = 0.5 # _not_ symbolic!
OLR3 = A2 + B2*T2

0.75

This OLR3 is not a symbolic expression and cannot be used to represent a process.

You can use your own variables in any predefined processes. You can define them by doing

@variables begin
    (T1_tropics(t) = 290.0), [bounds = (200.0, 350.0), description = "temperature in tropical box 1, in Kelvin"]
end

\[ \begin{equation} \left[ \begin{array}{c} \mathrm{T1}_{tropics}\left( t \right) \\ \end{array} \right] \end{equation} \]

and then assign them to the corresponding keyword argument

process = BasicRadiationBalance(T = T1_tropics)

BasicRadiationBalance(5.0e8, f(t), T1_tropics(t), ASR(t), OLR(t))

Defining variables with the extra bounds, description annotations is useful for integrating with the rest of the functionality of the library, and therefore it is strongly recommended.

Default values, limits, etc.

All predefined variables that could be dynamic variables (i.e., could have a time derivative applied to them) have a default value, a description, and plausible physical bounds.

To obtain the default value, use default_value(x). For the description, use getdescription(x). For the bounds, see plausible_limits.

Automatically named parameters

The majority of predefined processes create symbolic parameters that are automatically named based on the variables governing the processes. This default behaviour can be altered in two ways.

For example, IceAlbedoFeedback adds named parameters to the equations whose name is derived from the name of the variable representing ice albedo:

@variables my_ice_α(t) = 0.1 # don't forget the `(t)`!
ice_process = IceAlbedoFeedback(; α_ice = my_ice_α)
processes = [ice_process]

mtk = processes_to_mtkmodel(processes)
equations(mtk)

\[ \begin{align} \mathrm{my}_{ice\_\alpha}\left( t \right) =& max_{my\_ice\_\alpha} + 0.5 \left( - max_{my\_ice\_\alpha} + min_{my\_ice\_\alpha} \right) \left( 1 + \tanh\left( \frac{2 \left( - Tfreeze_{my\_ice\_\alpha} + \frac{1}{2} Tscale_{my\_ice\_\alpha} + T\left( t \right) \right)}{Tscale_{my\_ice\_\alpha}} \right) \right) \\ T\left( t \right) =& T_{0} \end{align} \]

parameters(mtk)

5-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 max_my_ice_α
 Tfreeze_my_ice_α
 Tscale_my_ice_α
 min_my_ice_α
 T_0

We can alter this behaviour by either providing our own named parameters to one of the keywords of the process, or wrapping a value around LiteralParameter to replace the parameter by a literal constant, like so:

@parameters myfreeze = 260.0
ice_process = IceAlbedoFeedback(;
    α_ice = my_ice_α,
    Tfreeze = myfreeze, # my custom parameter
    max = LiteralParameter(0.9) # don't make a parameter
)

mtk = processes_to_mtkmodel([ice_process])
equations(mtk)

\[ \begin{align} \mathrm{my}_{ice\_\alpha}\left( t \right) =& 0.9 + 0.5 \left( -0.9 + min_{my\_ice\_\alpha} \right) \left( 1 + \tanh\left( \frac{2 \left( \frac{1}{2} Tscale_{my\_ice\_\alpha} - myfreeze + T\left( t \right) \right)}{Tscale_{my\_ice\_\alpha}} \right) \right) \\ T\left( t \right) =& T_{0} \end{align} \]

parameters(mtk)

4-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 myfreeze
 Tscale_my_ice_α
 min_my_ice_α
 T_0

Integration with DynamicalSystems.jl

ConceptualClimateModels.jl integrates with DynamicalSystems.jl. It provides a summary function dynamical_system_summary, and also provides initial condition sampling to use when e.g., finding attractors and their basin fractions with DynamicalSystems.basins_fractions, via the function plausible_limits. Moreover, since all dynamical systems generated by ConceptualClimateModels.jl have symbolic bindings, one can use the symbolic variables in e.g., interactive GUI exploration or to access or set the parameters of the system.

plausible_limits(x::Num)

plausible_limits(ds::DynamicalSystem [, idxs])

plausible_ic_sampler(ds::DynamicalSystem)

plausible_grid(ds::DynamicalSystem, n = 101)

dynamical_system_summary(ds::DynamicalSystem)

API Reference

processes_to_coupledodes(processes [, default]; kw...)

processes_to_mtkmodel(processes::Vector [, default]; kw...)

Making new processes

To make a new processes you can:

1. Create a function that given some keyword arguments uses one of the existing generic processes to make and return a process instance. Or, it can return an equation directly, provided it satisfies the format of processes_to_mtkmodel. For an example of this, see the source code of SeparatedClearAllSkyAlbedo or EmissivityFeedbackTanh.
2. Create a new Process subtype. This is preferred, because it leads to much better printing/display of the list of processes. For an example of this, see the source code of IceAlbedoFeedback. To create a new Process see the API of ProcessBasedModelling.jl or read the documentation string of Process below.

Process

register_default_process!(process, m::Module; warn = true)

module MyProcesses
# ...

function __init__()
    register_default_process!.(
        [
            process1,
            process2,
            # ...
        ],
        Ref(MyProcesses)
    )
end

end # module

Generic processes

Processes that do not depend on any particular physical concept and instead provide a simple way to create new processes for a given climate variable:

ParameterProcess(variable, value = default_value(variable)) <: Process

@variables T(t) = 0.5
proc = ParameterProcess(T)

TimeDerivative(variable, expression [, τ])

ExpRelaxation(variable, expression [, τ]) <: Process

τn*Differential(t)(variable) ~ expression - variable

ExpRelaxation(process, τ)

AdditionProcess(process, added...)

SigmoidProcess <: Process
SigmoidProcess(variable, driver; left, right, scale, reference)

variable ~ left + (right - left)*(1 + tanh(2(driver - reference)/scale))*0.5