# Nonnormal Multidecadal Response of the Thermohaline Circulation Induced by Optimal Surface Salinity Perturbations

Florian Sévellec, Thierry Huck, Mahdi Ben Jelloul, Jérôme Vialard

► To cite this version:

Florian Sévellec, Thierry Huck, Mahdi Ben Jelloul, Jérôme Vialard. Nonnormal Multidecadal Response of the Thermohaline Circulation Induced by Optimal Surface Salinity Perturbations. Journal of Physical Oceanography, 2009, 39 (4), pp.852-872. 10.1175/2008JPO3998.1. hal-00384666

HAL Id: hal-00384666

https://hal.science/hal-00384666v1

Submitted on 22 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Nonnormal Multidecadal Response of the Thermohaline Circulation Induced by Optimal Surface Salinity Perturbations

FLORIAN SÉVELLEC

LOCEAN, Paris, France

THIERRY HUCK AND MAHDI BEN JELLOUL

Laboratoire de Physique des Océans, UMR 6523 CNRS/IFREMER/UBO, Brest, France

JÉRÔME VIALARD

IRD, LOCEAN, Paris, France

(Manuscript received 14 February 2008, in final form 8 October 2008)

## ABSTRACT

Optimal perturbations of sea surface salinity are obtained for an idealized North Atlantic basin using a 3D planetary geostrophic model—optimality is defined with respect to the intensity of the meridional overturning circulation. Both optimal initial and stochastic perturbations are computed in two experiments corresponding to two different formulations of the surface boundary conditions: the first experiment uses mixed boundary conditions (i.e., restoring surface temperature and prescribed freshwater flux), whereas the second experiment uses flux boundary conditions for both temperature and salinity. The latter reveals greater responses to both initial and stochastic perturbations that are related to the existence of a weakly damped oscillatory eigenmode of the Jacobian matrix, the optimal perturbations being closely related to its bi-orthogonal. The optimal initial perturbation induces a transient modification of the circulation after 24 yr. The spectral response to the optimal stochastic perturbation reveals a strong peak at 35 yr, corresponding to the period of this oscillatory eigenmode. This study provides an upper bound of the meridional overturning response at multidecadal time scales to freshwater flux perturbation: for typical amplitudes of Great Salinity Anomalies, initial perturbations can alter the circulation by +2.25 Sv (1 Sv ≡ 10⁶ m³ s⁻¹; i.e., 12.5% of the mean circulation) at most; stochastic perturbations with amplitudes typical of the interannual variability of the freshwater flux in midlatitudes induce a circulation variability with a standard deviation of 1 Sv (i.e., 5.5% of the mean circulation) at most.

## 1. Introduction

A strong modification of surface air temperature in the North Atlantic during the past century has been established (Mann et al. 1999)—in the context of global warming. This temperature modification is concomitant with a modification of sea surface salinity (SSS) in the same region noted since the mid-1970s and is related to an increase of precipitation in the North Atlantic sub-polar gyre (Josey and Marsh 2005). A similar salinity modification has also been measured in the deeper ocean

by Curry et al. (2003) and Curry and Mauritzen (2005), who point out the lack of information about the evaporation and precipitation over the ocean and the water cycle, as a result of the measurement uncertainties. In this paper, we will investigate the influence of such water cycle modifications on the ocean dynamics. The climate has many time scales of variability (Mitchell 1976), but the slow dynamics is largely controlled by the ocean: the slow part of the ocean circulation, called the thermohaline circulation, is associated with a net northward heat transport in the Atlantic and thus is able to generate strong modification of the current climate (Vellinga and Wood 2002). Given these considerations, we will focus on the influence of surface salinity and freshwater flux perturbations on the thermohaline circulation.

<footnote>*Corresponding author address:* Florian Sévellec, Dept. of Geology and Geophysics, Yale University, P.O. Box 208109, New Haven, CT 06520-8109.
E-mail: florian.sevellec@yale.edu</footnote>
Generalized stability analysis (Farrell and Ioannou 1996a,b) provides a powerful tool for studying the ocean response to perturbations in a linear framework. Because of advective–diffusive processes, the ocean dynamics is “nonnormal” (Farrell and Moore 1992)—that is, the ocean response does not have the same structure as the ocean sensitivity. This means that the most efficient perturbation to excite the fastest-growing (or least damped) mode, obtained by a classical stability analysis, is not necessarily this fastest-growing mode structure and can in fact have a very different spatial structure. To understand the growth of a perturbation over finite time rather than the asymptotic limit, we must thus use generalized stability analysis that takes into account the nonnormality of the dynamics.

In the context of nonnormal dynamics and generalized stability analysis, Farrell (1988) determined the optimal perturbations in an atmospheric model. The method was used in an oceanic context by Farrell and Moore (1992) and Moore and Farrell (1993) to find the most rapidly growing perturbation of a quasigeostrophic ocean model representing the Gulf Stream. Farrell and Ioannou (1993) then extended the approach to determine the optimal stochastic perturbation in an atmospheric model—that is, the pattern that, if modulated by a stochastic amplitude, induces the maximum variability. These methods were fully described by Farrell and Ioannou (1996a) for the autonomous operators (the background state is independent of time) and by Farrell and Ioannou (1996b) for nonautonomous cases (the background state evolves in time). These theories were applied by Moore et al. (2002) in an ocean circulation model and particularly in a tropical configuration to find the optimal perturbation of ENSO (Moore et al. 2003).

Regarding the thermohaline circulation variability, several studies were undertaken. Lohmann and Schneider (1999) investigated initial error growth and predictability in the Stommel (1961) two-box model by finding optimal initial and stochastic perturbations. More recently, Tziperman and Ioannou (2002) searched the optimal initial and stochastic patterns in a three-box model. The mechanism leading to the transient growth is due to the rapid decay of the temperature anomaly and thus the enhancement of a salinity-dominated density anomaly. The spectral response to the optimal stochastic perturbation did not reveal any peak despite the existence of a damped oscillatory eigenmode. The same kind of analysis was performed in a two-box model but within a nonlinear framework (Mu et al. 2004): the authors show the limit of the linear approach and the asymmetry of the advection term in the nonlinear equations. Sun et al. (2005) extended this work in a coupled ocean–atmosphere box model and suggested that for weak amplitude perturbations,

the linear approach remains valid. In this linear context, Zanna and Tziperman (2005) studied the transient growth at 40 yr of an optimal initial perturbation in a simple ocean–atmosphere coupled model (latitude–depth ocean model with only two levels on the vertical and a one-layer atmosphere): the growth mechanism is due to the advection by the circulation anomaly, which increases both temperature and salinity anomalies. Recently, this work has been extended to look at the optimal stochastic perturbation in the same model (Zanna and Tziperman 2008) where, despite the existence of oscillatory eigenmodes, the spectrum of the variability is mostly red without any peak. In a global ocean general circulation model, Sévellec et al. (2008) found two different optimal initial perturbations of the surface salinity linked to two different measures of the circulation [the intensity of the meridional overturning circulation (hereafter MOC) and the meridional heat transport]. These two different optimal perturbations induce transient growth at 10.5 and 2.2 yr, respectively. This study also provided the max bound of variability as a result of surface salinity perturbations. In a reduced space based on the main empirical orthogonal functions of a coupled ocean–atmosphere general circulation model, Tziperman et al. (2008) show the existence of transient growth both for an energy norm and a quadratic norm of the circulation. The same kind of theoretical tools are also used for sensitivity studies of the thermohaline circulation. For instance, Sirkes and Tziperman (2001) have studied the sensitivity of the meridional heat transport at 24°N and found an oscillatory mode with a centennial time scale. In the same way, Bugnion et al. (2006a,b) have illustrated the sensitivity of the ocean circulation to different forcing terms, such as the heat flux, the surface freshwater flux, or the wind stress.

Sévellec et al. (2007) found the optimal initial, constant, and stochastic surface salinity perturbations in a latitude–depth ocean model and provided upper bounds of the associated variability. Their study shows the importance of a centennial oscillatory weakly damped eigenmode of the linear tangent model and its biorthogonal in the ocean response and optimal patterns, respectively. Here we propose to perform the same analysis in a 3D ocean model, where the most weakly damped modes of variability are on decadal-to-multidecadal time scales. The use of a planetary geostrophic model, where the state vector is reduced to temperature and salinity fields, allows us to represent a realistic 3D dynamics but also to perform both linear and generalized stability analysis with respect to all eigenmodes and their biorthogonal. Our approach uses a linear measure (the intensity of the MOC) and constraints (namely, salt conservative surface salinity perturbations) to provide the answer to two
specific physical questions, in our case: what are the initial and stochastic surface salinity perturbations that have the greatest effect on the MOC? Computing the initial optimal perturbations can be useful to understanding the effect of a sudden modification of the sea surface salinity on the thermohaline circulation, like, for example, a sudden release of melting continental ice into the ocean. Stochastic optimal perturbations are relevant to understanding how high-frequency atmospheric forcing related to midlatitude weather regimes [the North Atlantic Oscillation (NAO), for instance] can induce long-term variability of the thermohaline circulation. In addition, these analyses provide the upper bounds on the modification of the MOC.

Some ocean models use constant fluxes for both heat and freshwater ("flux" boundary conditions), while others use fixed freshwater forcing and sea surface temperature (SST) relaxation ("mixed" boundary conditions, justified by the feedback of SST on air-sea fluxes). Weaver and Sarachik (1991) and Huck and Vallis (2001) have pointed out that the stability of the thermohaline circulation is modified by the use of one or the other conditions. In the same way, Arzel et al. (2006) have shown that the variability occurring under each condition strongly differs. We will thus compute the optimal forcing patterns for both flux and mixed boundary conditions in this study.

We first introduce the ocean model and configuration (section 2) and describe the steady state, which is valid for both experiments. Asymptotic linear stability analyses are then performed for each experiment (section 3). In section 4, we investigate the optimal initial perturbation and describe the induced transient growth for both experiments. In section 5, we perform the same analysis for the optimal stochastic perturbation and compute the spectral response. Finally, we discuss our results and propose some outlook (section 6).

