Climate

Climate: energy balance, variability, glacial cycles, low-order atmospheres, net energy, tipping and its warnings.

CarbonCycleElement(label, visual=None, *, c=83.58, w=2050.0, mu=250.0, b=4.0, theta=5.0, c_x=55.0, c_p=110.0, nu=0.0, w0=2000.0, gamma=4.0, f0=0.694, c_f=43.9, beta=1.7)

Bases: Element

Ocean of Rothman's excitable carbon cycle, its carbonate ion concentration and dissolved inorganic carbon.

The carbonate ion concentration c and the dissolved inorganic carbon w, in μmol kg⁻¹, follow dc/dt = f(c) (mu (1 - b s(c, c_p) - theta (1 - s(c, c_x)) - nu) + w - w0) and dw/dt = mu (1 - b s(c, c_p) + theta (1 - s(c, c_x)) + nu) - w + w0, with the Hill functions s(c, c*) = c^gamma / (c^gamma + c*^gamma) of burial and respiration and the buffer f(c) = f0 c^beta / (c^beta + c_f^beta).

Parameters:

Name	Type	Description	Default
label	str	The label of the ocean.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
c	Annotated[float, StateSpec(domain=open(0, oo))]	The initial carbonate ion concentration.	83.58
w	Annotated[float, StateSpec(domain=open(0, oo))]	The initial dissolved inorganic carbon.	2050.0
mu	Annotated[float, ParameterSpec()]	The characteristic rate of the carbon fluxes.	250.0
b	Annotated[float, ParameterSpec()]	The maximal burial of organic carbon.	4.0
theta	Annotated[float, ParameterSpec()]	The maximal respiration of organic carbon.	5.0
c_x	Annotated[float, ParameterSpec()]	The crossover carbonate concentration of respiration, which sets the excitability.	55.0
c_p	Annotated[float, ParameterSpec()]	The crossover carbonate concentration of burial.	110.0
nu	Annotated[float, ParameterSpec()]	The rate of external carbon injection.	0.0
w0	Annotated[float, ParameterSpec()]	The input of dissolved inorganic carbon by weathering and volcanism.	2000.0
gamma	Annotated[float, ParameterSpec()]	The Hill coefficient of burial and respiration.	4.0
f0	Annotated[float, ParameterSpec()]	The maximal buffer of the carbonate ion concentration.	0.694
c_f	Annotated[float, ParameterSpec()]	The carbonate concentration where the buffer is half its maximum.	43.9
beta	Annotated[float, ParameterSpec()]	The Hill coefficient of the buffer.	1.7

Examples:

>>> from skcomplex.models.physics.climate import CarbonCycleElement
>>> CarbonCycleElement('ocean').c_x
55.0

Source code in src/skcomplex/models/physics/climate/_carbon_cycle.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    c: Annotated[float, StateSpec(domain=RealSet.open(0, oo))] = 83.58,
    w: Annotated[float, StateSpec(domain=RealSet.open(0, oo))] = 2050.0,
    mu: Annotated[float, ParameterSpec()] = 250.0,
    b: Annotated[float, ParameterSpec()] = 4.0,
    theta: Annotated[float, ParameterSpec()] = 5.0,
    c_x: Annotated[float, ParameterSpec()] = 55.0,
    c_p: Annotated[float, ParameterSpec()] = 110.0,
    nu: Annotated[float, ParameterSpec()] = 0.0,
    w0: Annotated[float, ParameterSpec()] = 2000.0,
    gamma: Annotated[float, ParameterSpec()] = 4.0,
    f0: Annotated[float, ParameterSpec()] = 0.694,
    c_f: Annotated[float, ParameterSpec()] = 43.9,
    beta: Annotated[float, ParameterSpec()] = 1.7,
) -> None:
    super().__init__(label, visual=visual)
    self.c = c
    self.w = w
    self.mu = mu
    self.b = b
    self.theta = theta
    self.c_x = c_x
    self.c_p = c_p
    self.nu = nu
    self.w0 = w0
    self.gamma = gamma
    self.f0 = f0
    self.c_f = c_f
    self.beta = beta

ClimateAnomalyElement(label, visual=None, *, T=0.0, lam=0.5, sigma=0.3)

Bases: Element

Slow temperature anomaly of Hasselmann's stochastic climate, with its feedback rate and weather noise.

The relaxation of hasselmann damps the anomaly at the rate lam and drives it with weather noise of amplitude sigma.

Parameters:

Name	Type	Description	Default
label	str	The label of the anomaly.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
T	Annotated[float, StateSpec()]	The initial temperature anomaly.	0.0
lam	Annotated[float, ParameterSpec(domain=open(0, oo))]	The feedback rate that damps the anomaly.	0.5
sigma	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The amplitude of the weather noise.	0.3

Examples:

>>> from skcomplex.models.physics.climate import ClimateAnomalyElement
>>> ClimateAnomalyElement('ocean').lam
0.5

Source code in src/skcomplex/models/physics/climate/_climate_anomaly.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    T: Annotated[float, StateSpec()] = 0.0,
    lam: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    sigma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.3,
) -> None:
    super().__init__(label, visual=visual)
    self.T = T
    self.lam = lam
    self.sigma = sigma

DaisyworldElement(label, visual=None, *, white=0.2, black=0.2, L=1.0, S=917.0, sigma=5.67e-08, q=2060000000.0, p=1.0, gamma=0.3, k=0.003265, optimum=295.5, albedo_ground=0.5, albedo_white=0.75, albedo_black=0.25)

Bases: Element

Planet of Daisyworld, whose white and black daisies regulate its temperature through its albedo.

The areas covered by white and black daisies grow on the bare ground at a rate parabolic in their local temperature, with an optimum at optimum, and die at the rate gamma. The effective temperature follows the planetary albedo, sigma T^4 = S L (1 - A), and the local temperature of a daisy is warmer or cooler than it by q times the difference of their albedos.

Parameters:

Name	Type	Description	Default
label	str	The label of the planet.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
white	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial area fraction of white daisies.	0.2
black	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial area fraction of black daisies.	0.2
L	Annotated[float, ParameterSpec()]	The luminosity of the sun relative to its present value.	1.0
S	Annotated[float, ParameterSpec()]	The flux of the present sun in W m⁻².	917.0
sigma	Annotated[float, ParameterSpec()]	The Stefan-Boltzmann constant in W m⁻² K⁻⁴.	5.67e-08
q	Annotated[float, ParameterSpec()]	The transfer of heat between areas of different albedo, in K⁴.	2060000000.0
p	Annotated[float, ParameterSpec()]	The fraction of the planet fertile for daisies.	1.0
gamma	Annotated[float, ParameterSpec()]	The death rate of daisies.	0.3
k	Annotated[float, ParameterSpec()]	The curvature of the parabola of growth in K⁻².	0.003265
optimum	Annotated[float, ParameterSpec()]	The local temperature of the fastest growth in kelvin.	295.5
albedo_ground	Annotated[float, ParameterSpec()]	The albedo of the bare ground.	0.5
albedo_white	Annotated[float, ParameterSpec()]	The albedo of white daisies.	0.75
albedo_black	Annotated[float, ParameterSpec()]	The albedo of black daisies.	0.25

Examples:

>>> from skcomplex.models.physics.climate import DaisyworldElement
>>> DaisyworldElement('planet').albedo_white
0.75

Source code in src/skcomplex/models/physics/climate/_daisyworld.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    white: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
    black: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
    L: Annotated[float, ParameterSpec()] = 1.0,
    S: Annotated[float, ParameterSpec()] = 917.0,
    sigma: Annotated[float, ParameterSpec()] = 5.67e-08,
    q: Annotated[float, ParameterSpec()] = 2060000000.0,
    p: Annotated[float, ParameterSpec()] = 1.0,
    gamma: Annotated[float, ParameterSpec()] = 0.3,
    k: Annotated[float, ParameterSpec()] = 0.003265,
    optimum: Annotated[float, ParameterSpec()] = 295.5,
    albedo_ground: Annotated[float, ParameterSpec()] = 0.5,
    albedo_white: Annotated[float, ParameterSpec()] = 0.75,
    albedo_black: Annotated[float, ParameterSpec()] = 0.25,
) -> None:
    super().__init__(label, visual=visual)
    self.white = white
    self.black = black
    self.L = L
    self.S = S
    self.sigma = sigma
    self.q = q
    self.p = p
    self.gamma = gamma
    self.k = k
    self.optimum = optimum
    self.albedo_ground = albedo_ground
    self.albedo_white = albedo_white
    self.albedo_black = albedo_black

DelayedOscillatorElement(label, visual=None, *, T=0.1, alpha=0.75, delta=6.0)

Bases: Element

Eastern Pacific of the delayed oscillator of El Niño, its sea surface temperature anomaly.

The anomaly follows dT/dt = T - T^3 - alpha T(t - delta). It grows by local feedback and is damped by the delayed return of oceanic waves reflected at the western boundary.

Parameters:

Name	Type	Description	Default
label	str	The label of the ocean.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
T	Annotated[float, StateSpec()]	The initial temperature anomaly, also its history.	0.1
alpha	Annotated[float, ParameterSpec()]	The strength of the delayed negative feedback.	0.75
delta	Annotated[float, ParameterSpec()]	The delay of the wave transit.	6.0

Examples:

>>> from skcomplex.models.physics.climate import DelayedOscillatorElement
>>> DelayedOscillatorElement('pacific').delta
6.0

Source code in src/skcomplex/models/physics/climate/_delayed_oscillator.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    T: Annotated[float, StateSpec()] = 0.1,
    alpha: Annotated[float, ParameterSpec()] = 0.75,
    delta: Annotated[float, ParameterSpec()] = 6.0,
) -> None:
    super().__init__(label, visual=visual)
    self.T = T
    self.alpha = alpha
    self.delta = delta

EarlyWarningAnalysis(variables=None, *, window=None, detrend=True)

Bases: RollingAnalysis

Rolling indicators of critical slowing down before a tipping point, for every series of a source.

It is the rolling analysis of kind 'autocorrelation' with detrended windows, whose statistic and trend are the lag-one autocorrelation, and it adds the other indicators of Dakos et al. (2012): the rolling variance and skewness, and the reddening, the ratio of the mean power of the lowest eighth of frequencies to that of the highest half, which rises as fluctuations slow down. The trend of every indicator is Kendall's tau against time.

Parameters:

Name	Type	Description	Default
variables	Sequence[str] | None	The variables of a trajectory to read, or None for every scalar series.	None
window	int | None	The number of values in every rolling window.	None
detrend	bool	Whether to remove a linear trend in every window first.	True

Attributes:

Name	Type	Description
variance_		The variance of every series in every window, nan before the first full window.
autocorrelation_		The lag-1 autocorrelation of every series in every window, nan before the first full window.
skewness_		The skewness of every series in every window, nan before the first full window.
reddening_		The reddening of every series in every window, nan before the first full window.
variance_trend_		Kendall's tau of the variance of every series against time, positive when it rises.
autocorrelation_trend_		Kendall's tau of the autocorrelation of every series against time, positive when it rises.
skewness_trend_		Kendall's tau of the skewness of every series against time.
reddening_trend_		Kendall's tau of the reddening of every series against time.

