Climate

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

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.systems.physics.climate import CarbonCycleElement
>>> CarbonCycleElement('ocean').c_x
55.0

Source code in src/skcomplex/systems/physics/climate/_ice_ages.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.systems.physics.climate import ClimateAnomalyElement
>>> ClimateAnomalyElement('ocean').lam
0.5

Source code in src/skcomplex/systems/physics/climate/_variability.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.systems.physics.climate import DaisyworldElement
>>> DaisyworldElement('planet').albedo_white
0.75

Source code in src/skcomplex/systems/physics/climate/_energy_balance.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.systems.physics.climate import DelayedOscillatorElement
>>> DelayedOscillatorElement('pacific').delta
6.0

Source code in src/skcomplex/systems/physics/climate/_variability.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.systems.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).prepare(x).run().autocorrelation_trend_[0] > 0.5)
True

Source code in src/skcomplex/systems/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.systems.physics.climate import EnergySectorElement
>>> EnergySectorElement('energy').initial_eroi
20.0

Source code in src/skcomplex/systems/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

FoldFollowerElement(label, visual=None, *, y=-1.0, c=0.0, sigma=0.02)

Bases: Element

Stochastic fold dy = (-y³ + y + c) dt + sigma dW, whose control a leading element shifts.

The fold has a lower and an upper stable branch while |c| < 2 / (3 √3), and a coupling adds to its control.

Parameters:

Name	Type	Description	Default
label	str	The label of the follower.	required
visual	VisualSpec | None	How the follower is drawn, or None for default_visual_.	None
y	Annotated[float, StateSpec()]	The initial state.	-1.0
c	Annotated[float, ParameterSpec()]	The control of the fold, without the coupling.	0.0
sigma	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The amplitude of the noise.	0.02

Attributes:

Name	Type	Description
default_visual_		The visual of the follower when it sets none.

Examples:

>>> from skcomplex.systems.physics.climate import FoldFollowerElement
>>> FoldFollowerElement('follower').c
0.0

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

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

ForestCoverElement(label, visual=None, *, v=0.8, Tv=30.0, g0=2.0, beta=10.0, Topt=28.0, a=5.0, gamma=0.2)

Bases: Element

Vegetation cover of a tropical forest, whose growth rate falls on both sides of an optimal temperature.

The cover v runs from 0 for bare soil to 1 for full forest. Bare soil is warmer than forest by a, so the local temperature is T = Tv + a (1 - v), and the growth rate g = g0 (1 - ((T - Topt) / beta)²) peaks at Topt. A rate interaction gives the logistic growth g v (1 - v) and the disturbance gamma v, so the element holds the state and the parameters it reads.

Parameters:

Name	Type	Description	Default
label	str	The label of the forest.	required
visual	VisualSpec | None	How the forest is drawn, or None for default_visual_.	None
v	Annotated[float, StateSpec(domain=closed(0, 1))]	The initial vegetation cover.	0.8
Tv	Annotated[float, ParameterSpec()]	The ambient temperature, the temperature under full forest cover, in degrees Celsius.	30.0
g0	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The largest growth rate.	2.0
beta	Annotated[float, ParameterSpec(domain=open(0, oo))]	The half width of the growth rate against the temperature.	10.0
Topt	Annotated[float, ParameterSpec()]	The optimal temperature of growth.	28.0
a	Annotated[float, ParameterSpec()]	The warming of bare soil over forest.	5.0
gamma	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The disturbance rate.	0.2

Attributes:

Name	Type	Description
default_visual_		The visual of the forest when it sets none.

Examples:

>>> from skcomplex.systems.physics.climate import ForestCoverElement
>>> ForestCoverElement('forest').Topt
28.0

Source code in src/skcomplex/systems/physics/climate/_forests.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    v: Annotated[float, StateSpec(domain=RealSet.closed(0, 1))] = 0.8,
    Tv: Annotated[float, ParameterSpec()] = 30.0,
    g0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2.0,
    beta: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
    Topt: Annotated[float, ParameterSpec()] = 28.0,
    a: Annotated[float, ParameterSpec()] = 5.0,
    gamma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
) -> None:
    super().__init__(label, visual=visual)
    self.v = v
    self.Tv = Tv
    self.g0 = g0
    self.beta = beta
    self.Topt = Topt
    self.a = a
    self.gamma = gamma

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.systems.physics.climate import GeodynamoElement
>>> GeodynamoElement('dynamo').Gamma
0.9

Source code in src/skcomplex/systems/physics/climate/_low_order.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

