## Abstract

A salient feature of paleorecords of the last glacial interval in the North Atlantic is pronounced millennial variability, commonly known as Dansgaard–Oeschger events. It is believed that these events are related to variations in the Atlantic meridional overturning circulation and heat transport. Here, the authors formulate a new low-order model, based on the Howard–Malkus loop representation of ocean circulation, capable of reproducing millennial variability and its chaotic dynamics realistically. It is shown that even in this chaotic model changes in the state of the meridional overturning circulation are predictable. Accordingly, the authors define two predictive indices which give accurate predictions for the time the circulation should remain in the on phase and then stay in the subsequent off phase. These indices depend mainly on ocean stratification and describe the linear growth of small perturbations in the system. Thus, monitoring particular indices of the ocean state could help predict a potential shutdown of the overturning circulation.

## 1. Introduction

Ice core and proxy temperature data from the Northern Hemisphere show strong variability on millennial time scales, but only during the (last) glacial interval, which sharply contrasts the relative stability of the subsequent Holocene epoch (GRIP Members 1993; Alley 2000; Andersen et al. 2004). Distinct episodes, called Dansgaard–Oeschger (D–O) events, make up a series of abrupt climate shifts with a presumed periodicity of about 1470 ± 500 yr especially pronounced in *δ*^{18}O records (Dansgaard et al. 1993; Bond et al. 1997; Grootes and Stuiver 1997; Bond et al. 1999). Typically, these episodes have two phases—a rapid warming of 5°–10°C occurring over a few decades followed by a slow cooling over several centuries. It has been conjectured very early (Broecker et al. 1990) that D–O events are related to variations in the ocean thermohaline circulation (THC), and the warm and cold phases of the oscillation correspond to different states of the Atlantic meridional overturning circulation (AMOC), with stronger or weaker circulation, and poleward heat transport, respectively. For simplicity, we will refer to these states as the overturning on and off phases. The goal of this study is to formulate a simple low-order model of the AMOC that can simulate such millennial variability realistically and then use this model to study predictability of the circulation shifts.

To explain the millennial signature in the paleorecords, one can invoke two different paradigms: exogenous or endogenous (i.e., extrinsic or intrinsic to the system). In the former paradigm, millennial variability arises as the system response to a forcing with a 1500-yr period (Alley et al. 2001; Rahmstorf 2003; Ganopolski and Rahmstorf 2001) or to random forcing (Ditlevsen et al. 2005; Ditlevsen et al. 2007). The second paradigm, supported by the present study, suggests that D–O events arise as a self-sustained oscillation whose irregularity is determined by deterministic chaos. The truth may lie somewhere in the middle as indicated by studies that show that millennial oscillations become irregular when noise is added to the system (Timmermann et al. 2003; Stocker and Johnsen 2003).

Note that one study questioned the very existence of the millennial oscillation and suggested instead that the millennial variability might be attributable to the aliasing of the seasonal cycle in the *δ*^{18}O records (Wunsch 2000). However, other authors emphasized the nonstationary character of the times series of oxygen isotopes as only a few D–O events contribute to the spectral peak of 1470 yr, which makes the aliasing hypothesis unlikely (Schulz 2002).

Even before the discovery of D–O events, Stommel found multiple equilibria in a simple box model of the THC (Stommel 1961). Subsequently, the existence of multiple equilibria for the ocean overturning circulation has been confirmed in zonally averaged ocean models of various complexity (Thual and McWilliams 1992; Quon and Ghil 1992; Cessi and Young 1992; Sévellec and Fedorov 2011), in forced ocean general circulation models (Rahmstorf et al. 2005), and in some coarse coupled models (Manabe and Stouffer 1988, 1999). At the same time, millennial oscillations have been replicated in idealized ocean models (Weaver and Sarachik 1991; Weaver et al. 1993; Winton and Sarachik 1993; Winton 1993, 1995) or either statistical model (Rial and Sasha 2011). Episodic but persistent reductions in the ocean overturning circulation have been found in the proxy data (McManus et al. 2004).

Studies using zonally averaged and box models (some exposed to gradually varying surface freshwater fluxes) have demonstrated the existence of millennial oscillations for intermediate values of the freshwater flux and proposed different low-order conceptual dynamic systems describing millennial oscillation (Cessi 1996; Sakai and Peltier 1999; Colin de Verdière 2007; Colin de Verdière et al. 2006; Sévellec et al. 2010). Different authors emphasize the role of salinity (Sakai and Peltier 1999; Sévellec et al. 2010), convection (Cessi 1996; Colin de Verdière et al. 2006; Colin de Verdière 2007), or the interhemispheric bipolar seesaw (Stocker and Johnsen 2003) for these oscillations. Apparently, the period of the oscillations in these models depends on how far the system is from the bifurcation (Colin de Verdière et al. 2006; Sévellec et al. 2010) and can go from about 1000 yr to infinity even for a relative, weak change in the oscillation amplitude, as expected for an infinite-period bifurcation (Strogatz 1994).

