Melcher25

Melcher25(
    var_name='bistable',
    sigma=0.2,
    gamma=1.5,
    alpha=0.0,
    b0=0.625,
    q0=-9.0,
    q1=12.0,
    tau=0.902,
    state_variables=None,
    diagnostic_variables=None,
    *args,
    **kwargs,
)

Two-equation bistable Itô SDE for stochastic Dansgaard-Oeschger transitions.

Models the meridional buoyancy gradient Δb and buoyancy flux B as a coupled stochastic system that switches between stadial (weak AMOC) and interstadial (strong AMOC) states:

d(db)/dt = -B - |q0 + q1*(db - b0)| * (db - b0)  +  sigma * dW
dB/dt    = (db + alpha*B - gamma) / tau             +  sigma * dW

Parameters

var_name : str = 'bistable'

Label for model output. Default is ‘bistable’.

sigma : float = 0.2

Noise amplitude applied to both state variables. Range [0.15, 0.4]. Default is 0.2.

gamma : float or callable or cc.core.Forcing = 1.5

External forcing proportional to atmospheric CO2. Bistable behaviour occurs near γ ≈ 1.5; range [0.8, 3.2]. If callable, must follow the signature contract in contracts/signal_model_contract.md: (t), (t, state), or (t, state, model) with the first argument named t or time. Default is 1.5.

alpha : float or callable or cc.core.Forcing = 0.0

Southern Ocean AABW coupling; controls the asymmetry between stadial and interstadial durations. Range [-1.0, 0.0]. If callable, must follow the same signature contract as gamma. Default is 0.0.

b0 : float = 0.625

Reference buoyancy gradient that sets the centre of the bistable potential well. Default is 0.625 (calibrated from NGRIP, Melcher et al. 2025). Changing this value shifts both fixed points and may destroy bistability if moved far from the calibrated value.

q0 : float = -9.0

Constant term in the nonlinear restoring force. Default is -9.0 (calibrated from NGRIP). Together with q1 it controls the depth and separation of the two potential wells; large departures from the default can eliminate one well and collapse the system to monostable behaviour.

q1 : float = 12.0

Linear coefficient in the nonlinear restoring force. Default is 12.0 (calibrated from NGRIP). Must be positive for the restoring force to be bounded; values far from the default alter the transition rates.

tau : float = 0.902

Relaxation time-scale of the buoyancy-flux equation (in model time units). Default is 0.902 (calibrated from NGRIP). Very small values make B relax almost instantaneously and can suppress oscillations; very large values slow the recovery from transitions.

Notes

The parameters b0, q0, q1, and tau are calibrated from NGRIP data (Melcher et al. 2025). Their defaults are the suggested values that reproduce realistic Dansgaard-Oeschger statistics. You may change them to explore sensitivity or to fit a different dataset, but be aware that combinations far from the defaults can suppress bistability entirely.

All seven parameters (sigma, gamma, alpha, b0, q0, q1, tau) are registered in param_values and are therefore inspectable via model.doc('parameters') and model.list('parameters'). sigma, gamma, and alpha additionally support time-varying values via callables or Forcing objects.

This model sets uses_post_history = True: state variables are derived from the full solved trajectory. After every integrate() call, the states diagnostic variable is auto-populated with the hysteresis classification, and model.stadial_threshold / model.interstadial_threshold are set to the Jacobian-derived thresholds. Use method='heun_maruyama' for strong order 1.0 accuracy (preferred over 'euler_maruyama' for this additive-noise model).

See also

cc.core.CCModel : Base class. classify_bistable_states : Reclassify a Δb signal post-hoc without re-running the SDE.

Examples

Generate a synthetic DO record and inspect states

import numpy as np
import matplotlib.pyplot as plt
from climatecritters.model_critters.melcher25 import (
    Melcher25, classify_bistable_states
)

model = Melcher25(sigma=0.2, gamma=1.5, alpha=-0.4)
output = model.integrate(
    t_span=(0, 599.88), y0=[1.0, 0.0],
    method='heun_maruyama', dt=0.012,
    kwargs={'random_seed': 42, 'si': 0.12},
)

