Theoretical Investigation of the Atlantic Multidecadal Oscillation

A weakly damped mode of variability, corresponding to the oceanic signature of the Atlantic multidecadal oscillation (AMO)was found through the linear stability analysis of a realistic ocean general circulation model. A simple two-level model was proposed to rationalize both its period and damping rate. This model is extended here to three levels to investigate how themode can draw energy from the mean flow, as found in various ocean and coupled models. A linear stability analysis in this three-level model shows that the positive growth rate of the oscillatory mode depends on the zonally averaged isopycnal slope. This mode corresponds to a westward propagation of density anomalies in the pycnocline, typical of large-scale baroclinic Rossby waves. The most unstablemode corresponds to the largest scale one (at least for low isopycnal slope). Themode can be described in four phases composing a full oscillation cycle: 1) basin-scale warming of the North Atlantic (AMO positive phase), 2) decrease in upper-ocean poleward transport [hence a reduction of the Atlantic meridional overturning circulation (AMOC)], 3) basin-scale cooling (negative AMO), and 4) AMOC intensification. A criterion is developed to test, in oceanic datasets or numerical models, whether this multidecadal oscillation is an unstable oceanic internal mode of variability or if it is stable and externally forced. Consistent with the classical theory of baroclinic instability, this criterion depends on the vertical structure of the mode. If the upper pycnocline signature is in advance of the deeper pycnocline signature with respect to the westward propagation, the mode is unstable and could be described as an oceanic internal mode of variability.