Although these models have improved our understanding of the THC millennial variability, they typically produce too regular D–O events unless forced with high levels of noise. Also, recent studies have argued that the bifurcation behavior and the millennial cycle in the paleorecords are more typical of a 3 degree of freedom system (Sévellec et al. 2010). his limits the validity of conceptual models with fewer degrees of freedom (most of previous idealized models had 1 or 2 degrees of freedom).

In parallel to these studies, the early warning signal in the climate system (EWS), that is, detecting climate tipping points by statistical methods, has become an area of extensive research (e.g., Livina and Lenton 2007; Dakos et al. 2008; Scheffer et al. 2009; Lenton 2011). The EWS approach assumes that in the close proximity to the tipping point, the attraction of the system statistical equilibrium becomes weaker, leading to a variance increase. In this forced paradigm (internal or external noise is needed to sustain the system variance), EWS has been detected within a range of idealized models of D–O events (Cimatoribus et al. 2013). However, these authors have also suggested that there is no evidence of EWS, as defined by Livina and Lenton (2007), in the nonstochastic system.

In the present study, we formulate a new low-order idealized model of the AMOC millennial variability, which can be derived from the equations of motion for ocean circulation in their zonally averaged form. The model has 3 degrees of freedom and generates irregular millennial variability within the endogenous paradigm. That is, ocean dynamics in the model without any stochastic forcing give rise to a deterministic chaotic behavior and induce variability on millennial time scales. We will describe the bifurcation diagrams of the model and its routes to chaos when surface freshwater flux increases. Further, we will demonstrate that even in this deterministic chaotic system, the shutdown of the AMOC can be accurately predicted by evaluating the growth of small perturbations. Following the approach of Evans et al. (2004), who have found that the growth of small disturbances depends on the imminence of the tipping point, we will show the existence of early warning indicators, or more generally of both qualitative and quantitative predictability, of the AMOC regime shifts.

## 2. The model, steady states, and their stability

The objectives of this work is to study the AMOC dynamics and its long-term variability (from centennial to millennial) and to develop and validate a methodology for predicting abrupt climate changes related to the reorganization of the AMOC. To that end, here we formulate the simplest deterministic model for the ocean meridional overturning circulation capable of chaotic behavior, which can be derived from the zonally averaged equations of motion on the depth-latitude plane [described, e.g., in Wright and Stocker (1991, 1992)]. In this model, we approximate the AMOC as a purely rotational circulation on the depth-latitude plan, for which deformations of the zonally averaged velocity field are neglected and the flow follows a single streamline (Fig. 1a). This last approximation implies that the model does not distinguish between shallow- and deep-ocean circulations.

This model is a modification of the Howard–Malkus loop (Fig. 1b), originally proposed by Welander (1957, 1965, 1967) and further developed by Howard (1971) and Malkus (1972); for a review, see Welander (1986) and Wunsch (2005). The loop model has been shown to describe centennial variability of the THC (Winton and Sarachik 1993; Sévellec et al. 2006). One of the versions of the loop model was formulated for three variables: two salt gradients and the overturning circulation strength (Dewar and Huang 1995, 1996; Huang and Dewar 1996). Being analogous to waterwheel models where water is replaced by salt, this model showed chaotic behavior (Strogatz 1994). Here, we will start with this version of the loop model but will add a new term representing a steady component of the overturning circulation driven by mean temperature gradients and surface winds.

### a. Model equations

The model three variables are the varying component of the ocean meridional overturning *ω* and the vertical (bottom–top) and meridional (north–south) salinity gradients, *S*_{BT} and *S*_{NS}, respectively, each depending on time *t*. We assume that the steady component of the circulation Ω_{0} is set by the mean oceanic temperature gradients and surface winds, and the variable part of circulation *ω* is controlled by variations in salinity gradients. The total overturning is Ω = Ω_{0} + *ω*, where Ω, Ω_{0}, and *ω* are measured in the units of yr^{−1}, which reflects the ocean overturning rates. (See appendix A for the model derivation starting from the zonally averaged equations of motion for the ocean circulation and the main approximations.) The model equations connecting these variables are as follows:

The first equation describes the momentum balance, in which the rate of change of the ocean overturning is set by the buoyancy torque with the coefficient *ϵ* and a linear friction with the coefficient *λ*, and *β* is the haline contraction coefficient. The two other equations describe the evolution of the salinity gradients driven by advection, linear damping (with the coefficient *K*), and the surface salt flux *F*_{0}*S*_{0}/*h*, where *F*_{0} is freshwater flux intensity, *S*_{0} is a reference salinity, and *h* is the loop depth or thickness (more exactly the depth of the level of no motion for the baroclinic flow). Particular values of the parameters are given in Table 1. For simplicity, we use virtual salt rather than freshwater fluxes.

The use of a constant Ω_{0} for the thermally and wind-driven circulation is a key approximation in our model, which assumes that, to the leading order, the atmospheric equator-to-pole temperature gradient is fixed (set by different insolation at the poles and at the equator), as is the meridional temperature gradient in the ocean. This approximation is based on the quick adjustment time scales of the atmospheric dynamics (a few months to a few years) as compared to the ocean overturning time scales (hundreds to thousands of years). The fixed oceanic temperature gradient, and the fixed winds, drives a constant overturning circulation in the ocean. (The assumption of steady thermal gradients could be relaxed by introducing two more equations for the ocean vertical and meridional temperature gradients and modifying the momentum balance appropriately.)

For a nonzero Ω_{0}, the model equations are now able to describe millennial variability with a deterministic chaotic behavior (Fig. 2). Understanding this variability and predicting regime shifts within the model is our next goal.

### b. Steady states

Using the set of equations in (1), we easily determine that the system’s steady states (*d*_{t} ≡ 0) are given by

where , , and represent the overturning circulation and the vertical and horizontal salinity gradients for the steady state. [Note that the existence of a nonzero steady-state circulation in the absence of vertical mixing does not contradict the ideas of Sandström (1908), since this circulation could be very shallow.] Applying a linearization yields a set of equations for small perturbations around these steady states as

where are anomalies of the state vector. The conventional linear stability analysis yields the following equation for the system’s eigenvalues *γ*:

Analyzing these equations while using the freshwater flux intensity *F*_{0} as a control parameter, we obtain several important results. First, for very low values of the freshwater flux (*F*_{0} < 0.6 cm yr^{−1}), there exists a single stable steady state (Fig. 3a). It corresponds to a solution largely dominated by the thermally and wind-driven circulation with a small contribution from salinity.

The freshwater flux intensity equal to approximately 0.6 cm yr^{−1} gives rise to an imperfect supercritical pitchfork bifurcation in the system. After that there exist two other steady states, which are both unstable. A second bifurcation occurs at approximately 1.6 cm yr^{−1}, that is, a subcritical Hopf bifurcation, transforming the stable steady state into an unstable one.

After this bifurcation (for higher values of the freshwater flux), the system has three unstable steady states: a state with zero , controlled by the balance between salt diffusion and the surface freshwater forcing, and two other states where the dominant balance is between advection and the freshwater forcing. The latter two (unstable) steady states are analogous, to some extent, to the thermal and haline (stable) steady states of the Stommel model (Stommel 1961). Depending on the sign of , the total flow can be greater or smaller than Ω_{0}. Setting Ω_{0} = 0 (i.e., neglecting both wind- and thermally driven circulations) would lead to the classical results of the rotating waterwheel (Strogatz 1994), that is, a perfect supercritical pitchfork bifurcation followed by two subcritical Hopf bifurcations.

### c. The model homoclinic orbit and route to chaos

In the regime with three unstable steady states, numerical integrations of (1) reveal the existence of a homoclinic orbit in the model phase space, that is, an orbit connecting an unstable steady state to itself. The periodicity of this orbit has a millennial time scale, and it appears through a global bifurcation leading to an infinite period at the bifurcation point (Strogatz 1994; Dijkstra 2000). The period of the oscillation corresponding to this orbit follows a −⅔ power law as a function of the distance to the bifurcation (in terms of the freshwater flux intensity; see Fig. 3), whereas the amplitude of the oscillation does not change. These characteristics of the millennial oscillation are consistent with the previous analysis of a zonally averaged latitude–depth model (Sévellec et al. 2010) and of another idealized model (Colin de Verdière 2007).

The millennial oscillation could be described as rapid switches between the on and off phases of the AMOC (negative and zero values of Ω after centennial variability has been filtered out, respectively), characteristic of a relaxation oscillation (Winton and Sarachik 1993; Winton 1993, 1995; Cessi 1996; Sakai and Peltier 1999; Colin de Verdière et al. 2006). During the on phase, a centennial oscillation with an advection-related period of 2*π*/Ω ~ 250 yr (which is an essential part of the model variability) acts to increase the meridional salinity gradient until the circulation stops, Ω ≃ 0 (Sévellec et al. 2006). During the off phase, the surface freshwater flux slowly reduces the vertical salinity gradient until the stratification becomes unstable and the circulation restarts. In contrast to previous studies (e.g., Winton and Sarachik 1993; Sévellec et al. 2010; Colin de Verdière 2007), the unrealistic overshoot of the circulation during its reactivation, on the order of 100 Sv (1 Sv ≡ 10^{6} m^{3} s^{−1}) over a few centuries during the so-called flushing events, is absent in the loop model. The off state is related to the weak circulation regime identified by Sévellec and Fedorov (2011), wherein heat and freshwater transports are primarily carried on by horizontal diffusion.