Examples:

>>> import numpy as np
>>> from skcomplex.models.physics.climate import EarlyWarningAnalysis
>>> rng = np.random.default_rng(0)
>>> x = np.zeros(2000)
>>> for t in range(1, 2000):
...     x[t] = (0.2 + 0.79 * t / 2000) * x[t - 1] + rng.standard_normal()
>>> bool(EarlyWarningAnalysis(window=200).check().prepare(x).run().autocorrelation_trend_[0] > 0.5)
True

Source code in src/skcomplex/models/physics/climate/_early_warning.py

def __init__(
    self,
    variables: Sequence[str] | None = None,
    *,
    window: int | None = None,
    detrend: bool = True,
) -> None:
    super().__init__('autocorrelation', variables, window=window, detrend=detrend)

EnergySectorElement(label, visual=None, *, reserve=1000.0, capacity=0.1, population=1.0, reserve_size=1000.0, initial_eroi=20.0, decline=3.0, reinvestment=0.1, lifetime=20.0, birth_rate=0.03, death_rate=0.015)

Bases: Element

Energy sector of a society, a finite reserve, the capacity that extracts it, and the population it supplies.

The capacity draws the gross energy capacity reserve / reserve_size from the reserve. Part of it pays for the extraction work, the gross energy divided by the energy return on investment 1 + (initial_eroi - 1) (reserve / reserve_size)^decline, and the rest is the net energy society receives. A share reinvestment of the net energy builds capacity, which depreciates over lifetime, and the population grows at the rate birth_rate - death_rate.

Parameters:

Name	Type	Description	Default
label	str	The label of the sector.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
reserve	Annotated[float, StateSpec(domain=closed_open(0, oo), quantity=energy)]	The initial energy reserve.	1000.0
capacity	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial extraction capacity, the gross energy per year from a full reserve.	0.1
population	Annotated[float, StateSpec(domain=closed_open(0, oo), quantity=people)]	The initial population.	1.0
reserve_size	Annotated[float, ParameterSpec(domain=open(0, oo))]	The full size of the reserve.	1000.0
initial_eroi	Annotated[float, ParameterSpec(domain=open(0, oo))]	The energy return on investment of a full reserve.	20.0
decline	Annotated[float, ParameterSpec(domain=open(0, oo))]	The power of the remaining fraction of the reserve in the energy return on investment.	3.0
reinvestment	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The share of the net energy reinvested in extraction capacity.	0.1
lifetime	Annotated[float, ParameterSpec(domain=open(0, oo))]	The lifetime of the extraction capacity in years.	20.0
birth_rate	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The relative birth rate of the population per year.	0.03
death_rate	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The relative death rate of the population per year.	0.015

Examples:

>>> from skcomplex.models.physics.climate import EnergySectorElement
>>> EnergySectorElement('energy').initial_eroi
20.0

Source code in src/skcomplex/models/physics/climate/_net_energy.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    reserve: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo), quantity='energy')] = 1000.0,
    capacity: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    population: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo), quantity='people')] = 1.0,
    reserve_size: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1000.0,
    initial_eroi: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 20.0,
    decline: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 3.0,
    reinvestment: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    lifetime: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 20.0,
    birth_rate: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.03,
    death_rate: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.015,
) -> None:
    super().__init__(label, visual=visual)
    self.reserve = reserve
    self.capacity = capacity
    self.population = population
    self.reserve_size = reserve_size
    self.initial_eroi = initial_eroi
    self.decline = decline
    self.reinvestment = reinvestment
    self.lifetime = lifetime
    self.birth_rate = birth_rate
    self.death_rate = death_rate

GeodynamoElement(label, visual=None, *, Q=0.1, D=0.1, V=0.1, mu=0.119, nu=0.1, Gamma=0.9)

Bases: Element

Dynamo of Gissinger's geodynamo, its quadrupolar and dipolar fields and the velocity of the flow.

The quadrupolar field Q, the dipolar field D, and the velocity V follow dQ/dt = mu Q - V D, dD/dt = -nu D + V Q, and dV/dt = Gamma - V + Q D.

Parameters:

Name	Type	Description	Default
label	str	The label of the dynamo.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
Q	Annotated[float, StateSpec()]	The initial quadrupolar field.	0.1
D	Annotated[float, StateSpec()]	The initial dipolar field.	0.1
V	Annotated[float, StateSpec()]	The initial velocity.	0.1
mu	Annotated[float, ParameterSpec()]	The growth rate of the quadrupole.	0.119
nu	Annotated[float, ParameterSpec()]	The damping rate of the dipole.	0.1
Gamma	Annotated[float, ParameterSpec()]	The forcing of the flow.	0.9

Examples:

>>> from skcomplex.models.physics.climate import GeodynamoElement
>>> GeodynamoElement('dynamo').Gamma
0.9

Source code in src/skcomplex/models/physics/climate/_geodynamo.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    Q: Annotated[float, StateSpec()] = 0.1,
    D: Annotated[float, StateSpec()] = 0.1,
    V: Annotated[float, StateSpec()] = 0.1,
    mu: Annotated[float, ParameterSpec()] = 0.119,
    nu: Annotated[float, ParameterSpec()] = 0.1,
    Gamma: Annotated[float, ParameterSpec()] = 0.9,
) -> None:
    super().__init__(label, visual=visual)
    self.Q = Q
    self.D = D
    self.V = V
    self.mu = mu
    self.nu = nu
    self.Gamma = Gamma

IceAlbedoClimateElement(label, visual=None, *, T=290.0, s=1.0, epsilon=0.65, sigma=1.6e-10, albedo=0.5, contrast=0.2, T_ice=263.0, width=4.0)

Bases: Element

Global climate of a zero-dimensional energy balance, bistable through the ice-albedo feedback.

In units of the mean insolation, the global temperature follows dT/dt = s (1 - a(T)) - sigma epsilon T^4, with the albedo a(T) = albedo - contrast tanh((T - T_ice) / width), high on an icy planet and low on a warm one.

Parameters:

Name	Type	Description	Default
label	str	The label of the climate.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
T	Annotated[float, StateSpec(domain=open(0, oo))]	The initial global mean temperature in kelvin.	290.0
s	Annotated[float, ParameterSpec()]	The insolation relative to the present one.	1.0
epsilon	Annotated[float, ParameterSpec()]	The effective emissivity, lowered by the greenhouse effect.	0.65
sigma	Annotated[float, ParameterSpec()]	The Stefan-Boltzmann constant divided by the mean insolation.	1.6e-10
albedo	Annotated[float, ParameterSpec()]	The albedo at the middle of the transition from ice to no ice.	0.5
contrast	Annotated[float, ParameterSpec()]	Half the difference between the albedos of a frozen and of a warm planet.	0.2
T_ice	Annotated[float, ParameterSpec()]	The temperature of the middle of the transition.	263.0
width	Annotated[float, ParameterSpec()]	The width of the transition in kelvin.	4.0

Examples:

>>> from skcomplex.models.physics.climate import IceAlbedoClimateElement
>>> IceAlbedoClimateElement('earth').T_ice
263.0

Source code in src/skcomplex/models/physics/climate/_ice_albedo_climate.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    T: Annotated[float, StateSpec(domain=RealSet.open(0, oo))] = 290.0,
    s: Annotated[float, ParameterSpec()] = 1.0,
    epsilon: Annotated[float, ParameterSpec()] = 0.65,
    sigma: Annotated[float, ParameterSpec()] = 1.6e-10,
    albedo: Annotated[float, ParameterSpec()] = 0.5,
    contrast: Annotated[float, ParameterSpec()] = 0.2,
    T_ice: Annotated[float, ParameterSpec()] = 263.0,
    width: Annotated[float, ParameterSpec()] = 4.0,
) -> None:
    super().__init__(label, visual=visual)
    self.T = T
    self.s = s
    self.epsilon = epsilon
    self.sigma = sigma
    self.albedo = albedo
    self.contrast = contrast
    self.T_ice = T_ice
    self.width = width

LatitudeBandElement(label, visual=None, *, T=15.0, x=0.5, Q=334.4, s2=-0.482, A=203.3, B=2.09, D=0.649, C=10.0, albedo_ice=0.62, albedo_free=0.319, albedo_slope=0.202, T_ice=-10.0, width=1.0)

Bases: Element

Latitude band of the Budyko-Sellers energy balance, heated by the sun and exchanging heat with its neighbours.

The temperature of the band at x = sin(latitude) follows C dT/dt = Q s(x) (1 - a(x, T)) - (A + B T) + D d/dx((1 - x^2) dT/dx), with the annual insolation s(x) = 1 + s2 P2(x), where P2(x) = (3x^2 - 1) / 2. The albedo is albedo_free + albedo_slope P2(x) without ice and rises to albedo_ice where the temperature falls below T_ice.

Parameters:

Name	Type	Description	Default
label	str	The label of the band or the field.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
T	Annotated[float, StateSpec()]	The initial temperature in °C.	15.0
x	Annotated[float, ParameterSpec()]	The sine of the latitude of the band.	0.5
Q	Annotated[float, ParameterSpec()]	The mean insolation in W m⁻².	334.4
s2	Annotated[float, ParameterSpec()]	The coefficient of the second Legendre polynomial of the insolation.	-0.482
A	Annotated[float, ParameterSpec()]	The outgoing radiation at 0 °C in W m⁻².	203.3
B	Annotated[float, ParameterSpec()]	The increase of outgoing radiation per degree in W m⁻² K⁻¹.	2.09
D	Annotated[float, ParameterSpec()]	The coefficient of meridional heat transport in W m⁻² K⁻¹.	0.649
C	Annotated[float, ParameterSpec()]	The heat capacity in W yr m⁻² K⁻¹, so time is in years.	10.0
albedo_ice	Annotated[float, ParameterSpec()]	The albedo of an icy band.	0.62
albedo_free	Annotated[float, ParameterSpec()]	The mean albedo of an ice free band.	0.319
albedo_slope	Annotated[float, ParameterSpec()]	The coefficient of the second Legendre polynomial of the albedo of an ice free band.	0.202
T_ice	Annotated[float, ParameterSpec()]	The temperature in °C below which the band is icy.	-10.0
width	Annotated[float, ParameterSpec()]	The width of the transition to ice in kelvin.	1.0

Examples:

>>> from skcomplex.models.physics.climate import LatitudeBandElement
>>> LatitudeBandElement('band').albedo_ice
0.62

Source code in src/skcomplex/models/physics/climate/_latitude_band.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    T: Annotated[float, StateSpec()] = 15.0,
    x: Annotated[float, ParameterSpec()] = 0.5,
    Q: Annotated[float, ParameterSpec()] = 334.4,
    s2: Annotated[float, ParameterSpec()] = -0.482,
    A: Annotated[float, ParameterSpec()] = 203.3,
    B: Annotated[float, ParameterSpec()] = 2.09,
    D: Annotated[float, ParameterSpec()] = 0.649,
    C: Annotated[float, ParameterSpec()] = 10.0,
    albedo_ice: Annotated[float, ParameterSpec()] = 0.62,
    albedo_free: Annotated[float, ParameterSpec()] = 0.319,
    albedo_slope: Annotated[float, ParameterSpec()] = 0.202,
    T_ice: Annotated[float, ParameterSpec()] = -10.0,
    width: Annotated[float, ParameterSpec()] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.T = T
    self.x = x
    self.Q = Q
    self.s2 = s2
    self.A = A
    self.B = B
    self.D = D
    self.C = C
    self.albedo_ice = albedo_ice
    self.albedo_free = albedo_free
    self.albedo_slope = albedo_slope
    self.T_ice = T_ice
    self.width = width