HopfFollowerElement(label, visual=None, *, y=0.1, z=0.0, mu=0.0, omega=1.0, sigma=0.02)

Bases: Element

Stochastic Hopf normal form dz = ((mu + i omega) z - |z|² z) dt + sigma dW, whose mu a leader shifts.

The complex amplitude z = y + i z_im is held as two real states. For mu < 0 it decays to a focus at the origin, and for mu > 0 it settles on a limit cycle of radius sqrt(mu), and a coupling adds to mu.

Parameters:

Name	Type	Description	Default
label	str	The label of the follower.	required
visual	VisualSpec | None	How the follower is drawn, or None for default_visual_.	None
y	Annotated[float, StateSpec()]	The initial real part of the amplitude.	0.1
z	Annotated[float, StateSpec()]	The initial imaginary part of the amplitude.	0.0
mu	Annotated[float, ParameterSpec()]	The distance from the Hopf bifurcation, without the coupling.	0.0
omega	Annotated[float, ParameterSpec()]	The frequency of the oscillation.	1.0
sigma	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The amplitude of the noise.	0.02

Attributes:

Name	Type	Description
default_visual_		The visual of the follower when it sets none.

Examples:

>>> from skcomplex.systems.physics.climate import HopfFollowerElement
>>> HopfFollowerElement('follower').omega
1.0

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

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    y: Annotated[float, StateSpec()] = 0.1,
    z: Annotated[float, StateSpec()] = 0.0,
    mu: Annotated[float, ParameterSpec()] = 0.0,
    omega: Annotated[float, ParameterSpec()] = 1.0,
    sigma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.02,
) -> None:
    super().__init__(label, visual=visual)
    self.y = y
    self.z = z
    self.mu = mu
    self.omega = omega
    self.sigma = sigma

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.systems.physics.climate import IceAlbedoClimateElement
>>> IceAlbedoClimateElement('earth').T_ice
263.0

Source code in src/skcomplex/systems/physics/climate/_energy_balance.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.systems.physics.climate import LatitudeBandElement
>>> LatitudeBandElement('band').albedo_ice
0.62

Source code in src/skcomplex/systems/physics/climate/_energy_balance.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.systems.physics.climate import Lorenz84Element
>>> Lorenz84Element('atmosphere').F
8.0

Source code in src/skcomplex/systems/physics/climate/_low_order.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.systems.physics.climate import Lorenz96Element
>>> Lorenz96Element('ring').F
8.0

Source code in src/skcomplex/systems/physics/climate/_low_order.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

MitigatorsElement(label, visual=None, *, x=0.5, alpha=1.0, delta=1.0, fmax=5.0, w=3.0, Tc=1.0, tp=10.0, tf=15.0)

Bases: Element

Share of mitigators in a population that weighs the cost of mitigation against the perceived cost of warming.

The share x follows the replicator equation x' = x (1 - x) (-alpha + f + delta (2 x - 1)), with the cost of mitigation alpha, the strength delta of the social norm, and the cost of warming f = fmax / (1 + exp(-w (Tf - Tc))). The perceived warming Tf = tf / tp (T(t) - T(t - tp)) projects the warming of the last tp years over tf years. A rate interaction reads the temperature it perceives, so the element holds the state and the parameters it reads.

Parameters:

Name	Type	Description	Default
label	str	The label of the population.	required
visual	VisualSpec | None	How the population is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec(domain=closed(0, 1))]	The initial share of mitigators.	0.5
alpha	Annotated[float, ParameterSpec()]	The cost of mitigation.	1.0
delta	Annotated[float, ParameterSpec()]	The strength of the social norm.	1.0
fmax	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The largest cost of warming.	5.0
w	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The steepness of the cost of warming.	3.0
Tc	Annotated[float, ParameterSpec()]	The critical warming.	1.0
tp	Annotated[float, ParameterSpec(domain=open(0, oo))]	The years of past warming that the projection reads.	10.0
tf	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The years ahead of the projection.	15.0

Attributes:

Name	Type	Description
default_visual_		The visual of the population when it sets none.

Examples:

>>> from skcomplex.systems.physics.climate import MitigatorsElement
>>> MitigatorsElement('people').tp
10.0