In the most realistic regime (Ω_{0} = −0.025 yr^{−1}, *F*_{0} = 1 m yr^{−1}) the loop model shows a chaotic behavior. In this regime, the system state clearly depends on the initial conditions (Fig. 2a). The model variability displays both centennial and millennial frequencies (Fig. 2b). However, the spectral peak corresponding to millennial variability is rather broad. This indicates that the models behavior follows not a limit cycle but a strange attractor. In such a regime, phase trajectories move like “a ball in a pinball machine” (Strogatz 1994). This chaotic behavior of the AMOC can be compared to the dynamics of atmospheric convection described by the Lorenz (1963) attractor; in fact, discussing the sensitivity of trajectories to initial conditions, we could even refer to a “shrimp effect” instead of the well-known butterfly effect.

Now that we have studied the bifurcation structure of the model behavior, we will use the model’s most realistic setting (Ω_{0} = −0.025 yr^{−1}, *F*_{0} = 1 m yr^{−1}) to investigate what controls the timing of the circulation shutdown. Two questions will be addressed: in the circulation on phase, could we predict (i) the time left before the shutdown and (ii) the duration of the incoming off phase?

## 3. Methodology and results: Predicting AMOC regime shifts

In a recent study, Evans et al. (2004) has demonstrated that a prediction of regime change is possible in the famous Lorenz model (Lorenz 1963; Saltzman 1962). Perturbing the model chaotic trajectories by instantaneous random perturbations, the authors of that study showed that the growth of these perturbations may indicate the impending regime change, which gives an example of a qualitative prediction in a chaotic system. They also showed that longer durations of the perturbation growth before the regime change are associated with longer durations of the next regime. However, their quantitative prediction had a strong uncertainty.

In the present paper, we will follow a similar philosophy but apply a less ad hoc methodology to perturb model trajectories, a methodology based on a generalized stability analysis (Farrell and Ioannou 1996a,b) as applied to our loop model. The generalized stability analysis allows finding perturbations with a structure optimal (or most efficient) for perturbing the AMOC in a linear framework ( appendix B). We will demonstrate that the growth of the optimal perturbations provides a precursor of the AMOC change (Palmer 1999), which can be used for both qualitative and quantitative prediction.

### a. Qualitative prediction

As a first step, we introduce the optimal perturbation growth parameter *g*, giving the growth rate of the perturbations that induce the maximum change of ocean meridional overturning (see appendix B for the exact mathematical definition). This parameter *g* is a local in-time diagnostic computed over time scales (*τ* = 20 yr) significantly shorter than the characteristic time scales of the loop model (~250 yr and ~2 kyr). It varies along the system trajectories in the phase space and helps identify the regions of the attractor where a perturbation grows (*g* > 0) or decays (*g* < 0) (see Fig. 4). In the latter regime, small disturbances have no significant impact on the system long-term evolution, while in the former regime, *g* sets the bounds on the error growth in the system.

Computing this growth parameter in the model phase space reveals a cylindrical symmetry, that is, a radial symmetry on the plane (*S*_{BT}, *S*_{NS}) for a constant Ω (Fig. 4d). Thus, despite the chaotic nature of the strange attractor, the growth parameter *g* shows a regular structure. The axis defining the cylindrical symmetry (Ω constant) crosses *S*_{BT} ≃ −0.5 psu and *S*_{NS} ≃ 0.5 psu (defining the center of stability, i.e., the minimum of perturbation growth). This offset between the center of stability and the steady state implies that, during the on phase, the perturbation growth is greater for high values of salinity gradients ( and, but less significantly, ) and low values of circulation intensity ().

In summary, a close qualitative relationship between the growth of perturbations and the time left before the shutdown of the AMOC (Fig. 4a) exists. In other words, the closer we are to an AMOC shutdown, the stronger the expected growth of optimal perturbations. This statement suggests that our methodology has general qualitative predictive skills even in the chaotic regime, which is consistent with the previous idealized study of Evans et al. (2004). The next step is to investigate whether our approach has quantitative predictive skills as well. The key question is whether one could use the optimal perturbation approach to define a quantitative measure to assess the temporal proximity and future duration of the AMOC shutdown.

