Shallow Lake Regime Shift: Sediment Phosphorus, Hysteresis and the Fold Bifurcation

One equation, two stable lakes. Raise the nutrient loading until the water flips from clear to green, then find out how far you must cut it to flip back.

Live measurements

Equilibria versus loading L: solid = stable, dashed = unstable (the threshold between the two lakes). White dot = this lake now; trail = last 60 years.

Net rate dP/dt = L + recycling − losses. Where it crosses zero going down the lake is stable.

■ phosphorus P (µg/L) ■ loading L × 4 (last 80 years)

Stable states at this loading- Fold: jump to turbid above L =- Fold: jump back to clear below L =- Distance to the nearest fold- Return time to equilibrium 1/|f′|- Light reaching the bed- Plant cover (diagnostic)- Sediment recycling share of inputs-

Lake controls

Loading L

Sediment recycling r

Flushing + settling b

Weather noise σ

Time speed

L is nutrient input from the catchment (µg/L per year of lake volume). r is the maximum phosphorus the sediment can release (µg/L per year); dredging or fish removal lowers it. b is lost by outflow and settling.

Experiment: the hysteresis loop

Raises L from 0 to 40 and back in steps of 0.1 (each one relaxed to equilibrium with RK4) and records where the lake jumps. The analytic folds come from the maximum and minimum of L(P) = bP − recycling(P).

Independent numerical checks

Check	Reference	Model

The science

Shallow lakes can sit in one of two self-stabilising states: clear water with rooted plants and many fish, or turbid water full of algae. The same nutrient loading can be compatible with both, and a lake that has flipped to turbid does not return when the loading is cut back to the level at which it flipped: that is hysteresis (Scheffer et al. 1993, 2001).

Model (Carpenter, Ludwig & Brock 1999): dP/dt = L - bP + rP^q/(m^q+P^q), with P the water phosphorus (µg/L), L the loading, b the combined loss rate, and a sigmoidal recycling term (anoxic sediment releases phosphorus faster as the water gets richer). Defaults: b = 1/yr, r = 90 µg/L/yr, m = 60 µg/L, q = 4 (illustrative round numbers, not fitted to a lake). If the maximum slope of the recycling curve exceeds b there is an S-shaped equilibrium curve L(P) = bP − rP^q/(m^q+P^q): two saddle-node folds at its maximum (clear → turbid) and minimum (turbid → clear). In between, three equilibria coexist; the middle one is unstable and is the tipping threshold.

Critical slowing down: close to a fold the slope f′(P*) → 0, so the return time 1/|f′| grows without bound, and with weather noise the variance σ²/(2|f′|) of P increases. Both are widely proposed early-warning signals (Scheffer et al. 2009); here you can read the return time live and watch the noisy P wander more as L approaches the fold.

Diagnostics (empirical, illustrative): chlorophyll relaxes to 0.0724 P^1.45 (a Dillon-Rigler-type power law) with a 2-month lag, light attenuation K_d = 0.3 + 0.03 Chl (1/m), Secchi depth = 1.7/K_d, the share of surface light reaching the 3 m bed is exp(−K_d·3). Plant cover and fish visibility are drawn from this light and are consequences here; in the full theory plants feed back on turbidity, which is what makes the clear state so stable, and that feedback is only represented through the recycling term.

Simplifications, stated honestly: one well-mixed lake box, no stratification, no nitrogen, no zooplankton or fish dynamics, constant loading (except what you change), and parameters chosen for a clean demonstration. Real lakes show hysteresis loops that are rarely this sharp, and the thresholds of a real lake have to be measured, not computed from this table. Noise is additive Gaussian (Euler-Maruyama), standing in for storms and seasons. The lake is drawn at 1 unit = 1 m with a mean depth of about 2 m; plants, fish and the jetty are generic models.

Models and baked textures: Blender 5.2 (Cycles); rendering: Three.js r160.