Source code in src/skcomplex/systems/physics/climate/_forests.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(domain=RealSet.closed(0, 1))] = 0.5,
    alpha: Annotated[float, ParameterSpec()] = 1.0,
    delta: Annotated[float, ParameterSpec()] = 1.0,
    fmax: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 5.0,
    w: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 3.0,
    Tc: Annotated[float, ParameterSpec()] = 1.0,
    tp: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
    tf: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 15.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.alpha = alpha
    self.delta = delta
    self.fmax = fmax
    self.w = w
    self.Tc = Tc
    self.tp = tp
    self.tf = tf

PollutedRegionElement(label, visual=None, *, population=100000000.0, pollution=10000000.0, pollution_normal=0.25, pollution_standard=15000000.0, birth_normal=0.04, death_normal=0.028, subsistence_food=180.0, cultivated_land=7400000.0, crop_yield=2200.0)

Bases: Element

Population of a region and the pollution it generates, which the environment absorbs ever more slowly.

The pollution grows by pollution_normal per person and is absorbed over a time that a table lengthens with the pollution ratio, the pollution over its standard. The births and the deaths are the normal rates times table multipliers of the food ratio, the food per capita over the subsistence food, and of the pollution ratio. The food is the cultivated land times the yield, a cropping intensity that a table raises over time, and a table multiplier of the pollution ratio.

Parameters:

Name	Type	Description	Default
label	str	The label of the region.	required
visual	VisualSpec | None	How the region is drawn, or None for default_visual_.	None
population	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial population.	100000000.0
pollution	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial pollution.	10000000.0
pollution_normal	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The pollution one person generates per year.	0.25
pollution_standard	Annotated[float, ParameterSpec(domain=open(0, oo))]	The pollution that the pollution ratio divides by.	15000000.0
birth_normal	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The normal birth rate per year.	0.04
death_normal	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The normal death rate per year.	0.028
subsistence_food	Annotated[float, ParameterSpec(domain=open(0, oo))]	The food per capita that the food ratio divides by.	180.0
cultivated_land	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The cultivated land.	7400000.0
crop_yield	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The yield of the land per crop.	2200.0

Attributes:

Name	Type	Description
default_visual_		The visual of the region when it sets none.

Examples:

>>> from skcomplex.systems.physics.climate import PollutedRegionElement
>>> PollutedRegionElement('region').birth_normal
0.04

Source code in src/skcomplex/systems/physics/climate/_pollution.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    population: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1e8,
    pollution: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1e7,
    pollution_normal: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.25,
    pollution_standard: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.5e7,
    birth_normal: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.04,
    death_normal: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.028,
    subsistence_food: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 180.0,
    cultivated_land: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 7.4e6,
    crop_yield: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2200.0,
) -> None:
    super().__init__(label, visual=visual)
    self.population = population
    self.pollution = pollution
    self.pollution_normal = pollution_normal
    self.pollution_standard = pollution_standard
    self.birth_normal = birth_normal
    self.death_normal = death_normal
    self.subsistence_food = subsistence_food
    self.cultivated_land = cultivated_land
    self.crop_yield = crop_yield

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.systems.physics.climate import RampedFoldElement
>>> RampedFoldElement('system').rate
0.004

Source code in src/skcomplex/systems/physics/climate/_tipping_types.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.systems.physics.climate import SaltzmanMaaschElement
>>> SaltzmanMaaschElement('climate').q
2.5

Source code in src/skcomplex/systems/physics/climate/_ice_ages.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.systems.physics.climate import SoilTemperatureElement
>>> SoilTemperatureElement('soil').r
3.8

Source code in src/skcomplex/systems/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.systems.physics.climate import StochasticEnergyBalanceElement
>>> StochasticEnergyBalanceElement('earth').nu
1.0

Source code in src/skcomplex/systems/physics/climate/_tipping_types.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.systems.physics.climate import TippingElementElement
>>> TippingElementElement('ice sheet').x
-1.0

Source code in src/skcomplex/systems/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.systems.physics.climate import VanDerPolIceAgesElement
>>> VanDerPolIceAgesElement('climate').alpha
30.0

Source code in src/skcomplex/systems/physics/climate/_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 system, with the temperatures bands.T from the equator to the pole.

Examples:

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

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