### b. Quantitative prediction

To obtain proper statistics we have performed several 2 × 10^{6}-yr-long simulations. Each of the simulations has about 1000 successions of on- to off-phase intervals (red and blue, respectively, Fig. 5a). Next, using these simulations, we will explore the quantitative predictability of the system regime shifts.

#### 1) Predicting the time left in the on phase

To investigate the predictability of the regime shift, we compare the time left in the on phase against the perturbation growth parameter *g*. This comparison reveals a tight correspondence between the two variables (Fig. 5b). Thus, we have arrived at a remarkable result—any value of the growth parameter *g* corresponds to a particular time left in the on phase, within an uncertainty range. This result implies quantitative predictive skills: the knowledge of the growth parameter *g* provides an accurate prediction of the time left in the on phase. Larger values of *g* suggest a quickly approaching AMOC regime shift.

Despite this prediction skill, there is small uncertainty that it depends on *g* (Fig. 5b). However, measured as the standard deviation, this uncertainty decreases for higher values of *g* (closer to the regime shift). Even the relative error (standard deviation divided by the mean) decreases with the growth parameter for *g* > 0.1 yr^{−1}. Consequently, our prediction becomes more and more accurate close to the regime change.

#### 2) Predicting the duration of the incoming off phase

To investigate the quantitative predictability of the time the system will stay in the off phase, we use a slightly different procedure. For each on-phase interval, we find the maximum growth parameter *G* along the phase–space trajectory (stars in Fig. 5a): *G* = max(*g*)|_{on-phase}. This parameter corresponds to the most unstable part of the trajectory and gives the maximum possible growth rates of optimal perturbations during each on-phase interval.

Plotting the time spent in each off phase as a function of *G* evaluated over the preceding on-phase interval reveals a very good correspondence between the two variables (Fig. 6a). In fact, each *G* is associated with a unique duration of the next off-phase interval, which suggests a high level of predictability. In other words, knowing *G* for an on-phase interval gives an accurate prediction for the total duration of the next off phase. The prediction accuracy exceeds 98% (this number is given by the fraction of predictions that fit into the general pattern of Fig. 6a).

To understand why the time spent in the next off phase and the maximum growth parameter *G* during the preceding on phase are related, we should first look at what controls *G*. Plotting the three model variables, Ω, *S*_{BT}, and *S*_{NS}, at the time of the maximum perturbation growth as a function of *G* reveals that this parameter is most closely related to *S*_{BT} (Fig. 6). The connection of *G* to Ω and *S*_{NS} is more uncertain and is affected by the chaotic behavior of the system. With a perfect knowledge of *G* we are able to obtain the value of *S*_{BT} with a less than 0.1% uncertainty, but we were able to obtain the values of Ω or *S*_{NS} with a 4% and 12% uncertainty, respectively. In general, the close relation between the maximum growth parameter *G* and each of the three variables at the time of the maximum perturbation growth is an important result.

Now we are ready to discuss the relation between *G* and the off-phase duration and its different regimes. First, for large values of *G*, ocean circulation is weak but vertical stratification is strong (Figs. 6b,d). According to Fig. 5b, large values of *G* should lead to a rapid shutdown of the circulation, during which the ocean can retain this strong stratification. Then, it will take a long time to reduce this stratification by the surface freshwater flux forcing before the circulation can resume, which makes the time spent in the off phase quite long.

For low values of *G*, ocean circulation is strong and stratification is weak (Figs. 6b,d). The time left before the shutdown is rather long (Fig. 5b). The circulation needs a sufficient time to slow down before the shutdown. However, this rather long delay allows the building of a strong ocean stratification before the shutdown of the circulation through the interaction between the overturning flow and freshwater flux, that is, through the positive salinity feedback (Marotzke 1996) and its oscillatory counterpart (Sévellec et al. 2006). Again, after the shutdown, it will take a long time to reduce this stratification before the circulation can resume.

For intermediate values of *G*, the shortest possible duration of the future off phase is achieved (dashed vertical lines in Fig. 6a). This occurs because of the saturation of the north–south salinity gradient (Fig. 6c), which sustains the circulation. Over this saturation value, the ocean circulation is too fast, and the salinity gradients cannot increase any longer through the positive salinity feedback (Marotzke 1996; Sévellec et al. 2006), which is consistent with Stommel’s model (Stommel 1961).

## 4. Conclusions