Lorenz84Element(label, visual=None, *, x=1.0, y=1.0, z=1.0, a=0.25, b=4.0, F=8.0, G=1.0)

Bases: Element

Atmosphere of Lorenz-84, the general circulation reduced to a westerly jet and a travelling wave of eddies.

The jet x and the cosine and sine phases y and z of the eddies follow dx/dt = -y^2 - z^2 - a x + a F, dy/dt = x y - b x z - y + G, and dz/dt = b x y + x z - z.

Parameters:

Name	Type	Description	Default
label	str	The label of the atmosphere.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec()]	The initial strength of the westerly jet.	1.0
y	Annotated[float, StateSpec()]	The initial cosine phase of the eddies.	1.0
z	Annotated[float, StateSpec()]	The initial sine phase of the eddies.	1.0
a	Annotated[float, ParameterSpec()]	The damping of the jet relative to the eddies.	0.25
b	Annotated[float, ParameterSpec()]	The displacement of the eddies by the jet.	4.0
F	Annotated[float, ParameterSpec()]	The forcing of the jet by the equator-to-pole temperature contrast.	8.0
G	Annotated[float, ParameterSpec()]	The forcing of the eddies by the land-sea contrast.	1.0

Examples:

>>> from skcomplex.models.physics.climate import Lorenz84Element
>>> Lorenz84Element('atmosphere').F
8.0

Source code in src/skcomplex/models/physics/climate/_lorenz84.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = 1.0,
    y: Annotated[float, StateSpec()] = 1.0,
    z: Annotated[float, StateSpec()] = 1.0,
    a: Annotated[float, ParameterSpec()] = 0.25,
    b: Annotated[float, ParameterSpec()] = 4.0,
    F: Annotated[float, ParameterSpec()] = 8.0,
    G: Annotated[float, ParameterSpec()] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.y = y
    self.z = z
    self.a = a
    self.b = b
    self.F = F
    self.G = G

Lorenz96Element(label, visual=None, *, x=8.0, F=8.0)

Bases: Element

Ring of Lorenz-96 sites around a latitude circle, advected by their neighbours, damped, and forced.

Every site follows dx_i/dt = (x_{i+1} - x_{i-2}) x_{i-1} - x_i + F, with the indices taken around the ring.

Parameters:

Name	Type	Description	Default
label	str	The label of the ring.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[ArrayLike, StateSpec(vector=True)]	The initial values of all sites.	8.0
F	Annotated[float, ParameterSpec()]	The constant forcing.	8.0

Examples:

>>> from skcomplex.models.physics.climate import Lorenz96Element
>>> Lorenz96Element('ring').F
8.0

Source code in src/skcomplex/models/physics/climate/_lorenz96.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True)] = 8.0,
    F: Annotated[float, ParameterSpec()] = 8.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.F = F

RampedFoldElement(label, visual=None, *, x=-1.3247, a=-1.0, b=1.0, sigma=0.02, phi0=-1.0, phi1=1.0, rate=0.004)

Bases: Element

Stochastic fold dx = (a x^3 + b x + phi(t)) dt + sigma dW whose control rises linearly from phi0 to phi1.

For a < 0 < b the fold has two stable branches for |phi| < phi_c = sqrt(-4 b^3 / (27 a)), and the lower branch vanishes at phi = phi_c. The control rises at the rate rate and stays at phi1 once it reaches it.

Parameters:

Name	Type	Description	Default
label	str	The label of the system.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec()]	The initial state.	-1.3247
a	Annotated[float, ParameterSpec()]	The cubic coefficient, negative.	-1.0
b	Annotated[float, ParameterSpec()]	The linear coefficient, positive.	1.0
sigma	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The amplitude of the noise.	0.02
phi0	float	The initial control.	-1.0
phi1	float	The final control.	1.0
rate	float	The rate at which the control rises.	0.004

Examples:

>>> from skcomplex.models.physics.climate import RampedFoldElement
>>> RampedFoldElement('system').rate
0.004

Source code in src/skcomplex/models/physics/climate/_ramped_fold.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = -1.3247,
    a: Annotated[float, ParameterSpec()] = -1.0,
    b: Annotated[float, ParameterSpec()] = 1.0,
    sigma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.02,
    phi0: float = -1.0,
    phi1: float = 1.0,
    rate: float = 0.004,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.a = a
    self.b = b
    self.sigma = sigma
    self.phi0 = phi0
    self.phi1 = phi1
    self.rate = rate

SaltzmanMaaschElement(label, visual=None, *, x=-0.2, y=1.0, z=0.0, p=1.0, q=2.5, r=1.3, s=0.6, u=0.2, nu=0.0, t0=10.0)

Bases: Element

Climate of the Saltzman-Maasch glacial cycles, the anomalies of ice, carbon dioxide, and deep ocean temperature.

The anomalies of global ice volume x, atmospheric carbon dioxide y, and deep ocean temperature z follow dx/dt = -x - y - u z - nu F(t), dy/dt = -p z + r y - s y^2 - y^3, and dz/dt = -q (x + z), where F(t) = cos(2π t0 t / 22.1) - cos(2π t0 t / 44.1) is the orbital forcing of precession and obliquity.

Parameters:

Name	Type	Description	Default
label	str	The label of the climate.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec()]	The initial ice volume anomaly.	-0.2
y	Annotated[float, StateSpec()]	The initial carbon dioxide anomaly.	1.0
z	Annotated[float, StateSpec()]	The initial deep ocean temperature anomaly.	0.0
p	Annotated[float, ParameterSpec()]	The effect of the deep ocean temperature on carbon dioxide.	1.0
q	Annotated[float, ParameterSpec()]	The relaxation rate of the deep ocean temperature to the ice volume.	2.5
r	Annotated[float, ParameterSpec()]	The growth rate of carbon dioxide by its own feedback.	1.3
s	Annotated[float, ParameterSpec()]	The quadratic saturation of carbon dioxide.	0.6
u	Annotated[float, ParameterSpec()]	The effect of the deep ocean temperature on ice volume.	0.2
nu	Annotated[float, ParameterSpec()]	The amplitude of the orbital forcing.	0.0
t0	Annotated[float, ParameterSpec()]	The time unit in thousands of years, which scales the orbital periods.	10.0

Examples:

>>> from skcomplex.models.physics.climate import SaltzmanMaaschElement
>>> SaltzmanMaaschElement('climate').q
2.5

Source code in src/skcomplex/models/physics/climate/_saltzman_maasch.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = -0.2,
    y: Annotated[float, StateSpec()] = 1.0,
    z: Annotated[float, StateSpec()] = 0.0,
    p: Annotated[float, ParameterSpec()] = 1.0,
    q: Annotated[float, ParameterSpec()] = 2.5,
    r: Annotated[float, ParameterSpec()] = 1.3,
    s: Annotated[float, ParameterSpec()] = 0.6,
    u: Annotated[float, ParameterSpec()] = 0.2,
    nu: Annotated[float, ParameterSpec()] = 0.0,
    t0: Annotated[float, ParameterSpec()] = 10.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.y = y
    self.z = z
    self.p = p
    self.q = q
    self.r = r
    self.s = s
    self.u = u
    self.nu = nu
    self.t0 = t0

SoilTemperatureElement(label, visual=None, *, x=0.3, y=0.1, r=3.8, epsilon=0.1)

Bases: Element

Soil of the coupled temperature maps, its dimensionless surface temperature and deeper soil temperature.

Two maps of soil_temperature iterate the states. The surface temperature follows the energy balance of the ground, reduced to a logistic map, and the deeper soil temperature relaxes towards it: x_{n+1} = r x_n (1 - x_n) + epsilon y_n and y_{n+1} = epsilon (x_n + y_n).

Parameters:

Name	Type	Description	Default
label	str	The label of the soil.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec()]	The initial dimensionless surface temperature.	0.3
y	Annotated[float, StateSpec()]	The initial dimensionless deeper soil temperature.	0.1
r	Annotated[float, ParameterSpec()]	The logistic parameter of the surface energy balance.	3.8
epsilon	Annotated[float, ParameterSpec()]	The coupling between the surface and the deeper soil.	0.1

Examples:

>>> from skcomplex.models.physics.climate import SoilTemperatureElement
>>> SoilTemperatureElement('soil').r
3.8

Source code in src/skcomplex/models/physics/climate/_soil.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = 0.3,
    y: Annotated[float, StateSpec()] = 0.1,
    r: Annotated[float, ParameterSpec()] = 3.8,
    epsilon: Annotated[float, ParameterSpec()] = 0.1,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.y = y
    self.r = r
    self.epsilon = epsilon

StochasticEnergyBalanceElement(label, visual=None, *, T=290.0, mu=1.0, lam=0.0, drift=0.0, rho=0.0, b_ini=1.69e-05, b_fin=1.69e-05, a_ini=0.6927, nu=1.0, c0=C0, c1=C1)

Bases: Element

Global mean surface temperature of the Sutera-Fraedrich energy balance, with drifting and ramped controls.

The temperature follows dT = c0 [-T^4 + c1 mu (b2 T^2 - a2)] dt + sqrt(nu) dW, in kelvin and years. The control mu drifts at the rate -drift, and b2 = b_ini + lam (b_fin - b_ini) follows the logistic ramp dlam/dt = rho lam (1 - lam). The control a2 = a_ini + c1 mu (b2^2 - b_ini^2) / 4 moves with b2 so that the discriminant (c1 mu b2)^2 - 4 c1 mu a2 of the equilibria stays constant during the ramp.

Parameters:

Name	Type	Description	Default
label	str	The label of the climate.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
T	Annotated[float, StateSpec()]	The initial temperature in kelvin.	290.0
mu	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial relative solar input.	1.0
lam	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial progress of the ramp, between 0 and 1.	0.0
drift	Annotated[float, ParameterSpec()]	The rate at which mu decreases, per year.	0.0
rho	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rate of the logistic ramp of b2, per year.	0.0
b_ini	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The initial value of b2, per square kelvin.	1.69e-05
b_fin	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The final value of b2, per square kelvin.	1.69e-05
a_ini	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The initial value of a2.	0.6927
nu	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The variance of the noise per year, in square kelvin.	1.0
c0	Annotated[float, ParameterSpec(domain=open(0, oo))]	The scale of the energy balance, per cubic kelvin and year.	C0
c1	Annotated[float, ParameterSpec(domain=open(0, oo))]	The scale of the absorbed radiation, in kelvin to the fourth.	C1

Examples:

>>> from skcomplex.models.physics.climate import StochasticEnergyBalanceElement
>>> StochasticEnergyBalanceElement('earth').nu
1.0

