11institutetext: Faculty of Information Technology, Brno University of
Technology, Czech Republic

# Static Deadlock Detection in Low-Level C Code††thanks: The work was
supported by the project 20-07487S of the Czech Science Foundation and the
Brno Ph.D. Talent Scholarship Programme.

Dominik Harmim  Vladimír Marcin  Lucie Svobodová  Tomáš Vojnar

###### Abstract

We present a novel scalable deadlock analyser L2D2 capable of handling C code
with low-level unstructured lock manipulation. L2D2 runs along the call tree
of a program, starting from its leaves, and analyses each function just once,
without any knowledge of the call context. L2D2 builds function summaries
recording information about locks that are assumed or known to be locked or
unlocked at the entry, inside, and at the exit of functions, together with
lock dependencies, and reports warnings about possible deadlocks when cycles
in the lock dependencies are detected. We implemented L2D2 as a plugin of the
Facebook/Meta Infer framework and report results of experiments on a large
body of C as well as C++ code illustrating the effectiveness and efficiency of
L2D2.

##  1 Introduction

Nowadays, programs often use _multi-threading_ to utilise the many processors
of current computers better. However, concurrency does bring not only speed-
ups but also a much larger space for nasty errors easy to cause but difficult
to find. The reason why finding errors in concurrent programs is particularly
hard is that concurrently running threads may _interleave_ in many different
ways, with bugs hiding in just a few of them. Such interleavings are hard to
discover by testing even if it is many times repeated.

Coverage of such rare behaviours can be improved using approaches such as
_systematic testing_ [22] and _noise-based testing_ [7, 9, 10]. Another way is
to use _extrapolating dynamic checkers_ , such as [11, 12], which can report
warnings about possible errors even if those are not seen in the testing runs,
based on spotting some of their symptoms. Unfortunately, even though such
checkers have proven quite useful in practice, they can, of course, still miss
errors. Moreover, monitoring a run of large software through such checkers may
also be quite expensive.

On the other hand, approaches based on _model checking_ , i.e., exhaustive
state-space exploration, can guarantee the discovery of all potentially
present errors — either in general or at least up to some bound, which is
usually given in the number of context switches. However, so far, the
scalability of these techniques is not sufficient to handle truly large
industrial code, even when combined with methods such as _sequentialisation_
[17, 19], which represents one of the most scalable approaches in the area.

An alternative to the above approaches, which can scale better than model
checking and can find bugs not found dynamically (though for the price of
potentially missing some errors and/or producing false alarms), is offered by
approaches based on _static analysis_ , e.g., in the form of _abstract
interpretation_ [5] or _data-flow analysis_ [15]. The former approach is
supported, e.g., in Facebook/Meta Infer — an open-source framework for
creating highly scalable, compositional, incremental, and interprocedural
static analysers based on abstract interpretation [4].

Infer provides several analysers that check for various types of bugs, such as
buffer overflows, null-dereferencing, or memory leaks. However, most
importantly, Infer is a  _framework_ for building new analysers quickly and
easily. As for _concurrency-related bugs_ , Infer provides support for finding
some forms of _data races_ and _deadlocks_ , but it is limited to _high-level_
Java and C++ programs only and fails for C programs, which use a  _lower-level
lock manipulation_ [1, 6].

In this paper, we propose a  _deadlock checker_ that fits the common
principles of analyses used in Infer and is applicable even to _C code_ with
_lower-level lock manipulation_. Our checker is called L2D2 for low-level
deadlock detector.

As is common in Infer, L2D2 computes function summaries _upwards_ along the
call tree, starting from its leaves, and analyses every function just once,
without knowing anything about its call contexts. The summaries contain
various pieces of information about locks that are assumed to be
locked/unlocked at the entry of a function, that may be locked/unlocked at the
end of the function, that may be both locked and unlocked inside a function,
as well as about lock dependencies (saying that some lock is locked while
another is still held). If L2D2 detects a loop in the lock dependencies, it
warns about possible deadlocks. L2D2 uses multiple heuristics to reduce the
number of false alarms, such as detection of locks serving as gate locks.

To show the effectiveness and efficiency of L2D2, we present experiments in
which we managed to apply it to 930 programs with 10.3 million lines of code
(MLoC) in total, out of which 8 contained known deadlocks. L2D2 rediscovered
all the deadlocks, and, out of the remaining 922 programs, it claimed 909
deadlock-free and reported false alarms for 13 of the programs only. The code
included benchmarks coming from the CProver tool derived from the Debian
GNU/Linux distribution, the code of the grep, sort, tgrep, and memcached
utilities, and the eProsima/Fast-DDS middleware.

#### Related Work

To the best of our knowledge, L2D2 is the only currently existing, publicly
available, _compositional static deadlock analyser_ for _low-level code_.
Below, we briefly discuss approaches that we consider to be the closest to it.

RacerX [8] is a top-down, non-compositional, flow-sensitive and context-
sensitive analysis for C programs based on computing so-called _lock sets_ ,
i.e., sets of currently held locks, constructing a _static lock-order graph_ ,
and reporting possible deadlocks in case of cycles in it. It employs various
heuristics to reduce false-positive reports. Some of the ideas concerning the
lock sets are similar to those used in L2D2, and some of the heuristics used
in RacerX inspired those used in L2D2.

The deadlock analyser implemented within the CProver framework [16] targets C
code with POSIX threads and uses a combination of multiple analyses to create
a context-sensitive and sound analysis. It also builds a lock-order graph and
searches for cycles to detect deadlocks. Its most costly phase is the pointer
analysis used. An experimental comparison with this tool is given in Section
4.

Starvation [2] is implemented in the Infer framework, and hence it is bottom-
up, context-insensitive, and compositional. It detects deadlocks by deriving
lock dependencies for each function and checking whether some other function
uses the locks in an inverse order. It is thus similar to L2D2, but Starvation
is limited to _high-level_ Java and C++ programs with _balanced locks_ only.
Moreover, it implements many heuristics explicitly tailored for Android Java
applications.

GoodLock [13] is a well-known _dynamic analysis_ for Java programs implemented
in Java PathFinder (JPF) [14]. As a representative of dynamic analysers, it
inherits their dependence on the concrete execution (or executions) of the
given software seen for detecting possible deadlocks. It monitors the lock
acquisition history by creating a  _dynamic lock-order graph_ , followed by
checking the graph for the existence of deadlock candidates by searching for
cycles in it. To increase chances of spotting even rarely occurring deadlocks,
not directly seen in the given execution, it makes _deadlock predictions_
based on an exponential number of permutations of a single execution. A
drawback of this approach is that it may produce a high rate of false
positives.

AirLock [3] is one of the state-of-the-art dynamic deadlock analysers. It
adopts and improves the basic approach from GoodLock by applying various
optimisations to the extracted lock-order graph. Moreover, AirLock, operating
on-the-fly, runs a polynomial-time algorithm on the lock graph to eliminate
parts without cycles, followed by running a higher-cost algorithm to detect
actual lock cycles.

##  2 Static Deadlock Detection in Low-Level Concurrent C Code

This section presents the design of the L2D2 analyser. We first introduce the
main ideas of the analysis, and then discuss it in more detail.

As already mentioned, L2D2 is designed to handle _C code_ with _low-level,
unstructured lock manipulation_. It does not start the analysis from the entry
code location as done in classical inter-procedural analyses based, e.g., on
[20]. Instead, it performs the analysis of a program function-by-function
_along the call tree_ , _starting from its leaves_. Therefore, each function
is analysed just once without any knowledge of its possible call contexts. For
each analysed function, L2D2 derives a _summary_ that consists of a  _pre-
condition_ and a _post-condition_. The summaries are then used when analysing
functions higher up in the call hierarchy. The obtained analysis is
_compositional_ on the level of functions, and when used in conjunction with
some version control system, it allows one to focus on _modified functions_
and their dependants only with no need to re-analyse the unchanged functions
(which is typically a vast majority of the code).

⬇

1void f(Lock *L3’) {

2 lock(&L4);

3 unlock(&L3’);

4 lock(&L2);

5 …

6 unlock(&L4); }

7void *t1(…) {

8 lock(&L1);

9 lock(&L3);

10 …

11 f(&L3);

12 unlock(&L1); }

13void *t2(…) {

14 lock(&L2);

15 …

16 lock(&L1); }

Listing 1: A sample low-level code causing a deadlock

L2D2 does not perform a classical _alias analysis_ , i.e., a precise analysis
for saying whether some pairs of accesses to locks may alias (such an analysis
is considered too expensive — no such sufficiently precise analysis works
compositionally and at scale). Instead, L2D2 uses _syntactic access paths_
[18], computed by the Infer framework, to represent lock objects. Access paths
represent heap locations via expressions used to access them. In particular,
an access path consists of a base variable followed by a sequence of field
selectors. According to [1], the access paths’ syntactic equality is a
reasonably efficient way to say (in an under-approximate fashion) that heap
accesses touch the same address. The mechanism is indeed successfully used,
e.g., in the production checker RacerD [1] to detect data races in real-world
programs.

We will use Listing 1 to illustrate some ideas behind L2D2. It works in two
phases. In the first phase, it computes a summary for each function by looking
for lock and unlock events (lock/unlock calls in the listing) in the function.
When a call of a user-defined function appears in the analysed function during
the analysis (like on line 11 in the listing), L2D2 uses a summary of the
function if available. Otherwise, the function is analysed on demand,
effectively analysing the code bottom-up (when a recursive call is
encountered, it is skipped). The summary is then applied to an _abstract
state_ at the call site. In the listing, the summary of f will be applied to
the abstract state of t1.