In this work we have developed an idealized model of ocean overturning circulation, which reproduces millennial variability as a relaxation oscillation controlled solely by internal oceanic processes. This oscillation is characterized by rapid transition (that occurs over just a few decades) between two different phases of the AMOC. Our model is based on modified rotating waterwheel equations (Strogatz 1994), representing the ocean overturning circulation as a sum of two components—a constant component representing the thermally and wind-driven circulation and a time-varying component driven by variations in the vertical and horizontal salinity gradients (Howard 1971; Malkus 1972). This representation allows our model to describe a multiequilibrium regime characteristic of the thermohaline circulation (Stommel 1961). Despite the model simplicity, it exhibits a rich chaotic behavior qualitatively consistent with the characteristics of the Dansgaard–Oeschger events in the paleorecords.

To understand the behavior of the model we have obtained its bifurcation diagram by varying the freshwater flux intensity *F*_{0}. As *F*_{0} increases, the bifurcation diagram reveals a pitchfork bifurcation (where a single, stable steady state splits into two stable states and one unstable state). This pitchfork bifurcation is perfect if and only if Ω_{0} = 0 and imperfect otherwise, which introduces asymmetry in the system behavior for Ω_{0} ≠ 0. For sufficiently strong *F*_{0} (>2 cm yr^{−1}) and any Ω_{0}, the steady states become unstable through a subcritical Hopf bifurcation. It is after this second bifurcation that the millennial oscillation becomes possible. In this regime, the freshwater flux induces a vertical salinity stratification, which reduces density stratification. This result is consistent with the work of Colin de Verdière and te Raa (2010) and Arzel et al. (2010), who stressed the importance of the weak ocean mean stratification for the existence of millennial variability in ocean–atmosphere coupled models.

This oscillation, which is a 3 degree of freedom oscillation, is described by a homoclinic orbit in the phase space emerging in the absence of any stable steady states. The orbit period depends on the distance to the bifurcation as a −⅔ power law (as a function of freshwater flux). This implies that exactly at the bifurcation point the period of the homoclinic orbit is infinitely long, which is consistent with previous studies (Colin de Verdière 2007; Sévellec et al. 2010).

For realistic parameters (*F*_{0} = 1 m yr^{−1} and 2*π*/Ω_{0} ≃ 250 yr), the model exhibits a chaotic behavior (corresponding to positive Lyapunov exponents), and the system trajectories in the phase space follow a strange attractor. In this regime, with a strong sensitivity to the initial conditions, small oceanic perturbations can potentially lead to a shutdown of the AMOC (the shrimp effect).

The power spectrum of the model time integration in this realistic parameters regime has a broad peak in the millennial band, consistent with the available data (Wunsch 2000), as well as weaker peaks at centennial frequencies. Furthermore, our idealized model is able to represent the irregular strengthening and weakening of the AMOC, typically associated with the D–O events, driven solely by internal oceanic processes (without external noise or any time-dependent forcing).

The chaotic properties of our idealized model make it a useful tool to test predictability of the AMOC collapse. To that end, we applied a generalized stability analysis, imposing instantaneous linear optimal perturbations along phase trajectories and evaluating the AMOC intensity after a given time delay (20 yr). This analysis shows higher growth rates of the optimal perturbations prior to the circulation collapse, indicative of qualitative skills to predict and track the AMOC collapse.

To make this analysis more quantitative, we have introduced a perturbation growth parameter *g* to show the existence of an empirical law relating, within a small uncertainty, the time left before the AMOC collapse to this parameter. We have also established another empirical law connecting the maximum amplitude of this parameter during the on phase, *G* = max(*g*)|_{on-phase}, to the total duration of the subsequent off phase (i.e., the time lapse between the collapse and recovering of the AMOC). As soon as these laws are established we do not need to run the model anymore to predict the timing of the next collapse and recovery of the AMOC. Computing these two parameters, *g* and *G*, will provide these predictions.

Our study has also shown that the maximum perturbation growth parameter is tightly linked to ocean vertical stratification (and not the intensity of the AMOC or the north–south salinity gradient). In other words, in a 3 degree of freedom model, only one variable of the model affects the prediction of the AMOC collapse. The predictive skill of the salinity stratification needs to be tested further in more realistic models exhibiting millennial oscillations, including zonally averaged ocean and coupled models (Sévellec et al. 2010; Colin de Verdière and te Raa 2010) and intermediate coupled models (Arzel et al. 2010).

