Tutorial

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, 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. 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.CCMV

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

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

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)
 CO2(t) ~ CO2_0
 S(t) ~ 1.0
 α(t) ~ α_0
 f(t) ~ 3.7log2((1//400)*CO2(t))
 ASR(t) ~ 340.25S(t)*(1 - α(t))

and parameter values:
 τ_T ~ 1.4695077149155033e6
 α_0 ~ 0.3
 CO2_0 ~ 400.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) \\ \mathrm{CO2}\left( t \right) =& CO2_{0} \\ S\left( t \right) =& 1 \\ \alpha\left( t \right) =& \alpha_{0} \\ f\left( t \right) =& 3.7 \log_{2}\left( \frac{1}{400} \mathrm{CO2}\left( t \right) \right) \\ \mathrm{ASR}\left( t \right) =& 340.25 S\left( t \right) \left( 1 - \alpha\left( t \right) \right) \end{align} \]

# show parameters
parameters(mtk)

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

The symbolic variables and parameters can also be used to query or alter the dynamical system. For example, we can obtain, or alter, any parameter by providing the symbolic parameter index. There are multiple ways to obtain the symbolic index provided we know its name. First, we can re-create the symbolic parameter:

index = first(@parameters CO2_0)

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

Second, we can use the retrieved MTK model and access its CO2_0 parameter:

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$ by:

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

290.0

or obtain the $OLR$ (outgoing longwave radiation)

observe_state(ds, :OLR)

245.0

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 composed 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.

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 the absorbed solar radiation $ASR$!

Well, ConceptualClimateModels.jl has a list of predefined processes that are automatically included in every call to processes_to_coupledodes. The default processes 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 attempt to assign one to these variables.

If there are no default processes, it makes the variables themselves parameters with the same name but with a subscript 0.

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(), # note that for default CO2 values this is zero forcing
    ASR ~ S*(1-α), # add the processes for ASR, but not for S or α
]

# note the empty list as 2nd argument, which is otherwise
# the default processes. Notice also that we make an MTK model
# (which is the step _before_ converting to a dynamical system)
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 α was introduced in process of variable ASR(t).
│ However, a process for α 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(α)`.
└ @ ProcessBasedModelling ...\ProcessBasedModelling\src\make.jl:65

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

This automation does not occur if there is no default value. For example, variables that can never be dynamic state variables, such as $ASR$, do 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 was introduced in process of
variable T(t). However, a process for ASR was not provided,
there is no default process for ASR, and ASR doesn't have a default value.
Please provide a process for variable ASR.

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 makes the modelling experience stress-free, especially when large and complex models are being created.

Adding your own processes

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.

Premade symbolic variable instances

You might be wondering, when we wrote the equation ASR ~ S*(1-α) for the $ASR$ process, or when we wrote x ~ 0.5 * T^2, where did the variable bindings ASR, S, α come from? For convenience, ConceptualClimateModels.jl defines some symbolic variables for typical climate quantities and assigns default processes to them. we brought all of these into scope when we did

using ConceptualClimateModels.CCMV

where CCMV (standing for ConceptualClimateModels Variables) is a submodule that defines and exports the variables. We list all of these below. These default variables are used throughout the library as the default variables in predefined processes. When going through documentation strings of predefined processes, such as BasicRadiationBalance, you will notice that the function call signatures are like:

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

There are keywords that do not have an assignment like T, f above. These always represent climate variables, never parameters, and for the variables they use the existing predefined climate 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,

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!
OLR2 = A2 + B2*T2

0.75

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

You can use your own variables for any of the 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.

List of premade variables

The premade variables are not exported by default. To bring them into global scope you have to do:

using ConceptualClimateModels
using ConceptualClimateModels.CCMV

and then use them,

T, q, OLR

(T(t), q(t), OLR(t))

All the premade variables are:

PREDEFINED_CCM_VARIABLES

\[ \begin{equation} \left[ \begin{array}{c} T\left( t \right) \\ S\left( t \right) \\ f\left( t \right) \\ \alpha\left( t \right) \\ \alpha_{ice}\left( t \right) \\ \alpha_{cloud}\left( t \right) \\ {\Delta}T\left( t \right) \\ {\Delta}S\left( t \right) \\ \varepsilon\left( t \right) \\ C\left( t \right) \\ \ell\left( t \right) \\ \mathrm{CO2}\left( t \right) \\ q\left( t \right) \\ \mathrm{OLR}\left( t \right) \\ \mathrm{ASR}\left( t \right) \\ \end{array} \right] \end{equation} \]

Default values, limits, etc.

All premade variables 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, use:

physically_plausible_limits(x)

e.g.,

physically_plausible_limits(T)

(200.0, 350.0)

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 by providing initial condition sampling to use when e.g., finding attractors and their basin fractions with DynamicalSystems.basins_fractions, and with the function dynamical_system_summary. 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.

physically_plausible_limits(ds::DynamicalSystem)

physically_plausible_ic_sampler(ds::DynamicalSystem)

physically_plausible_grid(ds::DynamicalSystem, n = 101)

dynamical_system_summary(ds::DynamicalSystem)

API Reference

processes_to_coupledodes(processes, default = DEFAULT_PROCESSES; kw...)

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

Making new processes

To make a new processes you can:

1. Create a function that given some keyword arguments (including which symbolic variables to use) 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