Can You Tell a System Is About to Tip Before It Does

The Number We Wanted, and the Number We Got

Critical slowing down works. That goes first. Everything after it reads like an attack, and it is not one. Take a system with an alternative state waiting for it, push the control parameter toward the point where its current state stops existing, and the system does exactly what the theory says: it gets sluggish, and its wobbles take longer to die away while its variance grows. We measured that in 3000 independent simulated runs and it is not subtle. Lag-1 autocorrelation goes from 0.331 to 0.691. Variance goes up by a factor of three. You can see it without squinting, and other people have seen it in real systems that were pushed on purpose: a whole-lake food web manipulated to a regime shift [9], and yeast populations in a dilution series driven to collapse [13].

Then we built the control. The control is a system driven at exactly the same rate, through exactly the same states, which stops ten percent short of the threshold and holds there. It never tips. Nothing happens to it, ever.

Our detector fired on it 20.53% of the time. On the runs that actually tipped, the same detector fired 21.07% of the time. Those two numbers differ by less than their combined standard error.

An AUC of 0.502 is a coin. Given one run that crossed and one run that came close and stopped, and given the entire early warning statistic computed from both, the statistic picks the right one 50.2% of the time. Flip a coin instead and you get 50.0%. The gap is 0.2 percentage points, inside the noise of our own estimate.

This article is the long version of that sentence, and we set out to write it in the opposite direction. Two of us had read the AMOC early warning papers [17] and wanted to build something that showed how the trick works. A third had read Ditlevsen and Johnsen on the Greenland ice cores, who went looking for the same signal ahead of the Dansgaard-Oeschger warmings and titled the result "early warning and wishful thinking" [8], and was less optimistic from the start. It does work. It works at the thing it does. The phrase "early warning" makes you think it does something else.

What We Simulated, and What We Did Not

The club has no laboratory. No lake to push toward a regime shift, and no ice core either. Our research is computation. Everything below came out of a Python program we wrote and ran on a school laptop, and that program simulates a mathematical model rather than anything that exists. Nothing here was observed. When we say we measured the autocorrelation, we mean we measured it in our own simulated output, the way you would measure the fairness of a die by rolling it rather than by reasoning about cubes.

The model is the cusp normal form, one state variable \(x\) pushed around by one control parameter \(\mu\):

and then, separately, we pretend to measure it badly:

Two stable states exist for \(|\mu| < \mu_c\), separated by a saddle, and we start on the lower one at \(\mu = 0\), where \(x^* = -1\) exactly, then raise \(\mu\) at a constant rate. At \(\mu_c = 2/(3\sqrt{3}) = 0.38490018\) the lower state and the saddle merge and annihilate, and the system has nowhere to go but up to the other branch near \(x = 1.1547\). That jump is the tipping point, a fold: the simplest and most-studied kind, and the kind standard early warning theory was built for [1][6], which matters for every number below.

Choosing that model is a choice, and it is a generous one. A fold is the case where critical slowing down is guaranteed to happen. Systems can also tip because the noise kicks them over a barrier that never went away, or because the forcing moved faster than the state could follow, and in those cases the recovery rate never approaches zero and there is no slowing down to detect [10]. Ashwin and colleagues sorted these into three categories and showed that only one of them has the property everybody advertises [10]. We modelled that one. Every number in this article is therefore a best case.

The Promise, Written Out

Wissel showed in 1984 that the characteristic return time of a system diverges as a threshold is approached, that the divergence has a universal form [1], and that the whole argument is short enough to follow in full, which is one reason the idea spread as fast as it did.

Solve \(-x^3 + x + \mu = 0\) for the lower branch. The trigonometric solution of the cubic gives

and the recovery rate is minus the slope of the flow at that point,

Write \(x^* = x_{\text{fold}} + d\) with \(x_{\text{fold}} = -1/\sqrt{3}\) and the whole thing collapses to \(\kappa = 2\sqrt{3}\,|d| + 3d^2\). The recovery rate vanishes linearly in the distance between the state and the fold. In terms of the control parameter it vanishes like \(\sqrt{\mu_c - \mu}\), which is the same statement in different clothing, and under a linear ramp in time it falls like the square root of the time remaining. Our output prints both limits beside the exact value at nine values of \(\mu\), and at \(\mu/\mu_c = 0.9999\) the linear law is accurate to 0.41%.

Linearise around \(x^*\) and the deviations are an Ornstein-Uhlenbeck process, which sampled at interval \(\Delta t_s\) is an AR(1) with coefficient \(\alpha = e^{-\kappa \Delta t_s}\) and stationary variance \(\sigma^2/2\kappa\). As \(\kappa \to 0\), \(\alpha \to 1\) and the variance diverges, and that divergence is the whole promise; the two indicators everybody computes, rising autocorrelation and rising variance, are those same two formulas read backwards [2][4][5].