Source code in src/skcomplex/models/physics/climate/_stochastic_energy_balance.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    T: Annotated[float, StateSpec()] = 290.0,
    mu: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    lam: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    drift: Annotated[float, ParameterSpec()] = 0.0,
    rho: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    b_ini: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.69e-05,
    b_fin: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.69e-05,
    a_ini: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.6927,
    nu: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    c0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = C0,
    c1: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = C1,
) -> None:
    super().__init__(label, visual=visual)
    self.T = T
    self.mu = mu
    self.lam = lam
    self.drift = drift
    self.rho = rho
    self.b_ini = b_ini
    self.b_fin = b_fin
    self.a_ini = a_ini
    self.nu = nu
    self.c0 = c0
    self.c1 = c1

TippingElementElement(label, visual=None, *, x=-1.0, c=0.0, d=0.2)

Bases: Element

Tipping element in the normal form of a fold, pushed toward its upper state by tipped neighbours.

The state follows dx_i/dt = -x_i^3 + x_i + sqrt(4/27) c_i + d Σ_j A_ij (x_j + 1). Alone, the element has a lower and an upper stable state for |c| < 1 and only the upper one beyond the critical value c = 1. A neighbour in its lower state x_j ≈ -1 adds almost nothing, while a tipped neighbour at x_j ≈ 1 adds about 2 d. The element carries the fold, and the coupling of tipping_cascade_network adds the pushes of the neighbours.

Parameters:

Name	Type	Description	Default
label	str	The label of the element or the field.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec()]	The initial state, near -1 for an element in its lower state.	-1.0
c	Annotated[float, ParameterSpec()]	The control parameter, critical at 1.	0.0
d	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The coupling strength of every link.	0.2

Examples:

>>> from skcomplex.models.physics.climate import TippingElementElement
>>> TippingElementElement('ice sheet').x
-1.0

Source code in src/skcomplex/models/physics/climate/_tipping.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = -1.0,
    c: Annotated[float, ParameterSpec()] = 0.0,
    d: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.c = c
    self.d = d

VanDerPolIceAgesElement(label, visual=None, *, x=0.0, y=1.0, alpha=30.0, beta=0.75, gamma=0.0, tau=35.09)

Bases: Element

Climate of glacial cycles as a Van der Pol relaxation oscillator, a slow ice volume and a fast climate state.

With time in thousands of years, tau x' = -(y + beta + gamma F(t)) and tau y' = -alpha (y^3/3 - y - x), where F(t) = cos(2 pi t/23) + cos(2 pi t/41) stands for precession and obliquity.

Parameters:

Name	Type	Description	Default
label	str	The label of the climate.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec()]	The initial ice volume.	0.0
y	Annotated[float, StateSpec()]	The initial climate state.	1.0
alpha	Annotated[float, ParameterSpec()]	The time-scale separation of the fast climate state.	30.0
beta	Annotated[float, ParameterSpec()]	The asymmetry of the cycle, between build-up and termination.	0.75
gamma	Annotated[float, ParameterSpec()]	The strength of the orbital forcing.	0.0
tau	Annotated[float, ParameterSpec()]	The time scale, in thousands of years.	35.09

Examples:

>>> from skcomplex.models.physics.climate import VanDerPolIceAgesElement
>>> VanDerPolIceAgesElement('climate').alpha
30.0

Source code in src/skcomplex/models/physics/climate/_van_der_pol_ice_ages.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = 0.0,
    y: Annotated[float, StateSpec()] = 1.0,
    alpha: Annotated[float, ParameterSpec()] = 30.0,
    beta: Annotated[float, ParameterSpec()] = 0.75,
    gamma: Annotated[float, ParameterSpec()] = 0.0,
    tau: Annotated[float, ParameterSpec()] = 35.09,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.y = y
    self.alpha = alpha
    self.beta = beta
    self.gamma = gamma
    self.tau = tau

budyko_sellers(n_bands=30, T=15.0, *, Q=334.4, A=203.3, B=2.09, D=0.649, C=10.0)

Return the Budyko-Sellers energy balance of one hemisphere, resolved in the sine of latitude.

The temperature of a band at x = sin(latitude) follows C dT/dt = Q s(x) (1 - a(x, T)) - (A + B T) + D d/dx((1 - x^2) dT/dx), with the annual insolation s(x) = 1 - 0.482 P2(x), where P2(x) = (3x^2 - 1) / 2. The albedo is 0.319 + 0.202 P2(x) without ice and rises to 0.62 where the temperature falls below -10 °C. Bands are centred at x = (i + 1/2) / n_bands, with no heat flux across the equator or the pole. Present insolation keeps both a climate with polar ice caps and a snowball Earth.

Parameters:

Name	Type	Description	Default
n_bands	int	The number of latitude bands.	30
T	float	The initial temperature of every band in °C.	15.0
Q	float	The mean insolation in W m⁻².	334.4
A	float	The outgoing radiation at 0 °C in W m⁻².	203.3
B	float	The increase of outgoing radiation per degree in W m⁻² K⁻¹.	2.09
D	float	The coefficient of meridional heat transport in W m⁻² K⁻¹.	0.649
C	float	The heat capacity in W yr m⁻² K⁻¹, so time is in years.	10.0

Returns:

Name	Type	Description
system	System	The model, with the temperatures bands.T from the equator to the pole.

Examples:

>>> from skcomplex.models.physics.climate import budyko_sellers
>>> budyko_sellers(10).elements[0].check().prepare().size_
10

Source code in src/skcomplex/models/physics/climate/_energy_balance.py

@register_model(
    field='physics/climate',
    tags=('field', 'pde', 'energy balance', 'ice-albedo feedback', 'snowball'),
    references=('Budyko 1969', 'Sellers 1969', 'North 1975', 'Ghil and Lucarini 2020'),
    defaults={'until': 100.0, 'dt': 2.0},
)
def budyko_sellers(
    n_bands: int = 30,
    T: float = 15.0,
    *,
    Q: float = 334.4,
    A: float = 203.3,
    B: float = 2.09,
    D: float = 0.649,
    C: float = 10.0,
) -> System:
    """Return the Budyko-Sellers energy balance of one hemisphere, resolved in the sine of latitude.

    The temperature of a band at `x = sin(latitude)` follows `C dT/dt = Q s(x) (1 - a(x, T)) - (A + B T) +
    D d/dx((1 - x^2) dT/dx)`, with the annual insolation `s(x) = 1 - 0.482 P2(x)`, where `P2(x) = (3x^2 - 1) / 2`.
    The albedo is `0.319 + 0.202 P2(x)` without ice and rises to 0.62 where the temperature falls below -10 °C.
    Bands are centred at `x = (i + 1/2) / n_bands`, with no heat flux across the equator or the pole. Present
    insolation keeps both a climate with polar ice caps and a snowball Earth.

    Args:
        n_bands:
            The number of latitude bands.

        T:
            The initial temperature of every band in °C.

        Q:
            The mean insolation in W m⁻².

        A:
            The outgoing radiation at 0 °C in W m⁻².

        B:
            The increase of outgoing radiation per degree in W m⁻² K⁻¹.

        D:
            The coefficient of meridional heat transport in W m⁻² K⁻¹.

        C:
            The heat capacity in W yr m⁻² K⁻¹, so time is in years.

    Returns:
        system:
            The model, with the temperatures `bands.T` from the equator to the pole.

    Examples:
        >>> from skcomplex.models.physics.climate import budyko_sellers
        >>> budyko_sellers(10).elements[0].check().prepare().size_
        10
    """
    x = (np.arange(n_bands) + 0.5) / n_bands
    bands = FieldElement(
        'bands',
        element=LatitudeBandElement,
        space=LatticeSpace(n_bands, spacing=1 / n_bands),
        boundary={'T': 'zero-flux'},
        values={'T': np.full(n_bands, float(T)), 'x': x, 'Q': Q, 'A': A, 'B': B, 'D': D, 'C': C},
    )
    return System(elements=[bands])

carbon_cycle(c=83.58, w=2050.0, *, c_x=55.0, nu=0.0)

Return Rothman's excitable carbon cycle, where a large enough perturbation of the ocean sets off an excursion.

The carbonate ion concentration c and the dissolved inorganic carbon w, in μmol kg⁻¹, follow dc/dt = f(c) (mu (1 - b s(c, c_p) - theta (1 - s(c, c_x)) - nu) + w - w0) and dw/dt = mu (1 - b s(c, c_p) + theta (1 - s(c, c_x)) + nu) - w + w0, with the sigmoids s(c, c*) = c^gamma / (c^gamma + c*^gamma) of burial and respiration and the buffer f(c) = f0 c^beta / (c^beta + c_f^beta). The steady state near c = 83.6 and w = 2197 is stable, and from a dissolved carbon far enough below it the ocean returns only after a large excursion of c.

Parameters:

Name	Type	Description	Default
c	float	The initial carbonate ion concentration.	83.58
w	float	The initial dissolved inorganic carbon.	2050.0
c_x	float	The crossover carbonate concentration of respiration, which sets the excitability.	55.0
nu	float	The rate of external carbon injection.	0.0

Returns:

Name	Type	Description
system	System	The model, with the concentrations ocean.c and ocean.w.

Examples:

>>> from skcomplex.models.physics.climate import carbon_cycle
>>> carbon_cycle(nu=0.1).elements[0].nu
0.1

Source code in src/skcomplex/models/physics/climate/_ice_ages.py

@register_model(
    field='physics/climate',
    tags=('ode', 'excitability', 'carbon cycle', 'mass extinctions'),
    references=('Rothman 2019', 'Datseris and Parlitz 2022'),
    defaults={'until': 20.0, 'dt': 0.05},
)
def carbon_cycle(
    c: float = 83.58,
    w: float = 2050.0,
    *,
    c_x: float = 55.0,
    nu: float = 0.0,
) -> System:
    """Return Rothman's excitable carbon cycle, where a large enough perturbation of the ocean sets off an excursion.

    The carbonate ion concentration `c` and the dissolved inorganic carbon `w`, in μmol kg⁻¹, follow
    `dc/dt = f(c) (mu (1 - b s(c, c_p) - theta (1 - s(c, c_x)) - nu) + w - w0)` and
    `dw/dt = mu (1 - b s(c, c_p) + theta (1 - s(c, c_x)) + nu) - w + w0`, with the sigmoids
    `s(c, c*) = c^gamma / (c^gamma + c*^gamma)` of burial and respiration and the buffer
    `f(c) = f0 c^beta / (c^beta + c_f^beta)`. The steady state near `c = 83.6` and `w = 2197` is stable, and from a
    dissolved carbon far enough below it the ocean returns only after a large excursion of `c`.

    Args:
        c:
            The initial carbonate ion concentration.

        w:
            The initial dissolved inorganic carbon.

        c_x:
            The crossover carbonate concentration of respiration, which sets the excitability.

        nu:
            The rate of external carbon injection.

    Returns:
        system:
            The model, with the concentrations `ocean.c` and `ocean.w`.

    Examples:
        >>> from skcomplex.models.physics.climate import carbon_cycle
        >>> carbon_cycle(nu=0.1).elements[0].nu
        0.1
    """
    return System(elements=[CarbonCycleElement('ocean', c=c, w=w, c_x=c_x, nu=nu)])