## 2. The ocean model and steady state

The model is a planetary geostrophic model in spherical coordinates: only the thermodynamic equations are prognostic, with the equations for dynamics being diagnostic. This approximation corresponds to the limit of the quasigeostrophic approximation for spatial scale of the order of the earth's radius (Colin de Verdière 1988). The dynamic equations are

$$f \mathbf{k} \times \mathbf{u} = -\rho_0^{-1} \nabla_H P - \varepsilon \mathbf{u}, \quad (1a)$$

$$-\partial_z P - \rho g = 0, \quad \text{and} \quad (1b)$$

$$\nabla_H \cdot \mathbf{u} + \partial_z w = 0, \quad (1c)$$

where $\nabla_H$ is the horizontal gradient, $z$ the vertical coordinates, $f$ the Coriolis parameter, $\mathbf{u}$ and $w$ the horizontal and vertical velocity fields, $P$ the pressure, $\rho$ ($\rho_0$) the (reference) density, $\varepsilon$ the linear friction coefficient, and $\mathbf{k}$ the unitary vertical vector. The use of linear friction instead of Laplacian is not fundamental and is not as important as the choice of lateral boundary conditions, for instance—at least after its numerical implementation (Huck et al. 1999b). Here no-slip is used along all boundaries.

The thermodynamic equations are similar to the primitive equations:

$$\partial_t T + \mathbf{u} \cdot \nabla_H T + w \partial_z T = K_H \nabla_H^2 T + K_V \partial_z^2 T + C_T + \mathcal{F}_T, \quad (2a)$$

$$\partial_t S + \mathbf{u} \cdot \nabla_H S + w \partial_z S = K_H \nabla_H^2 S + K_V \partial_z^2 S + C_S + \mathcal{F}_S, \quad (2b)$$

$$\rho = \rho_0(1 - \alpha T + \beta S), \quad (2c)$$

where $t$ is the time, $T$ the temperature, $S$ the salinity, $K_H$ ($K_V$) the horizontal (vertical) tracer diffusion coefficient, $\alpha$ ($\beta$) the thermal expansion (haline contraction) coefficient, and $C_T$ ($C_S$) the convective adjustment effect on temperature (salinity). This convective adjustment sets in when $\partial_z \rho > 0$ and mixes $T$ and $S$ instantaneously (conserving salt and heat content) until the density profile is stable. Here, $\mathcal{F}_T$ ($\mathcal{F}_S$) is the atmospheric forcing applied to temperature (salinity) in the surface mixed layer $h_s$. The other boundary conditions do not allow heat or salinity flux through the boundaries.

As discussed in section 1, we will use both restoring and flux boundary conditions for temperature, which distinguish the two experiments presented in the following (sections 4 and 5). Obviously, the absence of a sea surface salinity feedback on atmospheric freshwater flux justifies the use of a prescribed flux boundary condition for salinity for both experiments, $\mathcal{F}_S = S_0$ SSF/$h_s$, where $S_0$ is the reference salinity and SSF is the freshwater flux (conserving the total salt content; Fig. 1) corresponding to the evaporation minus precipitation flux SSF = $F_0 \sin[2\pi(\theta - \theta_0)/(\theta_1 - \theta_0)]$. The first experiment corresponds to the restoring condition (mixed boundary conditions, hereafter MBC) $\mathcal{F}_T = \tau_T^{-1}$ (SST* - SST), where $\tau_T$ is the restoring time scale and SST(*) the (restoring) surface temperature (Fig. 1), according to SST* = 0.5($T_1 + T_0$) - 0.5($T_1 - T_0$) $\cos[\pi(\theta - \theta_0)/(\theta_1 - \theta_0)]$. The direct numerical time integration of (1) and (2) for this experiment leads to a steady state (Fig. 1). The second experiment corresponds to a prescribed heat flux boundary condition (FBC), $\mathcal{F}_T = (\rho_0 C_{pw} h_s)^{-1}$ SHF, where $C_{pw}$ is the seawater heat capacity and SHF the surface heat flux
<figure><img src="image_2.png" /><figcaption>FIG. 1. Steady state of the planetary geostrophic model. Surface (a) restoring temperature, (b) SHF, and (c) SSF. MOC intensity during the time integrations under (d) MBC and (e) FBC. (f) Surface temperature and (g) salinity. (h) ZOC and (i) MOC. Contours (b) from −80 to 50 W m⁻² by 10 W m⁻²; (f) from 4° to 24° by 2°C; (g) from 34.8 to 36.2 by 0.2 psu; (h) from −2 to 10 by 1 Sv; and (i) from 0 to 18 by 2 Sv. In all the figures, solid, dashed, and dotted lines correspond to positive, negative, and zero values, respectively.</figcaption></figure>

(its basin integral is zero to conserve the total heat content; Fig. 1). This flux is computed from the temperature-restoring term in the MBC experiment steady state. The time integration of (2) under FBC shows that this modification of the boundary condition does not affect the steady state, which remains the same (Fig. 1); nevertheless, the stability of the steady state might be modified, as shown in the next section.

The domain of the model corresponds to a flat-bottomed pie-shaped representation of the North Atlantic basin, from 10° to 60°N with 64° zonal extension. The uniform basin depth is 3400 m, discretized on 10 vertical levels with thicknesses of 100, 150, 200, 250, 300, 400 m, and four deep levels of 500 m—the first level corresponds to the surface mixed layer $h_s$. The other parameters used for the numerical time integrations have been chosen such that the steady state remains stable for both experiments (Table 1): 10 000 yr of numerical integration provides the steady state for the MBC experiment, and a restart for a 100-yr-long time integration of the FBC experiment, after the prescription of the equilibrium surface heat flux, confirms that the steady state is unchanged (Fig. 1).

The sea surface temperature decreases from 24°C to the south of the domain to 4°C in the north. There is also a gradient between salty water around 36 psu at ~25°N and fresher water around 34.8 psu to the north of the basin. This steady state induces a realistic MOC of 18 Sv (1 Sv ≡ 10⁶ m³ s⁻¹), intensified in the north and at the surface. This meridional circulation corresponds to a

northward surface transport and a downwelling along the northern boundary associated with deep convection; upwelling with a southward transport occurs in the rest of the basin (Fig. 1), as theoretically anticipated by Stommel and Arons (1960). The stratification mainly occurs in the first 1000 m (Figs. 2a,b) and induces a first

<table><caption>TABLE 1. Parameters used for the time integrations of the planetary geostrophic model.</caption><tr><td>n<sub>x</sub></td><td>16</td><td>Number of grid points in longitude</td></tr><tr><td>n<sub>y</sub></td><td>14</td><td>Number of grid points in latitude</td></tr><tr><td>n<sub>z</sub></td><td>10</td><td>Number of grid points on the vertical</td></tr><tr><td>dt</td><td>1 day</td><td>Numerical time step</td></tr><tr><td>H</td><td>3400 m</td><td>Uniform ocean depth</td></tr><tr><td>h<sub>s</sub></td><td>100 m</td><td>Surface mixed layer depth</td></tr><tr><td>W</td><td>64°</td><td>Zonal basin extent</td></tr><tr><td>θ<sub>0</sub></td><td>10°N</td><td>Southern boundary position</td></tr><tr><td>θ<sub>1</sub></td><td>60°N</td><td>Northern boundary position</td></tr><tr><td>K<sub>H</sub></td><td>1800 m<sup>2</sup> s<sup>-1</sup></td><td>Horizontal tracer diffusion</td></tr><tr><td>K<sub>V</sub></td><td>10<sup>-4</sup> m<sup>2</sup> s<sup>-1</sup></td><td>Vertical tracer diffusion</td></tr><tr><td>g</td><td>9.8 m s<sup>-2</sup></td><td>Gravity acceleration</td></tr><tr><td>ρ<sub>0</sub></td><td>1000 kg m<sup>-3</sup></td><td>Reference density</td></tr><tr><td>C<sub>ρω</sub></td><td>4000 J kg<sup>-1</sup> K<sup>-1</sup></td><td>Seawater heat capacity</td></tr><tr><td>a</td><td>6.4 × 10<sup>6</sup> m</td><td>Earth radius</td></tr><tr><td>S<sub>0</sub></td><td>35 psu</td><td>Reference salinity</td></tr><tr><td>α</td><td>2 × 10<sup>-4</sup> K<sup>-1</sup></td><td>Thermal expansion coefficient</td></tr><tr><td>β</td><td>8 × 10<sup>-4</sup> psu<sup>-1</sup></td><td>Haline contraction coefficient</td></tr><tr><td>τ<sub>T</sub></td><td>132 days</td><td>Temperature restoring time scale</td></tr><tr><td>F<sub>0</sub></td><td>50 cm yr<sup>-1</sup></td><td>Freshwater flux intensity</td></tr><tr><td>T<sub>0</sub></td><td>26°C</td><td>Restoring temperature at θ<sub>0</sub></td></tr><tr><td>T<sub>1</sub></td><td>3°C</td><td>Restoring temperature at θ<sub>1</sub></td></tr><tr><td>ε</td><td>4.4 × 10<sup>-6</sup> s<sup>-1</sup></td><td>Linear friction coefficient</td></tr></table>
<figure><img src="image_3.png" /><figcaption>FIG. 2. Stratification and circulation of the model steady state. Zonal mean (a) temperature and (b) salinity, (c) surface velocities averaged over the upper 1000 m, and (d) deep velocities averaged from 1000 m to the bottom. (c),(d) The contours represent the vertical velocities and the quivers, the horizontal velocities. Contours from (a) 6° to 24° by 2°C; (b) from 35 to 35.9 by 0.1 psu; (c) from -3.36 × 10⁻⁵ to 1.39 × 10⁻⁵ m s⁻¹ by 10⁻⁷ m s⁻¹; and (d) from -2.34 × 10⁻⁵ to 1.19 × 10⁻⁵ m s⁻¹ by 10⁻⁷ m s⁻¹.</figcaption></figure>