Observation noise complicates this, and the complication deserves the space, because it turns out to be the single most damaging entry on the sensitivity list in §11. If you measure \(z = y + \eta\) with \(\eta\) independent from sample to sample, then

Every lag is deflated by the same factor. So the lag-1 autocorrelation is biased toward zero by an amount that depends on how noisy your instrument is, and you do not know that amount. But the deflation cancels in a ratio of successive autocovariances, which gives an estimator of \(\kappa\) that does not care about observation noise at all. We built that as a third indicator and it behaves exactly as advertised, and it is useless, and the reason it is useless is the most instructive small result in this study. We come back to it in §5.

Checking Ourselves Against the Algebra

The code had to earn trust first, so we made it reproduce things already known. Four checks, each printing our value beside the accepted one.

The recovery rate, measured by kicking it. Hold \(\mu\) fixed, kill the noise, displace by 0.01, fit the decay. At \(\mu = 0\) we get 2.004874 against the analytic 2.000000, high by 0.24%. That residual is not a mystery. Euler-Maruyama shrinks a displacement by \((1 - \kappa\,\mathrm{d}t)\) per step rather than \(e^{-\kappa\,\mathrm{d}t}\), so the exact prediction for the scheme we actually integrated is \(-\ln(1-\kappa\,\mathrm{d}t)/\mathrm{d}t = 2.020271\), and our fit sits 0.76% below that. The two predictions bracket us. We print the discretisation gap in every table, so nobody mistakes it for physics.

Autocorrelation with no observation noise. Stationary runs, 400 replicates, 4000 samples each. At \(\mu = 0\) we measure 0.604244 against an exact 0.603465, z = 1.26. As \(\mu\) rises the measurement drifts above the linear prediction. At \(\mu/\mu_c = 0.85\) it runs 1.43% high, and the variance 9.19% high. The drift is real, and the model explains it. The cubic makes the well asymmetric: on the side facing the fold the restoring force is weaker than the linear term claims, so excursions in that direction decay more slowly than the linearisation allows. The bias flatters the early warning signal rather than hiding it.

The same thing once the instrument is imperfect. Nine combinations of \(\mu\) and \(\sigma_{\text{obs}}\). At \(\mu = 0\), \(\sigma_{\text{obs}} = 0.02\): 0.370194 measured, 0.369408 predicted, z = 0.95. The deflation formula holds everywhere we tested it.

The observation-noise-proof recovery rate. At \(\mu/\mu_c = 0.75\) the ratio estimator returns 0.891827, 0.891827 and 0.892 at observation noise levels of 0.01, 0.02 and 0.04 respectively, against a true 0.920698. The estimator runs 3.1% low, for the same asymmetry reason, and it runs the same 3.1% low at every noise level, which is the part that matters. The bias correction works. Its standard error is 38 times larger than the lag-1 estimator's at the same settings.

One of those does not pass silently. The largest disagreement, 23.88 standard errors at \(\mu/\mu_c = 0.85\), is a genuine failure of the linearised theory and not of the code, and we know it is the nonlinearity because it grows monotonically with \(\mu\), runs in the direction the cubic predicts, and appears identically in the variance. The linear formula near a fold is optimistic by a couple of percent, and anybody quoting it should say so.

The Signal Exists

Now the ramp. Three thousand replicates, \(\mu\) raised from zero to \(\mu_c\) over 250 time units and then past it, process noise 0.05, observation noise 0.02, indicators computed in sliding windows of 100 samples stepped 10 samples at a time, each window linearly detrended first.

Before reading the trajectory we measured what the window itself does to the answer, because a 100-sample linearly detrended window is not a clean view of anything. Fitting and removing a mean and a slope removes real low-frequency power, and short-record correlation estimates are biased toward zero anyway. Our flat control is stationary at \(\mu = 0\) with an exactly known answer, so the bias is directly measurable: the windowed estimator returns 0.322181 where the process has 0.369408. The windowed answer comes to a factor of 0.8722 below the truth, which is 12.8% of the autocorrelation thrown away by the act of measuring it. Variance loses 3.3%.

With that one number in hand the trajectory is a clean validation; raw, the measured autocorrelation runs 12.8% below the quasi-static prediction over the whole run. Divide by 0.8722 and the ratio is 1.002 at the start and 0.980 at the fold. The 2% that does not come back is the quasi-static approximation failing, and it has to fail: the approximation assumes the system equilibrates faster than \(\mu\) moves, and near the fold, by construction, it does not. The theory promises very slightly more signal than the run delivers, and it does so in exactly the region where the warning is supposed to arrive.