In the second phase, L2D2 looks through all computed summaries of the analysed
program and focuses on so-called _dependencies_ that are a part of the
summaries. These dependencies represent possible locking sequences of the
analysed program. The obtained set of dependencies is interpreted as a
relation. L2D2 computes the transitive closure of this relation and reports a
deadlock if some lock depends on itself in the closure. If we run L2D2 on the
code in Listing 1, it will report a potential deadlock due to the cyclic
dependency between the locks L1 and L2 that arises when the thread t1 holds L1
and waits on L2 and the thread t2 holds L2 and waits on L1.

###  2.1 Computing Function Summaries

This section outlines the structure and computation of the summaries used by
L2D2 when analysing some function $f$. Intuitively, the pre-condition
expresses what states of locks $f$ expects from its callers, and the post-
condition reflects the effect of $f$ on the locks. More precisely, the post-
condition includes the lockset and unlockset sets, holding information about
which locks _may be locked_ and _unlocked_ , resp., at the exit of $f$. The
pre-condition consists of the locked and unlocked sets, stating which locks
are _expected to be locked_ and _unlocked_ , resp., upon a call of $f$. Note
that the locked/unlocked sets are maintained but not used in the basic
algorithm introduced later in this section. They are used to detect possible
_double-locking_ /_unlocking_ , see Section 3. Next, the summary’s post-
condition contains the so-called _lock dependencies_ (deps) in the form of
pairs of locks $(\text{{{\color[rgb]{1,0,1}L2}}, {{\color[rgb]{0,.5,.5}L1}}})$
where locking of L1 was observed while L2 was locked. This exact situation can
be seen in Listing 1 on line 16.

⬇

f: PRE-CONDITION

locked={L3’}

unlocked={L2, L4}

POST-CONDITION

lockset={L2}

unlockset={L3’, L4}

wereLocked={L2, L4}

deps={(L4, L2)}

order={(L3’, L2)}

t1: PRE-CONDITION

unlocked={L1, L2, L3, L4}

POST-CONDITION

lockset={L2}

unlockset={L1, L3, L4}

wereLocked={L1, L2, L3, L4}

deps={(L1, L2), (L1, L3),

(L1, L4), (L3, L4)}

t2: PRE-CONDITION

unlocked={L1, L2}

POST-CONDITION

lockset={L1, L2}

wereLocked={L1, L2}

deps={(L2, L1)}

Listing 2: Summaries for the functions from Listing 1

Two more sets are a part of the summary’s post-condition. First, the
wereLocked set contains information on which _locks may be locked and then
again unlocked_ within $f$. This is needed to detect lock dependencies with
such locks in functions higher up in the call hierarchy. Such a situation can
be seen in Listing 1. The lock L4 is locked and then unlocked again within the
function f. In this case, the lock will not be in lockset, and we would have
no information that it was locked there. Consequently, we would not create any
lock dependencies w.r.t. this lock. However, this lock will appear in
wereLocked, so we can create dependencies with it (like the dependency
$(\text{{{\color[rgb]{0,.5,.5}L1}}, {{\color[rgb]{.75,.5,.25}L4}}})$ in the
function t1 when calling f on line 11, which could not be created otherwise).