# Δb time series as a pyleoclim Series — ready to plot or analyse
ts = output.to_pyleo('db')
ts.plot()
plt.savefig('docs/reference/figures/Melcher25_example.png',
        dpi=150, bbox_inches='tight')

# Stadial/interstadial classification is computed automatically
states = output.diagnostic_variables['states']   # 1 = cold, 0 = warm
print('stadial threshold:',      model.stadial_threshold)
print('interstadial threshold:', model.interstadial_threshold)

# Count DO events (interstadial onsets)
n_events = int(np.sum(np.diff(states) == -1))
print('number of DO events:', n_events)

Melcher25 example output

Add a time-varying forcing (e.g. a slow CO2 ramp across a glacial cycle):

import numpy as np
import matplotlib.pyplot as plt
from climatecritters.model_critters.melcher25 import Melcher25

t_arr = np.linspace(0, 599.88, 5000)
gamma_ramp = np.linspace(0.8, 3.2, 5000)   # γ increases over time

model = Melcher25(
    sigma=0.2,
    gamma=lambda t: float(np.interp(t, t_arr, gamma_ramp)),
    alpha=0.0,
)
output = model.integrate(
    t_span=(t_arr[0], t_arr[-1]), y0=[1.0, 0.0],
    method='heun_maruyama', dt=0.012,
    kwargs={'random_seed': 0, 'si': 0.12},
)
ts = output.to_pyleo('db')
ts.plot()
plt.savefig('docs/reference/figures/Melcher25_ramp_example.png',
        dpi=150, bbox_inches='tight')

Melcher25_ramp example output

References

Melcher et al. (2025), Clim. Past, 21, 115-132. https://cp.copernicus.org/articles/21/115/2025/

Methods

Name	Description
compute_stability_thresholds	Compute stadial and interstadial thresholds via Jacobian stability analysis.
dydt	Evaluate the deterministic drift of the bistable SDE at time t and state y.
populate_diagnostics_from_history	Classify states and compute thresholds from the full solved trajectory.
sde_noise	Return the diffusion coefficient vector at time t and state y.

compute_stability_thresholds

Melcher25.compute_stability_thresholds(alpha)

Compute stadial and interstadial thresholds via Jacobian stability analysis.

Parameters

alpha : float or array - like

Southern Ocean AABW coupling parameter. The mean is taken if an array is passed (the stadial threshold shifts by less than 0.05 units across the full alpha range, so the mean is a good approximation).

Returns

stadial_threshold : float

db drops below this → stadial state (1).

interstadial_threshold : float

db rises above this → interstadial state (0).

Notes

The lower threshold is the bifurcation point where trace(J) = 0 (Eq. 7, Melcher et al. 2025):

b_low = b0 + (-q0 - alpha/tau) / (2*q1)

The upper threshold is the non-differentiable point of the drift:

b_nd = b0 - q0/q1

which is independent of alpha.

dydt

Melcher25.dydt(t, y)

Evaluate the deterministic drift of the bistable SDE at time t and state y.

Called by the solver at each timestep. Because this model sets uses_post_history = True, state is not accumulated here as a side effect.

Parameters

t : float

Current time.

y : array - like

Current state vector [db, B].

Returns

dydt : list of float

Time-derivatives [d(db)/dt, dB/dt].

populate_diagnostics_from_history

Melcher25.populate_diagnostics_from_history(time, history)

Classify states and compute thresholds from the full solved trajectory.

Called automatically by CCModel.post_integrate after every integrate() call. Populates diagnostic_variables['states'] with the hysteresis classification and sets self.stadial_threshold / self.interstadial_threshold.

Parameters

time : array - like

Time axis of the solution.

history : ndarray of shape (n_times, 2)

State trajectory [db, B] at each time step.

sde_noise

Melcher25.sde_noise(t, y)

Return the diffusion coefficient vector at time t and state y.

Called by SDE solvers (euler_maruyama, heun_maruyama, milstein) to scale the Wiener increments. Noise is additive (state-independent), which means heun_maruyama achieves strong order 1.0 for this model.

Parameters

t : float

Current time.

y : array - like

Current state vector [db, B].

Returns

diffusion : ndarray of shape (2,)

Diffusion coefficient for each state variable.