So much for the physics. The indicator rises, it rises for the reason the textbook says, and it rises by an amount the textbook predicts to within a couple of percent once you account for your own instrument. If the question were "does critical slowing down happen", the answer is yes and we would stop here.

A word about the third indicator, the ratio estimator that corrects for observation noise. The ratio estimator is strictly better. It corrects exactly the bias that spoils the other one, and its detection power at 20% of the run before the fold is 0.164, against 0.428 for plain lag-1 autocorrelation and 0.505 for plain variance. The estimator that knows more about the problem loses to the estimator that knows less, by a factor of 2.6, because correcting the bias costs variance and the trend test is paying for variance. That trade is not special to our model and it is worth remembering the next time somebody proposes a cleverer indicator.

Notes From the Table

Thirty Percent

Here is the check that decides whether anything above is worth anything.

Take the flat control. Its \(\mu\) never moves. 3000 runs of a system sitting in a stable state with noise on it, and nothing in its future differs in any way from its past. Compute the indicator in sliding windows exactly as before, run a one-sided Mann-Kendall test for an increasing trend at nominal \(\alpha = 0.05\), and count how often it rejects.

0.3017, with a standard error of 0.0084.

Thirty standard errors from the nominal rate, and six times too large. The test does not have a 5% false alarm rate. It has a 30% false alarm rate, and it holds at roughly 30% at every lead time we evaluated, from 150 time units out to the fold itself. Point this analysis at a system going nowhere. One time in three, you publish a warning.

We want to be careful about the blame here. The indicator is fine. Critical slowing down is fine, and so is the physics in §5. The broken part is a statistical test applied to a dependent series, and the problem has a name and a literature of its own [3]. The mistake is also very easy to walk into, if you compute a rolling window and hand the output to a routine that expects independent samples.

The fix that works is boring. Stop overlapping the windows. With stride equal to window length our false alarm rate across ten evaluation points ranges from 0.0430 to 0.0707, and at the largest sample size, ten independent windows, it is 0.0520 against 0.0500, half a standard error out. The two values above 0.06 occur at six and seven windows, where the normal approximation to the Mann-Kendall null is being asked to do something it is not good at, and that is a small-sample problem rather than a dependence problem.

The other fix keeps the overlap and calibrates the threshold against the control instead of a table, and that is what we do for the rest of this article. The critical value comes out at \(z = 4.654\) at a lead of 100 rather than the tabulated 1.645, which tells you how far from valid the tabulated value was. With that threshold, the flat control fires on 4.97% of runs, and the test is exactly what it claims to be.

And if this article stopped at the end of this section it would be a useful methods note with a tidy ending. A common analysis is invalid, here is the size of the error, here are two ways to fix it, both of which work, and we would rather it had stopped there.

The Near Miss

A flat control is a soft control. It asks whether the detector can tell a system approaching a threshold from a system doing nothing at all, and of course it can, because one of them is changing and the other is not. The question a person actually has is harder: given a system that is visibly getting more sluggish, will it cross?

So we built a second control. Same noise, same sampling, same windows, same everything, driven at exactly the same rate toward exactly the same fold, and stopped at \(\mu = 0.9\mu_c\), where it holds for the rest of the record. The recovery rate falls from 2.0 to 0.556 and stays there. Of 3000 such runs, 42 crossed anyway on noise alone, which is a real fact about near-misses and we kept those runs in.

Look at what the figure is showing. Before \(t = 225\) the tipping and near-miss ensembles are not two ensembles. They are the same experiment, run twice, with different random seeds, and their indicator curves lie on top of each other because the systems generating them are identical. Nothing in the observed data before \(t = 225\) contains any information about which of the two you are looking at, because there is no such information to contain. The difference between them lives entirely in what the driver does next, and the driver is not in the data.

The numbers follow from that. At a lead time of 100, detection power on the tipping runs is 0.2107 and the alarm rate on near-miss runs is 0.2053, both with standard errors of 0.0074. At a lead of 50 it is 0.4280 against 0.4380, and at a lead of 25, 0.6103 against 0.6277. In several rows the near-miss control fires slightly more often than the real thing, which is noise, and which we report because rounding it away would be dishonest.

At lead zero, standing at the fold, having watched the two systems diverge for a full window length, the AUC climbs all the way to 0.5081. The statistic is dominated by 225 time units of shared history and barely notices the 25 units of difference.

