Stochastic energy balance with bifurcation-, noise-, and rate-induced tipping
The Sutera-Fraedrich stochastic energy balance set up for bifurcation-, noise-, or rate-induced tipping.
- Name:
sutera_fraedrich_tipping - Field:
physics/climate - Formalism: stochastic differential equations
- References: Fang, Wang, Kurths, and Fan 2025; Sutera 1981; Fraedrich 1978; Ashwin et al. 2012
Factory
@register_model(
field='physics/climate',
tags=('tipping points', 'energy balance', 'sde', 'B-tipping', 'N-tipping', 'R-tipping'),
references=('Fang, Wang, Kurths, and Fan 2025', 'Sutera 1981', 'Fraedrich 1978', 'Ashwin et al. 2012'),
defaults={'until': 300.0, 'dt': 0.5, 'solver': EulerMaruyamaSolver(dt=0.05)},
)
def sutera_fraedrich_tipping(tipping: str = 'noise', *, rho: float = 0.18, nu: float | None = None) -> System:
"""Return the Sutera-Fraedrich stochastic energy balance set up for bifurcation-, noise-, or rate-induced tipping.
The temperature follows `dT = c0 [-T^4 + c1 mu (b2 T^2 - a2)] dt + sqrt(nu) dW` with the parameters of Table 1 of
Fang et al. and starts at its warm stable equilibrium. In bifurcation-induced tipping, `mu` decreases at
`4e-4` per year until the warm equilibrium meets the unstable one in a fold near `mu = 0.76`. In noise-induced
tipping, the controls are fixed and noise of variance `nu = 1` carries the state over the unstable equilibrium. In
rate-induced tipping, `b2` and `a2` ramp logistically at the rate `rho` with the equilibria preserved, and the state
fails to track the moving warm equilibrium when `rho` exceeds about 0.175. A tipped state leaves the warm branch and
falls to the other stable root of the quartic.
Args:
tipping:
The type of tipping, `'bifurcation'`, `'noise'`, or `'rate'`.
rho:
The rate of the logistic ramp of rate-induced tipping, per year.
nu:
The variance of the noise per year, or `None` for the value of the type, 1 for noise-induced tipping and
0 otherwise.
Returns:
system:
The model, with the temperature `earth.T`, the control `earth.mu`, and the progress `earth.lam` of the ramp.
Raises:
ValidationError:
The type of tipping is unknown.
Examples:
>>> from skcomplex.models.physics.climate import sutera_fraedrich_tipping
>>> round(sutera_fraedrich_tipping('rate').elements[0].T, 1)
291.6
"""
if tipping not in TIPPING_TYPES:
error_msg = f'Unknown tipping type {tipping!r}. Use one of {", ".join(TIPPING_TYPES)}.'
raise ValidationError(error_msg)
values = TIPPING_TYPES[tipping]
earth = StochasticEnergyBalanceElement(
'earth',
T=_warm_state(1.0, values['b_ini'], values['a_ini']),
lam=values['lam'],
drift=values['drift'],
rho=rho if tipping == 'rate' else 0.0,
b_ini=values['b_ini'],
b_fin=values['b_fin'],
a_ini=values['a_ini'],
nu=values['nu'] if nu is None else nu,
)
return System(elements=[earth])
Components
class StochasticEnergyBalanceElement(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.
Args:
label:
The label of the climate.
visual:
How the element is drawn, or `None` for `default_visual_`.
T:
The initial temperature in kelvin.
mu:
The initial relative solar input.
lam:
The initial progress of the ramp, between 0 and 1.
drift:
The rate at which `mu` decreases, per year.
rho:
The rate of the logistic ramp of `b2`, per year.
b_ini:
The initial value of `b2`, per square kelvin.
b_fin:
The final value of `b2`, per square kelvin.
a_ini:
The initial value of `a2`.
nu:
The variance of the noise per year, in square kelvin.
c0:
The scale of the energy balance, per cubic kelvin and year.
c1:
The scale of the absorbed radiation, in kelvin to the fourth.
Examples:
>>> from skcomplex.models.physics.climate import StochasticEnergyBalanceElement
>>> StochasticEnergyBalanceElement('earth').nu
1.0
"""
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
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the energy balance with its noise, the drift of `mu`, and the logistic ramp of `b2` and `a2`.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
b2 = s.b_ini + s.lam * (s.b_fin - s.b_ini)
a2 = s.a_ini + s.c1 * s.mu * (b2**2 - s.b_ini**2) / 4
return [
RateLaw(s.T, s.c0 * (-(s.T**4) + s.c1 * s.mu * (b2 * s.T**2 - a2))),
NoiseLaw(s.T, s.nu**0.5),
RateLaw(s.mu, -s.drift),
RateLaw(s.lam, s.rho * s.lam * (1 - s.lam)),
]
Equations
earth.T: diff(T_earth(t), t) == -1/4*(4*T_earth(t)^4 - (4*((b_fin_earth - b_ini_earth)*lam_earth(t) + b_ini_earth)*T_earth(t)^2 - (((b_fin_earth - b_ini_earth)*lam_earth(t) + b_ini_earth)^2 - b_ini_earth^2)*c1_earth*mu_earth(t) - 4*a_ini_earth)*c1_earth*mu_earth(t))*c0_earth
earth.mu: diff(mu_earth(t), t) == -drift_earth
earth.lam: diff(lam_earth(t), t) == -rho_earth*(lam_earth(t) - 1)*lam_earth(t)