daisyworld(white=0.2, black=0.2, *, L=1.0, gamma=0.3)

Return Daisyworld, where white and black daisies regulate the temperature of their planet through its albedo.

The areas covered by white and black daisies grow on the bare ground at a rate parabolic in their local temperature, with an optimum at 22.5 °C, and die at the rate gamma. The effective temperature follows the planetary albedo, sigma T^4 = S L (1 - A), and is recorded as the observable temperature.

Parameters:

Name	Type	Description	Default
white	float	The initial area fraction of white daisies.	0.2
black	float	The initial area fraction of black daisies.	0.2
L	float	The luminosity of the sun relative to its present value.	1.0
gamma	float	The death rate of daisies.	0.3

Returns:

Name	Type	Description
system	System	The model, with the areas planet.white and planet.black and the observable temperature.

Examples:

>>> from skcomplex.models.physics.climate import daisyworld
>>> daisyworld(L=1.2).elements[0].L
1.2

Source code in src/skcomplex/models/physics/climate/_energy_balance.py

@register_model(
    field='physics/climate',
    tags=('ode', 'regulation', 'Gaia', 'vegetation'),
    references=('Watson and Lovelock 1983', 'Jensen 2023'),
    defaults={'until': 100.0, 'dt': 0.5},
    observables=lambda system: [ObservableSpec('temperature', _daisyworld_temperature)],
)
def daisyworld(
    white: float = 0.2,
    black: float = 0.2,
    *,
    L: float = 1.0,
    gamma: float = 0.3,
) -> System:
    """Return Daisyworld, where white and black daisies regulate the temperature of their planet through its albedo.

    The areas covered by white and black daisies grow on the bare ground at a rate parabolic in their local
    temperature, with an optimum at 22.5 °C, and die at the rate `gamma`. The effective temperature follows the
    planetary albedo, `sigma T^4 = S L (1 - A)`, and is recorded as the observable `temperature`.

    Args:
        white:
            The initial area fraction of white daisies.

        black:
            The initial area fraction of black daisies.

        L:
            The luminosity of the sun relative to its present value.

        gamma:
            The death rate of daisies.

    Returns:
        system:
            The model, with the areas `planet.white` and `planet.black` and the observable `temperature`.

    Examples:
        >>> from skcomplex.models.physics.climate import daisyworld
        >>> daisyworld(L=1.2).elements[0].L
        1.2
    """
    return System(
        elements=[DaisyworldElement('planet', white=white, black=black, L=L, gamma=gamma)],
    )

energy_balance(T=290.0, *, s=1.0, epsilon=0.65, sigma=1.6e-10, albedo=0.5, contrast=0.2, T_ice=263.0, width=4.0)

Return the zero-dimensional energy balance of the Earth, bistable through the ice-albedo feedback.

In units of the mean insolation, the global temperature follows dT/dt = s (1 - a(T)) - sigma epsilon T^4, with the albedo a(T) = albedo - contrast tanh((T - T_ice) / width), high on an icy planet and low on a warm one. For the default parameters a warm climate near 286 K and a frozen climate near 232 K are both stable.

Parameters:

Name	Type	Description	Default
T	float	The initial global mean temperature in kelvin.	290.0
s	float	The insolation relative to the present one.	1.0
epsilon	float	The effective emissivity, lowered by the greenhouse effect.	0.65
sigma	float	The Stefan-Boltzmann constant divided by the mean insolation.	1.6e-10
albedo	float	The albedo at the middle of the transition from ice to no ice.	0.5
contrast	float	Half the difference between the albedos of a frozen and of a warm planet.	0.2
T_ice	float	The temperature of the middle of the transition.	263.0
width	float	The width of the transition in kelvin.	4.0

Returns:

Name	Type	Description
system	System	The model, with the temperature earth.T.

Examples:

>>> from skcomplex.models.physics.climate import energy_balance
>>> energy_balance(T=230.0).elements[0].T
230.0

Source code in src/skcomplex/models/physics/climate/_energy_balance.py

@register_model(
    field='physics/climate',
    tags=('ode', 'bistability', 'ice-albedo feedback', 'energy balance'),
    references=('Budyko 1969', 'Sellers 1969', 'Datseris and Parlitz 2022', 'Ghil and Lucarini 2020'),
    defaults={'until': 500.0, 'dt': 5.0},
)
def energy_balance(
    T: float = 290.0,
    *,
    s: float = 1.0,
    epsilon: float = 0.65,
    sigma: float = 1.6e-10,
    albedo: float = 0.5,
    contrast: float = 0.2,
    T_ice: float = 263.0,
    width: float = 4.0,
) -> System:
    """Return the zero-dimensional energy balance of the Earth, bistable through the ice-albedo feedback.

    In units of the mean insolation, the global temperature follows `dT/dt = s (1 - a(T)) - sigma epsilon T^4`, with
    the albedo `a(T) = albedo - contrast tanh((T - T_ice) / width)`, high on an icy planet and low on a warm one. For
    the default parameters a warm climate near 286 K and a frozen climate near 232 K are both stable.

    Args:
        T:
            The initial global mean temperature in kelvin.

        s:
            The insolation relative to the present one.

        epsilon:
            The effective emissivity, lowered by the greenhouse effect.

        sigma:
            The Stefan-Boltzmann constant divided by the mean insolation.

        albedo:
            The albedo at the middle of the transition from ice to no ice.

        contrast:
            Half the difference between the albedos of a frozen and of a warm planet.

        T_ice:
            The temperature of the middle of the transition.

        width:
            The width of the transition in kelvin.

    Returns:
        system:
            The model, with the temperature `earth.T`.

    Examples:
        >>> from skcomplex.models.physics.climate import energy_balance
        >>> energy_balance(T=230.0).elements[0].T
        230.0
    """
    climate = IceAlbedoClimateElement(
        'earth',
        T=T,
        s=s,
        epsilon=epsilon,
        sigma=sigma,
        albedo=albedo,
        contrast=contrast,
        T_ice=T_ice,
        width=width,
    )
    return System(elements=[climate])

energy_return(reserve=1000.0, capacity=0.1, *, initial_eroi=20.0, reinvestment=0.1, birth_rate=0.03)

Return the extraction of a finite energy reserve whose energy return on investment declines as it depletes.

The extraction capacity draws the gross energy from the reserve at a rate that falls with the remaining fraction, since the easy sources go first. Part of the gross energy pays for the extraction work, the gross energy divided by the energy return on investment (EROI), which falls from initial_eroi to one with the cube of the remaining fraction as the reserve empties. The rest is the net energy that society receives. Society reinvests a share of the net energy in extraction capacity, which depreciates, and the population grows exponentially with no feedback from the energy it receives. Energy is conserved, every unit leaving the reserve is either cost or net energy, so the gross energy rises exponentially and peaks, the net energy peaks before it, and the net energy per person peaks earlier still. The structure follows Mobus (2022, Section 9.6) and its figures of peak net energy, with illustrative parameters in energy units, persons, and years.

Parameters:

Name	Type	Description	Default
reserve	float	The initial energy reserve, also its full size.	1000.0
capacity	float	The initial extraction capacity, the gross energy per year from a full reserve.	0.1
initial_eroi	float	The energy return on investment of a full reserve.	20.0
reinvestment	float	The share of the net energy reinvested in extraction capacity.	0.1
birth_rate	float	The relative birth rate of the population per year.	0.03

Returns:

Name	Type	Description
system	System	The model, with the reserve energy.reserve, the capacity energy.capacity, the population energy.population, and the flows energy.energy_cost, energy.net_energy, energy.births, and energy.deaths per year, over years. The observables of the model are energy.eroi, energy.gross_energy, energy.extraction_cost, and energy.net_per_capita.

Examples:

>>> from skcomplex.models.physics.climate import energy_return
>>> energy_return(initial_eroi=20.0).elements[0].initial_eroi
20.0

Source code in src/skcomplex/models/physics/climate/_net_energy.py

@register_model(
    field='physics/climate',
    tags=('energy', 'net energy', 'stocks and flows', 'system dynamics', 'ode'),
    references=('Mobus 2022',),
    defaults={'until': 250.0, 'dt': 1.0},
    observables=lambda system: [
        ObservableSpec('eroi', _eroi, element='energy'),
        ObservableSpec('gross_energy', _gross_energy, element='energy'),
        ObservableSpec('extraction_cost', _extraction_cost, element='energy'),
        ObservableSpec('net_per_capita', _net_per_capita, element='energy'),
    ],
)
def energy_return(
    reserve: float = 1000.0,
    capacity: float = 0.1,
    *,
    initial_eroi: float = 20.0,
    reinvestment: float = 0.1,
    birth_rate: float = 0.03,
) -> System:
    """Return the extraction of a finite energy reserve whose energy return on investment declines as it depletes.

    The extraction capacity draws the gross energy from the reserve at a rate that falls with the remaining fraction,
    since the easy sources go first. Part of the gross energy pays for the extraction work, the gross energy divided
    by the energy return on investment (EROI), which falls from `initial_eroi` to one with the cube of the remaining
    fraction as the reserve empties. The rest is the net energy that society receives. Society reinvests a share of
    the net energy in extraction capacity, which depreciates, and the population grows exponentially with no feedback
    from the energy it receives. Energy is conserved, every unit leaving the reserve is either cost or net energy, so
    the gross energy rises exponentially and peaks, the net energy peaks before it, and the net energy per person
    peaks earlier still. The structure follows Mobus (2022, Section 9.6) and its figures of peak net energy, with
    illustrative parameters in energy units, persons, and years.

    Args:
        reserve:
            The initial energy reserve, also its full size.

        capacity:
            The initial extraction capacity, the gross energy per year from a full reserve.

        initial_eroi:
            The energy return on investment of a full reserve.

        reinvestment:
            The share of the net energy reinvested in extraction capacity.

        birth_rate:
            The relative birth rate of the population per year.

    Returns:
        system:
            The model, with the reserve `energy.reserve`, the capacity `energy.capacity`, the population
            `energy.population`, and the flows `energy.energy_cost`, `energy.net_energy`, `energy.births`, and
            `energy.deaths` per year, over years. The observables of the model are `energy.eroi`,
            `energy.gross_energy`, `energy.extraction_cost`, and `energy.net_per_capita`.

    Examples:
        >>> from skcomplex.models.physics.climate import energy_return
        >>> energy_return(initial_eroi=20.0).elements[0].initial_eroi
        20.0
    """
    sector = EnergySectorElement(
        'energy',
        reserve=reserve,
        capacity=capacity,
        reserve_size=reserve,
        initial_eroi=initial_eroi,
        reinvestment=reinvestment,
        birth_rate=birth_rate,
    )
    return System(elements=[sector])

enso(T=0.1, *, alpha=0.75, delta=6.0)

Return the delayed oscillator of El Niño of Suarez and Schopf, dT/dt = T - T^3 - alpha T(t - delta).

The sea surface temperature anomaly of the eastern Pacific grows by local feedback and is damped by the delayed return of oceanic waves reflected at the western boundary. For alpha < 1 and a long enough delay the warm and cold states lose stability to an oscillation between El Niño and La Niña. The history before the start is the initial anomaly.

Parameters:

