`ConceptualClimateModels.ConceptualClimateModels`

— Module**ConceptualClimateModels.jl**

ConceptualClimateModels.jl is a Julia package for creating and analysing conceptual models of climate, such as energy balance models, glaciation cycle models, or climate tipping models. Such conceptual models are simplified representation of basic climate components, and the processes that connect them, such as flows of energy or mass. Within this context such models are typically coupled ordinary differential equations (with partial or stochastic DEs also being possible).

ConceptualClimateModels.jl accelerates both modelling and analysis aspects of working with such models by:

- Building upon ModelingToolkit.jl for creating equations from
*symbolic expressions*. - Utilizing ProcessBasedModelling.jl to provide a field-specific framework that allows easily testing different physical hypotheses regarding how climate variables couple to each other, or how climate processes are represented by equations.
- Offering many predefined processes from current literature and ongoing research. All predefined processes cite the literature rigorously using BiBTeX.
- Being easy to extend with more climate variables or physical processes.
- Allowing the straightforward coupling of different conceptual models with each other.
- Automating the naming of custom parameters relating to existing climate processes.
- Integrating with DynamicalSystems.jl to automate the start-up phase of typical nonlinear dynamics based workflows.

with other features coming soon, such as:

- Support for latitudinal models (where each variable is vector-valued over latitude circles)
- Support for stochastic ordinary differential equations

To install it, run `import Pkg; Pkg.add("ConceptualClimateModels")`

.

All further information is provided in the documentation, which you can either find online or build locally by running the `docs/make.jl`

file.

ConceptualClimateModels.jl development is funded by UKRI's Engineering and Physical Sciences Research Council, grant no. EP/Y01653X/1 (grant agreement for a EU Marie Sklodowska-Curie Postdoctoral Fellowship).

`ConceptualClimateModels.BasicRadiationBalance`

— Type`BasicRadiationBalance(; T, f, ASR, OLR, c_T = 5e8)`

Create the equation

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

representing the most basic radiative energy balance at the top of the atmosphere setting a global mean temperature, see e.g., any introductory article North1981, Ghil1981. `ASR`

is the absorbed solar radiation, which defaults to `S*(1 - α)`

in the default processes. `S`

is the received insolation, by default equal to `solar_constant`

, but could e.g., be any astronomical forcing such as `AstronomicalForcingDeSaedeleer`

. `α`

is the albedo and `f`

any additional forcing such as `CO2Forcing`

. `OLR`

defaults to `A + B*T`

. `c_T`

is the heat capacity of the system in J/K/m². However, for convenience, the parameter added to the final equation is `τ_T`

which is the timescale in seconds, i.e., `c_T/solar_constant`

.

`ConceptualClimateModels.CO2Forcing`

— Type`CO2Forcing(; f, CO2, CO2f = 3.7)`

Create the equation $f ~ CO2f \log_2(CO2/400)$ which describes the forcing added to the TOA energy balance due to CO2 concentrations, assumming the `OLR`

expression is calibrated for 0 added forcing at 400 ppm which is the default for OLR expressions provided by ConceptualClimateModels.jl.

The default value of $f$ comes from Eq. (3.2) of Bastiaansen2023 which cites IPCC-5, while Etminan2016 report practically the same value assuming a constant $f$ (note here the log is base 2). In reality $f$ depends on $CO2$ and other greenhouse gases concentrations due to spectral overlaps, see Etminan2016 Sec. 4.

`ConceptualClimateModels.CoAlbedoProduct`

— Type`CoAlbedoProduct(; α, albedo_variables = (α_ice, α_cloud))`

Create the equation `1 - α ~ prod(a -> (1 - a), albedo_variables)`

meaning that the co-albedo is the product of the co-albedos of all albedo variables. This would be e.g., the planetary albedo if all components were uniform layers, while the bottom-most layer (surface) had perfect absorption and all other layers had 0 absorption and finite reflection.

`ConceptualClimateModels.DirectAlbedoAddition`

— Type`DirectAlbedoAddition(; α, α_bg = 0.1, other_albedo_variables = (α_ice, α_clouds))`