Introduction
Through its poleward heat transport, the Atlantic meridional overturning circulation (AMOC) is a major actor in setting the climate in the Atlantic (Srokosz et al. 2012). In particular, the AMOC variations are a potential explanation for climate variability in the Atlantic on time scales from decadal to centennial and longer. However, studies from observation, modeling, and theory suggest numerous and potentially inconsistent mechanisms that can lead to temporal variations of the AMOC [as reviewed by Yoshimori et al. (2010)]. Even if only focusing on comprehensive climate models, studies show that the AMOC can vary on a broad range of time scales with proposed mechanisms greatly varying from one model to the other [e.g., IPCC's Fifth Assessment Report (AR5); Stocker et al. 2014]. Thus, the mechanisms of AMOC variability remain a subject of continuing debates.
On multidecadal time scales, the Atlantic multidecadal oscillation/variability (AMO/AMV), with typical periods between 50 to 70 yr (Kushnir 1994;Delworth and Mann 2000), is characterized by large-scale warming and cooling at the surface of the Atlantic (leading to spatial average changes of the order of 0.1 K at the surface). This variability is associated with variations in the AMOC intensity of several Sverdrups (Sv; 1 Sv [ 10 6 m 3 s 21 ) (Knight et al. 2005). It is suggested that the AMO exerts significant impacts on climate (Sutton and Hodson 2003), affecting precipitation over Europe and North America (Sutton and Hodson 2005) and hurricane activity over the North Atlantic (Goldenberg et al. 2001), for instance. Given the length of the historical record compared to the relatively long period of the AMO, observational studies are strongly limited. Thus, AMOlike multidecadal variability has been studied in climate models of different complexity (Yoshimori et al. 2010;Latif 1998). Apart from the already mentioned AMOC variations, several other physical mechanisms have been proposed, including changes in the subtropical gyre (Dong and Sutton 2005;D'Orgeville and Peltier 2009;Cheng et al. 2004) or changes in the subpolar gyre that alter the heat budget of the Labrador Sea (Danabasoglu 2008), just to give a few examples.
Together with the 50-70-yr variability, Frankcombe et al. (2008) and Chylek et al. (2011) have suggested that the North Atlantic also varies with a 20-30-yr time scale by studying observations of temperature and sea level. Frankcombe and Dijkstra (2009) confirmed the existence of such sea level variation in both observations and climate models. This variability is characterized by a westward propagation of subsurface temperature anomalies (Frankcombe et al. 2008), making large-scale baroclinic Rossby waves the ideal explanation. Our study will focus on this particular 20-30-yr time scale.
Studies with idealized ocean models under prescribed surface buoyancy fluxes (as opposed to restoring surface boundary conditions) also show spontaneous modes of multidecadal variability characterized by the westward propagation of buoyancy anomalies (Greatbatch and Zhang 1995;Huck et al. 1999). Colin de Verdière and  first suggested the oscillations are sustained by large-scale baroclinic instability, mostly on the basis of the vertical structure of the perturbations, and are supported by density variance budgets (te Raa and Dijkstra 2002;Arzel et al. 2006). Huck et al. (2001) performed a local stability analysis for quasigeostrophic waves and showed that the mode could be associated with growing large-scale Rossby waves propagating westward in the presence of the meridional thermal gradient. Through a global linear stability analysis, Huck and Vallis (2001) proved that the oscillation results from an unstable linear mode, which can be strongly damped by surface restoring boundary conditions. In parallel to these studies, te Raa and Dijkstra (2002) interpreted the mode as resulting from temperature anomalies propagating westward along the polar boundary, inducing geostrophically alternating changes in zonal and meridional overturning and growing through the phase difference between temperature and vertical velocity anomalies (referred to as ''generalized'' baroclinic instability). Moving away from idealized configurations, te Raa et al. (2004) and Dijkstra et al. (2006) have shown that the mode is not modified when the realistic North Atlantic coastline is used instead of an idealized rectangular basin.
These westward-propagating buoyancy anomalies have been recently invoked to explain oscillations in an idealized geometry coupled model (Buckley et al. 2012). The authors point out potential regions of variance growth on the eastern and western boundaries of the subtropical gyre, suggesting Rossby waves are unstable in these regions. In contrast with previous studies that may have linked too tightly the oscillation mechanism to AMOC variations, they rather suggest the AMOC variations may passively respond to the Rossby waves through the thermal wind relation and not play an active role in modifying the sea surface temperatures (SSTs) and upper-ocean heat content budget in the northern regions. This analysis has been extended to the oscillations in more realistic coupled simulations, namely, with NCAR CCSM3 and GFDL CM2.1, highlighting the role of upper-ocean density anomalies propagating around the subpolar gyre and geostrophically affecting the AMOC upon their arrival in the northwest region .
More recently, Sévellec and Fedorov (2013) rigorously demonstrated the existence of a multidecadal natural mode of oscillation solely controlled by ocean dynamics in a realistic ocean general circulation model (OGCM). Using a linear stability analysis, they found the least damped oscillatory eigenmode of the system. The mode period is about 24 yr, and its damping time scale is close to 40 yr. This oscillation corresponds to a large-scale baroclinic Rossby wave in the presence of mean zonal flow and meridional temperature gradient. Its westward propagation comes from the competition between three factors (e.g., Tulloch et al. 2009): (i) the mean eastward zonal advection, (ii) the equivalent westward advection or geostrophic self-advection that depends on the mean meridional thermal gradient in the ocean, and (iii) the westward advection typical of largescale baroclinic Rossby waves (which is related to the b effect). The damping term is controlled by horizontal diffusion, and no unstable mode could be found. These last conclusions also apply in a two-layer, shallow-water model of the large-scale circulation forced either by wind or buoyancy forcing, with and without bottom topography (Ferjani et al. 2013(Ferjani et al. , 2014. What salient feature is missing from these different numerical model simulations to allow unstable modes leading to self-sustained decadal oscillations is what we investigate further. To allow baroclinic instability on decadal time scales, the interactions of at least two baroclinic modes must be permitted (Liu 1999). Therefore, we extend the theoretical model of Sévellec and Fedorov (2013) by adding a third level to their idealized two-level model. This level is included by splitting the first level representing the pycnocline in two equal parts. As it will be demonstrated, this allows baroclinic motion inside the pycnocline and the instability of the multidecadal mode. This is equivalent to the difference between shallowwater models with 1.5 and 2.5 layers, the latter allowing baroclinic instability (Simonnet et al. 2003).
We show that this newly developed idealized model has a multidecadal mode controlled by westwardpropagating temperature anomalies typical of largescale baroclinic Rossby waves, leading to the following sequence of 1) warming (inducing a positive AMO), 2) decrease in upper-ocean poleward transport (corresponding to a decreasing AMOC), 3) cooling (negative AMO), and 4) increase in upper-ocean poleward transport (increasing AMOC), and so on. This mode is unstable if the zonally averaged pycnocline slope is large enough. In our experiment, the largest-scale mode is most unstable for weak slope of the pycnocline. The instability occurs through a positive feedback when the upper pycnocline anomaly propagates a quarter period in advance of the deeper pycnocline anomaly. In this configuration, the upper and deeper pycnocline anomalies reinforce each other through the advection by anomalous meridional velocity of the mean temperature gradient (which is typical of baroclinic instability; Pedlosky 1996).
The structure of the paper is as follows: In section 2, we formulate the idealized ocean model used to analytically obtain the multidecadal mode of variability. In section 3, we discuss the properties of this mode of variability and its potential instability. In section 4, we summarize the implications of this study, especially in the context of the AMO.

The set of equations and model configuration
In this section, we will derive an idealized model of the ocean circulation and perform its linear stability analysis. The latter consists of computing the eigenvalues and eigenmodes of the linearized ocean dynamics, that is, the typical time scales and shapes associated with intrinsic ocean variability. We will demonstrate that the leading eigenmodes correspond to a multidecadal largescale baroclinic Rossby wave, which can be interpreted as the AMO.
The theoretical model consists of a flat-bottom rectangular basin (Fig. 1) representing the North Atlantic (from y 0 5 108N to y 1 5 708N). The zonal extent of the ocean is W 5 608, and its depth is H 5 4500 m (Table 1). The rotation rate varies to represent the curvature of Earth (b plane).
The set of equations for this theoretical model derives from the horizontal momentum equations, hydrostatic balance, nondivergence, a linear equation of state, and the evolution of temperature by advective-diffusive FIG. 1. A schematic of the idealized three-level model comprising two upper-ocean levels (levels 1 and 2) and a deep ocean one (level 3). The model three prognostic variables are the temperatures at the three levels (T 0 1,2,3 ). The six diagnostic variables are meridional and vertical velocities at each level (y 0 1,2,3 and w 0 1,2,3 ). The main parameters are the ocean layer thicknesses h 1 , h 2 , and h 3 (with the total ocean depth H 5 h 1 1 h 2 1 h 3 ), the zonal extent of the Atlantic basin W, and the mean zonal flow u and temperature T. The intensity of the shading (lighter to darker) represents the spatial variation of the mean temperature (cooler to warmer), varying linearly with latitude y in the upper ocean (levels 1 and 2). In the deep ocean, we use constant values equal to temperature at the northern basin boundary of the upper ocean. Also, we implicitly assume that there is a nonzero vertical temperature gradient in the upper layers that can support baroclinic Rossby waves due to the b effect. The dependency of the model variables on space coordinates (zonal x, meridional y, and vertical z) and time t is shown in parentheses. processes (calculations are restricted to temperature for simplicity, but numerical applications are done for equivalent temperature, i.e., density, to implicitly account for salinity). These equations can be simplified by assuming low viscosity (Re 1, i.e., viscous forces are negligible) and small inertial terms (Ro 1, i.e., horizontal equations are in geostrophic balance), where Re is the Reynolds number and measures the ratio of inertial forces to viscous forces, and Ro is the Rossby number. These assumptions are valid for the largescale ocean circulation and lead to the geostrophic balance. Thus, we obtain the noninertial set of equations in Cartesian coordinates described by Salmon (1998) in the absence of friction or viscosity: r 5 r 0 [1 2 a T (T 2 T 0 )], and (1d) where t is time; x, y, and z, are the zonal, meridional, and vertical coordinates; u, y, and w are the zonal, meridional, and vertical velocities; P is the pressure; r (r 0 ) is the (reference) density; T (T 0 ) is the (reference) temperature; f is the Coriolis parameter; g is Earth's gravity acceleration; a T is the thermal expansion coefficient; › t , › x , › y , and › z are the time, zonal, meridional, and vertical partial derivatives; D t is the material derivative (5› t 1 u› x 1 y› y 1 w› z ); and k and k y are the horizontal and vertical eddy diffusivity coefficients, respectively. Note that this set of equations corresponds to the firstorder balance of the widely used quasigeostrophic dynamics (Pedlosky 1979). It is also often referred to as planetary geostrophic equations (geostrophic regime of type 2 in Phillips 1963;Colin de Verdière 1988;Salmon 1998), stressing its validity to the study of basin-scale ocean dynamics.
To further simplify the system, the equation of evolution of temperature is linearized. For that purpose, we split the temperature field in a time-mean term and an anomalous term. Then, by assuming that the latter is always significantly smaller than the former, we neglect second-order terms. Note that the zero-order term describes the balance of the mean state, which will not be explicitly solved but rather used as a control parameter of the linearized dynamics. A scaling argument using the basin horizontal and vertical scale shows that vertical diffusion is much smaller than horizontal diffusion and thus can be neglected. The mean meridional and vertical velocities are neglected (y 5 w 5 0, respectively), restricting the mean state to a zonal flow at leading order, which is particularly accurate in the outcropping region. This assumption does not imply the absence of meridional overturning but its overall weakness compared to the zonal jet. Finally, in the linearized advective-diffusive equation for temperature anomaly, we consider an infinitely large meridional extent (› y 5 0, for the temperature anomaly). This last assumption is not needed, but it allows for an easier mathematical treatment of the problem. It implies that the horizontal deformation and propagation are primarily zonal (the latter being also consistent with a weak zonal gradient of mean temperature).
Following these assumptions, we obtain a new set of equations: where u 0 , y 0 , and w 0 are the zonal, meridional, and vertical velocity anomalies; T 0 is the temperature anomaly; u is the mean zonal velocity; and › y T and › z T are the mean meridional and vertical gradients of temperature, respectively. The thermal wind balance [(2b)] derives from the geostrophic balance together with the hydrostatic equilibrium [(1)]. The baroclinic condition [(2c)] derives from the vertical integration of the nondivergence equation with surface and bottom rigid conditions (w 0 j z50 5 w 0 j z52H 5 0) and u 0 5 0 (since › y 5 0 for the temperature anomaly). After discretization over three levels (following Fig. 1) and assuming no mean terms in the deeper level, the advection diffusion [(2a)] becomes whereas the thermal wind balance [(2b)] and the baroclinic condition [(2c)] are h 1 y 0 2 1 h 2 y 0 2 1 h 3 y 0 where the indices correspond to the three levels from the surface to the deep ocean, and h 1,2,3 are the constant thicknesses of each level related by h 1 1 h 2 1 h 3 5 H. To obtain this last set of equations, the vertical discretization follows the Arakawa C grid (together with simple linear interpolations, if needed). Using the three last equations, we can express the meridional velocity in levels 1 and 2 of the idealized model: In principle, the mean zonal velocity u and the mean meridional and vertical gradients › y T and › z T, respectively, can be different between the two upper levels. However, here we choose to use the same values over the two upper levels for a more straightforward comparison with the two-level approach of Sévellec and Fedorov (2013) and to reduce the number of parameters in the linearized dynamical system.
Finally, using the nondivergence together with the thermal wind, we can express the vertical velocities at each level using Sverdrup (1947) balance. Thus, assuming w 0 j z50 5 0, the vertical motion is obtained by vertically integrating (2d): where b 5 › y f is the gradient of planetary vorticity; w 0 1/2 is the anomalous vertical flow between layers 1 and 2; and w 0 1 and w 0 2 are the anomalous vertical flow at levels 1 and 2, respectively.
Since we know the meridional velocity anomalies as function of temperature anomalies, we also have the expression of the anomalous vertical motion in the two upper levels: We have obtained both the expression of the vertical and meridional velocity in the first and second level, making the system fully closed. Now using (5) and (7) in (3), we obtain the closed set of partial differential equations: Further, we expand temperature anomalies into Fourier harmonics in the zonal direction: where T cn k and T sn k are the Fourier coefficients, n (51, 2, 3. . .) is the wavenumber, k (51, 2, 3) is the level index, and c and s refer to the cosine and sine function coefficients, respectively. Mathematically there is no limit to n. However, since the geostrophic assumption in (1a) and (1b) breaks down for scale smaller than the first baroclinic Rossby radius of deformation ('20 km at midlatitude), we set the physical limit to n max 5 200 (i.e., a length scale of ;25 km).
Lateral boundary conditions are needed to constrain the advection-diffusion equation and to fully solve the problem. To obtain such conditions, we define two boundary layers (at the west and the east of the basin). These boundary layers connect the geostrophic interior to the basin boundaries. In the appendix, we consider in details the ocean dynamics in the boundary layers to constrain the Fourier expansion, which only describes the interior solution of the problem where geostrophic balance applies (geostrophic balance filters Kelvin waves as well as short Rossby waves).
On long time scales, such as decadal, we assume that the basin boundary layers can be treated as being continuously adjusted. This follows the study of Johnson and Marshall (2002), for example, who showed that the basin boundary adjustment in the Atlantic occurs over 2-3 months and is largely achieved by Kelvin waves propagating along the basin boundaries including the equator. For the idealized model, this assumption corresponds to the adjustment at all times of the east and the west boundary layers together: › t hTi EB 5 › t hTi WB , where h.i EB and h.i WB denote 3D spatial averages for the east and west boundary layers, respectively. This means that warming (or cooling) in both east and west boundary layers occurs at the same time. We also assume that the Kelvin waves, involved in this westeast adjustment, conserve their vertical structure. This allows us to apply independently the boundary layer condition for each level (e.g., a warming of the west boundary layer at level 1 leads to a warming of the east boundary layer at level 1 but does not affect level 2 and level 3 east boundary layers).
Applying these assumptions together with the temperature advection-diffusion equation, we demonstrate in the appendix that the appropriate boundary conditions for the interior solution become › x T 0 j East 5 2› x T 0 j West and T 0 East 5 2T 0 West (for details, see also Sévellec and Fedorov 2013). We apply these conditions for the three levels. In summary, they are a consequence of two factors: 1) the basin boundary is set by Kelvin waves' dynamics, and 2) the boundary layers' dynamics are necessary to connect the geostrophic interior solution with the basin boundary. In the context of the Fourier expansion, these boundary conditions restrict the solution to odd wavenumbers (n 5 1, 3, 5. . .) and allow for continuous oscillations.
Using (8) leads to a linear dynamical system (with solely ordinary differential equations that make the system easier to solve) with 6 degrees of freedom corresponding to the six Fourier amplitudes: The off-diagonal nonzero terms in (9) are calculated as where terms including u correspond to mean advection; terms including › y T correspond to geostrophic selfadvection (i.e., advection linked to anomalous meridional velocities of the background meridional temperature gradient); and terms including b correspond to long (nondispersive) baroclinic Rossby wave propagation. In very idealized settings, such as 1.5-layer shallow-water equations (where the dynamics are solely controlled by the first baroclinic mode), the geostrophic self-advection and the mean advection exactly cancel each other; this is known as the non-Doppler effect (Rossby et al. 1939;Held 1983;Killworth et al. 1997;Liu 1999).

a. Linear stability analysis: A multidecadal oscillation
To investigate the stability and exhibit the preferred modes of variability of this model, we perform a linear stability analysis that consists of computing the eigenvalues and eigenmodes of the Jacobian matrix of the linearized dynamical system (9). Unlike correlationbased analyses (e.g., empirical orthogonal functions analysis), linear stability analysis ensures the dynamical coherence of the time evolution of the modes (Monahan et al. 2009). There exists as many eigenvalues as degrees of freedom in the system (in our case n max , where n max is the highest wavenumber n) that provide information on the stability of the system; if the real part of a single eigenvalue is strictly positive, the system is unstable. More specifically, each eigenvalue sets if perturbation in the ''natural'' direction of the system defined by the associated eigenmode is stable or not (negative or positive real part, respectively) and is oscillatory or not (nonzero or zero imaginary part, respectively). In the rest of the study, particular attention will be given to describe the changes in the stability (eigenvalues) and the natural directions of the system (eigenmodes) with modifications of the model parameters (namely, stratification and horizontal diffusivity).
As will be demonstrated later, the eigenvalues show the existence of a large-scale oscillatory mode at multidecadal time scales. The eigenmode associated with these time and spatial scales is characteristic of a largescale baroclinic Rossby wave. This oscillation can be unstable for relatively low stratification. When unstable, the eigenmodes are consistent with baroclinic instability (Pedlosky 1996). Overall, we will demonstrate, through this linear stability analysis, that ocean dynamics develops a mode of variability with period and pattern similar to the AMO, which is potentially unstable (depending on the stratification).
The linear system described by (9) reveals six eigenvalues (l 1,2,3,4,5,6 ), two of which are simply set by the diffusive time scale: l 5,6 5 2k(np/W) 2 . These two eigenvalues are associated with purely decaying standing modes (i.e., nonpropagating) with signatures in all three layers, associated with the horizontal diffusion. These standing modes allow anomalies in the deep ocean to slowly vanish. The other four eigenvalues follow To test the growth or decay of these eigenmodes (i.e., stability or instability of the system), we build a control parameter, measuring the intensity of the mean vertical gradient of density c such that where DT and DS are the mean temperature and salinity contrast across the meridional direction, respectively, and a S is the haline contraction coefficient. In the same way, we assume a mean meridional gradient of density: where L is the meridional extent of the basin. Mean vertical and meridional temperature gradients are computed in terms of equivalent temperature that takes into account salinity through a linear equation of state. Computing the solution explicitly for salinity would lead to exactly the same solution in terms of spatial shape and time scales, where salinity anomaly would partially compensate temperature anomaly (equivalently to the mean gradients).
By combining (10) and (11), c, which measures the intensity of the stratification, can be expressed as c 5 (h 1 1 h 2 )/(uL), where u 5 2› y T/› z T. Thus, c is also a measure of the isopycnal slope of the mean state u. More precisely, the coefficient c is the inverse of the isopycnal slope normalized by the aspect ratio of the pycnocline. In general, we expect c ' 1, whereas c 1 corresponds to flat isopycnals and c 1 corresponds to important outcropping of the isopycnals. Next, we will test the influence of this parameter on the linear stability analysis, that is, how the stratification, or the ispoycnal slope, modifies the system stability. Results from the linear stability (i.e., the eigenvalues l 1,2,3,4 ) with varying c are summarized in Fig. 2 using parameters values from Table 1.
The stability analysis reveals the existence of eigenmodes with decadal to multidecadal periods (Fig. 2a). For values of c larger than approximately two, the twolevel model of Sévellec and Fedorov (2013) captures well the lower branch, corresponding to a damped eigenmode with a vertical structure close to a first baroclinic mode. The addition of the third level leads to a major difference: the emergence of an unstable eigenmode when the meridional density gradient is strong compared to the vertical gradient (for c # 0.33). This eigenmode represents an oscillation with a period around 23 yr. This bifurcation appears through the interactions of the lower branch with the upper branch (for c . 2 in Fig. 2a). The latter branch corresponds to eigenmodes with a vertical structure closer to a second baroclinic mode. When unstable, the mode period is particularly robust, varying between 22 to 33 yr when c is modified by more than one order of magnitude (Fig. 2a). As will be demonstrated later, the instability mechanism has the same behavior and properties as classical baroclinic instability but is acting on a basin scale. Following Colin de Verdière and Huck (1999), we call it a large-scale baroclinic instability.
The modes can be described by their signature on large-scale ocean metrics such as mean sea surface temperature (i.e., the AMO index) and AMOC (Fig. 3). The AMOC anomaly is simply the zonal integration of FIG. 2. Characteristic time scales of the two pairs of oscillatory eigenvalues for n 5 1 mode as a function of c measuring the mean vertical density gradient coefficient (or the inverse of the isopycnal slope), given by (10). (a) Period and (b) real part of the eigenvalue (inverse of the growth rate). The period corresponds to the time for an anomaly to cross the basin twice (once for the positive anomaly to go from west to east and once for the following negative anomaly to go from west to east; describing a full oscillation). Gray lines follow the results including solely two levels, that is, a single upper and a deep level [following Sévellec and Fedorov (2013), modified for the parameters of Table 1], in which case the real part of the eigenvalue is the inverse of the diffusivity time scale. the meridional transport in the two upper layers [W(h 1 y 0 1 1 h 2 y 0 2 )]. All the modes (independent of c, n, or the particular pair of eigenmodes) show an AMOC anomaly that leads the zonally averaged SST anomaly by a quarter cycle. Therefore, the multidecadal mode of variability can be described in four phases: (i) an increase of AMOC corresponding to an increase in poleward heat transport, (ii) leading to a warming of SST. The associated geostrophic flow extracts from the mean meridional temperature gradient a positive and negative temperature anomaly along the western and eastern basin boundaries, respectively. (iii) This east-west density difference causes a negative AMOC anomaly [(5)], leading to a decrease of meridional heat transport and so (iv) a cooling of SST. This oscillation results from the westward propagation of generalized large-scale baroclinic Rossby waves, where the meridional gradient of potential vorticity includes both the classical planetary b effect and the meridional density gradient, referred to as geostrophic self-advection (Sévellec and Fedorov 2013), reinforcing each other (note the latter is often larger than the former) and leading to the exact same four phases of oscillation. This description is consistent with both OGCM and idealized analysis of Sévellec and Fedorov (2013) and also the mechanism described by te Raa and Dijkstra (2002).
Even if the three processes contributing to the advection of temperature (i.e., the mean advection, the self-advection, and b-effect propagation) can be academically separated by their mathematical expression [cf. Fig. 4, following (3)], they are still physically linked.
For example, the mean zonal advection is partially controlled by the mean meridional gradient of density, linking the mean flow to the self-advection and explaining the non-Doppler effects (Rossby et al. 1939;Held 1983;Killworth et al. 1997). More fundamentally for the description of the above paragraph, the selfadvection and b-effect propagation are tightly linked. Indeed, it is the divergence of the meridional flow or AMOC [through the set of (6)] that induces the thermocline undulation related to the baroclinic Rossby waves. This made the oscillatory sequence described by the four phases between AMOC and SST (see paragraph above) also fundamentally linked to the classical large-scale baroclinic Rossby wave propagation induced by the b effect. However, it is important to note that, except for low isopycnal slope and away from the bifurcation, we found that the self-advection (i.e., advection of mean meridional temperature gradient by meridional velocity anomalies) term dominates over the mean zonal advection and the b-effect propagation (Fig. 4).
This stability analysis demonstrates that the covariability of AMOC and SST is an intrinsic property of the ocean dynamics, and their relation is tied to a quarterphase delay. This relation derives from the existence of large-scale baroclinic Rossby waves. However, this result should be carefully interpreted. Actually, despite good qualitative relation between the AMOC variability and SST variability, the quantitative link is more intricate. Given a normalized temperature anomaly, the AMOC anomaly is independent of the wavenumber n, The amplitude is arbitrary. The exponential growth/decay is suppressed. The spatial patterns and temporal evolution associated with these eigenmodes are shown in Fig. 7.

SEPTEMBER 2015
S É V E L L E C A N D H U C K whereas zonally averaged SST anomaly decreases as n increases. This implies that the quantitative relation between AMOC and zonally averaged SST anomalies depends on the wavenumber. To summarize, it is crucial to know the horizontal spatial scale of the anomaly (i.e., the wavenumber n) to make any quantitative prediction on the relation between AMOC and SST anomalies.

b. Instability threshold of the multidecadal oscillation
For all n modes, for high value of c, the two pairs of eigenmodes are oscillatory decaying (Fig. 2b); the decay time scale is controlled solely by the horizontal diffusion process. Decreasing c changes the period and shape of the two pairs of eigenmodes but does not modify the decay time scale (Figs. 2, 5) until the imaginary part of the eigenvalue of the two pairs of oscillatory eigenmodes are strictly identical (so that both growth rate, period, and shape are identical). At this stage, the two pairs of eigenmodes are merging, and there is a nonunique way to define the eigenmodes, since the two shapes cannot be separated by their intrinsic time scales. After this point, the pairs of eigenmodes remain oscillatory decaying (for n 5 1 mode, the period of both pairs of oscillatory eigenmodes is 23 yr), but one is more stable and the other less stable than before the bifurcation, when the damping time scale is solely set by the horizontal diffusion (Fig. 2b, where the gray line indicates the damping set by the horizontal diffusion). It is only when the instability of the less unstable pair of oscillatory eigenmodes overcomes the horizontal diffusion process that the steady state exhibits an unstable multidecadal mode of variability (Fig. 2b).
The merging of the eigenmodes is not straightforward to describe from classic theories of large-scale Rossby waves (Pedlosky 1996). For the low value of the isopycnal slope (c . 2), the propagation is mainly controlled by the baroclinic Rossby wave term (in comparison to the propagation by geostrophic selfadvection). In this regime, as expected from classic theory, the vertical structure of the eigenmodes corresponds almost exactly to a first and second baroclinic mode (Fig. 5). The propagation speed increases linearly with the vertical stratification intensity (so that the period of the modes decreases linearly with c). In this regime, the solution of Sévellec and Fedorov (2013) captures the solution of the three-level model (gray line vs lower branch in Fig. 2a). Decreasing the stratification further (0.33 , c # 2) results in an increase of the propagation speed of the slow pair of eigenmodes (formerly second baroclinic mode) and a further decrease of the propagation speed of the fast pair of eigenmodes (formerly first baroclinic mode), the latter diverging from the solution of Sévellec and Fedorov (2013). During this phase both pairs of eigenmodes become more and more surface intensified (Fig. 5). This surface intensification is obtained through a mixture of first and second baroclinic modes. This implies that the faster (slower) pair of eigenmodes is partially controlled by the second (first) baroclinic mode, making it slow down (accelerate) and get closer to the propagation speed of second (first) baroclinic mode. Here, the propagation speed of two pairs of eigenmodes converges. Also, the mixture of the first and second baroclinic modes makes the two pairs of eigenmodes nonorthogonal and interaction becomes possible. The surface intensification FIG. 4. As in Fig. 3, but for the three advective terms of (3): mean zonal advection (black solid lines), meridional self-advection (black dashed lines), and vertical advection related to the b effect (gray solid lines). The advective terms are zonally and vertically averaged over the two first levels. In all cases, the self-advection dominates over the two other terms.
increases when the stratification is decreased until the signature of both pairs of eigenmodes are almost zero in the second level (Fig. 5). At the bifurcation, the two pairs of eigenmodes are virtually identical with the same surface-intensified pattern, the same period, and the same damping time scale.
The fork bifurcation, that is, the location of the critical value of c when l r is split from one branch to two branches (Fig. 2b), is independent of n. However, since the diffusive damping time scale is shorter for the highest zonal wavenumbers (with a n 22 dependency), the largest modes are unstable for a higher value of c. This suggests that the largest modes require lower isopycnal slope to become unstable (i.e., to overcome the horizontal diffusion). Only the three largest zonal modes (associated with n 5 1, 3, 5) are able to become unstable (Table 2), whereas higher zonal wavenumbers never overcome the horizontal diffusion.
Since several modes of different wavenumber can be unstable for the same value of c, we defined the leading unstable mode as the one with the highest growth rate. We obtain that for c # 0.05, the leading unstable mode is n 5 5; for 0.05 , c # 0.27, it is n 5 3, and for 0.27 , c # 0.33, it is n 5 1. In summary, we can identify four regimes: 1) for high isopycnal slope (c # 0.05), we expect the ocean to be first unstable through the n 5 5 mode, corresponding to an oscillatory period of ;6 yr; FIG. 5. Shape of the two pairs of oscillatory eigenmodes (each column) for n 5 1 mode as a function of c, given by (10). Each row corresponds to a variable of the two first levels of the idealized model (i.e., the Fourier coefficients for n 5 1): (top to bottom) T c1 1 , T s1 1 , T c1 2 , and T s1 2 . The eigenmodes have no signature in the deep ocean, the third level of the idealized model. For each variable, plus symbols (1) represent the real part and crosses (3) represent the imaginary part. For high values of c, the shape of the two pairs of eigenmodes (1 and 2 vs 3 and 4) are different but slowly converge when c is decreased until they are strictly identical at the bifurcation and for lower values of c. 2) for intermediate isopycnal slope such as 0.05 # c # 0.27, we expect the ocean to be unstable through the n 5 3 mode, corresponding to an oscillatory period between 10.5 and 7.5 yr; 3) for intermediate isopycnal slope such as 0.27 , c # 0.33, we expect the ocean to be unstable through the n 5 1 mode, corresponding to an oscillatory period between 23 and 21.5 yr; and 4) for low isopycnal slope (c . 0.33), we expect the ocean to remain stable.
Note that in the first case described above, for high isopycnal slope (c # 0.05), the n 5 1 mode and n 5 3 mode are also both unstable, but their growth rate is weaker than n 5 5 mode. Based on the OGCM study of Sévellec and Fedorov (2013), we expect c ' 1, that is, between cases 3 and 4 above.
To extensively map the instability, we compute the period of the leading unstable eigenmode as a function of the vertical density gradient and horizontal diffusion (Fig. 6). In the absence, or for a low value (,100 m 2 s 21 ), of horizontal diffusion, the period is controlled by the small spatial-scale Fourier mode, which corresponds to high-frequency oscillation. Since our model assumes geostrophy (Ro 1), this limit reaches the validity of our approach (Colin de Verdière 1986). However, as soon as horizontal diffusion is not negligible, the leading unstable mode corresponds to larger spatial-scale Fourier modes. The vertical density gradient influence on the period of the oscillatory eigenmodes remains limited except close to the bifurcation. The horizontal diffusion selects the spatial scale of the instability and thus the period (i.e., different Fourier modes have very different oscillatory periods), whereas the vertical density gradient (variation of c from 0 to 0.3) controls weaker variation of the period by roughly 10%. Indeed, the vertical stratification influences the wave speed through the classical b effect, and hence the period of the oscillatory eigenmodes, but its impact remains limited because of the larger geostrophic self-advection through the mean meridional temperature gradient (Fig. 4;Sévellec and Fedorov 2013). Since in this model the horizontal diffusion represents effective diffusion induced by the eddy field, one can expect a nonnegligible horizontal diffusion and so a selection of the largest Fourier modes (as shown in Fig. 6). Following our calculation, we also show that the largest spatialscale Fourier modes are more likely to be destabilized FIG. 6. Period of the leading unstable eigenmodes as a function of the mean vertical density gradient coefficient c and the horizontal diffusion coefficient. Each value was evaluated over the 21 largest Fourier modes of the idealized model. Contour interval is 1 yr. The dotted region corresponds to the absence of the unstable mode. Large change in the period is concomitant with the change in the Fourier mode associated with the leading unstable eigenmode. The five main regions are labeled by their respective Fourier mode number n, and hTi represents the averaged period for these regions.
for weak isopycnal slope (i.e., c close to the instability) than higher ones (Fig. 6).

c. Instability mechanism of the multidecadal oscillation
To further understand the instability, we focus on the eigenmodes computed before and after both the ''fork'' and the instability (c 5 0.35 and 0.30, respectively). As expected from the eigenvalues analysis, the two pairs of eigenmodes for both cases are oscillatory waves propagating westward (Fig. 7). However, the structures of the two pairs of eigenmodes differ. In one case, before the instability (c 5 0.35), a pair of eigenmodes corresponds to a level-1 temperature anomaly in phase with the level-2 temperature anomaly [Figs. 7b(1)-(2), close to a first vertical baroclinic mode], whereas the other pair of eigenmodes corresponds to temperature anomalies out of phase between level 1 and level 2 [Figs. 7b(3)-(4), close to a second vertical baroclinic mode]. The in-phase pair of eigenmodes propagates faster (period of 18 yr) than the out of phase one (period of 26 yr), and they are both decaying because of the horizontal diffusion process (with an e-folding time scale of 36 yr). Right before the instability, both pairs of eigenmodes become more and more surface intensified as the stratification decreases (c going from 2 to 0.35), so that at the bifurcation they have only a surface signature (with virtually no level-2 signature). This convergence in pattern is accompanied by a convergence in period, so that the two pairs of eigenmodes are indistinguishable. In the other case, after the instability (c 5 0.30), one pair corresponds to a level-1 temperature anomaly in advance by a quarter phase to the level-2 anomaly [Figs. 7a(1)-(2)], whereas the other pair corresponds to a level-1 anomaly delayed by a quarter phase to the level-2 anomaly [Figs. 7a(3)-(4)]. Both pairs correspond to an oscillation with a period of 22 yr. However, the former is unstable (with a growth time scale of 19 yr), and the latter is stable (with decay time scale of 9.2 yr). The enhanced instability of the eigenmodes leaning toward the mean flow (i.e., level-1 temperature anomaly in advance of level-2 temperature anomaly with respect to the westward propagation) is consistent with classic theory of baroclinic instability (Pedlosky 1996). Unlike in the stable regime, the phase mismatch between the two upper levels can be seen in the slight lag between zonally averaged SST and ocean heat content [Figs. 3b(1)-(2)].
Note that the period and growth of the pair of eigenmodes corresponding to level-1 and level-2 temperature anomalies in phase can accurately be approximated by the pair of eigenmodes obtained with a model having a single upper layer (corresponding to the average of levels 1 and 2; Sévellec and Fedorov 2013) only at low isopycnal slope (c . 2; Fig. 2). For high isopycnal slope, the solutions of the two models diverge, since the single upper-level model cannot represent, by construction, the second baroclinic mode, and so temperature FIG. 7. Eigenmodes for the Fourier n 5 1 mode (a) after and (b) before both the fork bifurcation and the instability. In both cases, there exists two pairs of eigenmodes (the two rows) with a signature in level 1 and 2 of the idealized model (the two columns). (a) and (b) correspond to c 5 0.3 and 0.35, respectively. Subpanels 1 and 3 denote the level 1 and subpanels 2 and 4 denote the level 2. The first and second pair of eigenmodes are displayed in 1-2 and 3-4, respectively. Thick solid line corresponds to zero value. Grayscale shading is used for positive values. Contours without shading indicate negative values. Contour intervals for level 1 and 2 are 0.06 K and 6 mK, respectively. The amplitude is arbitrary. The exponential growth/decay of anomalies is suppressed.
anomalies intensified at the surface or the quarter phase delay between the two upper levels.
Before the instability, one can understand the period difference between the in-phase eigenmode and the outof-phase eigenmode (the first and second baroclinic modes, respectively) by looking at the impact of temperature anomaly at levels 1 and 2 on the vertical and meridional velocities, following the two sets of (5) and (7). Actually, the anomalous meridional and vertical velocities result in the propagation of nondispersive large-scale baroclinic Rossby waves [for further comments see Sévellec and Fedorov (2013)]. Applying (5), which is derived from the thermal wind balance and the baroclinic condition, reveals a striking difference between meridional velocities induced by temperature anomalies at levels 1 or 2. A temperature anomaly at level 1 leads to meridional velocities at levels 1 and 2 of opposite sign, whereas the same anomaly at level 2 leads to meridional velocities of the same sign at levels 1 and 2 ( Fig. 8). On the other hand, for vertical velocities [(7)], temperature anomalies at level 1 and level 2 lead to velocities of the same direction (Fig. 8).
Thus, applying (5) and (7) suggests that in-phase anomalies in levels 1 and 2 have a constructive impact on velocities (Fig. 8), whereas out-of-phase anomalies have a compensated effect (except for the meridional velocities at level 2). This suggests that in-phase anomalies, that is, the first baroclinic mode, induce stronger velocities and hence propagate faster (top-left panel of Fig. 9) than out-of-phase anomalies, that is, the second baroclinic mode (top-right panel of Fig. 9), consistently with classical theory of baroclinic Rossby waves (Pedlosky 1996). This leads to the shorter period of the in-phase eigenmode compared to the out-of-phase eigenmode. Sévellec and Fedorov (2013) assumed a single upper-oceanic layer, corresponding to the average of layers 1 and 2 of this study; consequently, their solution corresponds to the in-phase eigenmode (i.e., the first FIG. 8. Impact on the meridional and vertical velocities of a temperature anomaly in the (a1-3) level 1 or (b1-3) 2 for the Fourier n 5 1 mode. (a1) and (b1) represent the temperature anomaly imposed in level 1 and 2, respectively. (a2) and (b2) represent the meridional and vertical velocities response in level 1 (black and gray lines, respectively). (a3) and (b3) represent the meridional and vertical velocities response in level 2 (black and gray lines, respectively). Solid and dashed lines are associated with the shape of the Fourier coefficient (cos and sin, respectively) and their impact on the meridional and vertical velocities in levels 1 and 2, so that solid and dashed temperature curves are associated with solid and dashed velocities curves, respectively. baroclinic mode), the fastest pair of eigenmodes (gray line vs the lower branch in Fig. 2a).
Following the same approach, one can understand how the structure of the two pairs of eigenmodes after the instability leads to an unstable and a stable eigenmode. On the one hand, we suppose a temperature anomaly such as the level-2 signature is in a quarterphase delay with the level 1 [e.g., the solid line in Fig. 8a(1) and the dashed line in Fig. 8b(1)]. The level-1 anomaly induces a northward flow in the center of the basin at level 2 [solid black line in Fig. 8a(3)]. Through advective interaction with the meridional temperature gradient, this induces a warming in the center of the basin at level 2 (bottom-left panel of Fig. 9). Since this warming corresponds to the shape of the initial level-2 temperature anomaly, the latter is reinforced. This level-2 temperature anomaly, with this quarter-phase delay, induces a northward flow in the west of the basin and southward flow in the east in level 1 [dashed black line in Fig. 8b(2)]. These velocities extract from the mean temperature gradient a warm anomaly in the west and a cold anomaly in the east in level 1 (bottom-left panel of Fig. 9), hence reinforcing the initial level-1 anomaly. Overall, the level-1 anomaly feeds the level-2 anomalies and vice versa (bottom-left panel of Fig. 9). In a linear framework, this positive feedback implies a possible instability or growth mechanism.
On the other hand, in a symmetric manner, if the level-2 anomaly is a quarter phase in advance of the level-1 anomaly, the anomalies cancel each other (i.e., negative feedback), which implies a damping of the initial anomalies in a linear framework (bottom-right panel of Fig. 9).
This change in the vertical structure of the multidecadal mode of variability can be used to determine if the mode is stable or unstable. If the multidecadal mode FIG. 9. Schematic of the geostrophic self-advection behavior for low and high isopycnal slope (leading to a stable and an unstable mode, respectively) in the two upper levels of the idealized model. The warmer to colder background temperature, equivalent in both levels, is indicated by the orange to blue color gradient. Cold and warm temperature anomalies are indicated by the blue and red circles, respectively. The solid and dashed arrows indicate the geostrophic flow induced by the first-and second-level anomalies, respectively. Red and blue arrows indicate warming and cooling induced by meridional transport by the anomalous velocity of the mean temperature gradient, respectively. (Propagation by mean advection of anomalous temperature and through the b effect exists but is not indicated in the figure.) Dark arrows indicate feedback, where the anomalous velocity directly impacts the core of the anomalous temperature, whereas light arrows indicate propagation through self-advection. For low isopycnal slope, the two pairs of eigenmodes are close to the first and second vertical baroclinic modes, where the constructive or destructive self-advection between the two levels modifies the propagation speed of the first-level anomaly. For high isopycnal slope, the first and second baroclinic modes are mixed leading to feedback (positive or negative depending on which level is in advance with respect to the westward propagation) between the two levels and thus to instability or damping. exhibits a deep signature in phase or out of phase with its surface signature (even with surface intensification), the mode is stable. If the mode has a deep signature delayed by a quarter phase from its surface signature, the mode is unstable. Thus, we have a criterion to distinguish between an exogenous and an endogenous mode of variability. From an oceanic perspective, the former requires external forcing to be stimulated (e.g., atmospheric synoptic noise), whereas the latter is self-sustained with no need of external forcing.
Since this instability process is based on meridional advection, its efficiency is proportional to the meridional density gradient (i.e., › y T in our idealized case). However, it is only possible if the b-effect propagation is not too fast (› z T is weak). This explains the overall structure of the bifurcation in the linear stability analyses; that is, stable for low isopycnal slope, unstable for high isopycnal slope.
The interaction between levels 1 and 2 for high isopycnal slope (i.e., in the unstable regime) is not represented in the Sévellec and Fedorov (2013) approach, which only represents anomaly behavior averaged over the pycnocline (first baroclinic modes). This explains the limit of validity of Sévellec and Fedorov (2013) study to low isopycnal slope, that is, away from instability (Fig. 6).
The description of the instability and the vertical tilt of anomaly, that is, the quarter-phase delay between level-1 and level-2 anomalies, is typical of baroclinic instability (Pedlosky 1996). In our case, on planetary scale, we also consider horizontal diffusion (by mesoscale eddies) that acts as a scale selection, allowing the gravest mode to be more sensitive to instability than short-scale mode.

Conclusions
In this study, we demonstrate the existence of an oceanic multidecadal mode of variability associated with the westward propagation of large-scale baroclinic Rossby waves. This mechanism of variability could be responsible for the observed North Atlantic variations on the multidecadal time scale, often referred to as the AMO, which shows similar westward propagation, period time scale, and basin spatial scale (Frankcombe et al. 2008;Frankcombe and Dijkstra 2009;Chylek et al. 2011). This multidecadal variability has significant impact on climate by modifying evaporation/precipitation of bordering regions of the North Atlantic (Sutton and Hodson 2005) and hurricane activity over the tropical North Atlantic (Goldenberg et al. 2001), for instance.
Our result on the importance of the large-scale baroclinic Rossby waves for the multidecadal variability of the AMOC is consistent with earlier studies of Huck et al. (1999), Colin de Verdière and , te Raa and Dijkstra (2002), and Sévellec and Fedorov (2013). This confirms that the large-scale baroclinic Rossby waves are a major candidate to explain the AMO. In particular, we generalize the result of Sévellec and Fedorov (2013) to include one more degree of freedom on the vertical (i.e., addition of a vertical level). Through this modification, we demonstrate that this multidecadal mode becomes unstable if the mean isopycnal slope is strong enough.
In our calculations, the largest-scale mode is the most unstable one (i.e., unstable for the weakest isopycnal slope) because of the large-scale selection due to horizontal eddy diffusion. The instability itself is induced by a feedback between anomalies in the upper pycnocline and in the deeper pycnocline. This process is controlled through the meridional advection by anomalous velocity of the mean meridional density gradient. It is only possible if the density anomaly in the upper pycnocline propagates in advance of the deeper pycnocline anomaly, with respect to the westward propagation. All these results are consistent with the classical theory of baroclinic instability (Charney 1947;Eady 1949) where, in our study, a scale selection process is added through horizontal eddy-induced diffusion and leads to a largescale baroclinic instability (Huck et al. 2001). In the classical framework of vertical normal modes of the stratification, our results suggest that the instability results from the interaction (coupling) of the first and second baroclinic modes, the only ones remaining in our simplified three-level model. This means that the stability of the mode can be determined through its vertical structure. This analytical result is consistent with numerical experiments in OGCM; in Sévellec and Fedorov (2013), the least damped oscillatory eigenmode is close to a first vertical baroclinic mode, whereas in Arzel et al. (2006), the sustained oscillatory mode leans toward the mean flow (encompassing the first and second vertical baroclinic modes). Sévellec and Fedorov (2013) have shown the existence of the large-scale mode in a low-resolution OGCM in realistic configuration with a coefficient of horizontal tracer diffusivity equivalent to the one used in our study. They also demonstrated that this mode was the leading oscillatory mode, suggesting that the scale selection by horizontal diffusion mentioned in our study is valid. However, since horizontal diffusion, represented by a Laplacian operator, is a parameterization of subgrid processes, this result is not easily transposable to the real ocean. It can be argued that this is the action of the operator itself, and therefore of the parameterization rather than the process, that selects the large scale.
Despite this caveat, the presence of this mode in a highresolution model (Huck et al. 2015) tends to confirm and potentially generalize our analytical result, where eddy turbulence selects the gravest scale. Despite being encouraging, the validity of our analytical result for numerical models with different resolutions, and that for dynamics covering a broad spectrum of processes of different scales, needs to be further clarified in the future.
In low-resolution models, fronts are usually far too diffusive, for instance along the western boundary and North Atlantic Current. This potentially explains the damping of the multidecadal oscillation in Sévellec and Fedorov (2013). Simply moving to higher-resolution models will tighten the main fronts and increase the amount of energy available for baroclinic instability, both at mesoscale and large scale. Buckley et al. (2012) suggests baroclinic instability of the western and/or eastern boundary currents as well could excite basin modes. Since buoyancy budgets may not allow a robust localization of the growth, it may be difficult to unambiguously identify which region contributes the most.
In a recent study, LaCasce and Pedlosky (2004) suggested that large-scale baroclinic Rossby waves could not be maintained because of the impact of small-scale baroclinic instability on the front of the large-scale wave. This result, a priori, contradicts our finding that largescale baroclinic Rossby waves are the most unstable in the presence of eddy-induced diffusion. To test the two views, we have performed high-resolution modeling with a primitive equations model in an idealized flatbottom rectangular basin configuration of the North Atlantic (Huck et al. 2015). This study demonstrates the existence of sustained large-scale baroclinic Rossby waves in the presence of eddy turbulence, although it remains unclear if the large-scale mode is stable and sustained by the eddies' synoptic noise or if the mode is unstable. Its mere presence supports the theoretical result presented here.
Large-scale baroclinic instability, as well as its subsequent adjustment through an oceanic multidecadal oscillation (controlled by large-scale Rossby wave dynamics), is a potential explanation for the 20-30-yr signature of the Atlantic multidecadal oscillation. Actually, this large-scale mode has an impact on the AMO index (defined as the spatial average of sea surface temperature in the North Atlantic) by modifying the temperature in the upper ocean on a basin scale, consistent with observations by Frankcombe et al. (2008), Frankcombe andDijkstra (2009), andChylek et al. (2011). In this context, the AMO can be explained by a 20-30-yr North Atlantic oscillation, corresponding to an intensification of the AMOC followed by a basin-scale warming, an AMOC decrease, a cooling, and so on. Thus, tracking this mode of variability (e.g., propagation of large-scale Rossby waves), as well as its possible instability (e.g., increase of meridional slope of the pycnocline as a precursor of the AMO), would increase decadal climate predictability.
On the other hand, following the hypothesis suggested by Griffies and Tziperman (1995), two studies demonstrate that this exact same multidecadal oscillation can be sustained by atmospheric noise (Frankcombe et al. 2009;Sévellec et al. 2009), even in its stable regime. Since we have obtained a criterion, consistent with the classic theory of the baroclinic instability, on the stability of the multidecadal mode of variability, in future work, we will apply this criterion to coupled GCMs [from the IPCC's Fourth Assessment Report (AR4) and AR5 of Stocker et al. (2014)] and in situ observations. Looking at the leaning angle of temperature anomalies will potentially answer if the multidecadal mode is an endogenous or exogenous oceanic mode of variability, that is, if it is a self-sustained oceanic mode or if it is sustained by external stimulation (e.g., from synoptic atmospheric dynamics), respectively. Sensitivity is certainly expected to the model resolution affecting the mean state, the subgrid-scale parameterizations (eddies), vertical mixing, and air-sea interactions.

Zonal Boundary Conditions
In this section, we derive the east and west boundary conditions needed for the idealized model. The idealized model, by assuming geostrophic balance, only describes the interior solution. The boundary layers are necessary to connect the interior flow to the basin boundaries.
To conform to the setup of the idealized model, we introduce two boundary layers (of width d) at the western and the eastern sides of the basin (Fig. A1). We assume that the geostrophic balance is achieved at the outer edge of the boundary layers, that is, at the interface between a boundary layer and the interior. We also assume that geostrophic flow slips freely along the outer edge, that is, at the interface between a boundary layer and the interior. Finally, and consistently with the results of Johnson and Marshall (2002) and with the behavior of the OGCM in Sévellec and Fedorov (2013), we assume that the two boundary layers are coupled at all times, that is, adjusting together to perturbations.
To formulate the boundary conditions, we start from the general advection-diffusion equation for temperature [independently of the set of (1)]: › t T 5 2› x (uT) 2 › y (yT) 2 › z (wT) 1 › x (k x › x T) 1 › y (k y › y T) 1 › z (k z › z T), where x, y, and z are the zonal, meridional, and vertical coordinates (Fig. A1); u, y, and w are the zonal, meridional, and vertical velocities; and k x , k y , and k z are the zonal, meridional, and vertical diffusivities. For the sake of simplicity, we make the horizontal diffusion coefficients constant and equal to k h .
After vertical and meridional integrations, with the assumption of zero heat flux at the ocean surface (k z › z Tj z50 5 0), the advection-diffusion equation reduces to › t hTi z,y 5 2› x huTi z,y 1 k h › 2 x hTi z,y , where hTi z,y 5 Ð 0

2H
Ð L 0 T dy dz. We also use the nondivergence of the velocity field: › x u 1 › y y 1 › z w 5 0.
After vertical and meridional integrations, this equation becomes › x hui z,y 5 0. Integrating from the outer edge of the western boundary layer to the outer edge of the eastern boundary layer, we obtain huj d i z,y 5 huj W2d i z,y .
Thermal wind balance together with the free-slip condition (› x yj d 5 › x yj W2d 5 0) yields › 2 x T 5 0 at the outer edge of the boundary layers; therefore, along the border between the boundary layers and the interior flow, we can neglect the second term on the right-hand side of (A1).
Adding the equations for integrated temperature at the outer edge of each boundary layer and applying the results from the flow nondivergence, we obtain › t (hTj d i z,y 1 hTj W2d i z,y ) 5 2huj d i z,y (h› x Tj d i z,y 2 h› x Tj W2d i z,y ) .
Returning to (A1) and integrating it zonally over the western and eastern boundary layers, with the assumption of zero heat flux at the solid boundaries, we obtain that the integrated boundary layer temperatures evolve as › t hTi WB 5 2h(uT)j d i z,y 1 k h h› x Tj d i z,y , › t hTi EB 5 h(uT)j W2d i z,y 2 k h h› x Tj W2d i z,y , where hTi WB 5 Ð d 0 hTi z,y dx and hTi EB 5 Ð W W2d hTi z,y dx. Now assuming that the two boundaries adjust together, through Kelvin waves on time scales much shorter than large-scale Rossby waves, › t hTi WB 5 › t hTi EB , and again applying the flow nondivergence, we obtain huj d i z,y (hTj d i z,y 1 hTj W2d i z,y ) 5 k h (h› x Tj d i z,y 1 h› x Tj W2d i z,y ) .
In the context of the idealized model, (A2) and (A3) become and u West (T 0 where West and East refer to the western and eastern boundaries of the interior, geostrophic region. Combining these two equations, we obtain Time integrating (A5) yields This means that for sufficiently long times (at which the two boundary layers are fully coupled, i.e., for t k h /u 2 West ' 50 days), we have T 0 West 1 T 0 East 5 0. Finally, using this equality and simplifying (A4b) yields a set of two boundary conditions for the interior solution in the idealized model: East 5 2T 0 West , and (A6a)