Transition Probabilities of Noise-induced Transitions of the Atlantic Ocean Circulation

The Atlantic Meridional Overturning Circulation (AMOC) is considered to be a tipping element of the climate system. As it cannot be excluded that the AMOC is in a multiple regime, transitions can occur due to atmospheric noise between the present-day state and a weaker AMOC state. For the first time, we here determine estimates of the transition probability of noise-induced transitions of the AMOC, within a certain time period, using a methodology from large deviation theory. We find that there are two types of transitions, with a partial or full collapse of the AMOC, having different transition probabilities. For the present-day state, we estimate the transition probability of the partial collapse over the next 100 years to be about 15%, with a high sensitivity of this probability to the surface freshwater noise amplitude.


Section A: Box model
The equations governing the evolution of the system read: d (V t S t ) dt = q S (θ (q S )S ts + θ (−q S )S t ) + q U S d − θ (q N )q N S t + r S (S ts − S t ) + r N (S n − S t ) + 2E s S 0 , d (V ts S ts ) dt = q Ek S s − q e S ts − q S (θ (q S )S ts + θ (−q S )S t ) + r S (S t − S ts ), where the function θ (x) is a step function, which is 1 for a positive argument and 0 otherwise. Through this function, we can represent the different circulation given by the different sign of q S and by the possibility that the downwelling (q N ) is zero. The last equation expresses the conservation of salinity in the basin. The transports depend, in turn, on the variables, via the following relations 1 : where the density of the generic box i is defined as The equation for the downwelling, introduced by Cimatoribus et al. 1 , is based on the northward geostrophic transport above the thermocline in the high latitudes. Thus, q N is determined by the density contrasts between the northern part and southern part of the Atlantic Ocean. The relation was found in several ocean GCM studies 2-4 , even with near eddy-resolving resolution. For instance, Fig. 8 in 2 shows the linear dependence of the NADW (North Atlantic Deep Water) flow on the density difference between two latitudinal strips (50-55 • N and 35-40 • S) at mid-depth (750 m). Moreover, this linear dependence was successfully tested in 5 .
Also the volume of the pycnocline, as well as the one of the deep box, depend on the state of the system, in particular on the variable D: The reference parameter values, together with their descriptions, are shown in table S1.

3/11
Section B: Details of TAMS Consider a system of SDAEs (Stochastic Differential Algebraic Equations), and assume that the corresponding deterministic system has two steady states, A and B, for a certain choice of the parameters. Several methods are available in order to find transition probabilities between A and B. In case the corresponding deterministic system can be described in terms of a potential, the transition rate for transitions induced by white additive noise is given by the Eyring-Kramers formula 6,7 . This result, later extended to non-gradient systems 8 , is derived under very strict assumptions on the noise, which are often not valid in real physical systems. Furthermore, the quantity we are interested in is the probability that the system undergoes a transition within a certain time, which is not always related to the transition rate, especially if the noise that forces the system is large. The Trajectory-Adaptive Multilevel Splitting algorithm 9 is based on the idea of simulating a large ensemble of trajectories, and discarding the ones that do not reach B and splitting (or branching) trajectories that are closer to B. As a consequence, the probability that a trajectory reaches B keeps increasing, which is why this method is more efficient than a brute-force method. At each time step, the trajectories are ranked according to a so-called reaction coordinate. The ones with the lowest value of the reaction coordinate are discarded, while new ones are generated, by randomly branching other trajectories. A certain weight w i , which is related to the number of discarded trajectories at each time iteration i, is a measure of the probability of a trajectory to reach the next step. An unbiased estimator of the transition probabilityp can be obtained multiplying these weights: where k is the number of iterations of the algorithm, N trans is the number of trajectories that eventually reached the destination equilibrium and N is the total number of trajectories in the ensemble. In our results in section 3 of the paper, we used 1000 trajectories for each repetition of the algorithm, which ensures an adequate precision of the results. A step-by-step description of the algorithm is given in 9 with the details on the statistical properties of the estimator of the probability presented in 10 .
The reaction coordinate gives a measure of how close to the state B -and far from the state A -each trajectory goes. Depending on the event one wants to study, an appropriate function has to be chosen.
For an F-type transition, we define the reaction coordinate in TAMS as a linear function of the downwelling (q N ), which assumes the value 0 when the system is in the 'on' state and 1 when q N = 0, hence where x represents the state vector of the system, x on the 'on' state, and q N (x) the value of the downwelling.
For an S-type transition, i.e., a full transition to the 'off' state, the reaction coordinate in TAMS is taken as where η = x u − x on 2 / x off − x on 2 is the normalized distance between the unstable steady state x u and the stable steady state x on .
As an example of the use of the reaction coordinates, consider the system of SDEs is the state vector, the constants σ 1 and σ 2 represent the amplitudes of the noise and W 1 ,W 2 two independent Wiener processes. The deterministic part of the equations describes a gradient system, where the double-well potential reads Among the three steady states of the system, x A = (−1, 0) and x B = (1, 0) are stable, while x C = (0, 0) is unstable. In order to find transition probabilities between the two equilibria with the TAMS algorithm, within a certain time T, the reaction coordinate is chosen according to equation (S7), where x on , x off and x u correspond respectively to x A , x B and x C . An example of iteration of the algorithm, based on the ranking given by the reaction coordinate, is shown in Fig. S1. The process is iterated, using a sufficient amount of trajectories (usually at least 1000), for a certain number of times, and the transition probability is calculated according to equation (S5). Figure S1. Example iteration of TAMS for the double-well gradient system (equations (S8)), using an ensemble of three trajectories. The trajectories (yellow, blue and purple) are initialised at x A , simulated for a certain amount of time T, and then ranked according to the reaction coordinate (S7). As the purple trajectory has the lowest maximum value of the reaction coordinate (Q 2 ), it is discarded and a new trajectory (green) is generated, branching randomly from one of the remaining trajectories (in this case, the blue one). This new trajectory has a new maximum value of the reaction coordinate (Q 2 ) which is, in this case, larger than Q 2 .