Create the equation `α ~ α_bg + other_albedo_variables...`

, meaning that planetary albedo `α`

is a direct sum of all specified albedos.

`ConceptualClimateModels.IceAlbedoFeedback`

— Type```
IceAlbedoFeedback(; T, α_ice,
max = 0.45, min = 0.1, Tscale = 10, Tfreeze = 275.15, τ = 0
)
```

Create an equation that assigns ice albedo `α_ice`

to a hyperbolic tangent of temperature `T`

. This represents an approximately linear decrease with `T`

, as ice melts over part of the earth, while it is constant for all `T`

for which the earth would be either entirely ice covered (`T < Tfreeze - scale`

) or ice free (`T > Tfreeze`

).

In essence this is a `TanhProcess`

with the given keywords as parameters with reference temperature `Tref = Tfreeze - scale/2`

.

This albedo is the most common used large-scale feedback in energy balance models, e.g., Ghil1981, although it is typically taken as a piece-wise linear function. There is little change with using a hyperbolic tangent instead, while the `tanh`

offers a differentiable flow.

The timescale `τ`

if not zero will make an `ExpRelaxation`

process relaxing to the hyperbolic tangent.

`ConceptualClimateModels.TanhProcess`

— Type```
TanhProcess(variable, driver, left, right, scale, reference) <: Process
TanhProcess(variable, driver; left, right, scale, reference) <: Process
```

A common process for when a `variable`

has a tanh-dependence on a `driver`

variable. The rest of the input arguments should be real numbers or `@parameter`

named parameters.

The process creates the expression:

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

i.e., a tanh formula that goes from value `left`

to value `right`

as a function of `driver`

over a range of `scale`

being centered at `reference`

.

If the values given to the parameters of the expression are real numbers, they become named parameters prefixed with the name of `variable`

, then the name of the `driver`

, and then `_tanh_left`

, `_tanh_right`

, `_tanh_rate`

and `_tanh_ref`

respectively. Use `LiteralParameter`

for parameters you do not wish to rename.

`ConceptualClimateModels.AbsoluteHumidityIGLCRH`

— Method`AbsoluteHumidityIGLCRH(; q = q, T = T, RH = 0.5)`

Create the equation

\[q = m_v \cdot RH \cdot e_s / (R\cdot T)\]

which approximates absolute humidity (water vapor concentration) in gr/m³ under two assumptions:

- ideal gas law
- constant relative humidity
`RH`

with warming, i.e.,`RH`

is a constant, so that`q = q(T)`

*only*.

Here `m_v = 18.016`

gr is the specific mass of water, `R`

the ideal gas constant, `e_s =`

`saturation_vapor_pressure`

`(T)`

in kPa.

Note that you could create a process for a *variable* `RH`

and pass the variable to the `RH`

keyword to skip assumption (2).

Use `absolute_humidity_igl(T, RH)`

to obtain `q`

at a given `T, RH`

.

`ConceptualClimateModels.AstronomicalForcingDeSaedeleer`

— Method`AstronomicalForcingDeSaedeleer(; S = S, extensive = false)`

Create the equation `S ~ astronomical_forcing_desaedeleer(t, extensive)`

which is Eq. (1) of DeSaedeleer2013:

`S = \sum_i s_i \sin(\omega_i t) + c_i \cos(\omega_i t)`

where the values of $\omega_i, s_i, c_i$ come from Berger1978 who performed a spectral expansion of the insolation. The validity range of this approximation is [-1, 0] Myr.

In the summation $i$ goes up to 35 if `extensive`

, otherwise up to 8. The components are sorted according to magnitude of the spectral line, so the default version has only the 8 most important spectra lines.

Note that in contrast to Eq. (1) of DeSaedeleer2013 we do not normalize $f$ and its value is in W/m² (the mean value is still deducted). Additionally, the values of $\omega_i$ have been adjusted to expect time in units of seconds.

`ConceptualClimateModels.BudykoOLR`

— Method```
BudykoOLR(; OLR=OLR, T=T, C=C,
BudykoOLR_A = -461.8068, BudykoOLR_B = 2.58978,
BudykoOLR_Ac = -377.22741, BudykoOLR_Bc = 1.536171
)
```

