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

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.

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.

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.

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.

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.

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.

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:

1. ideal gas law
2. 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.

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.

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.

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.

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.

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.

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.

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

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.

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.

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.

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.

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.

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.

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.

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.

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.

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

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.

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.

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

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