5/11
Section C: Additional Bifurcation diagrams The amplitude of the noise forcing our system was estimated from observational data relative to the freshwater forcing in the Atlantic Ocean. The values the P-E (precipitation minus evaporation) were year-averaged and integrated over two basins, corresponding to the northern and southern boxes in the model. The time series thus obtained were summed and subtracted from each other, resulting respectively in the symmetric and antisymmetric component of the freshwater forcing, namely E s and E a in the model (see Fig. S3). The value of f σ is then computed as the ratio between the standard deviation and the mean value of E a . Such value is representative of the interannual time scale variability of the forcing. If more high-frequency noise is used, ocean mixed layer processes (not represented in the box model) will integrate this noise to give interannual variability in the surface salinity affecting the north-south density gradient. Hence, in our box model these interannual variations are an adequate description of unresolved processes. As a first order approximation of the process, we assumed the noise to be white and normally distributed. In 11 , a similar choice was made, given that the decorrelation time of the stochastic freshwater forcing, computed by their model, is one year. We used f σ as the lower bound of the noise amplitude, as, if smaller time scales are considered, the associated variability increases.  Figure S3. From the observed yearly averaged P-E in the Atlantic Ocean, from ERA-Interim Archive at ECMWF 12 , its antisymmetric part is shown, respect to the North (50 • N -70 • N and South (South of 40 • S) Atlantic. The black horizontal lines indicate the mean value (solid line) and one standard deviation confidence interval (dashed lines). The value of f σ , that is the ratio between the standard deviation and the mean value, is about 0.1.

Section E: Sensitivity Analysis for F-type transitions
In order to test the robustness of our results, we computed transition probabilities for slightly different versions of our model, characterised by different values of certain parameters. In all the simulations, we kept the value of M ov and the noise constant. Note that, as the configuration of the system varies, a certain value of M ov does not correspond to a unique value ofĒ a . Based on the results of the analysis, shown in Fig. S4, we can conclude that the method seems to be robust under uncertainties in the estimation of the parameters.  Figure S5. Trajectory of the model (S1) on a long time scale (15,000 years), with reference parameters as in Table S1 and with The downwelling had already stopped after about 1000 years, indicating that an F-type transition had occurred much earlier.

9/11
Section G: F-type transitions in a control simulation of a GCM The AMOC varies over a wide range of time scales and such variations can be driven by buoyancy or wind forcing anomalies. In 13 a detailed analysis of the observed AMOC variability is performed, with particular attention to the anomaly in 2009/2010. It is argued that the extreme event was caused by westerly wind anomalies associated with the North Atlantic Oscillation (NAO) negative phase, which was found to be correlated with the surface wind-driven Ekman transport. When looking at the strength of the AMOC in pre-industrial control simulations of CMIP5 models, we can observe several extreme events in the circulation. In this section we show results from the CMIP5 model MIROC5 14 . The value of M 35 • S ov for this model was calculated as -0.036 Sv 15 , which indicates that the AMOC is in a multiple equilibria regime. Therefore, buoyancy-induced extreme events might occur in the circulation: with the terminology introduced in this paper, such events may be associated with the occurrence of F-type transitions. In Fig. S6, the time series of the maximum value of the AMOC, the zonally integrated Ekman transport and the maximum Sverdrup transport in the western part of the basin, at 26.5 • N, are shown. The comparison between the time series shows that some dips in strength of the AMOC are most likely related to the occurrence of wind anomalies, while others may be the consequence of changes in the buoyancy forcing.  Figure S6. 20 month Butterworth low-pass filtered time series of the maximum value of the AMOC at 26.5 • N (black), the zonally integrated Ekman transport (red) and the maximum Sverdrup transport at the same latitude (blue), in the pre-industrial control simulation of the CMIP5 model MIROC5. The dashed horizontal lines indicate, for each time series, the lower bound of the 2 standard deviations confidence interval around the mean. The circles identify the time points at which extreme events occur, defined as the values extending below the confidence interval. The dotted magenta vertical lines are drawn in correspondence with the extreme events of the AMOC on top of all the time series. In this way, it is possible to assess whether a certain event is related with anomalies in the wind circulation. Some of the events seem to occur independently of changes in both Ekman and Sverdrup transports (e.g. the one found between time 80000 and 100000).