Name	Type	Description	Default
T	float	The initial temperature anomaly, also its history.	0.1
alpha	float	The strength of the delayed negative feedback.	0.75
delta	float	The delay of the wave transit.	6.0

Returns:

Name	Type	Description
system	System	The model, with the anomaly pacific.T.

Examples:

>>> from skcomplex.models.physics.climate import enso
>>> enso(delta=4.0).elements[0].delta
4.0

Source code in src/skcomplex/models/physics/climate/_variability.py

@register_model(
    field='physics/climate',
    tags=('delay', 'oscillation', 'El Niño', 'ENSO'),
    references=('Suarez and Schopf 1988', 'Battisti and Hirst 1989', 'Ghil and Lucarini 2020'),
    defaults={'until': 200.0, 'dt': 0.5},
)
def enso(T: float = 0.1, *, alpha: float = 0.75, delta: float = 6.0) -> System:
    """Return the delayed oscillator of El Niño of Suarez and Schopf, `dT/dt = T - T^3 - alpha T(t - delta)`.

    The sea surface temperature anomaly of the eastern Pacific grows by local feedback and is damped by the delayed
    return of oceanic waves reflected at the western boundary. For `alpha < 1` and a long enough delay the warm and cold
    states lose stability to an oscillation between El Niño and La Niña. The history before the start is the initial
    anomaly.

    Args:
        T:
            The initial temperature anomaly, also its history.

        alpha:
            The strength of the delayed negative feedback.

        delta:
            The delay of the wave transit.

    Returns:
        system:
            The model, with the anomaly `pacific.T`.

    Examples:
        >>> from skcomplex.models.physics.climate import enso
        >>> enso(delta=4.0).elements[0].delta
        4.0
    """
    return System(elements=[DelayedOscillatorElement('pacific', T=T, alpha=alpha, delta=delta)])

geodynamo(Q=0.1, D=0.1, V=0.1, *, mu=0.119, nu=0.1, Gamma=0.9)

Return Gissinger's geodynamo, whose dipole field reverses at irregular times.

The quadrupolar field Q, the dipolar field D, and the velocity V follow dQ/dt = mu Q - V D, dD/dt = -nu D + V Q, and dV/dt = Gamma - V + Q D. The equations are symmetric under (Q, D) -> (-Q, -D), so the chaotic flow visits both polarities of the dipole.

Parameters:

Name	Type	Description	Default
Q	float	The initial quadrupolar field.	0.1
D	float	The initial dipolar field.	0.1
V	float	The initial velocity.	0.1
mu	float	The growth rate of the quadrupole.	0.119
nu	float	The damping rate of the dipole.	0.1
Gamma	float	The forcing of the flow.	0.9

Returns:

Name	Type	Description
system	System	The model, with the states dynamo.Q, dynamo.D, and dynamo.V.

Examples:

>>> from skcomplex.models.physics.climate import geodynamo
>>> geodynamo(mu=0.12).elements[0].mu
0.12

Source code in src/skcomplex/models/physics/climate/_low_order.py

@register_model(
    field='physics/climate',
    tags=('chaos', 'ode', 'geomagnetism', 'reversals'),
    references=('Gissinger 2012', 'Datseris and Parlitz 2022'),
    defaults={'until': 500.0, 'dt': 0.5},
)
def geodynamo(
    Q: float = 0.1,
    D: float = 0.1,
    V: float = 0.1,
    *,
    mu: float = 0.119,
    nu: float = 0.1,
    Gamma: float = 0.9,
) -> System:
    """Return Gissinger's geodynamo, whose dipole field reverses at irregular times.

    The quadrupolar field `Q`, the dipolar field `D`, and the velocity `V` follow `dQ/dt = mu Q - V D`,
    `dD/dt = -nu D + V Q`, and `dV/dt = Gamma - V + Q D`. The equations are symmetric under `(Q, D) -> (-Q, -D)`, so
    the chaotic flow visits both polarities of the dipole.

    Args:
        Q:
            The initial quadrupolar field.

        D:
            The initial dipolar field.

        V:
            The initial velocity.

        mu:
            The growth rate of the quadrupole.

        nu:
            The damping rate of the dipole.

        Gamma:
            The forcing of the flow.

    Returns:
        system:
            The model, with the states `dynamo.Q`, `dynamo.D`, and `dynamo.V`.

    Examples:
        >>> from skcomplex.models.physics.climate import geodynamo
        >>> geodynamo(mu=0.12).elements[0].mu
        0.12
    """
    return System(elements=[GeodynamoElement('dynamo', Q=Q, D=D, V=V, mu=mu, nu=nu, Gamma=Gamma)])

hasselmann(T=0.0, *, lam=0.5, sigma=0.3)

Return Hasselmann's stochastic climate, where fast weather noise drives a slow, damped temperature anomaly.

The anomaly follows dT = -lam T dt + sigma dW, an Ornstein-Uhlenbeck process whose red spectrum integrates white weather noise, with the stationary variance sigma^2 / (2 lam).

Parameters:

Name	Type	Description	Default
T	float	The initial temperature anomaly.	0.0
lam	float	The feedback rate that damps the anomaly.	0.5
sigma	float	The amplitude of the weather noise.	0.3

Returns:

Name	Type	Description
system	System	The model, with the anomaly ocean.T.

Examples:

>>> from skcomplex.models.physics.climate import hasselmann
>>> hasselmann(lam=0.1).elements[0].lam
0.1

Source code in src/skcomplex/models/physics/climate/_variability.py