Create the equation `OLR ~ A + B*T - C*(Ac + Bc*T)`

for the dependence of OLR on both temperature and cloud fraction (in 0-1). This is the same as Eq. (1) of Budyko1969. However, here `T`

is expected in Kelvin, and the coefficients have been extracted by fitting into CERES data in the same way as in `LinearOLR`

.

`ConceptualClimateModels.CloudAlbedoExponential`

— Method```
CloudAlbedoExponential(
α_cloud, C, a = 2.499219232848238, b = 17.596369331717433
)
```

Create the equation `α_cloud ~ sinh(a*C)/b`

relating cloud albedo to cloud fraction `C`

. This equation is exponential and not linear, as in observations. Engstrom2015 (and also Bender2017) discuss this exponential relation in detail, and provide as explanation that cloud effective albedo increases with latitude (due to solar zenith changes) while cloud fraction also increases with latitude.

Note that here however we modify the equation `α_cloud ~ exp(a*C - b)`

of Engstrom2015 to utilize the hyperbolic sine, so that `α_cloud = 0`

when `C = 0`

as is physically necessary. Then, `a, b`

are extracted by fitting CERES data, using as `α_cloud`

the energetically consistent cloud albedo as defined by Datseris2021, further yearly averaged and within latitudes (-60, 60) as in Bender2017. This albedo can be directly added to the clear sky albedo to produce the planetary albedo.

`ConceptualClimateModels.CloudAlbedoLinear`

— Method`CloudAlbedoLinear(; α_cloud, C, a = 0.2252861764703435)`

Same as in `CloudAlbedoExponential`

, but now the linear form `α_cloud ~ a*C`

is returned, with `a`

fitted from CERES data in the same way.

`ConceptualClimateModels.EmissivityFeedbackTanh`

— Method`EmissivityFeedbackTanh(; ε, Τ, left = 0.5, right = 0.4, rate = 0.5, Tref = 288.0)`

Create an equation that assigns emissivity `ε`

to hyperbolic tangent of temperature `T`

. This is an ad-hoc feedback that was used in Bastiaansen2023, similar to `EmissivitySellers1969`

. In essence this is a `TanhProcess`

with the given keywords as parameters.

`ConceptualClimateModels.EmissivityLogConcentration`

— Method```
EmissivityLogConcentration(; ε, T, CO2, q,
ε_ref = 0.6, ε_CO2 = 0.1, ε_H20 = 0.01
)
```

Create the equation

\[ε ~ ε_{ref} - ε_{CO2}*log2(CO2/400) - ε_{H2O}*log2(q/1e-3)\]

This equation is inspired by Eq. A.9 of VonderHeydt2016, which had only the first two terms. Here we added the last term to represent feedback due to water vapor. `q`

is absolute humidity (or, water vapor concentration), whose default process is `AbsoluteHumidityIGLCRH`

.

The logarithm of either water vapor or CO2 concentration is used because literature suggests that the radiative forcing of greenhouse gases scales with the logarithm of their concentration Huang2014. As `ε`

is a first order component of the OLR, it needs to scale linearly with the logarithm.

`ConceptualClimateModels.EmissivitySellers1969`

— Method`EmissivitySellers1969(; ε, T, m = 0.5)`

Create the equation `ε ~ 1 - m*tanh(19*T^6*1e-16)`

, which was used originally in Sellers1969 to represent the effect of "water vapor, carbon dioxide, dust and clouds on terrestrial radiation".

`ConceptualClimateModels.EmissivityStefanBoltzmanOLR`

— Method`EmissivityStefanBoltzmanOLR(; ε, T)`

Create the equation `OLR ~ ε*σ*T^4`

where `σ`

is the Stefan Boltzmann constant and `ε`

the effective emissivity, also known as the "grayness" of the system, or the deviation it has from being a perfect black body Ghil1981. `ε`

then needs to be parameterized itself to include greenhouse or other climate effects.

`ConceptualClimateModels.LinearClearSkyOLR`

— Method`LinearClearSkyOLR(; kw...)`