The crucial role of the salinity stratification, in the context of the AMOC collapse, has been investigated in a range of coupled models by Huisman et al. (2010) and Cimatoribus et al. (2013). These studies suggested that the meridional freshwater transport is an indicator of the AMOC bistable regime, so that a decrease of the ocean salinity stratification increases the possibility of an AMOC collapse. Unlike our study, these studies are exclusively in a steady state context, so direct comparisons are not possible.

Since it is nearly impossible to monitor the full trajectory of the actual ocean, the generalization of our results could be particularly significant: to monitor a potential collapse of the AMOC, we might not need to know the entire state of the ocean, but rather its specific indices, even though the choice of these indices is not straightforward. Even in our idealized model, before the full analysis of the model sensitivity to optimal perturbations, it was not obvious if one parameter was sufficient to anticipate the AMOC collapse or which particular parameter. Nor is it clear how many indices have to be monitored in a more realistic system. This highlights the need to perform such an analysis in more complex cases. Its implementation in ocean general circulation models will be a goal of our future work.

## Acknowledgments

This research was supported by grants from DOE Office of Science (DE-SC0007037), the David and Lucile Packard Foundation, and the Ti Ammo project funded through the French CNRS/INSU/LEFE program. We thank Inez Fung for helpful comments and suggestions.

### APPENDIX A

#### Derivation of the Idealized Model

To derive the idealized ocean model, we start with a description of the AMOC dynamics commonly used in 2D (latitude–depth) zonally averaged models (e.g., Sévellec and Fedorov 2011). This approximation invokes a linear frictional balance that relates the overturning mass transport to the meridional gradient of baroclinic pressure (e.g., Wright and Stocker 1991, 1992). In this paper, we simplify the dynamics further going from 2D to one dimension (1D). To do so, we neglect the deformation of the overturning flow and only consider its rotational part. This can be achieved by following a single streamline of the AMOC. In such an approach, shallow and deep circulations become virtually indistinguishable. Nevertheless, this simplification still allows an accurate representation of both the positive salinity feedback (Marotzke 1996) and the centennial variability of the THC, as demonstrated by Sévellec et al. (2006). These two processes are crucial for the millennial AMOC variability as discussed in section 2c and in Sévellec et al. (2010).

To express these ideas mathematically, it is convenient to place ourselves in the center of the ocean basin, so that we are inside the simple, closed curve formed by the streamline (Fig. 1). We introduce the angle *θ*, positive for counterclockwise motion and measured from the bottom of the ocean, as the new spatial coordinate. Then, the density gradient can be expressed as a buoyancy torque, while the flow intensity will not vary with *θ* (since the flow is incompressible). Following a linear frictional balance, this leads to an equation describing the changes in the overturning circulation:

where *θ* is the angle around the loop, Ω_{THC} is the thermohaline circulation intensity, *λ* is a friction coefficient, and *ϵ* is a proportionality factor between the overturning circulation and the buoyancy torque.

Retaining an inertial term in the momentum balance follows two arguments. (i) As shown by numerical calculations of Johnson and Marshall (2002), the ocean circulation does not adjust immediately to meridional pressure gradient. Consequently, incorporating inertial terms in zonally averaged model leads to a more realistic variability (Drbohlav and Jin 1998). (ii) Mathematically, there is no reason to exclude a long-term adjustment due to the inertial term. One could expand , where *σ* is a small parameter (≪1), using a scale separation between different ocean processes (*t* = *t*_{0} + *σt*_{1} + *σ*^{2}*t*_{2} + …). For example, even when the initial east–west adjustment is completed, the system can still undergo further adjustment on longer time scales.

Using a linear equation of state of seawater and assuming that thermally driven component of the ocean circulation is always in a steady state (controlled by fast atmospheric processes) we obtain

where *ω* and Ω_{T} are the haline and thermally driven circulations (such that Ω_{THC} = *ω* + Ω_{T}), and *α* and *β* are the thermal expansion and haline contraction coefficients, respectively. In the loop model we can also add a third term for the circulation representing the effect of the winds Ω_{W}. This term derives from the surface boundary conditions. We thus obtain that the constant overturning flow is Ω_{0} = Ω_{W} + Ω_{T}, and the total circulation is Ω = Ω_{0} + *ω*. This separation, which is used operationally to monitor the AMOC [Rapid Climate Change (RAPID); Hirschi and Marotzke 2007], is possible because of the linearity of the momentum equation. Further, we use this total circulation in the advective–diffusive equation controlling the evolution of salinity in the loop:

where *K* is the Laplacian eddy diffusion coefficient, is the freshwater flux, and *h* is the loop thickness or depth.

The periodicity of the loop implies that the spectrum of the solutions is discrete in terms of *θ* and allows the following Fourier decomposition:

where ℜ indicates the real part, with a reverse projection