The last sets that are a part of L2D2’s post-conditions are denoted as the
order sets. They comprise pairs of locks $(\text{{{\color[rgb]{0,1,1}L3'}},
{{\color[rgb]{1,0,1}L2}}})$ where locking of L2 was seen when L3’ was unlocked
before within the same function. Such a pair is produced, e.g., on line 4 in
Listing 1. These sets help L2D2 to better determine the order of operations in
functions. Without it, we would create, e.g., the non-existent dependency
$(\text{{{\color[rgb]{0,1,1}L3}}, {{\color[rgb]{1,0,1}L2}}})$ in the function
t1 when calling f on line 11. It should not be created because L3 is unlocked
in f on line 3 before L2 is locked on line 4. Note that the lock L3 from the
function t1 is passed to f as L3’. We resolve such situations by replacing the
function’s formal parameters with the actual ones at the concrete call site.

Listing 2 gives the summaries for the functions in Listing 1, omitting the
empty sets.

Data: lock $L$ being locked; abstract state $S$

1 def  _lock( $L$, $S$)_:

2 if  _$L\notin S.locked\cup S.unlocked$_ then

3 $S.unlocked\leftarrow S.unlocked\cup\\{L\\}$;

4

5 $S.lockset\leftarrow S.lockset\cup\\{L\\}$;

6 $S.unlockset\leftarrow S.unlockset\setminus\\{L\\}$;

7 $S.wereLocked\leftarrow S.wereLocked\cup\\{L\\}$;

8 $S.deps\leftarrow S.deps\cup(S.lockset\times\\{L\\})$;

9 $S.order\leftarrow S.order\cup(S.unlockset\times\\{L\\})$;

10

Algorithm 1 Lock acquisition

Data: lock $L$ being unlocked; abstract state $S$

1 def  _unlock( $L$, $S$)_:

2 if  _$L\notin S.locked\cup S.unlocked$_ then

3 $S.locked\leftarrow S.locked\cup\\{L\\}$;

4

5 $S.unlockset\leftarrow S.unlockset\cup\\{L\\}$;

6 $S.lockset\leftarrow S.lockset\setminus\\{L\\}$;

7

Algorithm 2 Lock release

The high-level algorithm for the summary’s computation is given in Algorithms
1–3. Algorithm 1 shows how the abstract state is updated whenever locking
occurs during the analysis. First, it updates the pre-condition by adding the
lock to the unlocked set if this locking is the first operation with that lock
in the given function $f$ (lines 2–3). Intuitively, this reflects that the
lock should be unlocked before calling $f$; otherwise, we would encounter
double-locking. Next, the lock acquisition takes place, meaning that the lock
is added to lockset and removed from unlockset (lines 4–5). Moreover, the lock
is added to wereLocked (line 6). Finally, we derive new dependencies and order
edges by considering all pairs $(L^{\prime},L)$ where $L^{\prime}$ is an
element of lockset and unlockset, resp., and $L$ is the acquired lock (lines
7–8). Algorithm 2 then updates the abstract state when some lock is released.
It is analogical to the algorithm for locking, but it does not update the
wereLocked, deps, and order sets.

Data: summary $\chi$ of a callee; abstract state $S$

1 def  _apply_summary( $\chi$, $S$)_:

2 $\chi\leftarrow\mathtt{replace\\_formals\\_with\\_actuals}(\chi)$;

3 if  _$\exists L:L\in\chi.unlocked\wedge L\notin S.unlockset$_ then
$S.unlocked\leftarrow S.unlocked\cup\\{L\\}$;

4 if  _$\exists L:L\in\chi.locked\wedge L\notin S.lockset$_ then
$S.locked\leftarrow S.locked\cup\\{L\\}$;

5 $S.lockset\leftarrow(S.lockset\cup\chi.lockset)\setminus\chi.unlockset$;

6 $S.unlockset\leftarrow(S.unlockset\setminus\chi.lockset)\cup\chi.unlockset$;

7 $S.wereLocked\leftarrow S.wereLocked\cup\chi.wereLocked$;

8 $S.deps\leftarrow
S.deps\cup((S.lockset\times\chi.wereLocked)\setminus\chi.order)$;

9

Algorithm 3 Integrating a summary of a callee

Algorithm 3 integrates a callee’s summary with the abstract state of an
analysed function. Initially, the summary is updated by replacing the formal
parameters with the actual ones (line 2). We also check that all the locks
that should be locked/unlocked before calling the callee are present in
lockset/unlockset, resp. If they are not, they must be locked/unlocked even
before the currently analysed function. Hence, we update the pre-condition
(lines 3–4). On lines 5–7, the lockset, unlockset, and wereLocked sets are
appropriately modified. At last, new dependencies between the currently held
locks and locks acquired in the callee are introduced (line 8). However, we
exclude all the dependencies from the order set to avoid adding such
$(L^{\prime},L)$ dependencies where $L^{\prime}$ was unlocked before locking
$L$ in the callee.

As L2D2 is based on abstract interpretation, we must further define the _join_
operator for combining states along _confluent program paths_ (e.g. in if
statements), the _entailment_ operator allowing the analysis to detect it has
reached a fixpoint and stop, and the _widening_ operator accelerating the
analysis of loops. Since we are interested in locking patterns along any
possible path, we define the join operator as the union of incoming states’
values for all the sets in the summaries. The entailment operator is defined
as testing for a subset on all the sets. The widening operator is made equal
to the join operator as we are working with summaries on finite and not too
large domains.

###  2.2 Reporting Deadlocks

Checking for deadlocks takes place after the summaries for all functions in
the analysed program are computed. L2D2 then merges all of the derived lock
dependencies into one set $R$. This set is interpreted as a relation, and its
transitive closure $R^{+}$ is computed. If any lock $L$ depends on itself in
the closure, i.e., $(L,L)\in R^{+}$, a potential for a deadlock has been
detected. For deadlocks using two locks, L2D2 then looks for dependencies that
cause the deadlock. In particular, it looks for a lock $L^{\prime}$ s.t.
$(L,L^{\prime})\in R^{+}\wedge(L^{\prime},L)\in R^{+}$ and reports the
dependencies (a generalisation to more locks is, of course, possible).

##  3 Increasing Analysis Accuracy

L2D2 further implements three heuristics intended to decrease the number of
possible false alarms. We now introduce the two most important (with the third
one being a simple support for recursive locks).

As _double-locking_ /_unlocking_ errors are quite rare in practice, the first
heuristic uses their detection as an indication that the analysis is over-
approximating too much. Instead of reporting such errors, L2D2 resets (some
of) the working sets. Namely, if a lock acquisition leads to double-locking,
it is assumed that L2D2 followed some non-existent path, and lockset is no
longer trustworthy. Therefore, it is erased, and the only lock left in it is
the currently acquired one as this is the only one about which we can safely
say it is locked. For that, the following statement is added to Algorithm 1:
$\mathbf{if}\ L\in S.lockset\ \mathbf{then}\ S.lockset\leftarrow\\{L\\};$.
When releasing a lock, we then check whether it may already be unlocked. If
so, lockset is erased, eliminating any dependencies that the locking error
would cause. For that, we add the following to Algorithm 2: $\mathbf{if}\ L\in
S.unlockset\ \mathbf{then}\ S.lockset\leftarrow\emptyset;$. Finally, we check
double-locking/unlocking when a function call is encountered. We ask whether
some lock that should be locked/unlocked in the callee is currently
released/held, resp. If such a lock is found, it is assumed that L2D2 used a
non-existent path to reach the function call, and so lockset is discarded, and
the lockset of the callee will be used instead. We implement this by adding
the following to Algorithm 3: $\mathbf{if}\
(S.lockset\cap\chi.unlocked\neq\emptyset)\vee$
$(S.unlockset\cap\chi.locked\neq\emptyset)\ \mathbf{then}\
S.lockset\leftarrow\chi.lockset;$.

The second heuristic used in L2D2 is the detection of so-called _gate locks_
[13], i.e., locks guarding other locks (upon which deadlocks on the nested
locks are not reported). Whenever we detect a possible deadlock — represented
by two reverse dependencies $d_{1}=(L,L^{\prime})$ and $d_{2}=(L^{\prime},L)$
— we check whether the same gate lock protects them. If so, we do not report a
deadlock. We check this by computing the intersection of the guards, i.e., all
locks locked before the program points where the dependencies $d_{1}$ and
$d_{2}$ were captured. In particular, we do not report a deadlock for
dependencies $d_{1}$ and $d_{2}$ if
$\text{{guards}}(d_{1})\cap\text{{guards}}(d_{2})\neq\emptyset$.

##  4 Experimental Evaluation

L2D2 has been implemented in OCaml as a plugin of Infer, and it is publicly
available111https://github.com/svobodovaLucie/infer. We now report on various
experiments we have performed with it. All of the experiments were run on a
machine with the AMD Ryzen 5 5500U CPU, 15 GiB of RAM, 64-bit Ubuntu 20.04.4
LTS, using Infer version v1.1.0-0e7270157.

Table 1: Results of L2D2 and CProver on non-deadlocking programs of the CProver test-suite | programs | programs | programs
---|---|---|---
checker | claimed | raising | failed
| safe | alarms | to analyse
CProver | 292 | 114 | 588
L2D2mode 1 | 906 | 11 | 77
L2D2mode 2 | 896 | 21 | 77

In our first set of experiments, we have applied L2D2 on a set of 1,002 C
programs with POSIX threads derived from a Debian GNU/Linux distribution,
originally prepared for evaluating the static deadlock analyser based on the
CProver framework proposed in [16]. The benchmark consists of 11.3 MLoC. Eight
of the programs contain a known deadlock. Like CProver, L2D2 was able to
detect all the deadlocks. The results for the remaining 994 programs are shown
in Table 1 (for L2D2, mode 1/mode 2 refer to using/not using the double-
locking-based heuristic), with some more details also in Table 2 discussed
below. We can see that, in mode 1, L2D2 produced 11 false alarms only (77
programs failed to compile since the Infer’s front-end did not support some of
the constructions used). We find this very encouraging, considering that the
CProver’s deadlock detector produced 114 false alarms. Moreover, L2D2 consumed
83 minutes only whereas CProver needed 4 hours to handle the programs it
correctly analysed, producing 453 timeouts (w.r.t. a 30-minute time limit),
and ran out of the available 24 GB of RAM in 135 cases (according to [16], the
results were obtained on Xeon X5667 at 3 GHz running Fedora 20 with 64-bit
binaries).

Table 2: Detailed results on eProsima/Fast-DDS, sort, grep, memcached, tgrep | kLoC |  alarms mode 1 |  alarms mode 2 |  dead- locks |  runtime (mm:ss)
---|---|---|---|---|---
Fast-DDS | 110 | 3 | 6 | 0 | 06:53
memcached | 31 | 6 | 7 | 0 | 00:08
sort | 7.2 | 0 | 0 | 0 | 00:02
grep | 8.7 | 0 | 0 | 0 | 00:03
tgrep | 2.4 | 0 | 0 | 0 | 00:01
CProver | 10,164 | 23 | 80 | 8 | 83:23

Table 2 provides our further experimental results. Unlike Table 1, the table
gives not only numbers of programs in which an alarm was raised, but it gives
concrete numbers of the alarms (more alarms can be raised in a single
program). Moreover, it shows how L2D2 behaved on multiple further real-life
programs. In particular, eProsima/Fast-DDS 2.6.1 is a C++ implementation of
the Data Distribution Service of the Object Management Group. For its
analysis, we replaced the C++ guard lock used, which is so far not supported
by L2D2, by a normal lock (exploiting the fact that Infer automatically adds
all needed unlock calls). Next, we analysed memcached version 1.6.10, a
distributed memory object caching system. The source code of this program was
pre-processed by Frama-C [21], and we report on the size of the pre-processed
code (likewise with all the further mentioned programs). Finally, we also
analysed grep 3.7, tgrep (a multi-threaded version of find combined with grep
by Ron Winacott), and GNU Coreutils sort 8.32. The alarms raised for Fast-DDS
are false alarms caused by some intricacy of C++ locks for which L2D2 was not
prepared. We were not able to check the status of the alarms raised for
memcached, but we consider them likely false alarms. However, we find the
results provided by L2D2 as quite encouraging since the numbers of false
alarms are low w.r.t. the number of programs and their extent, and, moreover,
we believe that there is space for further improvements (especially, but not
only for C++ locks).

#### Remark

This preprint was submitted for publication in the post-proceedings of the
18th International Conference on Computer Aided Systems Theory — Eurocast’22.
The preprint has not undergone any post-submission improvements or
corrections. The Version of Record of this contribution will be published in
LNCS, Springer.

## References

  * [1] S. Blackshear, N. Gorogiannis, P. O’Hearn, and I. Sergey.  RacerD: Compositional Static Race Detection.  Proc. of ACMPL, 2(OOPSLA):144:1–144:28, 2018.
  * [2] J. Brotherston, P. Brunet, N. Gorogiannis, and M. Kanovich.  A Compositional Deadlock Detector for Android Java.  In Proc. of ASE’21. IEEE, 2021.
  * [3] Y. Cai, R. Meng, and J. Palsberg.  Low-Overhead Deadlock Prediction.  In Proc. of ICSE’20. ACM, 2020.
  * [4] C. Calcagno, D. Distefano, J. Dubreil, D. Gabi, P. Hooimeijer, M. Luca, et al.  Moving Fast with Software Verification.  In Proc. of NFM’15, volume 9058 of LNCS. Springer, 2015.
  * [5] P. Cousot and R. Cousot.  Abstract Interpretation: A Unified Lattice Model for Static Analysis of Programs by Construction or Approx. of Fixpoints.  In Proc. of POPL’77. ACM, 1977.
  * [6] D. Distefano, M. Fähndrich, F. Logozzo, and P. O’Hearn.  Scaling Static Analyses at Facebook.  Commun. ACM, 62(8):62–70, 2019.
  * [7] O. Edelstein, E. Farchi, E. Goldin, Y. Nir, G. Ratsaby, and S. Ur.  Framework for Testing Multi-Threaded Java Programs.  Concur. Computat.: Pract. Exper., 15(3–5):485–499, 2003.
  * [8] D. Engler and K. Ashcraft.  RacerX: Effective, Static Detection of Race Conditions and Deadlocks.  In Proc. of SOSP’03. ACM, 2003.
  * [9] J. Fiedor, V. Hrubá, B. Křena, Z. Letko, S. Ur, and T. Vojnar.  Advances in Noise-Based Testing of Concurrent Software.  Softw. Test. Verif. Reliab., 25(3):272–309, 2015.
  * [10] J. Fiedor, M. Mužikovská, A. Smrčka, O. Vašíček, and T. Vojnar.  Advances in the ANaConDA Framework for Dynamic Analysis….  In Proc. of ISSTA’18. ACM, 2018.
  * [11] C. Flanagan and S. Freund.  FastTrack: Efficient and Precise Dynamic Race Detection.  In Proc. of PLDI’09. ACM, 2009.
  * [12] C. Flanagan, S. Freund, and J. Yi.  Velodrome: A Sound and Complete Dynamic Atomicity Checker for Multithreaded Programs.  In Proc. of PLDI’08. ACM, 2008.
  * [13] K. Havelund.  Using Runtime Analysis to Guide Model Checking of Java Programs.  In SPIN Model Checking and Software Verification, volume 1885 of LNCS. Springer, 2000.
  * [14] K. Havelund and T. Pressburger.  Model Checking Java Programs Using Java PathFinder.  Inter. Jour. on STTT, 2(4):366–381, 2000.
  * [15] G. Kildall.  A Unified Approach to Global Program Optimization.  In Proc. of POPL’73. ACM, 1973.
  * [16] D. Kroening, D. Poetzl, P. Schrammel, and B. Wachter.  Sound Static Deadlock Analysis for C/Pthreads.  In Proc. of ASE’16. ACM, 2016.
  * [17] A. Lal and T. Reps.  Reducing Concurrent Analysis Under a Context Bound to Sequential Analysis.  In Proc. of CAV’08. Springer, 2008.
  * [18] J. Lerch, J. Späth, E. Bodden, and M. Mezini.  Access-Path Abstraction: Scaling Field-Sensitive Data-Flow Analysis with Unbound. Access Paths.  In Proc. of ASE’15. IEEE, 2015.
  * [19] T. Nguyen, B. Fischer, S. Torre, and G. Parlato.  Lazy Sequentialization for the Safety Verif. of Unbound. Concur. Programs.  In Proc. of ATVA’16, volume 9938 of LNCS. Springer, 2016.
  * [20] T. Reps, S. Horwitz, and M. Sagiv.  Precise Interprocedural Dataflow Analysis via Graph Reachability.  In Proc. of POPL’95. ACM, 1995.
  * [21] J. Signoles, P. Cuoq, F. Kirchner, N. Kosmatov, V. Prevosto, and B. Yakobowski.  Frama-C: A Software Analysis Perspective.  Formal Aspects of Computing, 27, 2012.
  * [22] J. Wu, Y. Tang, H. Hu, H. Cui, and J. Yang.  Sound and Precise Analysis of Parallel Programs through Schedule Specialization.  In Proc. of PLDI’12. ACM, 2012.


# New Physics-Based Early Warning Signal shows AMOC is on Tipping Course
_First simulated AMOC tipping event in a state-of-the-art climate model_

René M. van Westen1∗, Michael Kliphuis1 and Henk A. Dijkstra1

1Institute for Marine and Atmospheric research Utrecht, Utrecht University,
Princetonplein 5, 3584 CC Utrecht, the Netherlands

∗To whom correspondence should be addressed; E-mail<EMAIL_ADDRESS>
> One of the most prominent climate tipping elements is the Atlantic
> Meridional Overturning Circulation (AMOC), which can potentially collapse
> due to the input of fresh water in the North Atlantic. Although AMOC
> collapses have been induced in state-of-the-art Global Climate Models (GCMs)
> by strong freshwater forcing, an AMOC tipping event has so far not been
> found. Here, we show results of the first AMOC tipping event in such a GCM,
> the Community Earth System Model, including the large climate impacts of the
> associated AMOC collapse. Using these results, we develop a new physics-
> based and observable early warning signal of AMOC tipping: the minimum of
> the AMOC induced freshwater transport at the southern boundary of the
> Atlantic. Available observational and reanalysis data indicate that the
> present-day AMOC is on route to tipping. The new early warning signal is a
> useful alternative to classical statistical ones which, when applied to our
> simulated AMOC tipping event, turn out to be sensitive to the analysed time
> interval before tipping.

## Introduction

The Atlantic Meridional Overturning Circulation (AMOC) effectively transports
heat and salt through the global ocean (?) and strongly modulates regional and
global climate. Continuous section measurements of the AMOC, available since
2004 at 26∘N from the RAPID-MOCHA array (?), have shown that the AMOC strength
has decreased by a few Sverdrups (1 Sv = $10^{6}$ m3s-1) from 2004 to 2012 and
thereafter it has strengthened (?) again. Longer time scale variability of the
AMOC strength, estimated by using sea surface temperature (SST) time series
based on ‘fingerprint’ patterns (?), indicates that the AMOC weakened by $3\pm
1$ Sv since about 1950. From proxy records, it has been suggested that the
AMOC is currently in its weakest state in over a millennium (?).

The AMOC has been labelled as one of the tipping elements in the climate
system (?, ?), indicating that it may undergo a relatively rapid change under
a slowly developing forcing. The AMOC is particularly sensitive to the ocean’s
freshwater forcing, either through the surface freshwater flux (e.g.,
precipitation) or by input of fresh water due to river run-off or ice melt
(e.g., from the Greenland Ice Sheet). Although no AMOC tipping has been found
in historical observations, there is much evidence from proxy records that
abrupt AMOC changes have occurred in the geological past during the so-called
Dansgaard-Oeschger events (?, ?, ?).

Classical early warning indicators, such as the increase in the variance
and/or the (lag-1) auto-correlation, when applied to SST-based time series,
suggest that the present-day AMOC approaches a tipping point before the end of
this century (?, ?). Apart from the fact that the SST-based AMOC fingerprints
may not represent the AMOC behaviour adequately, many (statistical)
assumptions are required to estimate the approaching AMOC tipping point (?, ?,
?, ?). Hence, there is strong need for a more physics-based, observable, and
reliable early warning indicator which characterises the AMOC tipping point.

## Results

## AMOC Collapse

To develop such an early warning indicator, we performed a targeted simulation
to find an AMOC tipping event in the Community Earth System Model (CESM).
Specifically, we use the f19_g16 configuration of CESM1.0.5, with horizontal
resolutions of 1∘ for the ocean/sea ice and 2∘ for the atmosphere/land
components, which has been used in the Coupled Model Intercomparison Project
(CMIP), phase 5. We start from a statistical equilibrium solution of a pre-
industrial control simulation and keep greenhouse gas, solar and aerosol
forcings constant to pre-industrial levels during the simulation. A quasi-
equilibrium approach (?, ?) is followed, by adding a slowly-varying freshwater
flux anomaly $F_{H}$ in the North Atlantic over the region between latitudes
20∘N and 50∘N. This freshwater flux anomaly is compensated over the rest of
the domain as shown in the inset of Figure 1a. We linearly increased the
freshwater flux forcing with a rate of $3\times 10^{-4}$ Sv yr-1 until model
year 2,200 where a maximum of $F_{H}=0.66$ Sv is reached. This is a
computationally highly expensive calculation and so cannot easily be repeated
for a suite of different GCMs.

Under increasing freshwater forcing, we find a gradual decrease (Figure 1a) in
the AMOC strength (see Methods). Natural variability dominates the AMOC
strength in the first 400 years but, after model year 800, a clear negative
trend appears due to the increasing freshwater forcing. Then, after 1,750
years of model integration, we find an abrupt AMOC collapse. The AMOC strength
is about 10 Sv in model year 1,750 and decreases to 2 Sv 100 years later
(model year 1,850) and eventually becomes slightly negative after model year
2,000. Such a transient AMOC response (model years 1,750 – 1,850) is
spectacular considering the slow change in the freshwater forcing (i.e.,
$\Delta F_{H}=0.03$ Sv). The characteristic meridional overturning circulation
and associated northward heat transport in the Atlantic Ocean have decreased
to nearly zero and by 75% (at 26∘N), respectively, after model year 2,000
(Figures 1b,c,d and S1a).

Figure 1: AMOC collapse. (a): The AMOC strength at 1,000 m and 26∘N, where the
yellow shading indicates observed ranges (?, ?). The cyan-coloured lines
indicate the magnitude of $F_{H}$. Inset: The hosing experiment where fresh
water is added to the ocean surface between 20∘N – 50∘N in the Atlantic Ocean
($+F_{H}$) and is compensated over the remaining ocean surface ($-F_{H}$). The
black sections indicate the 26∘N and 34∘S latitudes over which the AMOC
strength and $F_{\mathrm{ovS}}$ are determined, respectively. (b – d): The
AMOC streamfunction ($\Psi$) and the Atlantic meridional heat transport (MHT,
see also Figure S1) for model years 1 – 50, 1,701 – 1,750 and 2,151 – 2,200.
The contours indicate the isolines of $\Psi$ for different values.

This result differs substantially from earlier model simulations with GCMs
which have employed extremely large freshwater forcing (e.g., 1 Sv per year
over 50∘N – 70∘N (?)) or large initial salinity perturbations (?). The AMOC in
these simulations is a direct response to the very strong forcing, whereas in
our model simulations, which are more akin the simulations in Earth System
Models of Intermediate Complexity (?), the collapse is purely a response due
to internal feedbacks. This can be quantified by looking at the AMOC change
per cumulative change in North Atlantic freshwater forcing which is about
$\frac{8~{}\mathrm{Sv}}{1.5~{}\mathrm{Sv}}=5.3$ in our hosing simulations
(model years 1,750 – 1,850) and about
$\frac{18~{}\mathrm{Sv}}{50~{}\mathrm{Sv}}=0.36$ for the strong forcing of 1
Sv per year over a 50-year period for the results in (?). Also based on the
change in the AMOC per forcing change (here about 8 Sv due to a change in 0.03
Sv surface freshwater flux), it is clear that we found an AMOC tipping event
(?) in the CESM simulation which is the first one found in a state-of-the-art
GCM.

The differences in important ocean observables between the two different AMOC
states (averages over model years 2,151 – 2,200 minus years 1 – 50) are
presented in Figure S2. Figure S2a shows a cooling of the Northern Hemispheric
SSTs when the AMOC collapses, with SST differences as large as 10∘C near
Western Europe. On the contrary, the SSTs in the Southern Hemisphere increase
due to the collapse resulting in a distinct see-saw pattern between the
hemispheres (?). This pattern arises from the reduced meridional heat exchange
between the hemispheres (Figure S1). The North Atlantic upper 100 m salinities
are also strongly influenced under the AMOC collapse (Figure S2b). Note that
salinities outside of the Atlantic have increased partly due to the freshwater
flux compensation used in the set-up of the quasi-equilibrium experiment. From
the changes in the annual maximum mixed-layer depth (Figure S2c), it can be
deduced that deep convection ceases in the North Atlantic (around Greenland)
and convection increases in the Southern Ocean, in accordance with the
reversed AMOC state (Figure 1d). The weakening of the AMOC results, via
geostrophic balance, in dynamic sea-level rise in the Atlantic Ocean (Figure
S2d) and some coastal regions experience more than 70 cm of dynamic sea-level
rise.

### Climate Impacts

The SST changes due to AMOC collapse also affect the atmosphere and global
sea-ice distribution. The atmospheric responses (Figure S3) consist of a see-
saw pattern in the 2-meter surface temperature, a southward Inter-Tropical
Convergence Zone (ITCZ) shift, and the strengthening of the Hadley Cell in the
Northern Hemisphere. The stronger meridional temperature gradient over the
Northern Hemisphere amplifies the subtropical jet, while the opposite happens
in the Southern Hemisphere. During the gradual AMOC weakening over the first
1,400 model years, there are no significant trends ($p>0.05$, two-sided t-test
(?)) in the global mean surface temperature nor in the global sea-ice area.
Under the AMOC collapse, the Arctic (March) sea-ice pack extends down to 50∘N
and there is a gradual retreat of the Antarctic (September) sea-ice pack
(Figure S4). The vast expansion of the Northern Hemispheric sea-ice pack
amplifies further Northern Hemispheric cooling via the ice-albedo feedback.
These findings are qualitatively similar to those in (?) in which AMOC is
strongly weakened to 3 – 4 Sv.

The aforementioned ocean, atmosphere and sea-ice responses strongly influence
the regional climates across the globe (Figure 2). The European climate is
significantly different after the AMOC collapse, whereas for other regions
only specific months undergo significant changes. The Amazon Rainforest also
shows a drastic change in their precipitation patterns due to ITCZ shifts and
the dry season becomes the wet season and vice versa. These AMOC-induced
precipitation changes could severely disrupt the ecosystem of the Amazon
Rainforest (?, ?, ?) and potentially lead to cascading tipping (?, ?, ?). The
Northern Hemisphere shows cooler temperatures after the AMOC collapse while
the opposite is true for the Southern Hemisphere, although not all changes are
significantly different (due to large inter-annual variability).

---
Figure 2: Climograph for different regions. The climograph for 6 different
regions (spatial average over the $10^{\circ}\times 10^{\circ}$ boxes), where
the bars indicate the monthly precipitation and the curves the monthly
temperatures. The climograph is determined over model years 1 – 50 (red bars
and curves) and model years 2,151 – 2,200 (blue bars and curves). Note the
different vertical ranges for each climograph. The letters P and T in the bars
indicate significant ($p<0.05$, two-sided Welch’s t-test) monthly differences
for precipitation and temperature, respectively.

The European climate is greatly affected (Figure 3a) under the AMOC collapse.
Note that the corresponding changes occur within a relatively short period
(model years 1,750 – 1,850) and under a very small change in surface
freshwater forcing. The yearly-averaged atmospheric surface temperature trend
exceeds 1∘C per decade over a broad region in northwestern Europe and for
several European cities, temperatures are found to drop by 5∘C to 15 ∘C
(Figure 3c). The trends are even more striking when considering particular
months (Figure 3b). As an example, February temperatures for Bergen (Norway)
will drop by about 3.5∘C per decade (Figure 3d). These relatively strong
temperature trends are associated with the sea-ice albedo feedback through the
vast expansion of the Arctic sea-ice pack.

Figure 3: Surface temperature response during AMOC collapse. (a): The yearly-
averaged 2-meter surface temperature trend (model years 1,750 – 1,850), the
markers indicate non-significant trends ($p\geq 0.05$, two-sided t-test (?)).
(b): Similar to panel a, but now for the February 2-meter surface temperature
trend. The red dots indicate 5 different cities used in panels c and d. Note
the different colourbar ranges between panels a and b. (c): The temperature
difference (w.r.t. model year 1,600) for 5 different cities, including the
AMOC strength. The trends are determined over model years 1,750 – 1,850
(yellow shading) during which the AMOC strength strongly decreases. (d): The
monthly temperature trends for the 5 different cities.

### Physics-based Early Warning Indicator

From idealised ocean-climate models it has been suggested that the freshwater
transport of the AMOC at 34∘S, indicated by $F_{\mathrm{ovS}}$ (see Methods),
is an important indicator of AMOC stability (?, ?, ?, ?, ?). The reason is
that this quantity is a measure of the salt-advection feedback strength,
thought crucial in AMOC tipping. This feedback describes the amplification of
a freshwater perturbation in the North Atlantic through a weakening of the
AMOC which leads to less northward salt transport and hence amplification of
the initial freshwater perturbation (?, ?).

In the CESM results here (Figure 4a), $F_{\mathrm{ovS}}$ is positive at the
beginning of the simulation, which indicates that the AMOC exports net
salinity (w.r.t. reference salinity of 35 g kg-1) out of the Atlantic. This is
not in agreement with observations (?, ?) which is a well-known bias in CMIP
phase 3 (?), phase 5 (?) and phase 6 (?) models. In the CMIP phase 6 (CMIP6)
models, this bias is mainly due to large biases (compared to observations) in
the freshwater flux over the Indian Ocean (?).

Figure 4: Freshwater transport by the AMOC and $F_{\mathrm{ovS}}$ minimum
estimate. (a): The freshwater transport by the AMOC at 34∘S,
$F_{\mathrm{ovS}}$. The cyan-coloured lines indicate the magnitude of $F_{H}$.
Inset: Same dimensions as the main figure, but now for the freshwater
transport at 34∘S ($F_{\mathrm{ovS}}$, black curve), 60∘N ($F_{\mathrm{ovN}}$,
blue curve) and the azonal (gyre) component at 34∘S ($F_{\mathrm{azS}}$, cyan
curve). The time series are displayed as 25-year averages (to reduce the
variability of the time series). (b): The $F_{\mathrm{ovS}}$ and AMOC strength
over time, the time series are displayed as 25-year averages (to reduce the
variability of the time series). The markers indicate the 50-year average over
a particular period. The cyan-coloured curve indicates the present-day (1994 –
2020) CMIP6 regression and 1 standard deviation (?). The yellow shading in
panels a and b indicates observed ranges (?, ?, ?, ?) for $F_{\mathrm{ovS}}$
and AMOC strength. (c): The $F_{\mathrm{ovS}}$ minimum estimate procedure
(here for 50-year averages, see Methods). The time derivative of
$F_{\mathrm{ovS}}$ ($\mathrm{d}_{t}F_{\mathrm{ovS}}$, red curve) is determined
over the available and future data (up to model year 1,695), then a linear
trend is determined (model years 1,500 – 1,695) which is extrapolated to find
the zero (diamond label). (d): The $F_{\mathrm{ovS}}$ minimum estimate for
varying model years (i.e., available and future data) and different averaging
periods, the dotted lines and diamond label are similar to the ones from panel
c.

With increasing anomalous North Atlantic freshwater forcing, a greater
salinity transport into (and/or greater fresh water export out of) the
Atlantic Ocean is needed to balance the Atlantic’s salinity budget (?)
resulting in a declining $F_{\mathrm{ovS}}$. Note that the AMOC becomes weaker
under the freshwater forcing, so the salinity fields at 34∘S, which partly
control $F_{\mathrm{ovS}}$, must change accordingly to balance the Atlantic’s
salinity budget (Figure S5). The range of $F_{\mathrm{ovS}}$ and AMOC changes
in the northward overturning regime (until model year 1,750) are within
(Figure 4b) that of present-day simulations of CMIP6 models (?). But the most
important result here is that $F_{\mathrm{ovS}}$ goes through a minimum very
close to the collapse (near model year 1,750). Conceptual AMOC models (?, ?)
clearly identify such a minimum with a saddle-node bifurcation (?) which is
the AMOC tipping point in these models.

At the $F_{\mathrm{ovS}}$ minimum, an increment in the anomalous surface
freshwater forcing weakens the AMOC further and this weakening can not induce
sufficient salinity transport into (and/or fresh water transport out of) the
Atlantic Ocean (hence $F_{\mathrm{ovS}}$ increases) resulting in the
$F_{\mathrm{ovS}}$ minimum. In other words, the dominant balance between
$F_{H}$ and $F_{\mathrm{ovS}}$ changes cannot be sustained under AMOC
weakening. The tipping point is at slightly higher values of $F_{H}$ (and
hence slightly later in the simulation) than the value at the
$F_{\mathrm{ovS}}$ minimum, so the latter is a lower bound for tipping. This
is because the (azonal) gyre component at the southern boundary (34∘S,
$F_{\mathrm{azS}}$) and the freshwater transport changes at the northern
boundary (60∘N, $F_{\mathrm{ovN}}$) also contribute to the Atlantic Ocean’s
salinity budget (?). Although these responses under the surface freshwater
forcing are fairly small compared to $F_{\mathrm{ovS}}$ changes (inset in
Figure 4a) they allow the existence of a near-equilibrium AMOC state for
larger $F_{H}$ values than that at the $F_{\mathrm{ovS}}$ minimum. Eventually
the ocean circulation and salinities adjust to the collapsed state (Figure S5)
and $F_{\mathrm{ovS}}$ becomes positive again, which is in line with the
analyses from idealised climate model studies (?, ?).

The $F_{\mathrm{ovS}}$ minimum in our CESM simulation can be reasonably
estimated almost 100 years in advance (Figure 4d). We use cubic splines to
find a point in the ‘future’ where the time derivative of $F_{\mathrm{ovS}}$
(indicated by $\mathrm{d}_{t}F_{\mathrm{ovS}}$) goes through zero, marking the
$F_{\mathrm{ovS}}$ minimum (see Methods). Using 100 years of data (i.e.,
Available data in Figure 4c) is insufficient for a reliable $F_{\mathrm{ovS}}$
minimum estimate, but including more data to the analysis (i.e, Future data in
Figure 4c) improves the estimate in the CESM substantially. The result is also
robust to the averaging interval, when longer than 35 years (Figure 4d).

The historical $F_{\mathrm{ovS}}$, derived from reanalysis and assimilation
products (Figure 5a), are consistent in the sign of $F_{\mathrm{ovS}}$ when
comparing those to direct observations (?, ?). The reanalysis product mean
shows a robust and significant negative $F_{\mathrm{ovS}}$ trend (of -1.20 mSv
yr-1) over the last 40 years (Figure 5b) and its magnitude is close to the
projected CMIP6 mean trend (of -1.06 mSv yr-1, 2000 – 2100 (?)) under a high-
end climate change scenario. This multi-reanalysis mean negative trend
suggests that the AMOC is on course to tipping. Although the reanalysis
products are known to have different biases (?), this trend estimate is the
best result which can be obtained at the moment. However, these products are
too short (maximum of $\sim$100 years) at the moment to adequately estimate
the distance to the $F_{\mathrm{ovS}}$ minimum.

Figure 5: Historical $F_{\mathrm{ovS}}$ and trend. (a): The historical
$F_{\mathrm{ovS}}$ for 5 different reanalysis and assimilation products. The
horizontal resolution for the ocean component and the time span are indicated
in the legend. The $F_{\mathrm{ovS}}$ multi-reanalysis mean (inset) is the
yearly average over all available products. The yellow shading indicates
observed ranges for $F_{\mathrm{ovS}}$. (b): The present-day (1980 – 2022)
$F_{\mathrm{ovS}}$ trend. We use a sliding window (varying length of 20 to 43
years) over the available time series (1980 – 2022) to determine all
$F_{\mathrm{ovS}}$ trends and then determine the mean trend, which is
displayed for each reanalysis product and the multi-reanalysis mean. The
markers indicate the least significant (two-sided t-test (?)) trend for a
given sliding window and the trends at the maximum sliding window length (only
one trend possible) are also displayed.

By analysing SST-based proxies of the AMOC strength (?), it was suggested that
the real present-day AMOC approaches a tipping point (?, ?). Following the
same procedure as outlined in (?) on our AMOC tipping event in the CESM, we
find no consistent increase in the classical early warning indicators (i.e.,
variance and lag-1 auto-correlation) prior to the AMOC collapse (Figures
6a,b). The early warning indicators (temporarily and non-significantly) can
increase within each 300-year period indicated (shaded intervals), but over
longer time scales these quantities decline when approaching the tipping
point. Conducting the same analysis over the last 150 years (model years 1,550
– 1,700, 50 years before AMOC tipping event) indicates that the inter-annual
variability, variance ($\approx$ 0.05 – 0.07∘C2) and auto-correlation
($\approx$ 0.4 – 0.6) are very similar between the CESM and the observed SST
time series (Figure 4 in (?)). The variance is also increasing ($p<0.1$) over
this 150-year period, but the auto-correlation is decreasing ($p>0.1$). For
model years 1,427 – 1,557 an accurate tipping point estimate (following the
procedure outlined in (?)) can be made which is consistent with the timing of
the $F_{\mathrm{ovS}}$ minimum (red curves in Figure S6). This estimate is
only accurate when both the variance and auto-correlation increase. For
various 150-year periods (and when $F_{\mathrm{ovS}}<0$) either variance or
auto-correlation does not increase, resulting in inaccurate (blue curves in
Figure S6) tipping point estimates (Figure S6d). Because both variance and
autocorrelation are increasing in the SST-based time series in (?), their
estimate of the tipping point (2025 – 2095, 95%-confidence level) could be
accurate. On the other hand, our results (Figure S6) also show that the
accuracy is sensitive to the time interval analysed, due to decadal
variability in the SST time series, and that most 150-year time intervals do
not provide an accurate estimate of the tipping point.

Figure 6: Classical early warning indicators. (a): The spatially-averaged SST
over the subpolar region (blue outlined region in inset) and AMOC strength
(similar to Figure 1a) prior to the AMOC collapse. The SST time series is
averaged over the months November to May, where first the global mean SST
(November to May (?)) and then the time mean over the first 50 years are
subtracted. (b): The variance and (lag-1) auto-correlation of the SST subpolar
region for three 300-year periods. For each 300-year period, the linear trend
is removed (red lines in panel a) before determining the variance and auto-
correlation. The variance and auto-correlation are determined over a 70-year
sliding window and the significance of their trends (dashed lines) are
indicated for each period (see (?) for more details). (c & d): Similar to
panel b, but now for the (c): AMOC strength (panel a) and (d): the
$F_{\mathrm{ovS}}$ (Figure 4a).

Using the classical early warning indicators directly on the AMOC strength
time series (Figure 6c) gives similar results as the SST time series. The
trends over the 300-year periods are significant but the overall variance and
auto-correlation are decreasing towards the tipping point. Only for the
$F_{\mathrm{ovS}}$ (Figure 6d), we find a consistent and significant increase
in its variance when approaching the tipping point and this can also be
directly observed from the full time series (Figure 4a). The quantity
$F_{\mathrm{ovS}}$, in particular its minimum in combination with its variance
increase, is hence a very promising early warning signal for a (future) AMOC
collapse.

## Discussion

The results in this paper give a clear answer to a long-standing problem
around in the climate research community concerning the existence of AMOC
tipping behaviour in state-of-the-art GCMs (?, ?, ?, ?, ?, ?). Yes, it does
occur in these models! This is bad news for the climate system and humanity as
up till now one could think that AMOC tipping was only a theoretical concept
and tipping would disappear as soon as the full climate system, with all its
additional feedbacks, was considered. On the other hand, the tipping is
consistent with the wealth on paleo-climate evidence that rapid changes have
occurred in the AMOC, in particular during Dansgaard-Oeschger events (?).

The AMOC collapse dramatically changes the redistribution of heat (and salt)
and results in a cooling of the Northern Hemisphere while the Southern
Hemisphere slightly warms. Atmospheric and sea-ice feedbacks, which were not
considered in idealised climate models studies (?, ?, ?, ?), further amplify
the AMOC-induced changes resulting in a very strong and rapid cooling of the
European climate with temperatures trends of more than 3∘C per decade. In
comparison with the present-day global mean surface temperature trend (due to
climate change) of about 0.2∘C per decade, no realistic adaptation measures
can deal with such rapid temperature changes under an AMOC collapse (?, ?).

We have developed a physics-based, and observable (?, ?), early warning signal
characterising the tipping point of the AMOC: the minimum of the AMOC induced
freshwater transport at 34∘S in the Atlantic, here indicated by
$F_{\mathrm{ovS}}$. The quantity $F_{\mathrm{ovS}}$ has a strong basis in
conceptual models, where it is an indicator of the salt-advection feedback
strength. Although $F_{\mathrm{ovS}}$ has been shown to be a useful measure of
AMOC stability in GCMs (?), the minimum feature has so far not been connected
to the tipping point because an AMOC tipping event had up to now not been
found in these models. The $F_{\mathrm{ovS}}$ indicator is observable and
reanalysis products show that its value and, more importantly its trend, are
negative at the moment. This provides a strong indication that the AMOC is on
route towards tipping.

In addition, with future observations, an estimate of the distance to the AMOC
tipping point can in principle be obtained. Deploying machine learning
techniques on $F_{\mathrm{ovS}}$, in combination with its variance, could also
help in estimating the distance to AMOC tipping. We have shown that current
reanalysis products provide insufficient information to adequately estimate
this distance. Sustained future section measurements (available since 2009) at
34∘S from the SAMBA array (?, ?, ?) are therefore of utmost importance and
will become crucial to estimate the distance to an AMOC collapse.

In the CESM simulation here, AMOC tipping occurs at relatively large values of
the freshwater forcing. This is due to biases in precipitation elsewhere in
the models and mainly over the Indian Ocean (?). Hence, we needed to integrate
the CESM to rather large values of the freshwater forcing ($\sim 0.6$ Sv,
about a factor 80 times larger than the present-day melt rate of the Greenland
Ice Sheet (?)) to find the AMOC tipping event. The effect of the biases can be
seen from the value of the AMOC induced freshwater transport at 34∘S,
$F_{\mathrm{ovS}}$, which is positive at the start of the simulation. When
biases are corrected in the CESM, it is expected that the AMOC tipping is
expected to occur at smaller values of the freshwater forcing. As also the
present-day background climate state and the climate change forcing are
different than in our simulations, the real present-day AMOC may be much
closer to its tipping point than in the simulations shown here. Note that the
analysis of the early warning signal is not affected by these biases as this
analysis is independent of the background state and precise forcing details.

## Methods

Climate Model Simulations The hosing experiment was branched off from the end
(model year 2,800) of the pre-industrial CESM control simulation from Baatsen
et al. (?). Here it is shown that the upper 1,000 m of the ocean is well
equilibrated after 2,800 years of model integration. The horizontal
resolutions are 1∘ for the ocean/sea ice and 2∘ for the atmosphere land model
components in CESM. We added a freshwater flux anomaly over the North Atlantic
ocean surface hosing region between 20∘N and 50∘N (inset in Figure 1a). To
conserve salinity, this freshwater flux anomaly was compensated over the ocean
surface outside the hosing region. We linearly increased the freshwater
forcing with a trend $3\times 10^{-4}$ Sv yr-1 to a maximum of $F_{H}=0.66$ Sv
in model year 2,200.

The Freshwater Transport The freshwater transport by the overturning component
($F_{\mathrm{ovS}}$) and azonal (gyre) component ($F_{\mathrm{azS}}$) at 34∘S
are determined as:

$\displaystyle F_{\mathrm{ovS}}=F_{\mathrm{ov}}(y=34^{\circ}\mathrm{S})$
$\displaystyle=-\frac{1}{S_{0}}\int_{-H}^{0}\left[\int_{x_{W}}^{x_{E}}v^{*}\mathrm{d}x\right]\left[\langle
S\rangle-S_{0}\right]\mathrm{d}z$ (1a) $\displaystyle
F_{\mathrm{azS}}=F_{\mathrm{az}}(y=34^{\circ}\mathrm{S})$
$\displaystyle=-\frac{1}{S_{0}}\int_{-H}^{0}\int_{x_{W}}^{x_{E}}v^{\prime}S^{\prime}\mathrm{d}x\mathrm{d}z$
(1b)

where $S_{0}=35$ g kg-1 is a reference salinity. The $v^{*}$ indicates the
baroclinic velocity and is defined as $v^{*}=v-\hat{v}$, where $v$ is the
meridional velocity and $\hat{v}$ the barotropic meridional velocity (i.e.,
the section spatially-averaged meridional velocity). The $\langle S\rangle$
indicate the zonally-averaged salinity and primed quantities ($v^{\prime}$ and
$S^{\prime}$) are deviations from their respective zonal means. For more
details we refer to Jüling et al. (2021) (?).

The AMOC strength The AMOC strength is defined as the total meridional volume
transport at 26∘N over the upper 1,000 m:

$\mathrm{AMOC}(y=26^{\circ}\mathrm{N})=\int_{-1000}^{0}\int_{x_{W}}^{x_{E}}v~{}\mathrm{d}x\mathrm{d}z$
(2)

$F_{\mathrm{ovS}}$ minimum estimate To estimate the $F_{\mathrm{ovS}}$
minimum, we use cubic splines that interpolate piecewise, between so-called
knots, cubic polynomials which are twice continuously differentiable and we
impose that the second derivate is zero at the first and last knot. The knots
are determined over $n$-year averages of the $F_{\mathrm{ovS}}$ time series
for different starting years (1, 2, …, $n-1$) and result in $n$ different
cubic splines and their respective derivatives
($\mathrm{d}_{t}F_{\mathrm{ovS}}$, grey curves in Figure 4c). Using a linear
fit over the cubic spline mean derivative, we estimate where the derivative
goes through zero (i.e., the $F_{\mathrm{ovS}}$ minimum). A minimum of $n\geq
35$ (year averages and cubic splines) is required to substantially reduce the
variability of the time series and find a consistent $F_{\mathrm{ovS}}$
minimum estimate.

## Software and Model Output

The (processed) model output and analysis scripts will be made available on
Zenodo upon publication. The reanalysis and assimilation products can be
accessed through: GLORYS12V1 (https://doi.org/10.48670/moi-00021), SODA3.15.2
(http://www.soda.umd.edu) ORAS5 (https://doi.org/10.24381/cds.67e8eeb7),
ORA-20C (https://icdc.cen.uni-
hamburg.de/thredds/catalog/ftpthredds/EASYInit/ora20c/opa0/catalog.html) and
ECCO-V4r4 (https://www.ecco-group.org/products-ECCO-V4r4.htm).

## Acknowledgements

The model simulation and the analysis of all the model output was conducted on
the Dutch National Supercomputer Snellius within NWO-SURF project 17239.
R.M.v.W. and H.A.D. are funded by the European Research Council through the
ERC-AdG project TAOC (project 101055096).

## Author Contributions Statement

R.M.v.W. and H.A.D. conceived the idea for this study. M.K. performed the
model simulation with the CESM. R.M.v.W. conducted the analysis and prepared
all figures. All authors were actively involved in the interpretation of the
analysis results and the writing process.

## References

  * 1. W. E. Johns, et al., Continuous, array-based estimates of Atlantic Ocean heat transport at 26.5 N, Journal of Climate 24, 2429 (2011).
  * 2. M. A. Srokosz, H. L. Bryden, Observing the Atlantic Meridional Overturning Circulation yields a decade of inevitable surprises, Science 348, 1255575 (2015).
  * 3. B. I. Moat, et al., Pending recovery in the strength of the meridional overturning circulation at 26∘N, Ocean Science 16, 863 (2020).
  * 4. L. Caesar, S. Rahmstorf, A. Robinson, G. Feulner, V. Saba, Observed fingerprint of a weakening Atlantic Ocean overturning circulation, Nature 556, 191 (2018).
  * 5. L. Caesar, G. D. McCarthy, D. J. R. Thornalley, N. Cahill, S. Rahmstorf, Current Atlantic meridional overturning circulation weakest in last millennium, Nature Geoscience 14, 118 (2021).
  * 6. T. M. Lenton, et al., Tipping elements in the Earth’s climate system, Proceedings of the National Academy of Sciences of the United States of America 105, 1786 (2008).
  * 7. D. I. Armstrong McKay, et al., Exceeding 1.5 C global warming could trigger multiple climate tipping points, Science 377, eabn7950 (2022).
  * 8. S. Rahmstorf, Ocean circulation and climate changes during the past 120,000 years, Nature 419, 207 (2002).
  * 9. L. G. Henry, et al., North Atlantic ocean circulation and abrupt climate change during the last glaciation, Science 353, 470 (2016).
  * 10. J. Lynch-Stieglitz, The Atlantic Meridional Overturning Circulation and Abrupt Climate Change, Annual Review of Marine Science 9, 83 (2016).
  * 11. N. Boers, Observation-based early-warning signals for a collapse of the Atlantic Meridional Overturning Circulation, Nature Climate Change 11, 680 (2021).
  * 12. P. Ditlevsen, S. Ditlevsen, Warning of a forthcoming collapse of the Atlantic meridional overturning circulation, Nature Communications 14, 4254 (2023).
  * 13. P. D. Ditlevsen, S. J. Johnsen, Tipping points: Early warning and wishful thinking, Geophysical Research Letters 37, L19703 (2010).
  * 14. C. Kuehn, A mathematical framework for critical transitions: Bifurcations, fast–slow systems and stochastic dynamics, Physica D-Nonlinear Phenomena 240, 1020 (2011).
  * 15. S. Qin, C. Tang, Early-warning signals of critical transition: Effect of extrinsic noise, Physical Review E 97, 032406 (2018).
  * 16. S. Rahmstorf, et al., Thermohaline circulation hysteresis: a model intercomparison, Geophysical Research Letters L23605, 1 (2005).
  * 17. A. Hu, et al., Role of the Bering Strait on the hysteresis of the ocean conveyor belt circulation and glacial climate stability, Proceedings of the National Academy of Sciences 109, 6417 (2012).
  * 18. D. A. Smeed, et al., The North Atlantic Ocean is in a state of reduced overturning, Geophysical Research Letters 45, 1527 (2018).
  * 19. E. L. Worthington, et al., A 30-year reconstruction of the Atlantic meridional overturning circulation shows no decline, Ocean Science 17, 285 (2021).
  * 20. B. Orihuela-Pinto, M. H. England, A. S. Taschetto, Interbasin and interhemispheric impacts of a collapsed Atlantic Overturning Circulation, Nature Climate Change 12, 558 (2022).
  * 21. J. Mecking, S. Drijfhout, L. Jackson, M. Andrews, The effect of model bias on Atlantic freshwater transport and implications for AMOC bi-stability, Tellus A: Dynamic Meteorology and Oceanography 69, 1299910 (2017).
  * 22. T. F. Stocker, The Seesaw Effect, Science 282, 61 (1998).
  * 23. B. D. Santer, Statistical significance of trends and trend differences in layer-average atmospheric temperature time series, et al., Journal of Geophysical Research: Atmospheres 105, 7337 (2000).
  * 24. M. Hirota, M. Holmgren, E. H. Van Nes, M. Scheffer, Global resilience of tropical forest and savanna to critical transitions, Science 334, 232 (2011).
  * 25. N. Boers, N. Marwan, H. M. Barbosa, J. Kurths, A deforestation-induced tipping point for the South American monsoon system, Scientific reports 7, 41489 (2017).
  * 26. M. M. Dekker, A. S. Von Der Heydt, H. A. Dijkstra, Cascading transitions in the climate system, Earth System Dynamics 9, 1243 (2018).
  * 27. N. Wunderling, J. F. Donges, J. Kurths, R. Winkelmann, Interacting tipping elements increase risk of climate domino effects under global warming, Earth System Dynamics 12, 601 (2021).
  * 28. A. K. Klose, N. Wunderling, R. Winkelmann, J. F. Donges, What do we mean, ‘tipping cascade’?, Environmental Research Letters 16, 125011 (2021).
  * 29. S. Rahmstorf, On the freshwater forcing and transport of the Atlantic thermohaline circulation, Climate Dynamics 12, 799 (1996).
  * 30. P. de Vries, S. L. Weber, The Atlantic freshwater budget as a diagnostic for the existence of a stable shut down of the meridional overturning circulation, Geophysical Research Letters 32 (2005).
  * 31. H. A. Dijkstra, Characterization of the multiple equilibria regime in a global ocean model, Tellus A: Dynamic Meteorology and Oceanography 59, 695 (2007).
  * 32. S. E. Huisman, M. Den Toom, H. A. Dijkstra, S. Drijfhout, An indicator of the multiple equilibria regime of the Atlantic meridional overturning circulation, Journal of Physical Oceanography 40, 551 (2010).
  * 33. W. Weijer, et al., Stability of the Atlantic Meridional Overturning Circulation: A review and synthesis, Journal of Geophysical Research: Oceans 124, 5336 (2019).
  * 34. J. Marotzke, Abrupt climate change and thermohaline circulation: Mechanisms and Predictability, Proc. Natl. Acad. Sci. 97, 1347 (2000).
  * 35. W. R. Peltier, G. Vettoretti, Dansgaard-Oeschger oscillations predicted in a comprehensive model of glacial climate: A ”kicked” salt oscillator in the Atlantic, Geophysical Research Letters 41, 7306 (2014).
  * 36. H. L. Bryden, B. A. King, G. D. McCarthy, South Atlantic overturning circulation at 24 S, Journal of Marine Research 69, 38 (2011).
  * 37. R. M. van Westen, H. A. Dijkstra, Persistent Climate Model Biases in the Atlantic Ocean’s Freshwater Transport, EGUsphere pp. 1–29 (2023), 10.5194/egusphere-2023-1502.
  * 38. S. S. Drijfhout, S. L. Weber, E. van der Swaluw, The stability of the MOC as diagnosed from model projections for pre-industrial, present and future climates, Climate Dynamics 37, 1575 (2011).
  * 39. S. L. Garzoli, M. O. Baringer, S. Dong, R. C. Perez, Q. Yao, South Atlantic meridional fluxes, Deep Sea Research Part I: Oceanographic Research Papers 71, 21 (2013).
  * 40. P. Cessi, A simple box model of stochastically forced thermohaline flow, Journal of Physical Oceanography 24, 1911 (1994).
  * 41. H. Stommel, Thermohaline convection with two stable regimes of flow, Tellus 13, 224 (1961).
  * 42. L. C. Jackson, et al., The evolution of the North Atlantic meridional overturning circulation since 1980, Nature Reviews Earth & Environment 3, 241 (2022).
  * 43. R. J. Stouffer, et al., Investigating the causes of the response of the thermohaline circulation to past and future climate changes, Journal of Climate 19, 1365 (2006).
  * 44. J. Mecking, S. S. Drijfhout, L. C. Jackson, T. Graham, Stable AMOC off state in an eddy-permitting coupled climate model, Climate Dynamics 47, 2455 (2016).
  * 45. L. C. Jackson, R. A. Wood, Timescales of AMOC decline in response to fresh water forcing, Climate Dynamics 51, 1333 (2018).
  * 46. W. Liu, Z. Liu, E. C. Brady, Why is the AMOC monostable in coupled general circulation models?, Journal of Climate 27, 2427 (2014).
  * 47. L. Jackson, R. Wood, Hysteresis and Resilience of the AMOC in an Eddy‐Permitting GCM, Geophysical Research Letters 45, 8547 (2018).
  * 48. C. Vogel, K. O’Brien, Vulnerability and global environmental change: rhetoric and reality, Aviso: an information bulletin on global environmental change and human security, issue no. 13, March 2004 (2004).
  * 49. J. Birkmann, M. Garschagen, F. Kraas, N. Quang, Adaptive urban governance: new challenges for the second generation of urban adaptation strategies to climate change, Sustainability Science 5, 185 (2010).
  * 50. L. C. Jackson, Shutdown and recovery of the AMOC in a coupled global climate model: The role of the advective feedback, Geophysical Research Letters 40, 1182 (2013).
  * 51. C. S. Meinen, et al., Meridional overturning circulation transport variability at 34.5S during 2009–2017: Baroclinic and barotropic flows and the dueling influence of the boundaries, Geophysical Research Letters 45, 4180 (2018).
  * 52. M. Kersalé, et al., Highly variable upper and abyssal overturning cells in the South Atlantic, Science Advances 6, eaba7573 (2020).
  * 53. M. Kersalé, et al., Multi‐Year Estimates of Daily Heat Transport by the Atlantic Meridional Overturning Circulation at 34.5S, Journal of Geophysical Research: Oceans 126 (2021).
  * 54. I. Sasgen, et al., Return to rapid ice loss in Greenland and record loss in 2019 detected by the GRACE-FO satellites, Communications Earth & Environment 1, 1 (2020).
  * 55. M. Baatsen, et al., The middle to late Eocene greenhouse climate modelled using the CESM 1.0.5, Climate of the Past 16, 2573 (2020).
  * 56. A. Jüling, X. Zhang, D. Castellana, A. S. Von Der Heydt, H. A. Dijkstra, The Atlantic’s freshwater budget under climate change in the Community Earth System Model with strongly eddying oceans, Ocean Science 17, 729 (2021).
  * 57. A. Mamalakis, et al., Zonally contrasting shifts of the tropical rain belt in response to climate change, Nature Climate Change 11, 143 (2021).

Figure S1: Meridional heat transport (a & b): The oceanic and atmospheric
meridional heat transport (MHT) for model years 1 – 50 (red curves) and 2,151
– 2,200 (blue curves) and differences (black curves) between the two periods.
In panel a the meridional heat transport is shown for the global ocean (solid
curves), Atlantic Ocean (dashed curves) and Indian-Pacific Ocean (dotted
curves). Note the different vertical ranges between the two panels. Figure S2:
Oceanic response. (a): The sea surface temperature differences between the two
AMOC states (model years 2,151 – 2,200 minus 1 – 50), the markers indicate
non-significant ($p\geq 0.05$, two-sided Welch’s t-test) differences. (b – d):
Similar to panel a, but now for the (b): vertically-averaged (0 – 100 m)
salinity, (c): yearly-maximum mixed layer depth and (d): dynamic sea level.
Figure S3: Atmospheric response. (a): The 2-meter surface temperature
differences between the two AMOC states (model years 2,151 – 2,200 minus 1 –
50), the markers indicate non-significant ($p\geq 0.05$, two-sided Welch’s
t-test) differences. The red (blue) curves show positive (negative) values of
sea-level pressure differences with magnitudes of (-)1 hPa and (-)2 hPa for
the dashed and solid curves, respectively. (b – f): Similar to panel a, but
now for the (b): ITCZ location, (c): 850 hPa geopotential height (shading) and
850 hPa horizontal velocities (quivers), (d): precipitation. (e): meridional
mass streamfunction and (f): zonally-averaged zonal velocity. The ITCZ
location is determined from the monthly-averaged joint distribution of
outgoing longwave radiation and precipitation and are then converted to yearly
probabilities (?). The curves in panel e show the meridional mass
streamfunction for model years 1 – 50, where the red (blue) curves are
positive (negative) value with magnitudes of $(-)1\times 10^{10}$ kg s-1 and
$(-)2\times 10^{10}$ kg s-1 for the dashed and solid curves, respectively. The
curves in panel f show the zonal velocity for model years 1 – 50, where the
black curve is 0 m s-1 and the coloured curves are spaced every +5 m s-1.
Figure S4: Sea-ice response. (a): The Arctic sea-ice fraction for March for
model years 2,151 – 2,200. The dark blue curve shows the Arctic sea-ice edge
(i.e, the 15% sea-ice fraction isoline) for March for model years 1 – 50. (b):
Similar to panel a, but now for the Antarctic sea-ice fractions for September.
(c): The Northern Hemispheric sea-ice area for March and September, including
the 2-meter surface temperature for the global mean and Northern Hemisphere.
The sea-ice area is based on all the grid cells with sea-ice fractions larger
than 15%. The 2-meter surface temperature time series are displayed as 5-year
averages (to reduce the variability of the time series). (d): Similar to panel
c, but now for the Southern Hemispheric sea-ice area and the Southern
Hemispheric 2-meter surface temperature.  Figure S5: Water mass properties at
34∘S. (Upper row): The zonally-averaged meridional velocity at 34∘S for three
periods (model years 1 – 50, 1,701 – 1,750 and 2,151 – 2,200.). (Middle row):
The salinity along 34∘S for the three periods. (Lower row): The freshwater
transport with depth at 34∘S for the three periods. The different water masses
are derived from the velocity profile (?) and is only applicable for the
northward overturning circulation (left and middle column) and the names are:
Atlantic Surface Water (ASW), Antarctic Intermediate Water (AAIW), North
Atlantic Deep Water (NADW) and Antarctic Bottom Water (AABW). Figure S6:
Tipping point estimate from SST subpolar time series. (a): The monthly SST
subpolar time series (solid curves) for two 150-year periods. For each
150-year period, the monthly SST over the subpolar region (blue outlined
region in inset in Figure 6a) and the monthly globally-averaged SST are
determined. Next, the monthly mean (i.e., January months, February months,
etc.) is removed for both time series and then the globally-averaged SST
anomaly time series is subtracted from the SST subpolar region anomaly time
series. We did not subtract twice the globally-averaged SST (?) as there is no
climate change scenario in the CESM. (b & c): The variance and lag-1 auto-
correlation (solid curves) of the monthly SST subpolar time series over a
50-year sliding window. A linear trend is removed over each sliding window
before determining the variance and auto-correlation. (d): The tipping point
estimate, $t_{c}$, from the procedure outlined in (?) over 150-year periods
and varying starting year. An example of the tipping point estimate (auto-
correlation approaches 1) is shown by the dashed curves in panels a – c. The
shading indicates the 95%-confidence interval and is determined by varying the
length of the stationary part (50 to 80 years).