Two of the AUC numbers are worth reading together. Against the flat control the detector is genuinely informative, 0.7247 at a lead of 100 and 0.9675 at the fold; against the near-miss it manages 0.5021 and 0.5081. Same statistic, same runs, same threshold, same 3000 replicates. What separates the two columns is the question you put to the detector, rather than the quality of the detector itself.

Kéfi and colleagues made a version of this point in 2013 in a paper we found late and wish we had found early: early warning signals also precede transitions that are not catastrophic at all [14]. A system that changes smoothly, with no fold anywhere, can produce rising autocorrelation. The indicator responds to the flattening of a basin. Whether that flattening ends in a jump or in a slow recovery is a separate question that the indicator was never designed to answer and cannot answer.

Arithmetic

The ramp lasts 250 time units. Samples arrive every 0.25 units, so there are 1000 of them before the fold. The window is 100 samples, the stride 10. 91 windows end at or before the fold.

At a lead time of 100, the test uses the first 51. Mann-Kendall on 51 points compares \(51 \times 50 / 2 = 1275\) pairs. Under independence the statistic has variance \(n(n-1)(2n+5)/18 = 51 \times 50 \times 107/18 = 15158.3\), so a standard deviation of 123.1, and the nominal one-sided 5% critical value is \(1.645 \times 123.1 = 202.5\) concordant pairs in excess. The calibrated critical value, taken from the control, is \(4.654 \times 123.1 = 573\). The tabulated threshold is 2.8 times too small.

Process standard deviation at the start of the run is \(\sqrt{0.000631} = 0.0251\). Observation noise is 0.02. The deflation factor is \(0.000631 / (0.000631 + 0.0004) = 0.612\). Predicted lag-1 autocorrelation: \(0.603465 \times 0.612 = 0.3694\). We measured 0.370194.

Replicate counts. Stationary validation, \(5 \times 400\) plus \(9 \times 400\), is 5600. The baseline and its two controls are \(3 \times 3000\), which is 9000. Four further configurations, each with a run and two controls. \(4 \times 3 \times 2000\), or 24,000. Total 38,600.

Tipping times. Mean 261.52 against a deterministic fold at 250.00, standard deviation 7.55. First percentile 242.75. One run in a hundred tips 7.25 units early. 6.90% tip before the deterministic fold.

Detection power at 20% of the run before the fold, by configuration, in order: 0.7535 at the slowest ramp, 0.5143 with no observation noise, 0.4879 at the loudest process noise, 0.4280 at baseline, 0.3940 at window 200, 0.3887 at window 50, 0.2657 at observation noise 0.04, 0.2495 at the fastest ramp, 0.2160 at the quietest process noise, 0.1243 at observation noise 0.08. Near-miss alarm rates run 0.7570, 0.5323, 0.5000, 0.4380, 0.3773, 0.4000, 0.2730, 0.2165, 0.2050, 0.1200.

Mean absolute difference between those two lists: 0.0127. Largest: 0.0330.

Runtime 192.2 seconds.

The Strongest Objection to All of This

M. put it at meeting six and it took us a while to answer properly. The objection is that the near-miss result is a tautology dressed as a finding.

Stated fully: before \(t = 225\) the two ensembles are generated by identical stochastic processes, so of course no statistic can distinguish them, and no future statistic will either. You have constructed two things that are the same, expressed surprise that they look the same, and what you have shown is that statistics cannot see the future, which nobody claimed and which the early warning literature certainly does not claim. What that literature claims is that the indicator estimates resilience, our own §5 shows it estimates resilience rather well, and resilience is a property of the system now. Whether the forcing continues is a question about the forcing, which is a matter for climate projections and policy, and not something an autocorrelation was ever supposed to deliver.

All of that is correct. We accept every clause of it.

Here is why we think the result is still worth publishing: the tautology is only visible once you have built the near-miss control, and building it is not standard practice. Read the applied early warning literature and the validation you are usually offered is a comparison against a surrogate or a null model that does not approach anything, which is our flat control, against which the detector scores an AUC of 0.7247 and looks respectable. The near-miss control is the one that converts "this statistic is informative" into "this statistic is informative about the question you asked", and the two are different claims that the word "warning" runs together.

The objection also has a practical edge. If you have decided in advance that the forcing will continue, then the indicator does carry real information: it tells you roughly how close you are, and §5 says it does that to within a couple of percent. But if you already know the forcing will continue, you can often estimate the distance to the threshold from the forcing directly, which is cheaper and does not need a noisy time series. The regime where the indicator is most valuable, where you genuinely do not know what the driver will do, is precisely the regime where our AUC is 0.502.