baroclinic mode circulation with strong northward surface velocities along the western boundary (Figs. 2c,d). The downwelling mainly occurs in the north of 50°N, with intensification along the northern boundary and the northern part of the eastern boundary. The rest of the basin corresponds to a large-scale upwelling, intensified in the western boundary.

### 3. Asymptotic stability analyses

In this section, we perform a linear stability analysis for both MBC and FBC experiments (which have exactly the same steady state). Such analyses allow us to examine the asymptotic evolution of a small perturbation around a steady state. The prognostic Eq. (2) can be rewritten as a general dynamical system,

$$d_t|U\rangle = \mathcal{N}(|U\rangle), \quad (3)$$

where $\mathcal{N}$ is a nonlinear operator and $|U\rangle$ the state vector consisting of temperature and salinity at every grid point, written as a ket; the associated bra $\langle U|$ is defined through the Euclidian scalar product $\langle U|U\rangle$. We decompose the state vector as $|U\rangle = |\bar{U}\rangle + |u\rangle$, where $|\bar{U}\rangle$

represents the steady state, such that $\mathcal{N}(|\bar{U}\rangle) = 0$, and $|u\rangle$, the perturbation. The time evolution of the perturbation reads

$$d_t|u\rangle = \mathbf{A}|u\rangle, \quad \mathbf{A} = \left.\frac{\partial \mathcal{N}}{\partial|U\rangle}\right|_{|\bar{U}\rangle}, \quad (4)$$

where $\mathbf{A}$ is the Jacobian matrix, which is a function of the steady state $|\bar{U}\rangle$ (autonomous system). We can integrate this last relation to obtain the perturbation as a function of time:

$$|u(t)\rangle = \exp(\mathbf{A}t)|u(0)\rangle = \mathbf{M}(t)|u(0)\rangle, \quad (5)$$

where $\mathbf{M}(t)$ is the propagator of the linearized dynamics.

The Jacobian matrix is numerically evaluated by introducing weak positive and negative perturbations in the full nonlinear model and then evaluating the linear response as the mean of these responses (Huck and Vallis 2001). Once the Jacobian matrix is computed, an eigenanalysis is performed. For both MBC and FBC experiments, all eigenvalues have a negative real part, confirming that the steady state is stable. The pairs of
<table><caption>TABLE 2. Decay time scales (yr) and periods (yr) of the leading eigenvalues ($\lambda_r + i\lambda_i$) for MBC and FBC experiments.</caption><thead><tr><th rowspan="2"></th><th colspan="2">MBC</th><th colspan="2">FBC</th></tr><tr><th>τ = 1/λ<sub>r</sub></th><th>T = 2π/λ<sub>i</sub></th><th>τ = 1/λ<sub>r</sub></th><th>T = 2π/λ<sub>i</sub></th></tr></thead><tbody><tr><td>Eigenmode 1</td><td>-163.3</td><td>—</td><td>-207.4</td><td>34.2</td></tr><tr><td>Eigenmode 2</td><td>-58.1</td><td>328.5</td><td>-74.1</td><td>—</td></tr><tr><td>Eigenmode 3</td><td>-24.2</td><td>—</td><td>-45.7</td><td>339.1</td></tr><tr><td>Eigenmode 4</td><td>-24.0</td><td>257.3</td><td>-22.6</td><td>284.8</td></tr><tr><td>Eigenmode 5</td><td>-21.0</td><td>—</td><td>-21.8</td><td>—</td></tr></tbody></table>

complex conjugate eigenvalues represent damped oscillatory modes.

Nevertheless, the two spectra are different (Table 2): the least damped eigenmode for the MBC experiment corresponds to a purely damped mode, whereas for FBC it corresponds to a damped oscillatory mode with a period of 34 yr; the second least-damped eigenmode is a damped oscillatory mode (328-yr period) for the MBC experiment and a purely damped mode for the FBC experiment. In agreement with previous works (Huck and Vallis 2001; Arzel et al. 2006), this result shows the crucial importance of the physics of the surface boundary conditions, since we have computed the Jacobian matrix with strictly the same steady state.

The optimal excitation of these eigenmodes can be obtained by computing their biorthogonal (Farrell and Ioannou 1996a). This is done by diagonalizing the adjoint of the Jacobian matrix for both experiments. In a nonnormal system, the biorthogonals differ from the eigenmodes of the Jacobian matrix. This analysis is thus useful here to understand what is the most efficient pattern for exciting the least damped eigenmode of the Jacobian matrix. As the theory predicts, we obtain the same eigenvalue spectrum as for the Jacobian matrix, but for both experiments the structure of the eigenvectors is different (the respective least damped eigenmodes are shown in Figs. 3–5). Hence, as expected from advective–diffusive processes, the linear problem is nonnormal ($AA^† - A^†A \neq 0$, where $^†$ denotes the adjoint matrix).

The eigenmode analysis of the Jacobian matrix for the FBC experiments leads to a damped interdecadal oscillatory mode, which is fully described in the literature (Greatbatch and Zhang 1995; Chen and Ghil 1996; Huck et al. 1999a; Colin de Verdière and Huck 1999; te Raa and Dijkstra 2002). Dijkstra et al. (2006) have tried to reconcile the interdecadal variability found in idealized models under FBC and the observed Atlantic multidecadal oscillation (AMO; Delworth and Mann 2000; Knight et al. 2005); they conclude that these two variables are intrinsically the same, and thus the AMO is a mode of ocean variability. They also show that the

atmospheric damping decreases the growth rate of this mode, in agreement with the absence of the mode in our MBC experiment.