@register_model(
    field='physics/climate',
    tags=('sde', 'Ornstein-Uhlenbeck', 'red noise', 'stochastic climate'),
    references=('Hasselmann 1976', 'Ghil and Lucarini 2020'),
    defaults={'until': 100.0, 'dt': 0.1, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def hasselmann(T: float = 0.0, *, lam: float = 0.5, sigma: float = 0.3) -> System:
    """Return Hasselmann's stochastic climate, where fast weather noise drives a slow, damped temperature anomaly.

    The anomaly follows `dT = -lam T dt + sigma dW`, an Ornstein-Uhlenbeck process whose red spectrum integrates white
    weather noise, with the stationary variance `sigma^2 / (2 lam)`.

    Args:
        T:
            The initial temperature anomaly.

        lam:
            The feedback rate that damps the anomaly.

        sigma:
            The amplitude of the weather noise.

    Returns:
        system:
            The model, with the anomaly `ocean.T`.

    Examples:
        >>> from skcomplex.models.physics.climate import hasselmann
        >>> hasselmann(lam=0.1).elements[0].lam
        0.1
    """
    feedback = RelaxationInteraction(
        'feedback',
        target='ocean',
        state='T',
        goal=0.0,
        time=_get_damping_time,
        noise=_get_weather,
    )
    return System(elements=[ClimateAnomalyElement('ocean', T=T, lam=lam, sigma=sigma)], interactions=[feedback])

lorenz_84(x=1.0, y=1.0, z=1.0, *, a=0.25, b=4.0, F=8.0, G=1.0)

Return Lorenz-84, the general circulation reduced to a westerly jet and a travelling wave, chaotic for F = 8.

The jet x and the cosine and sine phases y and z of the eddies follow dx/dt = -y^2 - z^2 - a x + a F, dy/dt = x y - b x z - y + G, and dz/dt = b x y + x z - z.

Parameters:

Name	Type	Description	Default
x	float	The initial strength of the westerly jet.	1.0
y	float	The initial cosine phase of the eddies.	1.0
z	float	The initial sine phase of the eddies.	1.0
a	float	The damping of the jet relative to the eddies.	0.25
b	float	The displacement of the eddies by the jet.	4.0
F	float	The forcing of the jet by the equator-to-pole temperature contrast.	8.0
G	float	The forcing of the eddies by the land-sea contrast.	1.0

Returns:

Name	Type	Description
system	System	The model, with the states atmosphere.x, atmosphere.y, and atmosphere.z.

Examples:

>>> from skcomplex.models.physics.climate import lorenz_84
>>> lorenz_84(G=1.25).elements[0].G
1.25

Source code in src/skcomplex/models/physics/climate/_low_order.py

@register_model(
    field='physics/climate',
    tags=('chaos', 'ode', 'atmosphere', 'general circulation'),
    references=('Lorenz 1984', 'Datseris and Parlitz 2022'),
    defaults={'until': 100.0, 'dt': 0.05},
)
def lorenz_84(
    x: float = 1.0,
    y: float = 1.0,
    z: float = 1.0,
    *,
    a: float = 0.25,
    b: float = 4.0,
    F: float = 8.0,
    G: float = 1.0,
) -> System:
    """Return Lorenz-84, the general circulation reduced to a westerly jet and a travelling wave, chaotic for `F = 8`.

    The jet `x` and the cosine and sine phases `y` and `z` of the eddies follow `dx/dt = -y^2 - z^2 - a x + a F`,
    `dy/dt = x y - b x z - y + G`, and `dz/dt = b x y + x z - z`.

    Args:
        x:
            The initial strength of the westerly jet.

        y:
            The initial cosine phase of the eddies.

        z:
            The initial sine phase of the eddies.

        a:
            The damping of the jet relative to the eddies.

        b:
            The displacement of the eddies by the jet.

        F:
            The forcing of the jet by the equator-to-pole temperature contrast.

        G:
            The forcing of the eddies by the land-sea contrast.

    Returns:
        system:
            The model, with the states `atmosphere.x`, `atmosphere.y`, and `atmosphere.z`.

    Examples:
        >>> from skcomplex.models.physics.climate import lorenz_84
        >>> lorenz_84(G=1.25).elements[0].G
        1.25
    """
    return System(elements=[Lorenz84Element('atmosphere', x=x, y=y, z=z, a=a, b=b, F=F, G=G)])

lorenz_96(n_sites=20, *, F=8.0, perturbation=0.01)

Return Lorenz-96, dx_i/dt = (x_{i+1} - x_{i-2}) x_{i-1} - x_i + F on a ring of sites around a latitude circle.

The sites start at the steady state F, and the first one is displaced by perturbation. For F = 8 the steady state is unstable and the ring is chaotic.

Parameters:

Name	Type	Description	Default
n_sites	int	The number of sites, at least 4, the dimension of the state x.	20
F	float	The constant forcing.	8.0
perturbation	float	The displacement of the first site from the steady state.	0.01

Returns:

Name	Type	Description
system	System	The model, with the values ring.x of all sites.

Examples:

>>> from skcomplex.models.physics.climate import lorenz_96
>>> lorenz_96(6).space.n_dim
6

Source code in src/skcomplex/models/physics/climate/_low_order.py

@register_model(
    field='physics/climate',
    tags=('chaos', 'ode', 'atmosphere', 'spatiotemporal chaos', 'data assimilation'),
    references=('Lorenz 1996', 'Datseris and Parlitz 2022'),
    defaults={'until': 20.0, 'dt': 0.05},
)
def lorenz_96(
    n_sites: int = 20,
    *,
    F: float = 8.0,
    perturbation: float = 0.01,
) -> System:
    """Return Lorenz-96, `dx_i/dt = (x_{i+1} - x_{i-2}) x_{i-1} - x_i + F` on a ring of sites around a latitude circle.

    The sites start at the steady state `F`, and the first one is displaced by `perturbation`. For `F = 8` the
    steady state is unstable and the ring is chaotic.

    Args:
        n_sites:
            The number of sites, at least 4, the dimension of the state `x`.

        F:
            The constant forcing.

        perturbation:
            The displacement of the first site from the steady state.

    Returns:
        system:
            The model, with the values `ring.x` of all sites.

    Examples:
        >>> from skcomplex.models.physics.climate import lorenz_96
        >>> lorenz_96(6).space.n_dim
        6
    """
    values = np.full(n_sites, float(F))
    values[0] += perturbation
    return System(elements=[Lorenz96Element('ring', x=values, F=F)], space=ContinuousSpace(n_dim=n_sites))

ramped_stochastic_fold(*, rate=0.004, sigma=0.02, a=-1.0, b=1.0, phi0=-1.0, phi1=1.0)

Return a stochastic fold whose control is ramped through the tipping point of its lower branch.

The state follows dx = (a x^3 + b x + phi(t)) dt + sigma dW with phi = phi0 + rate t until it reaches phi1, Eq. 4 of Fang et al. with phi = -1 + 0.004 t. It starts on the lower branch, which vanishes at the static fold phi_c = sqrt(-4 b^3 / (27 a)), 2 / (3 sqrt(3)) for a = -1 and b = 1. Near the fold, recovery slows down, so the variance and the lag-1 autocorrelation of the state rise. A ramp faster than the relaxation carries the control past phi_c before the state tips, by a lag that grows with the rate.

Parameters:

Name	Type	Description	Default
rate	float	The rate at which the control rises.	0.004
sigma	float	The amplitude of the noise.	0.02
a	float	The cubic coefficient, negative.	-1.0
b	float	The linear coefficient, positive.	1.0
phi0	float	The initial control.	-1.0
phi1	float	The final control.	1.0

Returns:

Name	Type	Description
system	System	The model, with the state system.x, starting on the lower branch at phi0.

Examples:

>>> from skcomplex.models.physics.climate import ramped_stochastic_fold
>>> ramped_stochastic_fold(rate=0.01).elements[0].rate
0.01

Source code in src/skcomplex/models/physics/climate/_tipping_types.py

@register_model(
    field='physics/climate',
    tags=('tipping points', 'fold', 'ramp', 'sde', 'early warning signals', 'critical slowing down'),
    references=('Fang, Wang, Kurths, and Fan 2025', 'Scheffer et al. 2009'),
    defaults={'until': 500.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def ramped_stochastic_fold(
    *,
    rate: float = 0.004,
    sigma: float = 0.02,
    a: float = -1.0,
    b: float = 1.0,
    phi0: float = -1.0,
    phi1: float = 1.0,
) -> System:
    """Return a stochastic fold whose control is ramped through the tipping point of its lower branch.

    The state follows `dx = (a x^3 + b x + phi(t)) dt + sigma dW` with `phi = phi0 + rate t` until it reaches `phi1`,
    Eq. 4 of Fang et al. with `phi = -1 + 0.004 t`. It starts on the lower branch, which vanishes at the static fold
    `phi_c = sqrt(-4 b^3 / (27 a))`, `2 / (3 sqrt(3))` for `a = -1` and `b = 1`. Near the fold, recovery slows down,
    so the variance and the lag-1 autocorrelation of the state rise. A ramp faster than the relaxation carries the
    control past `phi_c` before the state tips, by a lag that grows with the rate.

    Args:
        rate:
            The rate at which the control rises.

        sigma:
            The amplitude of the noise.

        a:
            The cubic coefficient, negative.

        b:
            The linear coefficient, positive.

        phi0:
            The initial control.

        phi1:
            The final control.

    Returns:
        system:
            The model, with the state `system.x`, starting on the lower branch at `phi0`.

    Examples:
        >>> from skcomplex.models.physics.climate import ramped_stochastic_fold
        >>> ramped_stochastic_fold(rate=0.01).elements[0].rate
        0.01
    """
    fold = RampedFoldElement(
        'system',
        x=_lower_branch(a, b, phi0),
        a=a,
        b=b,
        sigma=sigma,
        phi0=phi0,
        phi1=phi1,
        rate=rate,
    )
    return System(elements=[fold])

saltzman_glaciation(x=-0.2, y=1.0, z=0.0, *, p=1.0, q=2.5, r=1.3, s=0.6, u=0.2, nu=0.0)

Return the glacial cycles of Saltzman and Maasch, a self-sustained oscillation of ice, carbon, and ocean.

The anomalies of global ice volume x, atmospheric carbon dioxide y, and deep ocean temperature z follow dx/dt = -x - y - u z - nu F(t), dy/dt = -p z + r y - s y^2 - y^3, and dz/dt = -q (x + z). The time unit is 10 kyr, and F(t) = cos(2π t / 2.21) - cos(2π t / 4.41) is a quasiperiodic orbital forcing with the periods of precession and obliquity. Without forcing the default parameters give a limit cycle of about 100 kyr.

Parameters:

Name	Type	Description	Default
x	float	The initial ice volume anomaly.	-0.2
y	float	The initial carbon dioxide anomaly.	1.0
z	float	The initial deep ocean temperature anomaly.	0.0
p	float	The effect of the deep ocean temperature on carbon dioxide.	1.0
q	float	The relaxation rate of the deep ocean temperature to the ice volume.	2.5
r	float	The growth rate of carbon dioxide by its own feedback.	1.3
s	float	The quadratic saturation of carbon dioxide.	0.6
u	float	The effect of the deep ocean temperature on ice volume.	0.2
nu	float	The amplitude of the orbital forcing.	0.0

Returns:

Name	Type	Description
system	System	The model, with the anomalies climate.x, climate.y, and climate.z.

Examples:

>>> from skcomplex.models.physics.climate import saltzman_glaciation
>>> saltzman_glaciation(nu=0.5).elements[0].nu
0.5

Source code in src/skcomplex/models/physics/climate/_ice_ages.py

@register_model(
    field='physics/climate',
    tags=('ode', 'limit cycle', 'ice ages', 'paleoclimate'),
    references=('Saltzman and Maasch 1991', 'Datseris and Parlitz 2022'),
    defaults={'until': 50.0, 'dt': 0.05},
)
def saltzman_glaciation(
    x: float = -0.2,
    y: float = 1.0,
    z: float = 0.0,
    *,
    p: float = 1.0,
    q: float = 2.5,
    r: float = 1.3,
    s: float = 0.6,
    u: float = 0.2,
    nu: float = 0.0,
) -> System:
    """Return the glacial cycles of Saltzman and Maasch, a self-sustained oscillation of ice, carbon, and ocean.

    The anomalies of global ice volume `x`, atmospheric carbon dioxide `y`, and deep ocean temperature `z` follow
    `dx/dt = -x - y - u z - nu F(t)`, `dy/dt = -p z + r y - s y^2 - y^3`, and `dz/dt = -q (x + z)`. The time unit is
    10 kyr, and `F(t) = cos(2π t / 2.21) - cos(2π t / 4.41)` is a quasiperiodic orbital forcing with the periods of
    precession and obliquity. Without forcing the default parameters give a limit cycle of about 100 kyr.

    Args:
        x:
            The initial ice volume anomaly.

        y:
            The initial carbon dioxide anomaly.

        z:
            The initial deep ocean temperature anomaly.

        p:
            The effect of the deep ocean temperature on carbon dioxide.

        q:
            The relaxation rate of the deep ocean temperature to the ice volume.

        r:
            The growth rate of carbon dioxide by its own feedback.

        s:
            The quadratic saturation of carbon dioxide.

        u:
            The effect of the deep ocean temperature on ice volume.

        nu:
            The amplitude of the orbital forcing.

    Returns:
        system:
            The model, with the anomalies `climate.x`, `climate.y`, and `climate.z`.

    Examples:
        >>> from skcomplex.models.physics.climate import saltzman_glaciation
        >>> saltzman_glaciation(nu=0.5).elements[0].nu
        0.5
    """
    return System(elements=[SaltzmanMaaschElement('climate', x=x, y=y, z=z, p=p, q=q, r=r, s=s, u=u, nu=nu)])

soil_temperature(x=0.3, y=0.1, *, r=3.8, epsilon=0.1)

Return the coupled maps of the soil surface temperature and the deeper soil temperature.

The dimensionless surface temperature x follows the energy balance of the ground, reduced to a logistic map, and the deeper soil temperature y relaxes towards it: x_{n+1} = r x_n (1 - x_n) + epsilon y_n and y_{n+1} = epsilon (x_n + y_n). Apart from the origin, the fixed point is x* = (r + epsilon^2 / (1 - epsilon) - 1) / r and y* = epsilon x* / (1 - epsilon).

Parameters:

Name	Type	Description	Default
x	float	The initial dimensionless surface temperature.	0.3
y	float	The initial dimensionless deeper soil temperature.	0.1
r	float	The logistic parameter of the surface energy balance.	3.8
epsilon	float	The coupling between the surface and the deeper soil.	0.1

Returns:

Name	Type	Description
system	System	The model, with the temperatures soil.x and soil.y.

Examples:

>>> from skcomplex.models.physics.climate import soil_temperature
>>> soil_temperature(r=2.5).elements[0].r
2.5

Source code in src/skcomplex/models/physics/climate/_soil.py

@register_model(
    field='physics/climate',
    tags=('discrete map', 'coupled maps', 'chaos', 'land surface'),
    references=('Mimić, Mihailović, and Budinčević 2013', 'Mihailović et al. 2023'),
    defaults={'until': 200, 'dt': 1},
)
def soil_temperature(x: float = 0.3, y: float = 0.1, *, r: float = 3.8, epsilon: float = 0.1) -> System:
    """Return the coupled maps of the soil surface temperature and the deeper soil temperature.

    The dimensionless surface temperature `x` follows the energy balance of the ground, reduced to a logistic map,
    and the deeper soil temperature `y` relaxes towards it: `x_{n+1} = r x_n (1 - x_n) + epsilon y_n` and
    `y_{n+1} = epsilon (x_n + y_n)`. Apart from the origin, the fixed point is
    `x* = (r + epsilon^2 / (1 - epsilon) - 1) / r` and `y* = epsilon x* / (1 - epsilon)`.

    Args:
        x:
            The initial dimensionless surface temperature.

        y:
            The initial dimensionless deeper soil temperature.

        r:
            The logistic parameter of the surface energy balance.

        epsilon:
            The coupling between the surface and the deeper soil.

    Returns:
        system:
            The model, with the temperatures `soil.x` and `soil.y`.

    Examples:
        >>> from skcomplex.models.physics.climate import soil_temperature
        >>> soil_temperature(r=2.5).elements[0].r
        2.5
    """
    surface = MapInteraction(
        'surface',
        target='soil',
        state='x',
        function=lambda s: s.target.r * s.target.x * (1 - s.target.x) + s.target.epsilon * s.target.y,
    )
    depth = MapInteraction(
        'depth',
        target='soil',
        state='y',
        function=lambda s: s.target.epsilon * (s.target.x + s.target.y),
    )
    return System(
        elements=[SoilTemperatureElement('soil', x=x, y=y, r=r, epsilon=epsilon)],
        interactions=[surface, depth],
    )

sutera_fraedrich_tipping(tipping='noise', *, rho=0.18, nu=None)

Return the Sutera-Fraedrich stochastic energy balance set up for bifurcation-, noise-, or rate-induced tipping.

The temperature follows dT = c0 [-T^4 + c1 mu (b2 T^2 - a2)] dt + sqrt(nu) dW with the parameters of Table 1 of Fang et al. and starts at its warm stable equilibrium. In bifurcation-induced tipping, mu decreases at 4e-4 per year until the warm equilibrium meets the unstable one in a fold near mu = 0.76. In noise-induced tipping, the controls are fixed and noise of variance nu = 1 carries the state over the unstable equilibrium. In rate-induced tipping, b2 and a2 ramp logistically at the rate rho with the equilibria preserved, and the state fails to track the moving warm equilibrium when rho exceeds about 0.175. A tipped state leaves the warm branch and falls to the other stable root of the quartic.

Parameters:

Name	Type	Description	Default
tipping	str	The type of tipping, 'bifurcation', 'noise', or 'rate'.	'noise'
rho	float	The rate of the logistic ramp of rate-induced tipping, per year.	0.18
nu	float | None	The variance of the noise per year, or None for the value of the type, 1 for noise-induced tipping and 0 otherwise.	None

Returns:

Name	Type	Description
system	System	The model, with the temperature earth.T, the control earth.mu, and the progress earth.lam of the ramp.

Raises:

Type	Description
ValidationError	The type of tipping is unknown.

Examples:

>>> from skcomplex.models.physics.climate import sutera_fraedrich_tipping
>>> round(sutera_fraedrich_tipping('rate').elements[0].T, 1)
291.6

Source code in src/skcomplex/models/physics/climate/_tipping_types.py

@register_model(
    field='physics/climate',
    tags=('tipping points', 'energy balance', 'sde', 'B-tipping', 'N-tipping', 'R-tipping'),
    references=('Fang, Wang, Kurths, and Fan 2025', 'Sutera 1981', 'Fraedrich 1978', 'Ashwin et al. 2012'),
    defaults={'until': 300.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=0.05)},
)
def sutera_fraedrich_tipping(tipping: str = 'noise', *, rho: float = 0.18, nu: float | None = None) -> System:
    """Return the Sutera-Fraedrich stochastic energy balance set up for bifurcation-, noise-, or rate-induced tipping.

    The temperature follows `dT = c0 [-T^4 + c1 mu (b2 T^2 - a2)] dt + sqrt(nu) dW` with the parameters of Table 1 of
    Fang et al. and starts at its warm stable equilibrium. In bifurcation-induced tipping, `mu` decreases at
    `4e-4` per year until the warm equilibrium meets the unstable one in a fold near `mu = 0.76`. In noise-induced
    tipping, the controls are fixed and noise of variance `nu = 1` carries the state over the unstable equilibrium. In
    rate-induced tipping, `b2` and `a2` ramp logistically at the rate `rho` with the equilibria preserved, and the state
    fails to track the moving warm equilibrium when `rho` exceeds about 0.175. A tipped state leaves the warm branch and
    falls to the other stable root of the quartic.

    Args:
        tipping:
            The type of tipping, `'bifurcation'`, `'noise'`, or `'rate'`.

        rho:
            The rate of the logistic ramp of rate-induced tipping, per year.

        nu:
            The variance of the noise per year, or `None` for the value of the type, 1 for noise-induced tipping and
            0 otherwise.

    Returns:
        system:
            The model, with the temperature `earth.T`, the control `earth.mu`, and the progress `earth.lam` of the ramp.

    Raises:
        ValidationError:
            The type of tipping is unknown.

    Examples:
        >>> from skcomplex.models.physics.climate import sutera_fraedrich_tipping
        >>> round(sutera_fraedrich_tipping('rate').elements[0].T, 1)
        291.6
    """
    if tipping not in TIPPING_TYPES:
        error_msg = f'Unknown tipping type {tipping!r}. Use one of {", ".join(TIPPING_TYPES)}.'
        raise ValidationError(error_msg)
    values = TIPPING_TYPES[tipping]
    earth = StochasticEnergyBalanceElement(
        'earth',
        T=_warm_state(1.0, values['b_ini'], values['a_ini']),
        lam=values['lam'],
        drift=values['drift'],
        rho=rho if tipping == 'rate' else 0.0,
        b_ini=values['b_ini'],
        b_fin=values['b_fin'],
        a_ini=values['a_ini'],
        nu=values['nu'] if nu is None else nu,
    )
    return System(elements=[earth])

tipping_cascade_network(graph=None, *, c=0.0, d=0.2, x=-1.0)

Return a network of tipping elements, where one element pushed past its tipping point can tip its neighbours.

Every node is a fold dx_i/dt = -x_i^3 + x_i + sqrt(4/27) c_i + d Σ_j A_ij (x_j + 1) that starts in its lower state. An element whose control passes the critical value 1 tips to its upper state, and its coupling then raises the effective control of each neighbour by 2 d / sqrt(4/27), which can tip neighbours that would not tip alone.

Parameters:

Name	Type	Description	Default
graph	Graph | None	The NetworkX graph of the interactions, or None for a path of four elements.	None
c	float | Sequence[float]	The control of every element, one value or one per node in node order.	0.0
d	float	The coupling strength of every link.	0.2
x	float | Sequence[float]	The initial states, one value or one per node in node order.	-1.0

Returns:

Name	Type	Description
system	System	The model, with the states elements.x, one value per node.

Examples:

>>> import networkx as nx
>>> from skcomplex.models.physics.climate import tipping_cascade_network
>>> tipping_cascade_network(nx.cycle_graph(5), c=[1.2, 0, 0, 0, 0]).elements[0].check().prepare().size_
5

Source code in src/skcomplex/models/physics/climate/_tipping.py

@register_model(
    field='physics/climate',
    tags=('tipping points', 'cascades', 'network', 'bistability', 'ode'),
    references=('Wunderling et al. 2021', 'Klose et al. 2020'),
    defaults={'until': 40.0, 'dt': 0.1},
)
def tipping_cascade_network(
    graph: nx.Graph | None = None,
    *,
    c: float | Sequence[float] = 0.0,
    d: float = 0.2,
    x: float | Sequence[float] = -1.0,
) -> System:
    """Return a network of tipping elements, where one element pushed past its tipping point can tip its neighbours.

    Every node is a fold `dx_i/dt = -x_i^3 + x_i + sqrt(4/27) c_i + d Σ_j A_ij (x_j + 1)` that starts in its lower
    state. An element whose control passes the critical value 1 tips to its upper state, and its coupling then raises
    the effective control of each neighbour by `2 d / sqrt(4/27)`, which can tip neighbours that would not tip alone.

    Args:
        graph:
            The NetworkX graph of the interactions, or `None` for a path of four elements.

        c:
            The control of every element, one value or one per node in node order.

        d:
            The coupling strength of every link.

        x:
            The initial states, one value or one per node in node order.

    Returns:
        system:
            The model, with the states `elements.x`, one value per node.

    Examples:
        >>> import networkx as nx
        >>> from skcomplex.models.physics.climate import tipping_cascade_network
        >>> tipping_cascade_network(nx.cycle_graph(5), c=[1.2, 0, 0, 0, 0]).elements[0].check().prepare().size_
        5
    """
    graph = nx.path_graph(4) if graph is None else graph
    size = graph.number_of_nodes()
    values = {
        'x': np.broadcast_to(np.asarray(x, dtype=float), size).copy(),
        'c': np.broadcast_to(np.asarray(c, dtype=float), size).copy(),
        'd': d,
    }
    elements = FieldElement('elements', element=TippingElementElement, space=NetworkSpace(graph), values=values)
    coupling = CouplingInteraction('coupling', target='elements', function=_push)
    return System(elements=[elements], interactions=[coupling], time=ContinuousTime())

van_der_pol_ice_ages(x=0.0, y=1.0, *, alpha=30.0, beta=0.75, gamma=0.0, tau=35.09)

Return glacial cycles as a Van der Pol relaxation oscillator, optionally paced by orbital forcing.

Time is in thousands of years. The slow variable x stands for ice volume and the fast variable y for the climate state, with tau x' = -(y + beta + gamma F(t)) and tau y' = -alpha (y^3/3 - y - x), where F(t) = cos(2 pi t/23) + cos(2 pi t/41) stands for precession and obliquity. For |beta| < 1 the unforced oscillator has a relaxation limit cycle of slow ice build-up and fast terminations, and gamma sets how strongly the orbit paces it.

Parameters:

Name	Type	Description	Default
x	float	The initial ice volume.	0.0
y	float	The initial climate state.	1.0
alpha	float	The time-scale separation of the fast climate state.	30.0
beta	float	The asymmetry of the cycle, between build-up and termination.	0.75
gamma	float	The strength of the orbital forcing.	0.0
tau	float	The time scale, in thousands of years.	35.09

Returns:

Name	Type	Description
system	System	The model, with the states climate.x and climate.y.

Examples:

>>> from skcomplex.models.physics.climate import van_der_pol_ice_ages
>>> van_der_pol_ice_ages(gamma=0.4).elements[0].gamma
0.4

Source code in src/skcomplex/models/physics/climate/_ice_ages.py

@register_model(
    field='physics/climate',
    tags=('ode', 'limit cycle', 'relaxation oscillator', 'ice ages', 'paleoclimate'),
    references=('De Saedeleer, Crucifix, and Wieczorek 2013', 'Datseris and Parlitz 2022'),
    defaults={'until': 400.0, 'dt': 0.5},
)
def van_der_pol_ice_ages(
    x: float = 0.0,
    y: float = 1.0,
    *,
    alpha: float = 30.0,
    beta: float = 0.75,
    gamma: float = 0.0,
    tau: float = 35.09,
) -> System:
    """Return glacial cycles as a Van der Pol relaxation oscillator, optionally paced by orbital forcing.

    Time is in thousands of years. The slow variable `x` stands for ice volume and the fast variable `y` for the climate
    state, with `tau x' = -(y + beta + gamma F(t))` and `tau y' = -alpha (y^3/3 - y - x)`, where
    `F(t) = cos(2 pi t/23) + cos(2 pi t/41)` stands for precession and obliquity. For `|beta| < 1` the unforced
    oscillator has a relaxation limit cycle of slow ice build-up and fast terminations, and `gamma` sets how strongly
    the orbit paces it.

    Args:
        x:
            The initial ice volume.

        y:
            The initial climate state.

        alpha:
            The time-scale separation of the fast climate state.

        beta:
            The asymmetry of the cycle, between build-up and termination.

        gamma:
            The strength of the orbital forcing.

        tau:
            The time scale, in thousands of years.

    Returns:
        system:
            The model, with the states `climate.x` and `climate.y`.

    Examples:
        >>> from skcomplex.models.physics.climate import van_der_pol_ice_ages
        >>> van_der_pol_ice_ages(gamma=0.4).elements[0].gamma
        0.4
    """
    return System(
        elements=[VanDerPolIceAgesElement('climate', x=x, y=y, alpha=alpha, beta=beta, gamma=gamma, tau=tau)],
    )