The counter-counter-argument, which we also accept, is that real systems are not driven by a knob somebody can stop: a lake receiving phosphorus, or an ice sheet receiving heat, has forcing with its own momentum, and a study of that system would be right to assume the driver continues. Boettiger and Hastings made a related point about the base rate: if tipping is rare, a test with decent power and a small false positive rate still produces mostly false positives [11]. Our contribution to that argument is a number for how small the false positive rate can be made when the alternative hypothesis is a near-miss rather than nothing. The answer: 0.2053 at a lead of 100, against a power of 0.2107.

One more thing M. was right about. We tested one fixed near-miss, stopping at 90% of the fold. A control that stops at 50% would be much easier to separate, and a real forecasting problem is a continuum of possible stopping points rather than the two cases we ran. Weighting that continuum requires a prior over what drivers do, which is exactly the thing the indicator was supposed to save you from needing.

Where a Different Choice Would Have Changed the Answer

Six places, worst first. Anybody transporting these numbers to a different field should be warned that the same indicators have been derived for epidemic thresholds, where the transition is a transcritical bifurcation rather than a fold and the algebra behind the slowing down is a different algebra [15].

The near-miss stopping point. Ours stops at \(0.9\mu_c\), and that single choice sets the headline: at \(0.5\mu_c\) the two ensembles diverge at \(t = 125\) instead of \(t = 225\), and half the record would carry real information. We chose 0.9 because a system that is not close to its threshold does not generate an early warning signal in the first place, so the interesting controls are the close ones. A reader who thinks 0.9 is unfairly hard has a point, and the interactive model lets you move it.

Detrending. We detrend each window linearly. Gaussian-kernel detrending of the whole series is the more common choice in the ecological literature [12] and removes different frequencies, which would change the window bias we measured as 0.8722 and could change the trend test's power in either direction. We did not test it and we do not know the sign.

The trend statistic. Mann-Kendall on the indicator series is the standard, and it is the thing that produced the 30% false alarm rate. A test built for autocorrelated series, or a likelihood-based approach fitted to the whole trajectory rather than to a summary of it, would behave differently against the flat control, and none of it would help against the near-miss, because the near-miss problem is an information problem rather than a statistical one.

Driving speed. The slowest ramp we ran gives 0.7535 power at 20% of the run out against 0.2495 for the fastest, a factor of three. Anybody quoting a warning time from a study like this is quoting a number that depends on how fast their particular system was being pushed, and there is no dimensionless form in which those numbers are comparable across systems.

The type of tipping. We simulated a fold, where slowing down is guaranteed. Rate-induced and noise-induced transitions produce no slowing down at all [10]. Structured ecological models can collapse with no warning whatsoever [7]. A study that sampled across tipping types, rather than assuming a fold the way we did, would report a lower detection rate than ours, never a higher one.

One variable. Our system has a single state. Hastings and Wysham showed that multi-species systems can shift without any of the classical precursors [7], and a fold in a high-dimensional system can be approached along a direction that the variable you happen to measure does not see.

What would not have changed the answer is the number of replicates: Figure 5A shows all three headline probabilities flat from about 500 replicates onward, and we ran 3000.

Where we would go next, if the club had another term. The honest missing piece is the case where the driver is itself uncertain and modelled, so that the question becomes a proper forecast with a probability attached rather than a yes or no. The useful version of this question takes that form, and it needs a prior over drivers, and priors over drivers are where the disagreement in the real debate actually lives. Dakos and colleagues reviewed what resilience indicators can and cannot deliver a decade after the field started [16], and the near-miss control is the piece of that argument we would want to put in front of any detector that comes next, learned or hand-built.

Readers who want the same style of question applied to the reverse problem, how many null observations it takes before absence means something, should see the club's occupancy detection study in Issue 1. The bistable system we are modelling here is the same object that gives the Sturtian glaciation its hysteresis loop in our Snowball Earth piece, seen from the other side: that article asks what a bistable system does when it flips, and this one asks whether you could have known it was about to.

Reproducing this

One file, one command. Python 3.12 and numpy, nothing else.

python threshold-detection.py > threshold-detection-output.txt

Expected runtime is about three and a half minutes on a modern laptop; ours took 192.2 s with numpy 2.4.2 on Python 3.12.3. Peak memory is modest, set by the slowest ramp cell. Its sample array is 2000 × 2400. The master seed is 20251217, hard-coded at the top of the file, and every batch draws an independent stream from it through numpy's SeedSequence.spawn, so the output is deterministic to the last digit. If your numbers differ from ours by more than the printed standard errors, something is wrong and we want to hear about it. The interactive model runs a smaller version of the same experiment in your browser, including the near-miss control, so you can watch the two systems fail to separate in real time.

References