Equivalent with `LinearOLR(; A = A = -326.0, B = 2.09, kw...)`

and provided as a convenience for the clear sky fit to CERES data.

`ConceptualClimateModels.LinearOLR`

— Method`LinearOLR(; OLR, T, A = -277.0, B = 1.8)`

Create the equation `OLR ~ A + B*T`

. This is a linearized outgoing longwave radiation (OLR), and is the same equation as (7) of North1981: $OLR = A + BT$ with $T$ temperature in Kelvin and $A, B$ constants. However, default $A, B$ are fitted from current CERES all sky OLR and using ERA5 data for the 2-meter temperature. We assume `T`

in Kelvin. This linear approximation is quite accurate for temporally averaged data $T \in (220, 280)$ however drops drastically in accuracy after that due to the nonlinear effects of clouds (as evident by observational data).

Koll2018 provide a "proof" of the linearity of the clear sky OLR due to spectral properties of water vapor.

We note a big difference between current CERES data and the values reported in North1981: here `A=214.67`

(assuming $T$ in Celcius) and `B=1.8`

versus the values `A=203.3`

and `B=2.09`

in North1981.

If instead of all sky, if we fit the clear sky CERES data, we get `A = -326.0, B = 2.09`

. Interestingly, coefficient `B`

here is the same as that reported by North1981, but `A=244.88`

(assuming `T`

in Celcius) is not.

`ConceptualClimateModels.SeparatedClearAllSkyAlbedo`

— Method`SeparatedClearAllSkyAlbedo(; α, α_cloud, C, α_clr = 0.15)`

Create the equation `α ~ α_cloud*C + α_clr*(1 - C)`

.

Bender2017 argue that one can assume a separation between clear-sky and cloud albedo, so that `α = α_cloud*C + α_clr*(1 - C)`

with `C`

the cloud fraction and `α_clr`

the clear sky albedo. They further cite Cess1976 to facilitate the claim Additionally, Eq. (20) of Barker1999 provides an identical expression.

In most cases you want to provide a variable with its own process for `α_clr`

.

`ConceptualClimateModels.SoedergrenClearSkyEmissivity`

— Method`SoedergrenClearSkyEmissivity(; ε, T, CO2, RH = 0.8, H_H20 = 2.0)`

Create Eq. 10 of Soedergren2018, which is the same as Eq. 21 of Barker1999 for the effective emissivity of clear sky atmosphere:

\[\varepsilon = 1 - \exp(0.082 - (2.38*0.1*e_s*RH*H_{H2O} + 40.3*CO2*1e-6)^{0.294})\]

with $e_s$ the `saturation_vapor_pressure`

. The equation assumes CO2 concentration is in ppm and vapor pressure in kPa hence the conversion factors 0.1 and 1e-6.

Be advised: this process is included for reference only. It should not be used because it is physically wrong. Emissivity increases with temperature, while it should decrease: higher temperature → stronger greenhouse effect → smaller effective emissivity required for higher temperature to have the same OLR as per the basic equation $OLR = ε σ Τ^4$.

`ConceptualClimateModels.dynamical_system_summary`

— Method`dynamical_system_summary(ds::DynamicalSystem)`

Return a printable/writable string containing a summary of `ds`

, which outlines its current status and lists all symbolic equations and parameters that constitute the system, if a referrence to a ModelingToolkit.jl exists in `ds`

.

`ConceptualClimateModels.physically_plausible_grid`

— Function`physically_plausible_grid(ds::DynamicalSystem, n = 101)`

Return a `grid`

that is a tuple of `range`

objects that each spans the `physically_plausible_limits`

of `ds`

. `n`

can be an integer or a vector of integers (for each dimension). The resulting `grid`

is useful to give to `DynamicalSystems.AttractorsViaRecurrences`

.

`ConceptualClimateModels.physically_plausible_ic_sampler`

— Method`physically_plausible_ic_sampler(ds::DynamicalSystem)`

Return a `sampler`

that can be called as a 0-argument function `sampler()`

, which yields random initial conditions within the hyperrectangle defined by the `physically_plausible_limits`

of `ds`

. The `sampler`