where the subscripts *r* and *i* stand for the real and imaginary parts, respectively. Since the freshwater forcing cannot impact modes with *n* ≠ 1 (no projection onto these modes), we restrict our consideration to the first mode (no sources of potential energy are available for other modes and hence they thus have only trivial solutions). Using notations *S*_{NS} = *S*_{i1} and *S*_{BT} = *S*_{r1}, to explicitly represent the north–south and bottom–top salinity gradients, we now arrive at the incompressible loop model equations rewritten as the following:

The time integration of the model equations follows a fourth-order Runge–Kutta method. Two initial conditions, randomly chosen in the phase space of the model, are used to plot Fig. 2 (black and red lines). This allows the testing of the system’s sensitivity to initial conditions and, at the same time, demonstrates the deterministic chaotic behavior in the model.

### APPENDIX B

#### Computing Optimal Perturbations

To compute the optimal perturbations we first define

where |**U〉** is model along a the trajectory. We also define 〈**U**| associated with the Euclidian norm: 〈**U** | **U**〉 (a scalar product).

The prognostic equations in (1) can be rewritten as a general autonomous dynamical system:

where is a nonlinear operator.

Next, we define the anomalous state vector as a perturbation to the model trajectory:

where |**u**〉 is the anomalous state vector, *ω*′ is the anomaly of circulation intensity, is the anomaly of the bottom–top salinity gradient, and is the anomaly of the north–south salinity gradient.

The time evolution of the perturbation is described by

where (*t*) is the Jacobian matrix evaluated for a particular value |**U**(*t*)〉. Integrating this equation from initial time *t*_{i} to the final time *t*_{f}, we obtain |**u**〉 as a function of time (Farrell and Ioannou 1996b) as

where (*t*_{f}, *t*_{i}) is the nonautonomous propagator of the linearized dynamics from time *t*_{i} to *t*_{f}.

To define linear optimal perturbations, we also need to define the optimality and choose a measure to maximize. Here, following Sévellec et al. (2007), as the measure we use the anomalous circulation intensity: *ω*′ = 〈**F** | **u**(*t*_{f})〉, with 〈**F**| = (1, 0, 0). Finally, we need to define a norm that will ensure nondegenerate solutions for optimization in a linear framework. For this norm we choose

where is the operator defining the inner product associated with this particular norm.

Subsequently, the linear optimal perturbations are obtained using a Lagrangian function that allows us to maximize the anomalous circulation intensity under the constraint of a normalized initial perturbation (Sévellec et al. 2007):

where *ν* is a Lagrangian multiplier and *τ* = *t*_{f} − *t*_{i} is the duration (i.e., the delay after which the impact of the initial perturbation is assessed). The solution of the maximization problem is given by the condition *d*(*t*_{i}, *τ*) = 0. Applying this condition, we find the linear optimal perturbations as

where † denotes the adjoint defined through the scalar product.

Using this expression, we now define the linear growth parameter of the optimal perturbation *g* at a particular time along the trajectory *t*_{i} as

This parameter shows how much the norm of the linear optimal perturbation will grow for the specified time delay *τ*. For our computations we choose *τ* = 20 yr; testing other values of *τ* shows no qualitative changes to the results. The only constraint we impose is that this time scale be significantly shorter than the characteristic time scales of the loop model dynamics (~250 yr and ~2 kyr), which makes *g* a local in-time diagnostic.

## REFERENCES

*Mechanisms of Global Climate Change at Millennial Time Scales, Geophys. Monogr.,*Vol. 112, Amer. Geophys. Union, 35–58.

^{−3}- to 10

^{−5}-year time resolution

*Geophys. Res. Lett.,*

**34,**L03712, doi:10.1029/2006GL028672.

*Mém. Soc. Roy. Sci. Liége,*

**4,**125–128.

*Decadal Climate Variability: Dynamics and Predictability,*D. L. T. Anderson and J. Willebrand, Eds., NATO ASI Series I, Vol. 44, Springer, 333–378, doi:10.1007/978-3-662-03291-6_8.

*Geophys. Res. Lett.,*

**30,**1510, doi:10.1029/2003GL017115.

*Geophys. Res. Lett.,*

**32,**L23605, doi:10.1029/2005GL023655.

*Large-Scale Transport Processes in Oceans and Atmosphere,*D. L. T. Anderson and J. Willebrand, Eds., NATO ASI Series C., Vol. 190, Springer, 163–200, doi:10.1007/978-94-009-4768-9_4.

*Ice in the Climate System,*W. R. Peltier, Ed., NATO ASI Series I, Vol. 12, Springer, 417–432, doi:10.1007/978-3-642-85016-5_24.