@register_system(
    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 system, with the temperatures `bands.T` from the equator to the pole.

    Examples:
        >>> from skcomplex.systems.physics.climate import budyko_sellers
        >>> budyko_sellers(10).elements[0].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 system, with the concentrations ocean.c and ocean.w.

Examples:

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

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

@register_system(
    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 system, with the concentrations `ocean.c` and `ocean.w`.

    Examples:
        >>> from skcomplex.systems.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 system, with the areas planet.white and planet.black and the observable temperature.

Examples:

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

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

@register_system(
    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 system, with the areas `planet.white` and `planet.black` and the observable `temperature`.

    Examples:
        >>> from skcomplex.systems.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 system, with the temperature earth.T.

Examples:

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

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

@register_system(
    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 system, with the temperature `earth.T`.

    Examples:
        >>> from skcomplex.systems.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 system, 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 system are energy.eroi, energy.gross_energy, energy.extraction_cost, and energy.net_per_capita.

Examples:

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

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

@register_system(
    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 system, 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 system are `energy.eroi`,
            `energy.gross_energy`, `energy.extraction_cost`, and `energy.net_per_capita`.

    Examples:
        >>> from skcomplex.systems.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 system, with the anomaly pacific.T.

Examples:

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

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

@register_system(
    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 system, with the anomaly `pacific.T`.

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

forest_dieback(v=0.8, *, Tv=30.0, g0=2.0, beta=10.0, Topt=28.0, a=5.0, gamma=0.2)

Return the forest dieback model of the Amazon, whose cover tips to bare soil when the ambient temperature rises.

The cover follows v' = g v (1 - v) - gamma v, with the growth rate g = g0 (1 - ((T - Topt) / beta)²) at the local temperature T = Tv + a (1 - v). Losing forest warms the soil, which lowers the growth and loses more forest, a positive feedback. Bare soil v = 0 is always a fixed point, and the forest survives at the root of (1 - v) g(v) = gamma. For the defaults the forest and bare soil are both stable for ambient temperatures between about 32.5 and 34.7 degrees, the bistable range of a hysteresis loop, and the forest dies beyond it.

Parameters:

Name	Type	Description	Default
v	float	The initial vegetation cover, between 0 and 1.	0.8
Tv	float	The ambient temperature under full forest cover, in degrees Celsius.	30.0
g0	float	The largest growth rate.	2.0
beta	float	The half width of the growth rate against the temperature.	10.0
Topt	float	The optimal temperature of growth.	28.0
a	float	The warming of bare soil over forest.	5.0
gamma	float	The disturbance rate.	0.2

Returns:

Name	Type	Description
system	System	The system, with the vegetation cover forest.v.

Examples:

>>> from skcomplex.systems.physics.climate import forest_dieback
>>> forest_dieback(Tv=33.0).elements[0].Tv
33.0

Source code in src/skcomplex/systems/physics/climate/_forests.py

@register_system(
    field='physics/climate',
    tags=('tipping points', 'vegetation', 'bistability', 'hysteresis', 'ode'),
    references=('Van Nes et al. 2014', 'Bury, Bauch, and Anand 2019', 'Liu and Bauch 2023'),
    defaults={'until': 200.0, 'dt': 0.5},
)
def forest_dieback(
    v: float = 0.8,
    *,
    Tv: float = 30.0,
    g0: float = 2.0,
    beta: float = 10.0,
    Topt: float = 28.0,
    a: float = 5.0,
    gamma: float = 0.2,
) -> System:
    """Return the forest dieback model of the Amazon, whose cover tips to bare soil when the ambient temperature rises.

    The cover follows `v' = g v (1 - v) - gamma v`, with the growth rate `g = g0 (1 - ((T - Topt) / beta)²)` at the
    local temperature `T = Tv + a (1 - v)`. Losing forest warms the soil, which lowers the growth and loses more forest,
    a positive feedback. Bare soil `v = 0` is always a fixed point, and the forest survives at the root of
    `(1 - v) g(v) = gamma`. For the defaults the forest and bare soil are both stable for ambient temperatures between
    about 32.5 and 34.7 degrees, the bistable range of a hysteresis loop, and the forest dies beyond it.

    Args:
        v:
            The initial vegetation cover, between 0 and 1.

        Tv:
            The ambient temperature under full forest cover, in degrees Celsius.

        g0:
            The largest growth rate.

        beta:
            The half width of the growth rate against the temperature.

        Topt:
            The optimal temperature of growth.

        a:
            The warming of bare soil over forest.

        gamma:
            The disturbance rate.

    Returns:
        system:
            The system, with the vegetation cover `forest.v`.

    Examples:
        >>> from skcomplex.systems.physics.climate import forest_dieback
        >>> forest_dieback(Tv=33.0).elements[0].Tv
        33.0
    """
    forest = ForestCoverElement('forest', v=v, Tv=Tv, g0=g0, beta=beta, Topt=Topt, a=a, gamma=gamma)
    growth = RateInteraction('growth', target='forest', state='v', rate=lambda s: _logistic(s.target, 1))
    return System(elements=[forest], interactions=[growth])

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 system, with the states dynamo.Q, dynamo.D, and dynamo.V.

Examples:

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

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

@register_system(
    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 system, with the states `dynamo.Q`, `dynamo.D`, and `dynamo.V`.

    Examples:
        >>> from skcomplex.systems.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 system, with the anomaly ocean.T.

Examples:

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

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

@register_system(
    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 system, with the anomaly `ocean.T`.

    Examples:
        >>> from skcomplex.systems.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])

leading_following_cascade(follower='fold', *, coupling=0.4, rate=0.004, sigma=0.02, phi0=-1.0, phi1=1.0, c=0.0, mu=0.0, omega=1.0)

Return a cascade of two tipping elements, a leading fold driven by a ramp and a follower forced by the leader.

The leading element follows dx = (-x³ + x + phi(t)) dt + sigma dW, whose control rises as phi = phi0 + rate t until phi1, and it tips from its lower to its upper branch when phi passes 2 / (3 √3). The following element is forced linearly by the state of the leader, Eq. 10 of Fang et al. after Dekker et al. A fold follower obeys dy = (-y³ + y + c + coupling x) dt + sigma dW, and it tips too, a domino cascade, when the tipped leader moves its control c + coupling x past 2 / (3 √3). A Hopf follower obeys dw = ((mu + coupling x + i omega) w - |w|² w) dt + sigma dW, and the tipped leader can carry it through its Hopf bifurcation onto an oscillation of radius sqrt(mu + coupling x). The follower can show rising variance and autocorrelation before the leader tips even when it stays far from its own tipping point.

Parameters:

Name	Type	Description	Default
follower	str	The normal form of the following element, 'fold' or 'hopf'.	'fold'
coupling	float	The strength of the linear forcing of the follower by the leader.	0.4
rate	float	The rate at which the control of the leader rises.	0.004
sigma	float	The amplitude of the noise of both elements.	0.02
phi0	float	The initial control of the leader.	-1.0
phi1	float	The final control of the leader.	1.0
c	float	The control of a fold follower, without the coupling.	0.0
mu	float	The distance of a Hopf follower from its bifurcation, without the coupling.	0.0
omega	float	The frequency of a Hopf follower.	1.0

Returns:

Name	Type	Description
system	System	The system, with the leader leader.x and the follower follower.y, and also follower.z for the Hopf follower.

Raises:

Type	Description
ValidationError	The follower is neither 'fold' nor 'hopf'.

Examples:

>>> from skcomplex.systems.physics.climate import leading_following_cascade
>>> leading_following_cascade('hopf').elements[1].omega
1.0

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

@register_system(
    field='physics/climate',
    tags=('tipping points', 'cascades', 'fold', 'hopf bifurcation', 'ramp', 'sde'),
    references=('Dekker, von der Heydt, and Dijkstra 2018', 'Fang, Wang, Kurths, and Fan 2025'),
    defaults={'until': 500.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def leading_following_cascade(
    follower: str = 'fold',
    *,
    coupling: float = 0.4,
    rate: float = 0.004,
    sigma: float = 0.02,
    phi0: float = -1.0,
    phi1: float = 1.0,
    c: float = 0.0,
    mu: float = 0.0,
    omega: float = 1.0,
) -> System:
    """Return a cascade of two tipping elements, a leading fold driven by a ramp and a follower forced by the leader.

    The leading element follows `dx = (-x³ + x + phi(t)) dt + sigma dW`, whose control rises as `phi = phi0 + rate t`
    until `phi1`, and it tips from its lower to its upper branch when `phi` passes `2 / (3 √3)`. The following element
    is forced linearly by the state of the leader, Eq. 10 of Fang et al. after Dekker et al. A fold follower obeys
    `dy = (-y³ + y + c + coupling x) dt + sigma dW`, and it tips too, a domino cascade, when the tipped leader moves
    its control `c + coupling x` past `2 / (3 √3)`. A Hopf follower obeys
    `dw = ((mu + coupling x + i omega) w - |w|² w) dt + sigma dW`, and the tipped leader can carry it through its Hopf
    bifurcation onto an oscillation of radius `sqrt(mu + coupling x)`. The follower can show rising variance and
    autocorrelation before the leader tips even when it stays far from its own tipping point.

    Args:
        follower:
            The normal form of the following element, `'fold'` or `'hopf'`.

        coupling:
            The strength of the linear forcing of the follower by the leader.

        rate:
            The rate at which the control of the leader rises.

        sigma:
            The amplitude of the noise of both elements.

        phi0:
            The initial control of the leader.

        phi1:
            The final control of the leader.

        c:
            The control of a fold follower, without the coupling.

        mu:
            The distance of a Hopf follower from its bifurcation, without the coupling.

        omega:
            The frequency of a Hopf follower.

    Returns:
        system:
            The system, with the leader `leader.x` and the follower `follower.y`, and also `follower.z` for the Hopf
            follower.

    Raises:
        ValidationError:
            The follower is neither `'fold'` nor `'hopf'`.

    Examples:
        >>> from skcomplex.systems.physics.climate import leading_following_cascade
        >>> leading_following_cascade('hopf').elements[1].omega
        1.0
    """
    if follower not in FOLLOWERS:
        error_msg = f'Unknown follower {follower!r}. Use one of {", ".join(FOLLOWERS)}.'
        raise ValidationError(error_msg)
    start = _lower_branch(-1.0, 1.0, phi0)
    leader = RampedFoldElement('leader', x=start, sigma=sigma, phi0=phi0, phi1=phi1, rate=rate)
    if follower == 'fold':
        follows: Element = FoldFollowerElement(
            'follower',
            y=_lower_branch(-1.0, 1.0, c + coupling * start),
            c=c,
            sigma=sigma,
        )
        forcing = [
            RateInteraction(
                'forcing',
                target='follower',
                state='y',
                sources=['leader'],
                rate=lambda s: coupling * s.sources[0].x,
            ),
        ]
    else:
        follows = HopfFollowerElement('follower', mu=mu, omega=omega, sigma=sigma)
        forcing = [
            RateInteraction(
                'forcing_y',
                target='follower',
                state='y',
                sources=['leader'],
                rate=lambda s: coupling * s.sources[0].x * s.target.y,
            ),
            RateInteraction(
                'forcing_z',
                target='follower',
                state='z',
                sources=['leader'],
                rate=lambda s: coupling * s.sources[0].x * s.target.z,
            ),
        ]
    return System(elements=[leader, follows], interactions=forcing)

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 system, with the states atmosphere.x, atmosphere.y, and atmosphere.z.

Examples:

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

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

@register_system(
    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 system, with the states `atmosphere.x`, `atmosphere.y`, and `atmosphere.z`.

    Examples:
        >>> from skcomplex.systems.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 system, with the values ring.x of all sites.

Examples:

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

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

@register_system(
    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 system, with the values `ring.x` of all sites.

    Examples:
        >>> from skcomplex.systems.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))

population_pollution(population=100000000.0, pollution=10000000.0, *, pollution_normal=0.25, pollution_standard=15000000.0, birth_normal=0.04, death_normal=0.028)

Return the pollution model of Bala, Arshad, and Noh, a population that degrades its environment.

The population generates pollution at pollution_normal per person per year, and the environment absorbs the pollution over a time that lengthens from 0.6 to 20 years as the pollution ratio, pollution over its standard, rises from 0 to 60. This makes the absorption saturate, so a growing population can drive the pollution up steeply. Pollution lowers births, food, and the quality of life, and raises deaths, through table multipliers, while food per capita over the subsistence food raises births and lowers deaths. The cropping intensity rises from 1.75 to 2.5 crops per year over 50 years. The model follows Section 4.11.9 of Bala et al. (2017), with its tables and parameters, and time in years.

Parameters:

Name	Type	Description	Default
population	float	The initial population.	100000000.0
pollution	float	The initial pollution.	10000000.0
pollution_normal	float	The pollution one person generates per year.	0.25
pollution_standard	float	The pollution that the pollution ratio divides by.	15000000.0
birth_normal	float	The normal birth rate per year.	0.04
death_normal	float	The normal death rate per year.	0.028

Returns:

Name	Type	Description
system	System	The system, with the stocks region.population and region.pollution, and the observables region.food_per_capita and region.quality_of_life.

Examples:

>>> from skcomplex.systems.physics.climate import population_pollution
>>> population_pollution(pollution_normal=0.5).elements[0].pollution_normal
0.5

Source code in src/skcomplex/systems/physics/climate/_pollution.py

@register_system(
    field='physics/climate',
    tags=('pollution', 'population', 'stocks and flows', 'system dynamics', 'lookup tables', 'ode'),
    references=('Bala et al. 2017',),
    defaults={'until': 50.0, 'dt': 0.25},
    observables=lambda system: [
        ObservableSpec('food_per_capita', PollutedRegionElement._food_per_capita, element='region'),
        ObservableSpec('quality_of_life', PollutedRegionElement._quality_of_life, element='region'),
    ],
)
def population_pollution(
    population: float = 1e8,
    pollution: float = 1e7,
    *,
    pollution_normal: float = 0.25,
    pollution_standard: float = 1.5e7,
    birth_normal: float = 0.04,
    death_normal: float = 0.028,
) -> System:
    """Return the pollution model of Bala, Arshad, and Noh, a population that degrades its environment.

    The population generates pollution at `pollution_normal` per person per year, and the environment absorbs the
    pollution over a time that lengthens from 0.6 to 20 years as the pollution ratio, pollution over its standard,
    rises from 0 to 60. This makes the absorption saturate, so a growing population can drive the pollution up
    steeply. Pollution lowers births, food, and the quality of life, and raises deaths, through table multipliers,
    while food per capita over the subsistence food raises births and lowers deaths. The cropping intensity rises
    from 1.75 to 2.5 crops per year over 50 years. The model follows Section 4.11.9 of Bala et al. (2017), with its
    tables and parameters, and time in years.

    Args:
        population:
            The initial population.

        pollution:
            The initial pollution.

        pollution_normal:
            The pollution one person generates per year.

        pollution_standard:
            The pollution that the pollution ratio divides by.

        birth_normal:
            The normal birth rate per year.

        death_normal:
            The normal death rate per year.

    Returns:
        system:
            The system, with the stocks `region.population` and `region.pollution`, and the observables
            `region.food_per_capita` and `region.quality_of_life`.

    Examples:
        >>> from skcomplex.systems.physics.climate import population_pollution
        >>> population_pollution(pollution_normal=0.5).elements[0].pollution_normal
        0.5
    """
    region = PollutedRegionElement(
        'region',
        population=population,
        pollution=pollution,
        pollution_normal=pollution_normal,
        pollution_standard=pollution_standard,
        birth_normal=birth_normal,
        death_normal=death_normal,
    )
    return System(elements=[region])

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 system, with the state system.x, starting on the lower branch at phi0.

Examples:

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

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

@register_system(
    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 system, with the state `system.x`, starting on the lower branch at `phi0`.

    Examples:
        >>> from skcomplex.systems.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 system, with the anomalies climate.x, climate.y, and climate.z.

Examples:

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

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

@register_system(
    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 system, with the anomalies `climate.x`, `climate.y`, and `climate.z`.

    Examples:
        >>> from skcomplex.systems.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)])

social_forest_dieback(v=0.8, x=0.5, *, Tv=30.0, Tc=1.0, alpha=1.0, delta=1.0, fmax=5.0, w=3.0, tp=10.0, tf=15.0)

Return the forest dieback model coupled with the replicator dynamics of mitigators who perceive its warming.

The forest grows as in forest_dieback, with its growth rate scaled by 0.2 + 0.4 x, so that mitigation by the share x of the population helps it. The share follows the replicator equation x' = x (1 - x) (-alpha + f(Tf) + delta (2 x - 1)), where the cost of warming f = fmax / (1 + exp(-w (Tf - Tc))) reads the warming Tf = tf / tp (T(t) - T(t - tp)) projected from the last tp years, a delay differential equation. Before the start the cover keeps its initial value. With no warming the shares 0 and 1 are stable and the share x* = 1 - (fmax / 2) / (1 + exp(w Tc)) between them is unstable, for the defaults of alpha, delta, and w. The population converges to all or no mitigators, and the forest to a positive cover or to bare soil, the green, brown, and barren outcomes, with long swings for some temperatures.

Parameters:

Name	Type	Description	Default
v	float	The initial vegetation cover, and its value before the start.	0.8
x	float	The initial share of mitigators.	0.5
Tv	float	The ambient temperature under full forest cover, in degrees Celsius.	30.0
Tc	float	The critical warming.	1.0
alpha	float	The cost of mitigation.	1.0
delta	float	The strength of the social norm.	1.0
fmax	float	The largest cost of warming.	5.0
w	float	The steepness of the cost of warming.	3.0
tp	float	The years of past warming that the projection reads.	10.0
tf	float	The years ahead of the projection.	15.0

Returns:

Name	Type	Description
system	System	The system, with the vegetation cover forest.v and the share of mitigators people.x.

Examples:

>>> from skcomplex.systems.physics.climate import social_forest_dieback
>>> social_forest_dieback(Tc=2.5).elements[1].Tc
2.5

Source code in src/skcomplex/systems/physics/climate/_forests.py

@register_system(
    field='physics/climate',
    tags=('tipping points', 'vegetation', 'social dynamics', 'replicator', 'delay', 'dde'),
    references=('Liu and Bauch 2023', 'Bury, Bauch, and Anand 2019'),
    defaults={'until': 200.0, 'dt': 0.5},
)
def social_forest_dieback(
    v: float = 0.8,
    x: float = 0.5,
    *,
    Tv: float = 30.0,
    Tc: float = 1.0,
    alpha: float = 1.0,
    delta: float = 1.0,
    fmax: float = 5.0,
    w: float = 3.0,
    tp: float = 10.0,
    tf: float = 15.0,
) -> System:
    """Return the forest dieback model coupled with the replicator dynamics of mitigators who perceive its warming.

    The forest grows as in `forest_dieback`, with its growth rate scaled by `0.2 + 0.4 x`, so that mitigation by the
    share `x` of the population helps it. The share follows the replicator equation
    `x' = x (1 - x) (-alpha + f(Tf) + delta (2 x - 1))`, where the cost of warming `f = fmax / (1 + exp(-w (Tf - Tc)))`
    reads the warming `Tf = tf / tp (T(t) - T(t - tp))` projected from the last `tp` years, a delay differential
    equation. Before the start the cover keeps its initial value. With no warming the shares 0 and 1 are stable and the
    share `x* = 1 - (fmax / 2) / (1 + exp(w Tc))` between them is unstable, for the defaults of `alpha`, `delta`, and
    `w`. The population converges to all or no mitigators, and the forest to a positive cover or to bare soil, the
    green, brown, and barren outcomes, with long swings for some temperatures.

    Args:
        v:
            The initial vegetation cover, and its value before the start.

        x:
            The initial share of mitigators.

        Tv:
            The ambient temperature under full forest cover, in degrees Celsius.

        Tc:
            The critical warming.

        alpha:
            The cost of mitigation.

        delta:
            The strength of the social norm.

        fmax:
            The largest cost of warming.

        w:
            The steepness of the cost of warming.

        tp:
            The years of past warming that the projection reads.

        tf:
            The years ahead of the projection.

    Returns:
        system:
            The system, with the vegetation cover `forest.v` and the share of mitigators `people.x`.

    Examples:
        >>> from skcomplex.systems.physics.climate import social_forest_dieback
        >>> social_forest_dieback(Tc=2.5).elements[1].Tc
        2.5
    """
    forest = ForestCoverElement('forest', v=v, Tv=Tv)
    people = MitigatorsElement('people', x=x, alpha=alpha, delta=delta, fmax=fmax, w=w, Tc=Tc, tp=tp, tf=tf)
    growth = RateInteraction(
        'growth',
        target='forest',
        state='v',
        sources=['people'],
        rate=lambda s: _logistic(s.target, 0.2 + 0.4 * s.sources[0].x),
    )
    mitigation = RateInteraction('mitigation', target='people', state='x', sources=['forest'], rate=_replicator)
    return System(elements=[forest, people], interactions=[growth, mitigation])

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 system, with the temperatures soil.x and soil.y.

Examples:

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

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

@register_system(
    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 system, with the temperatures `soil.x` and `soil.y`.

    Examples:
        >>> from skcomplex.systems.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 system, with the temperature earth.T, the control earth.mu, and the ramp progress earth.lam.

Raises:

Type	Description
ValidationError	The type of tipping is unknown.

Examples:

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

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

@register_system(
    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 system, with the temperature `earth.T`, the control `earth.mu`, and the ramp progress `earth.lam`.

    Raises:
        ValidationError:
            The type of tipping is unknown.

    Examples:
        >>> from skcomplex.systems.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 system, with the states elements.x, one value per node.

Examples:

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

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

@register_system(
    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 system, with the states `elements.x`, one value per node.

    Examples:
        >>> import networkx as nx
        >>> from skcomplex.systems.physics.climate import tipping_cascade_network
        >>> tipping_cascade_network(nx.cycle_graph(5), c=[1.2, 0, 0, 0, 0]).elements[0].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 system, with the states climate.x and climate.y.

Examples:

>>> from skcomplex.systems.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/systems/physics/climate/_ice_ages.py

@register_system(
    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 system, with the states `climate.x` and `climate.y`.

    Examples:
        >>> from skcomplex.systems.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)],
    )