The adjoint eigenmode analysis reveals one property of this interdecadal oscillation: the least damped eigenmode of $A^†$ is strongly dominated by the salinity at the surface (Figs. 5g–j). This structure is the most efficient for stimulating the least damped eigenmode of $A$ (Farrell and Ioannou 1996a). That is, although this latter oscillatory eigenmode of $A$ is temperature dominated (Figs. 5a–d), it is more efficiently stimulated by a perturbation in surface salinity than in surface temperature. This first result shows the importance of the surface salinity perturbations for stimulating the multidecadal oscillatory mode in our model. This view is in agreement with previous studies showing the importance of the freshwater flux for the production of variability by the system itself (Weaver et al. 1991, 1993) and particularly for North Atlantic multidecadal variability (Chen and Ghil 1995; te Raa and Dijkstra 2003. Moreover, this result validates a posteriori the choice of our analysis to look at the optimal surface salinity perturbation.

The existence of a stable steady state with nonnormal dynamics motivates the following generalized stability analysis. We will now investigate the initial and stochastic surface salinity perturbations that induce the largest transient modification and the largest variance of the circulation, respectively.

## 4. Optimal initial surface salinity perturbation

We address here the first question: what initial perturbation in surface salinity influences the thermohaline circulation the most? As fully discussed in Sévellec et al. (2007), we choose a linear measure of the intensity of the MOC for this purpose. Following appendix A, we express the Lagrangian function as

$$ \mathcal{L}_{\text{ini}}(\tau) = \langle F | u(\tau) \rangle - \gamma [\langle u(0) | \mathbf{S} | u(0) \rangle - 1], \quad (6) $$

where $\langle F | u(\tau) \rangle$ is the intensity of the MOC at time $\tau$, $\mathbf{S}$ the matrix defining the quadratic norm of the state vector, and $\gamma$ the Lagrange multiplier associated with the quadratic norm constraint. The optimal initial surface salinity perturbation [using $|u(0)\rangle = \mathbf{P}|u_0\rangle$] follows as

$$ |u_0\rangle = \pm \frac{(\mathbf{P}^\dagger \mathbf{S} \mathbf{P})^{-1} \mathbf{P}^\dagger \mathbf{M}^\dagger(\tau) |F\rangle}{\sqrt{\langle F | \mathbf{M}(\tau) \mathbf{P} (\mathbf{P}^\dagger \mathbf{S} \mathbf{P})^{-1} \mathbf{P}^\dagger \mathbf{M}^\dagger(\tau) |F\rangle}}. \quad (7) $$

These optimal perturbations are computed as a function of the time delay $\tau$ for which the Lagrangian and $\langle F | u(\tau) \rangle$
<figure><img src="image_4.png" /><figcaption>FIG. 3. Linear stability analysis for MBC experiment. The least damped eigenmode of the (a),(c),(e) Jacobian matrix and its (b),(d),(f) adjoint: (a),(b) surface temperature and (c),(d) salinity, both in terms of density; (e),(f) associated MOC; and (g) eigenvalue spectrum. Contour intervals (CIs) are 0.1, 0.1, and $5 \times 10^3$ Sv for (a),(c),(e) and 0.005, 0.1, and $5 \times 10^4$ Sv for (b),(d),(f), respectively.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4. Same as in Fig. 3, but for the second least-damped eigenmodes of the (left) linear tangent matrix and its (right) adjoint. The eigenvectors are a complex pair. The real part of the eigenvectors of the linear tangent (adjoint) matrix is chosen to maximize the MOC (minimize the density norm); the imaginary part is one quarter period later (82 yr). CIs are (a)–(f) 0.01, 0.02, and $10^3$ Sv and (g)–(l) 0.001, 0.1, and $10^5$ Sv, respectively.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIG. 5. Same as in Fig. 4, but for the least damped eigenvectors of the FBC experiment. The imaginary part is one quarter period later (8.5 yr). CIs are (a)–(f) 0.1, 0.02, and $10^3$ Sv and (g)–(l) 0.005, 0.1, and $10^5$ Sv, respectively.</figcaption></figure>

(cost function) are maximized, for both experiments MBC and FBC (Figs. 6, 7).

The maximum response of the MOC intensity perturbation for a normalized perturbation, $\langle F|u(\tau) \rangle / \langle u(0)|S|u(0) \rangle$, as a function of the time delay are displayed in Figs. 6b, 7b. The comparison of the two experiments reveals a sharp difference. For the MBC experiment, the cost function strongly decays for delays increasing from 1 to 50 yr, and a weak maximum appears for a delay of 111 yr. For the FBC experiment, the cost function decreases slowly with the delay but shows decadal oscillations. The order of magnitude of the MOC perturbations is significantly different too for delays longer than 10 yr. A local maximum exists for both experiments. Because of the small spatial scale of the measure ($|F\rangle$ is a Dirac function in latitude), the system needs some time to propagate the information and select a large spatial scale (climatically relevant) through diffusion. Note that, as tested in Sévellec et al. (2007), the use of a smoother function in latitude also shows a large response on rapid time scale. To study the climatically relevant optimal perturbations, we choose to study the local maxima corresponding to a modification

after a delay of 111 yr for the MBC and 24 yr for the FBC. The nonlinear direct time integrations perturbed by the two corresponding optimal patterns show transient growth of the MOC 111 and 24 yr, respectively, after the perturbation (Figs. 6c, 7c).

### a. Mixed boundary conditions

For the MBC experiment, the optimal perturbation corresponds to a meridional gradient of the surface salinity, with a weak incursion of freshwater along the western boundary (Fig. 6a). Because of the thermal wind relation, this perturbation, corresponding to a weak zonal gradient of salinity, is associated with a negative anomaly of the MOC intensity, corresponding to a southward surface flow. This structure has a weak correspondence to the surface salinity part of the least damped eigenmodes of $A^\dagger$ but a strong correspondence with the one of the second least-damped eigenmodes (Fig. 3d and Figs. 4i,j). Indeed, the normalized projections [Proj($|X_1\rangle$, $|X_2\rangle$) = $\langle X_1|X_2 \rangle / \sqrt{\langle X_1|X_1 \rangle \langle X_2|X_2 \rangle}$] between the two SSS are 0.21 and 0.94, respectively (Table 3). The strong projection shows the importance of the second least-damped eigenvectors of $A^\dagger$ to understand
<figure><img src="image_7.png" /><figcaption>FIG. 6. (a) Optimal initial SSS perturbation inducing the largest MOC response after 111 yr for the MBC experiment. The pattern is rescaled to 1 psu maximum and induces a MOC response of 0.13 Sv. CI is 0.05 psu. (b) Maximum response of the MOC intensity [⟨F|u(τ)⟩ / ⟨u(0)|S|u(0)⟩]: a local maximum appears for a delay of 111 yr. The apparently large response is due to the choice of the norm S. (c) Time evolution of the MOC intensity for the nonlinear model initialized by the optimal initial SSS perturbation: a local maximum of the MOC appears after 111 yr (vertical line). The time decaying envelope (dashed) corresponds to the single-mode approximation (12). The integration is made with a +0.01 psu perturbation, but the plot is rescaled to an equivalent perturbation of 1 psu.</figcaption></figure>

the optimal stimulation of the MOC intensity. A strong
projection (0.89) exists also with the fourth eigenmodes
of A†. The optimal surface salinity perturbation evolves
during 111 yr, and the snapshots of SST, SSS, zonal
overturning circulation (ZOC), and MOC at 111 yr (Fig.
8) show that the initial SSS structure has disappeared
and a strong meridional gradient of SST anomaly has
taken its place. This new density structure corresponds
to a positive meridional circulation that reinforces the

circulation of the steady state. The structure of the
anomaly has a strong normalized projection on one phase
of the second least-damped linear eigenmode (0.85; Fig.
8 vs. Figs. 4a,c,e; Table 3). The modification of the MOC
intensity induced by this phase corresponds to 87% of
the maximum value. The projection with the fourth
least-damped eigenmode is significantly less important
(0.42; Table 3). This analysis shows that several eigen-
modes are stimulated by the optimal initial SSS. But the

<figure><img src="image_8.png" /><figcaption>FIG. 7. Same as in Fig. 6, but for the FBC experiment. (b) A local maximum appears for a delay of 24 yr (vertical line), which we choose as the optimal time. (c) The time decaying envelope (dashed) corresponds to the single-mode approximation (14).</figcaption></figure>
<table><caption>TABLE 3. Normalized projection [Proj(|X₁⟩, |X₂⟩) = ⟨X₁|X₂>/√⟨X₁|X₁⟩⟨X₂|X₂⟩] between the optimal initial SSS perturbations or their optimal response and the eigenmodes of A† and A, respectively (in case of complex eigenmodes, the phase leading to the strongest projection is used). Strong projections (≥0.8) are boldfaced.</caption><thead><tr><th>MBC</th><th>Proj[|SSS(0)⟩,|SSS†⟩]</th><th>Proj[|u(111 yr)⟩,|u<sub>i</sub>⟩]</th></tr></thead><tbody><tr><td>Eigenmode 1</td><td>0.21</td><td>0.51</td></tr><tr><td>Eigenmode 2</td><td><b>0.94</b></td><td><b>0.85</b></td></tr><tr><td>Eigenmode 3</td><td>-0.29</td><td>0.32</td></tr><tr><td>Eigenmode 4</td><td><b>0.89</b></td><td>0.42</td></tr><tr><td>Eigenmode 5</td><td>-0.10</td><td>-0.30</td></tr><tr><td>FBC</td><td>Proj[|SSS(0)⟩,|SSS†⟩]</td><td>Proj[|u(24 yr)⟩,|u<sub>i</sub>⟩]</td></tr><tr><td>Eigenmode 1</td><td><b>0.96</b></td><td><b>0.95</b></td></tr><tr><td>Eigenmode 2</td><td>0.01</td><td>-0.46</td></tr><tr><td>Eigenmode 3</td><td>0.64</td><td>0.49</td></tr><tr><td>Eigenmode 4</td><td>0.27</td><td>0.53</td></tr><tr><td>Eigenmode 5</td><td>0.01</td><td>0.44</td></tr></tbody></table>

structure of the perturbation, when the intensity of the MOC reaches its maximum, is mainly controlled by the second least-damped eigenmodes of **A**.

To evaluate the importance of the least damped eigenmodes (of **A** and of **A†**), we can expand **A** as in Sévellec et al. (2007):

$$ \mathbf{A} = \sum_k |u_k\rangle \lambda_k \langle u_k^\dagger|, \quad (8) $$

where $|u_k\rangle$ and $|u_k^\dagger\rangle$ are the eigenvectors of **A** and **A†**, and $\lambda_k$ the eigenvalues. We can write the propagator as

$$ \mathbf{M}(t) = \sum_k |u_k\rangle e^{\lambda_k t} \langle u_k^\dagger|. \quad (9) $$

To isolate the influence of the pair of the second least-damped eigenmodes, we approximate the propagator as

$$ \mathbf{M}(t) \approx |u_2\rangle e^{\lambda_2 t} \langle u_2^\dagger| + \text{c.c.} \quad (10) $$

Using $|u_0\rangle$ as the optimal initial SSS, the intensity of the MOC then reads as

$$ \langle F | \mathbf{M}(t) | u_0 \rangle \approx \langle F | (|u_2\rangle e^{\lambda_2 t} \langle u_2^\dagger| + \text{c.c.}) | u_0 \rangle. \quad (11) $$

Decomposing the conjugate eigenvalues as the sum of their real and imaginary parts ($\lambda_2 = \lambda_{2r} + i\lambda_{2i}$), we obtain

<figure><img src="image_9.png" /><figcaption>FIG. 8. Snapshots at 111 yr of the nonlinear time integration initialized by the optimal SSS perturbation for MBC: anomalies of surface (a) temperature and (b) salinity, and (c) ZOC and (d) MOC. The initial perturbation is scaled for a +0.01-psu maximum amplitude; however, the figures are rescaled for a +1-psu maximum amplitude perturbation. CIs are 5 × 10⁻⁴ °C, 5 × 10⁻⁴ psu, 0.02 Sv, and 0.02 Sv, respectively.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 9. Transient MOC modification for the MBC experiment initialized by the optimal SSS perturbation. (a) North–south vertically averaged density difference (dashed and dotted lines represent the thermal and salinity differences in terms of density: $\rho_0\alpha T_{N-s}$, $\rho_0\beta S_{N-s}$, respectively). (b) MOC anomaly. The vertical solid lines show the time when the MOC intensity reaches its maximum. The vertical dashed lines separate the phases dominated by temperature and salinity. The initial perturbation was scaled for a +0.01 psu maximum amplitude; however, the figures are rescaled for a perturbation amplitude of +1 psu.</figcaption></figure>

$$ \langle F | \mathbf{M}(t) | u_0 \rangle \approx [\langle F | (|u_2| e^{i\lambda_2 t} \langle u_2^\dagger | + \text{c.c.}) | u_0 \rangle] e^{\lambda_2 t} \quad \text{and} \\ |\langle F | \mathbf{M}(t) | u_0 \rangle| \le 2 |\langle F | u_2 \rangle| |\langle u_2^\dagger | u_0 \rangle| e^{\lambda_2 t}. \qquad (12) $$

This last relation measures the influence of the second eigenmode exponential decay on the intensity of the MOC. This approximation bounds the variations of the circulation and stresses the predominance of the second least-damped eigenmode at 111 yr (Fig. 6). This analysis shows that using only the second eigenmodes of **A** and **A**† is a good approximation in this experiment for two different reasons. On the one hand, some eigenmodes are not stimulated by the optimal SSS perturbation (eigenmode 1, 3, or 5; Table 3). On the other hand, some eigenmodes (such as eigenmode 4), even though they are stimulated, are too strongly damped on the time scale of the maximization (111 yr). Because of these eigenmode, the bound (12) is incorrect for short time scales (Fig. 6).

The time evolution of the optimal initial SSS perturbation can be conceptually separated into three phases (Fig. 9). The first phase (less than 1 yr) corresponds to

the rapid adjustment of the circulation to the initial SSS perturbation. During the second phase (around the first 20 yr) the evolution of the large-scale density gradient is due to the evolution of the temperature. Finally, during the third phase, the advection of the salinity anomaly by the mean flow sets the large-scale density gradient and thus the optimal circulation response. The two last phases could be approximated as a large-scale modification of the north–south density gradient. We thus define

$$ X_{N-s} = \bar{X}^{\text{North}} - \bar{X}^{\text{South}}, $$

where $\bar{X}^{\text{South}}$ and $\bar{X}^{\text{North}}$ denote the average value of the scalar field $X$ on the south and north part of the basin (from 10° to 35°N and from 35° to 60°N over the whole depth and longitudes), respectively. We further assume that on the decadal time scale, the north–south density difference controls the MOC intensity (Stommel 1961). However, as a result of the 3D dynamics of our model, a delay will appear between the modifications of $\rho_{N-s}$ and
the MOC intensity. This delay corresponds to the east-west density adjustment by Rossby waves. During the first phase, the circulation is adjusted on a short time scale (<1 yr) to the initial salinity perturbation ($S'_{N-S} < 0$) and leads to a negative anomaly of the MOC (MOC' <0). After this initial shock, the temperature anomaly and the circulation evolve during ~20 yr as predicted by the temperature negative feedback on the circulation (Marotzke 1996): the negative anomaly of the MOC induces a negative anomaly of the meridional heat transport. The latter produces a positive anomaly of temperature in the south and a negative anomaly in the north ($T'_{N-S} < 0$). These temperature anomalies progressively reduce the negative MOC anomaly intensity until the meridional density difference is close to zero. After this temperature feedback, which corresponds to the large-scale adjustment of temperature anomaly on salinity, the third phase takes place. It consists in the advection of the salinity dominated anomaly by the mean flow. Actually, the mean flow advection moves the north-south negative density anomaly around the MOC cell and induces a north-south positive density anomaly. This last mechanism induces the modification of the MOC intensity at 111 yr. In this experiment, the optimal response mechanism is similar to the one found in Sévellec et al.'s (2007) 2D latitude-depth model under MBC, where both the optimal time scale (67 yr) and linear oscillatory mode (150 yr) periods were half as long as here.

## b. Flux boundary conditions

For the FBC experiment, we analyze the optimal surface salinity inducing the maximal modification after 24 yr. Its structure corresponds to a salinity minimum intensified in the north of the basin, with an incursion in the west of the basin. The zero anomaly is oriented from the southwest to the northeast (Fig. 7a). This structure is strongly linked to the one of the least damped eigenvectors of $A^\dagger$ (Figs. 5i,j). We obtain a normalized projection of the surface salinity component of 0.96 (Table 3). This projection shows the importance of the least-damped oscillatory eigenmode of $A^\dagger$ in the determination of the optimal surface salinity perturbation (the other eigenmodes are also stimulated, as shown by their nonzero normalized projection on the optimal SSS but not as efficiently as the least damped eigenmodes; Table 3). The nonlinear time integration perturbed by this optimal pattern shows that the MOC intensity has a maximum at 24 yr. Indeed, the initial perturbation induces a weak negative MOC anomaly that evolves into a positive anomaly, reaching a maximum after 24 yr (Fig. 7c). The associated patterns of SST, SSS, ZOC, and MOC anomalies (Fig. 10) are very close to those of

the least-damped oscillatory eigenmodes of $A$ at its phase, inducing the maximum intensity of the MOC (Figs. 5a,c,e). Indeed, we obtain a strong normalized projection for one phase (0.95; Table 3). This phase induces a MOC intensity modification that corresponds to 99% of the maximal modification value. This result is coherent with a decrease of the other eigenmodes of $A$ in 24 yr and the appearance of a transient response through the least damped eigenmode. Note that the choice of the phase is only due to the choice of the function we maximize, but it is definitely not due to the choice of the quadratic norm; the latter is only introduced for normalization purposes (which is needed to remove the degeneracy of a maximization problem in a linear framework).

To isolate the influence of the pair of the least damped eigenmodes, we approximate the propagator by the relation (which is asymptotically correct)

$$M(t) \approx |u_1\rangle e^{\lambda_1 t} \langle u_1^\dagger| + \text{c.c.} \quad (13)$$

We thus obtain

$$ \begin{aligned} \langle F | M(t) | u_0 \rangle &\approx [\langle F | (|u_1\rangle e^{i\lambda_1 t} \langle u_1^\dagger| + \text{c.c.}) |u_0\rangle] e^{\lambda_1 t} \quad \text{and} \\ |\langle F | M(t) | u_0 \rangle| &\le 2 |\langle F | u_1 \rangle| |\langle u_1^\dagger | u_0 \rangle| e^{\lambda_1 t}. \end{aligned} \quad (14) $$

This last relation measures the influence of the optimal initial perturbation on the intensity of the MOC modulated by the exponential decay as a result of the oscillatory eigenmode. This approximation bounds the variations of the circulation and stresses the predominance of the least damped eigenmode at 24 yr (Fig. 7). Despite the fact that the intensity of the response is controlled by the least damped eigenmode of $A^\dagger$, this analysis, compared to the perturbation propagated by the full equations (Fig. 7c), also shows that the MOC modification is mainly the result of the change of phase of the oscillatory eigenmode of $A$ between the initial time and after 24 yr. Indeed, each phase of the least damped eigenmodes does not have the same influence on the circulation intensity.

In this experiment, the response of the MOC intensity is strongly linked to the weakly damped oscillatory eigenmode, as discussed previously and shown in Fig. 11. Because of the FBC, the model state, which is defined by both temperature and salinity variables, can be described by the single-density variable. The two equations of evolution for temperature and salinity [Eqs. (2a,b)] can be unified in only one for the density [by application of (2c)]. This last point differs from the previous experiment, where temperature and salinity act separately to induce the transient modification. Moreover, the
<figure><img src="image_11.png" /><figcaption>FIG. 10. Same as in Fig. 8, but for after 24 yr of the nonlinear time integration perturbed by the optimal SSS for FBC. CIs are 0.05°C, $5 \times 10^{-3}$ psu, 0.2 Sv, and 0.2 Sv, respectively.</figcaption></figure>

mechanism here involves the full 3D dynamics in contrast with the analysis conducted above. The initial perturbation induces a density anomaly such that $\rho'_{N-S} < 0$. By the thermal wind relation, it is associated with a reduction of the ZOC ($ZOC' < 0$, where the prime denotes the anomaly). In presence of a stable vertical gradient of the mean density, the negative anomaly of the ZOC induces a decrease of the absolute value of the zonal density transport ($ZDT' > 0$, $\overline{ZDT} < 0$). This affects the zonal gradient of density and induces $\rho'_{W-E} < 0$ (this latter notation defines the difference between the density anomalies averaged over the western and eastern halves of the basin). By the thermal wind relation, the zonal gradient of the density anomaly $\rho'_{W-E} < 0$ is associated with a reduction of the MOC ($MOC' < 0$) and an increase of the meridional density transport ($MDT' > 0$). The latter affects the meridional density gradient and leads to the opposite phase from the initial perturbation $\rho'_{N-S} > 0$. This chain of events repeats itself, leading to an oscillation, where the delay between $\rho'_{N-S}$ (and its associated ZOC anomaly) and $\rho'_{W-E}$ (and its associated MOC anomaly) sustains this oscillatory mode, as fully described in the Raa and Dijkstra (2002). Through this oscillatory mechanism, the initial $\rho'_{N-S} < 0$ leads to an intensification of the MOC after 24 yr, as expected by the adjoint analysis.

### c. Discussion

We now discuss the range of validity of the linear approximation by computing the nonlinear time integrations perturbed by the optimal SSS pattern for a range of perturbation amplitudes from -2 psu to +2 psu for both experiments (Fig. 12). The optimal MOC responses for the linear model overestimate the nonlinear model response by 15% for the MBC experiment and only 30% for the FBC experiment. The larger difference in the former case may be simply due to the longer time integration for the optimal response (111 yr) compared to the latter case (24 yr), since the differences between the linear and nonlinear model will cumulate with time. Note that this does not validate the optimality of the initial SSS perturbations but only the weak sensitivity of their evolution to neglected nonlinear terms. Indeed, we cannot verify if the optimal perturbation in the linear framework remains optimal for the full nonlinear problem (Mu and Zhang 2006).

Finally, we can compute an upper bound of variability for both experiments from the optimal initial perturbations. Indeed, since we obtained the optimal SSS perturbation, we know that no other perturbation can lead to stronger modification of the MOC intensity in the linear limit. Using a typical amplitude for the Great
<figure><img src="image_12.png" /></figure>

FIG. 11. Transient MOC modification for the FBC experiment initialized by the optimal SSS perturbation. (a) North-south and (b) west-east vertically averaged density differences. (c) ZOC and (d) MOC anomalies. The vertical solid lines show the time when the MOC intensity reaches its maximum. The initial perturbation is scaled for a +0.01-psu maximum amplitude; however, the figures are rescaled for an amplitude of +1-psu perturbation.

Salinity Anomalies (GSA) of -0.5 psu over a 250-m-thick layer located in the north region (Belkin et al. 1998) to set the intensity of the SSS perturbation, we obtain an upper bound of the MOC intensity modification of +0.16 Sv—that is, 1% of the mean circulation—for the MBC experiment and +2.25 Sv—that is, 12.5% of the mean circulation—for the FBC experiment.

For these two experiments, the optimal perturbations show a large difference in terms of both intensity and time-scale responses and thus of the transient response mechanism. These differences appear despite identical steady states. To develop more ideas along this line, we propose to focus now on stochastic perturbations. We will thus be able to quantify how the surface temperature boundary conditions can alter the response of a stable oceanic circulation not only in the transient regime but also in an asymptotic regime.

## 5. Optimal stochastic surface salinity flux

In this section, we will consider the variability of the oceanic circulation induced by a stochastic perturbation

of the surface freshwater flux (SSF). Since the oceanic adjustment time scale is much longer than the synoptic atmospheric forcing time scales, the stochastic noise is approximated by a white noise (expected value is $E[a(t)] = 0 \text{ yr}^{-1/2}$ and variance is $E[a^2(t)] = 1 \text{ yr}^{-1}$; appendix B). As in the previous section, we look for the optimal SSF in two experiments corresponding to different surface temperature boundary conditions. The optimal stochastic perturbations are obtained following the methodology summarized in appendix B.

### a. Mixed boundary conditions

For the MBC experiment, the optimal stochastic perturbation corresponds to a zonal dipole, with the negative part strongly intensified in the northwest of the basin and a positive part in the high latitudes but in the east of the basin (Fig. 13a). This perturbation is mainly located in the north region, extending over 20° around the location of the meridional streamfunction maximum. The spectrum of the response to the stochastic forcing corresponds to a red noise, as expected for the ocean acting as an integrator of an atmospheric
<figure><img src="image_13.png" /><figcaption>FIG. 12. Validation of the linear approximation for the (a) MBC and (b) FBC experiments: linear theoretical values of the response of the MOC intensity (solid line) are compared to the nonlinear model time integrations (+) for perturbation amplitudes from -2 to +2 psu. The relative error for GSA is less than 15% (30%) for the MBC (FBC) experiment. Note the much larger response in the latter case.</figcaption></figure>

white-noise stimulation (Frankignoul and Hasselmann 1977) but with a small bump around 300 yr (Fig. 14a). This bump is the signature of the least-damped oscillatory eigenmode of the Jacobian matrix and its adjoint, which period is about 329 yr (Table 2). However, its strong damping (with a 58-yr e-folding time scale much shorter than its period) only induces a weak signature on the spectrum.

## b. *Flux boundary conditions*

For the FBC experiment, the optimal stochastic perturbation is a large-scale zonal dipole (Fig. 13b). The

shape of this optimal stochastic perturbation is very close to the one of the optimal initial perturbation for the same experiment (Fig. 7a). The perturbation is intensified in the northwest of the basin, with a weak incursion along the western boundary. The perturbation in the rest of the basin has an opposite sign to conserve the total salt amount. The standard deviation of the MOC in this experiment is 4 times as large as in the MBC experiment, and the spectrum of the MOC intensity response to this stochastic forcing is significantly different (Fig. 14b). Moreover, on top of the red-noise shape, a strong peak now appears at 35 yr. The peak

<figure><img src="image_14.png" /><figcaption>FIG. 13. Optimal stochastic perturbations for the (a) MBC and (b) FBC experiments. The maximum of the SSF structure is normalized to 1 psu yr⁻¹/². Perturbations induce a variability of the circulation intensity with a std dev of 3.6 and 14.5 Sv, respectively; CI is 0.05 psu yr⁻¹/².</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>FIG. 14. Spectrum of the response to the optimal stochastic perturbations for the (a),(b) MBC and (c)–(d) FBC experiments. (a),(c) Response of the MOC intensity during a nonlinear time integration perturbed by the stochastic optimal. The integrations were made with a std dev of +0.01 psu yr⁻¹ anomalies, but the plots are rescaled for 1 psu yr⁻¹. (b),(d) Comparison between the theoretical spectra (thick line) and those from the nonlinear time integrations (gray line). The vertical gray lines represent the CI of the time integration spectra. The shapes are in good agreement and a weak underestimation of the power spectral intensity also appears.</figcaption></figure>

corresponds to the least-damped oscillatory eigenmodes of the Jacobian matrix and its adjoint, which have a period of 34 yr and a decaying time scale one order of magnitude longer (207 yr). This mode can have several oscillation cycles before disappearing and thus generates a strong spectral signature. The increase of the power spectral density between the two experiments appears through the multidecadal weakly damped eigenmode, which is not significant in the MBC experiment but is preponderant in the FBC.

### c. Discussion

As in the previous section, we try to estimate the range of validity of our linear approximation. Several time integrations of the nonlinear model perturbed by the optimal stochastic perturbations modulated by different intensity of the standard deviation are performed (Fig. 15). For a range of standard deviations from 0 to 0.1 psu yr⁻¹ (equivalent to 28.6 cm yr⁻¹ in freshwater flux), we obtain a maximum error of 13% for the MBC experiment and 10% for the FBC experiment. Since the linearization of the convection, a strongly nonlinear process, is only approximated, it is a likely candidate for producing these errors. The optimal stochastic perturbation for the MBC experiment is more localized in the north of the basin than for the FBC experiment, that is, in the weak stratification zones where convection is active. This difference in the location where the two optimal perturbations are maximum may be an explanation for the difference in the error between the two experiments.

Now that the error is estimated, we can compute an upper bound of the multidecadal time-scale variability induced by stochastic surface salinity flux perturbations from these optimal stochastic perturbations. A typical value of the standard deviation of the stochastic noise at midlatitude is of the order of 20 cm yr⁻¹ (estimated from the interannual standard deviation at midlatitudes in National Centers for Environmental Prediction (NCEP) and 40-yr European Centre for Medium-Range Weather Forecasts Re-Analysis (ERA-40) re-analyses. Moreover the variance induced by the NAO, in terms of freshwater flux, is of the same order of magnitude (Mignot and Frankignoul 2003). Such noise induces a maximum standard deviation of the variability of the MOC intensity of 0.25 Sv (1.5% of the mean circulation) for the MBC experiment and 1 Sv (5.5% of the mean circulation) for the FBC experiment.

As in the optimal initial perturbation experiments, both experiments show a large difference in the intensity of their responses, although their optimal perturbation patterns are really close. To clarify this point, we have used the optimal stochastic of MBC experiment to force the FBC model, and vice versa. These changes of perturbation pattern do not at all affect the shape of the spectrum response, only its intensity (the relative change is ~50%). This last experiment confirms that the choice of the surface boundary condition for temperature is more crucial for the amplitude of the stochastic variability than the details of the stochastic pattern.
<figure><img src="image_16.png" /><figcaption>FIG. 15. Same as in Fig. 12, but for the variance of the MOC intensity. The relative error for typical NAO noise is less than 13% (10%) for the MBC (FBC) experiment. Note again the much larger variance in the latter case.</figcaption></figure>

## 6. Conclusions

Global warming is likely to induce changes in the global water cycle (Held and Soden 2006), such as massive melting of glaciers at high latitudes (i.e., Greenland), or a modification of the weather regime characteristics (i.e., NAO), and thus of high-frequency freshwater forcing. We have investigated the influence of surface salinity (flux) modifications on the ocean circulation, paying special attention to the thermohaline circulation, the slow component of the ocean and thus of the climate system, transporting heat poleward. In this study, we have addressed two types of optimal sea surface salinity perturbations: initial perturbations, inducing the largest change in the MOC intensity, and stochastic perturbations, inducing the largest variability of the MOC. In this framework, we have focused our attention on the importance of the surface boundary conditions for temperature (i.e., restoring term vs. flux forcing) and how it could modify the optimal solutions and the associated ocean response.

In agreement with previous studies (Huck and Vallis 2001; Arzel et al. 2006), our results stress the importance of the surface boundary conditions on the pattern of the optimal forcing and the intensity of the MOC response. We have compared two experiments corresponding to mixed boundary conditions (MBC), where surface temperature is restored to prescribed fields, and constant flux boundary conditions (FBC); in both cases, the surface salinity is forced by a constant freshwater flux. Although we have applied our computation on strictly the same stable steady state for both experiments and the optimal initial or stochastic perturbation patterns are qualitatively similar, the response of the ocean model to the

different perturbations show large differences. For the optimal initial perturbation, the response of the MOC intensity for MBC is 14 times weaker than for FBC. For the optimal stochastic perturbation, the standard deviation of the MOC for MBC is 4 times weaker than for FBC. These results can be interpreted in terms of the eigenmodes provided by linear stability analysis. In the FBC experiment, the most weakly damped eigenmode is a 34-yr period oscillatory mode with a 207-yr decaying time scale; this mode (and its biorthogonal) mainly explains the variability triggered by initial and stochastic perturbations, but such a strong variability does not occur in the MBC experiment, where the most weakly damped eigenmode of the Jacobian matrix is purely real. This dynamical difference is fundamental for the results of the two experiments, although the steady state is strictly the same: it raises the importance of the good representation of ocean-atmosphere interactions and feedbacks to reproduce realistically the variability and sensitivity in ocean numerical models.

The second result is the computation of an upper bound on the ocean response to perturbations, here in terms of the MOC. We have found that this upper bound is crucially dependent on the surface boundary conditions. We use the FBC to compute an upper bound, because this experiment induces the largest responses. For an initial perturbation of the scale of the Great Salinity Anomalies (–0.5 psu in the upper 250 m located in the northern region), we obtain a circulation intensity increase of +2.25 Sv (i.e., 12.5% of the mean circulation) after 24 yr. In a realistic global configuration of an ocean general circulation model under FBC, Sévellec et al. (2008) found a –0.8-Sv response of the
overturning 10.5 yr after a typical GSA anomaly amplitude optimal perturbation was introduced, corresponding to 11% of the mean overturning. Although the optimal pattern and relative response amplitude are in good agreement, the time scale and the response sign suggest that the multidecadal linear mode, if it exists, is more heavily damped in the global model. Consequently, the response to the SSS perturbation would correspond better to the first 10 yr of the present model response, before the MOC changes sign (as in the MBC experiment). For a stochastic perturbation of the amplitude of the estimated standard deviation of interannual freshwater flux at midlatitudes (20 cm yr⁻¹), we obtain a standard deviation of the circulation intensity of 1 Sv (5.5% of the mean circulation).

The third result is the importance of the multidecadal oscillation in the optimal nonnormal response, as in the optimal stochastic variability. These results reveal the nonnormal behavior of the multidecadal mode, which could be responsible for the Atlantic multidecadal oscillation according to Dijkstra et al. (2006) and its effect on the MOC. As shown in the analysis of the FBC experiment, the nonnormal response of the circulation intensity as a result of the optimal initial perturbation can be mainly explained (87%) by the existence of the linear weakly damped multidecadal oscillation and its biorthogonal. In the same way, the optimal stochastic perturbation induces a peak at the period of this multidecadal oscillation. Moreover, this 35-yr peak dominates the spectrum, which corresponds to a red noise for other frequencies. Our study highlights the nonnormal response of a stable ocean circulation to atmospheric white noise-inducing Atlantic multidecadal variability.

One methodological outlook of this work is the application of the full nonlinear equation model in the search of optimal perturbations. Such methods, known as conditional nonlinear optimal perturbations, were performed in several recent papers (Mu et al. 2004; Sun et al. 2005; Mu and Zhang 2006) and show that the linear behavior is a strong limitation in the computation of optimal perturbations. Given the influence of surface boundary conditions on the optimal patterns, time scales, and oceanic response, we have found it would be natural to extend these analyses in coupled ocean-atmosphere models, for instance, as performed recently by Tziperman et al. (2008). Another continuation of the present work could be the study of the optimal perturbation but with a mean state that evolves through the seasonal cycle (nonautonomous system). This would show the eventual differences in the sensitivity of the circulation, depending on the season at which perturbations are introduced and point out which processes control this sensitivity (like the mixed layer depth). And

finally, a strong forcing of the ocean circulation that has not been considered here is the wind—looking at its optimal patterns in the same manner is underway.

*Acknowledgments.* The lead author performed part of this work while funded by the ENSEMBLE-FP6 European project. We gratefully acknowledge comments on a preliminary version of the paper (part of F. Sévellec's Ph.D. thesis) by Andrew M. Moore and Michael Ghil. Some of the computations have been performed on the NEC vectorial machine at the CNRS national computing center (IDRIS, Orsay, France).

## APPENDIX A

### Constrained Maximization Method

The constrained maximization method, using Lagrange multipliers, which is applied in this study, has been fully described in Sévellec et al. (2007). Following the same notations, we briefly develop the solution. We are looking for the state vector ($|u\rangle$), which under $n$ constraints [$C_i(|u\rangle)$, for $i = 1, \dots, n$] maximizes a scalar quantity [$G(|u\rangle)$]. We then introduce the Lagrangian:

$$ \mathcal{L}(|u\rangle, \gamma_i) = G(|u\rangle) - \sum_{i=0}^{n} \gamma_i C_i(|u\rangle), \quad (A1) $$

where $\gamma_i$ are the Lagrange multipliers associated with the constraints. The state vector maximizing $G$ under the constraints thus verifies

$$ \frac{dG(|u\rangle)}{d|u\rangle} - \sum_{i=0}^{n} \gamma_i \frac{dC_i(|u\rangle)}{d|u\rangle} = 0 \quad \text{and} \quad (A2) $$

$$ C_i(|u\rangle) = 0, \quad \forall i \in [1, n]. \quad (A3) $$

We choose to maximize the intensity of the MOC, because it is a good indicator for the ocean poleward mass and heat transport, with a potential effect on North Atlantic climate. This scalar index can be written as a projection of the state vector on the bra $\langle F|$ . More precisely, we define $\langle F|U\rangle$ as the value of the MOC streamfunction at the latitude and depth of its maximum. The linearization induces $\langle F|U\rangle = \langle F|\bar{U}\rangle + \langle F|u\rangle$ and then $G(|u\rangle) = \langle F|u\rangle$ defines the function of the state vector, leading to the scalar we want to maximize (also called the cost function).

In this study we use three constraints. First, the degeneracy of a maximization in a linear context is required to normalize the initial perturbation. We choose the quadratic density norm expressed as $\langle u(0)|S|u(0)\rangle = \Sigma_i[(\alpha^2 T'_i (0)^2 + \beta^2 S'_i(0)^2)v_i] / \Sigma_i(v_i) = 1$, where $T'_i$ and $S'_i$ are the first two components of the eigenvector of $S$ corresponding to the $i$-th mode.
APPENDIX B

Optimal Stochastic Perturbation

First, we write down the equation of the perturbation evolution under stochastic forcing:

$$d_t|u(t)\rangle = \mathbf{A}|u(t)\rangle + a(t)|g\rangle, \quad (B1)$$

where $a(t)$ is the stochastic temporal part of the forcing and $|g\rangle$ is the time-independent spatial pattern. Its time integration leads to

$$|u(\tau)\rangle = \mathbf{M}(\tau)|u(0)\rangle + \int_{0}^{\tau} ds \mathbf{M}(\tau-s)a(s)|g\rangle. \quad (B2)$$

Without loss of generality for asymptotic results, we assume $|u(0)\rangle = 0$. The variation of the forcing is defined as a white noise (expected value is $E[a(t)] = 0 \text{ yr}^{-1/2}$ and variance is $E[a^2(t)] = 1 \text{ yr}^{-1}$). Thus, the autocorrelation is defined by the classical Dirac delta function $E[a(t)a(t')] = \delta(t - t')$, where E denotes the expected value. Thus, we can rewrite the perturbation time dependency as

$$|u(\tau)\rangle = \int_{0}^{\tau} ds \mathbf{M}(\tau - s)a(s)|g\rangle. \quad (B3)$$

We now seek the optimal pattern inducing the maximal variance of MOC intensity; thus, the Lagrangian function becomes

$$\mathcal{L}_{\text{sto}}(\tau) = \text{Var}[\langle F | u(\tau) \rangle] - \gamma(\langle g | \mathbf{S} | g \rangle - 1). \quad (B4)$$

This variance can be written as

$$ \begin{aligned} \operatorname{Var}[\langle F | u(\tau) \rangle] &= \operatorname{E}[\langle F | u(\tau) \rangle^2] - \langle F | \operatorname{E}[u(\tau)] \rangle^2 = \operatorname{E}[\langle F | u(\tau) \rangle^2], \\ &= \langle g | \int_0^\tau ds \int_0^\tau ds' \operatorname{E}[a(s)a(s')] \mathbf{M}^\dagger(\tau - s) |F\rangle \langle F | \mathbf{M}(\tau - s') |g\rangle, \quad \text{and} \\ &= \langle g | \int_0^\tau ds \mathbf{M}^\dagger(\tau - s) |F\rangle \langle F | \mathbf{M}(\tau - s) |g\rangle. \end{aligned} \quad (B5) $$

are the temperature and salinity anomaly at the grid points $i$ and $v_i$ their relative volume. This constraint is written explicitly, that is, associated with a Lagrange multiplier, $C(|u\rangle) = \langle u(0)|\mathbf{S}|u(0)\rangle - 1 = 0$. Note that a quadratic norm is used to measure the amplitude of the initial perturbation. This is needed, since a linear measure of the initial perturbation would result in a degenerate problem. As a result, the quantity we measure is not MOC growth [$\langle F|u(t) \rangle / \langle F|u(0) \rangle$] but the quantity [$\langle F|u(t) \rangle / \langle u(0)|\mathbf{S}|u(0) \rangle$] that we will name MOC response. Since the numerator and denominator do not have the same units, the value of this quantity does not denote an actual growth of perturbation; thus one should not be surprised by the large values of this quantity found in this paper. The second constraint ensures the salt conservation, as in the full model equation. Indeed, since we want to see a transient effect, we need to verify strictly the conservation property of the set of equations. Otherwise, we would allow a drift in the total salt content that would not verify anymore the same steady state. The last and third constraint allows perturbations in surface salinity only, which remains coherent with our scientific aim (surface salinity perturbation influencing the ocean circulation). These two last constraints can be written nonexplicitly, that is, not associated with a Lagrange multiplier: $|u\rangle = \mathbf{P}|u'\rangle$, where $|u'\rangle$ is the surface salinity vector conserving the salt (dimension $n_x n_y - 1$) and $\mathbf{P}$ is an operator going from the subspace of surface salinity conserving the salt to the one of the full state vector. For instance, $\mathbf{P}$ is a matrix (dimension $2n_x n_y n_z$, $n_x n_y - 1$) composed of zero everywhere except for one (+1) on each line of the $(n_x n_y - 1)$ surface salinity points and a full line of -1 on the remaining surface salinity point.

Applying the same nonexplicit constraints as before—that is, the perturbation modifies only the surface salinity and conserves the total salt content—we write $|g\rangle = \mathbf{P}|g'\rangle$. By maximizing the Lagrangian function (B4) in permanent regime (infinite time), we obtain the eigenvalue problem

$$\mathbf{N}^{-1} \mathbf{H}(\infty) |g'\rangle = \gamma |g'\rangle, \quad \text{where} \quad (B6)$$

$$ \mathbf{H}(\infty) = \lim_{\tau \to \infty} \int_0^\tau ds \mathbf{P}^\dagger \mathbf{M}^\dagger(\tau - s) |F\rangle \times \langle F | \mathbf{M}(\tau - s) \mathbf{P} \rangle \quad \text{and} \quad \mathbf{N} = \mathbf{P}^\dagger \mathbf{S} \mathbf{P}. $$

Following the decomposition on the eigenvectors of $\mathbf{A}$ ($|u_k\rangle$ and $\mathbf{A}^\dagger(|u_k^\dagger\rangle)$, we can express $\mathbf{H}(\infty)$ as

$$\mathbf{H}(\infty) = \mathbf{P}^\dagger \sum_{ij} \frac{-1}{\lambda_i^* + \lambda_j} |u_i^\dagger\rangle \langle u_i|F\rangle \langle F|u_j\rangle \langle u_j^\dagger|\mathbf{P}, \quad (B7)$$
where $\lambda_k^*$ denotes the complex conjugates of $\lambda_k$. We know that the optimal pattern corresponds to the eigenmode of (B6) with the largest real part.

By using the same notation, we can obtain the power spectral density of the intensity of the MOC in terms of eigenvectors and eigenvalues of $\mathbf{A}$ and $\mathbf{A}^\dagger$ following Ioannou (1995), who has shown its enhancement in nonnormal dynamics context:

$$ 
\begin{aligned}
\text{psd}(\omega) &= \langle \hat{u}(\omega) | F \rangle \langle F | \hat{u}(\omega) \rangle, \\
&= \langle g | (-i\omega\mathbf{I} - \mathbf{A}^\dagger)^{-1} | F \rangle \langle F | (i\omega\mathbf{I} - \mathbf{A})^{-1} | g \rangle, \\
&= \sum_{ij} \langle g | u_i^\dagger \rangle \frac{1}{-i\omega - \lambda_i^*} \langle u_i | F \rangle \langle F | u_j \rangle \frac{1}{i\omega - \lambda_j} \langle u_j^\dagger | g \rangle.
\end{aligned} 
\quad (B8) $$

REFERENCES

Arzel, O., T. Huck, and A. Colin de Verdière, 2006: The different nature of the interdecadal variability of the thermohaline circulation under mixed and flux boundary conditions. *J. Phys. Oceanogr.*, **36**, 1703–1718.

Belkin, I. S., S. Levitus, J. Antonov, and S.-A. Malmberg, 1998: “Great salinity anomalies” in the North Atlantic. *Prog. Oceanogr.*, **41**, 1–68.

Bugnion, V., C. Hill, and P. H. Stone, 2006a: An adjoint analysis of the meridional overturning circulation in a hybrid coupled model. *J. Climate*, **19**, 3751–3767.

—, —, and —, 2006b: An adjoint analysis of the meridional overturning circulation in an ocean model. *J. Climate*, **19**, 3732–3750.

Chen, F., and M. Ghil, 1995: Interdecadal variability of the thermohaline circulation and high-latitude surface flux. *J. Phys. Oceanogr.*, **25**, 2547–2568.

—, —, and —, 1996: Interdecadal variability in a hybrid coupled ocean-atmosphere model. *J. Phys. Oceanogr.*, **26**, 1561–1578.

Colin de Verdière, A., 1988: Buoyancy driven planetary flows. *J. Mar. Res.*, **46**, 215–265.

—, and T. Huck, 1999: Baroclinic instability: An oceanic wave-maker for interdecadal variability. *J. Phys. Oceanogr.*, **29**, 893–910.

Curry, R., and C. Mauritzen, 2005: Dilution of the northern North Atlantic Ocean in recent decades. *Science*, **308**, 1772–1774.

—, B. Dickson, and I. Yashayaev, 2003: A change in freshwater balance of the Atlantic Ocean over the past four decades. *Nature*, **426**, 826–829.

Delworth, T. L., and M. E. Mann, 2000: Observed and simulated multidecadal variability in the Northern Hemisphere. *Climate Dyn.*, **16**, 661–676.

Dijkstra, H. A., and Coauthors, 2006: On the physics of the Atlantic Multidecadal Oscillation. *Ocean Dyn.*, **56**, 36–50.

Farrell, B. F., 1988: Optimal excitation of neutral Rossby waves. *J. Atmos. Sci.*, **45**, 163–172.

—, and A. M. Moore, 1992: An adjoint method for obtaining the most rapidly growing perturbation to oceanic flows. *J. Phys. Oceanogr.*, **22**, 338–349.

—, and P. J. Ioannou, 1993: Stochastic forcing of perturbation variance in unbounded shear and deformation flows. *J. Atmos. Sci.*, **50**, 200–211.

—, and —, 1996a: Generalized stability theory. Part I: Autonomous operators. *J. Atmos. Sci.*, **53**, 2025–2040.

—, and —, 1996b: Generalized stability theory. Part II: Non-autonomous operators. *J. Atmos. Sci.*, **53**, 2041–2053.

Frankignoul, C., and K. Hasselmann, 1977: Stochastic climate models. Part II: Application to sea-surface temperature anomalies and thermocline variability. *Tellus*, **29**, 289–305.

Greatbatch, R. J., and S. Zhang, 1995: An interdecadal oscillation in an idealized ocean basin forced by constant heat flux. *J. Climate*, **8**, 81–91.

Held, I. M., and B. J. Soden, 2006: Robust responses of the hydrological cycle to global warming. *J. Climate*, **19**, 5686–5699.

Huck, T., and G. K. Vallis, 2001: Linear stability analysis of three-dimensional thermally-driven ocean circulation: Application to interdecadal oscillations. *Tellus*, **53A**, 526–545.

—, A. Colin de Verdière, and A. J. Weaver, 1999a: Interdecadal variability of the thermohaline circulation in box-ocean models forced by fixed surface fluxes. *J. Phys. Oceanogr.*, **29**, 865–892.

—, A. J. Weaver, and A. Colin de Verdière, 1999b: On the influence of the parameterization of lateral boundary layers on the thermohaline circulation in coarse-resolution ocean models. *J. Mar. Res.*, **57**, 387–426.

Ioannou, P. J., 1995: Nonnormality increases variance. *J. Atmos. Sci.*, **52**, 1155–1158.

Josey, S. A., and R. Marsh, 2005: Surface freshwater flux variability and recent freshening of the North Atlantic in the eastern subpolar gyre. *J. Geophys. Res.*, **110**, C05008, doi:10.1029/2004JC002521.

Knight, J. R., R. J. Allan, C. K. Folland, M. Vellinga, and M. E. Mann, 2005: A signature of persistent natural thermohaline circulation cycles in observed climate. *Geophys. Res. Lett.*, **32**, L20708, doi:10.1029/2005GL024233.

Lohmann, G., and J. Schneider, 1999: Dynamics and predictability of Stommel’s box model: A phase space perspective with implications for decadal climate variability. *Tellus*, **51A**, 326–336.

Mann, M. E., R. S. Bradley, and M. K. Hughes, 1999: Northern Hemisphere temperature during the past millenium: Inferences, uncertainties, and limitations. *Geophys. Res. Lett.*, **26**, L759–L762.

Marotzke, J., 1996: Analysis of thermohaline feedbacks. *Decadal Climate Variability: Dynamics and Predictability*, D. L. T. Anderson and J. Willebrand, Eds., NATO ASI Series, Vol. 44, Springer, 333–378.

Mignot, J., and C. Frankignoul, 2003: On the interannual variability of surface salinity in the Atlantic. *Climate Dyn.*, **20**, 555–565.

Mitchell, J. M., 1976: An overview of climatic variability and its causal mechanisms. *Quat. Res.*, **6**, 481–493.

Moore, A. M., and B. F. Farrell, 1993: Rapid perturbation growth on spatially and temporally varying oceanic flows determined using an adjoint method: Application to the Gulf Stream. *J. Phys. Oceanogr.*, **23**, 1682–1702.

—, C. L. Perez, and J. Zavala-Garay, 2002: A non-normal view of the wind-driven ocean circulation. *J. Phys. Oceanogr.*, **32**, 2681–2705.

—, J. Vialard, A. T. Weaver, D. L. T. Anderson, R. Kleeman, and J. R. Johnson, 2003: The role of air-sea interaction in controlling the optimal perturbations of low-frequency tropical coupled ocean-atmosphere modes. *J. Climate*, **16**, 951–968.

Mu, M., and Z. Zhang, 2006: Conditional nonlinear optimal perturbations of a two-dimensional quasigeostrophic model. *J. Atmos. Sci.*, **63**, 1587–1604.
—, L. Sun, and H. A. Dijkstra, 2004: The sensitivity and stability of the ocean's thermohaline circulation to finite-amplitude perturbations. *J. Phys. Oceanogr.*, **34**, 2305–2315.

Sévellec, F., M. Ben Jelloul, and T. Huck, 2007: Optimal surface salinity perturbations influencing the thermohaline circulation. *J. Phys. Oceanogr.*, **37**, 2789–2808.

—, T. Huck, M. Ben Jelloul, N. Grima, J. Vialard, and A. Weaver, 2008: Optimal surface salinity perturbations of the meridional overturning and heat transport in a global ocean general circulation model. *J. Phys. Oceanogr.*, **38**, 2739–2754.

Sirkes, Z., and E. Tziperman, 2001: Identifying a damped oscillatory thermohaline mode in a general circulation model using an adjoint model. *J. Phys. Oceanogr.*, **31**, 2297–2305.

Stommel, H., 1961: Thermohaline convection with stable regimes flow. *Tellus*, **13**, 224–230.

—, and A. B. Arons, 1960: On the abyssal circulation of the world ocean. I. Stationary planetary flow patterns on a sphere. *Deep-Sea Res.*, **6**, 140–154.

Sun, L., M. Mu, D.-J. Sun, and X.-Y. Yin, 2005: Passive mechanism of decadal variation of thermohaline circulation. *J. Geophys. Res.*, **110**, C07025, doi:10.1029/2005JC002897.

te Raa, L. A., and H. A. Dijkstra, 2002: Instability of the thermohaline ocean circulation on interdecadal timescales. *J. Phys. Oceanogr.*, **32**, 138–160.

—, and —, 2003: Sensitivity of North Atlantic multidecadal variability to freshwater flux forcing. *J. Climate*, **16**, 2586–2601.

Tziperman, E., and P. J. Ioannou, 2002: Transient growth and optimal excitation of thermohaline variability. *J. Phys. Oceanogr.*, **32**, 3427–3435.

—, L. Zanna, and C. Penland, 2008: Nonnormal thermohaline circulation dynamics in a coupled ocean–atmosphere GCM. *J. Phys. Oceanogr.*, **38**, 588–604.

Vellinga, M., and R. A. Wood, 2002: Global climatic impacts of a collapse of the Atlantic thermohaline circulation. *Climatic Change*, **54**, 251–267.

Weaver, A. J., and E. S. Sarachik, 1991: The role of mixed boundary conditions in numerical models of ocean's climate. *J. Phys. Oceanogr.*, **21**, 1470–1493.

—, —, and J. Marotzke, 1991: Freshwater flux forcing of decadal and interdecadal oceanic variability. *Nature*, **353**, 836–838.

—, J. Marotzke, P. F. Cummins, and E. S. Sarachik, 1993: Stability and variability of the thermohaline circulation. *J. Phys. Oceanogr.*, **23**, 39–60.

Zanna, L., and E. Tziperman, 2005: Nonnormal amplification of the thermohaline circulation. *J. Phys. Oceanogr.*, **35**, 1593–1605.

—, and —, 2008: Optimal surface excitation of the thermohaline circulation. *J. Phys. Oceanogr.*, **38**, 1820–1830.