is useful to give to e.g., `DynamicalSystems.basins_fractions`

.

`ConceptualClimateModels.physically_plausible_limits`

— Method`physically_plausible_limits(x)`

Return a tuple (min, max) of plausible limiting values for the variable `x`

. If the variable does not have defined `bounds`

metadata, then the default value ± 20% is used. If there is no default value, a heuristic is tried, and an error is thrown if it fails.

`ConceptualClimateModels.physically_plausible_limits`

— Method`physically_plausible_limits(ds::DynamicalSystem)`

Return a vector of limits (min, max) for each dynamic state variable in `ds`

, assumming `ds`

has been made using variables with bounds (all default symbolic variables of ConceptualClimateModels.jl satisfy this).

`ConceptualClimateModels.processes_to_coupledodes`

— Function`processes_to_coupledodes(processes, default = DEFAULT_PROCESSES; kw...)`

Convert a given `Vector`

of processes to a `DynamicalSystem`

, in particular `CoupledODEs`

. All processes represent symbolic equations managed by ModelingToolkit.jl. `default`

is a vector for default processes that "process-less" variables introduced in `processes`

will obtain. Use `processes_to_mtkmodel`

to obtain the MTK model before it is structurally simplified and converted to a `DynamicalSystem`

. See also `processes_to_mtkmodel`

for more details on what `processes`

is, or see the online Tutorial.

**Keyword arguments**

`diffeq`

: options passed to DifferentialEquations.jl ODE solving when constructing the`CoupledODEs`

.`inplace`

: whether the dynamical system will be in place or not. Defaults to`true`

if the system dimension is ≤ 5.`split = false`

: whether to split parameters as per ModelingToolkit.jl. Note the default is not ModelingToolkit's default, i.e., no splitting occurs. This accelerates parameter access, assuming all parameters are of the same type.`kw...`

: all other keywords are propagated to`processes_to_mtkmodel`

.

`ConceptualClimateModels.saturation_vapor_pressure`

— Method`saturation_vapor_pressure(T)`

Given surface temperature (in K) return saturation pressure for water vapor (in kPa) using the Tetens-Murray formula from TetensEquation2023, which is `A*exp(B*T/(C+T))`

. Different formula is used for when `T`

is less than the freezing point.

`ConceptualClimateModels.ΔTLinearRelaxation`

— Method`ΔTLinearRelaxation(; ΔT, T, τ = 5e6, A = 36.53, B = 0.658)`

Create the equation

\[\tau_{\Delta T}\frac{\Delta T}{dt} = \Delta T_{ref}(T) - \Delta T\]

which exponentially relaxes the equator-to-pole temperature difference `ΔT`

to its reference value $\Delta T_{ref}(T) = A - B*(T - 275.15)$, i.e., it decreases linearly with global mean temperature $T$ (in Kelvin). The default values for `A, B`

are obtained from Equation (2) of Gaskell2022. We also fitted paleoclimate data of Osman2021 and found very similar results, `A = 35.8, B = -1.11`

for north hemisphere and `A = 27.4, b = -0.513`

for south.

Here `ΔT`

is defined as the temperature difference between average temperatures at (0, 30) and (60, 90) latitudes. The timescale is taken as 2 months, although if `τ = 0`

is given, the equation $\Delta T ~ \Delta T_{ref}(T)$ is created instead.

`ConceptualClimateModels.ΔTStommelModel`

— Method`ΔTStommelModel(; ΔT=ΔT, ΔS=ΔS, η1 = 2, η2 = 1, η3 = 0.3)`

Create the equations

\[\dot{\Delta T} = \eta_1 - \Delta T - |\Delta T - \Delta S| \Delta T \dot{\Delta S} = \eta_2 - \eta_3\Delta S - |\Delta T - \Delta S| \Delta S\]

which are the two equations of the Stommel box model for Atlantic thermohaline circulation Stommel1961, here presented in nondimensionalized form Lohmann2021, so that temperature and sality are normalized by their coefficients $a_T, a_S$ relating them to the density of water

\[\rho = \rho_0 [1 - a_T(T - T_0) + a_S(S-S_0)]\]

for some reference values.