## Abstract

The optimal excitation of Atlantic meridional overturning circulation (MOC) anomalies is investigated in an ocean general circulation model with an idealized configuration. The optimal three-dimensional spatial structure of temperature and salinity perturbations, defined as the leading singular vector and generating the maximum amplification of MOC anomalies, is evaluated by solving a generalized eigenvalue problem using tangent linear and adjoint models. Despite the stable linearized dynamics, a large amplification of MOC anomalies, mostly due to the interference of nonnormal modes, is initiated by the optimal perturbations.

The largest amplification of MOC anomalies, found to be excited by high-latitude deep density perturbations in the northern part of the basin, is achieved after about 7.5 years. The anomalies grow as a result of a conversion of mean available potential energy into potential and kinetic energy of the perturbations, reminiscent of baroclinic instability. The time scale of growth of MOC anomalies can be understood by examining the time evolution of deep zonal density gradients, which are related to the MOC via the thermal wind relation. The velocity of propagation of the density anomalies, found to depend on the horizontal component of the mean flow velocity and the mean density gradient, determines the growth time scale of the MOC anomalies and therefore provides an upper bound on the MOC predictability time.

The results suggest that the nonnormal linearized ocean dynamics can give rise to enhanced MOC variability if, for instance, overflows, eddies, and/or deep convection can excite high-latitude density anomalies in the ocean interior with a structure resembling that of the optimal perturbations found in this study. The findings also indicate that errors in ocean initial conditions or in model parameterizations or processes, particularly at depth, may significantly reduce the Atlantic MOC predictability time to less than a decade.

## 1. Introduction

The present-day North Atlantic Ocean exhibits variability on a wide range of time scales, from days to several decades (e.g., Marshall et al. 2001b; Hurrell et al. 2006; Vallis 2009). Interannual to multidecadal fluctuations of North Atlantic sea surface temperature (SST) and corresponding atmospheric surface pressures are often associated with the variability of the ocean large-scale meridional overturning circulation (MOC). The Atlantic MOC, defined as the zonally averaged meridional flow and forced by both wind and buoyancy fluxes, transports about 15 Sv (1 Sv = 10^{6} m^{3} s^{−1}) of water and 1 PW (1 PW = 10^{15} W) of heat from low to high latitudes (estimated at 40°N; Ganachaud and Wunsch 2000). Cunningham et al. (2007) estimated a year-long-average MOC of 18.7 ± 5.6 Sv from the transatlantic Rapid Climate Change (RAPID) array at 26°N. Consequently, the MOC is expected to affect high-latitude climate by, for example, maintaining the Arctic sea ice edge (Poulsen et al. 2001; Winton 2003). Furthermore, the large localized sinking at high latitudes is associated with high CO_{2} uptake in these regions. It is therefore essential to understand the response of the circulation to perturbations to explain present-day Atlantic climate variability and to potentially anticipate future changes.

The variability and stability of the Atlantic MOC have been extensively studied (e.g., Delworth et al. 1993; Weaver et al. 1993; Griffies and Tziperman 1995; Delworth and Greatbatch 2000; Bryan et al. 2006). It is widely believed that the MOC responds to various forcing due to external factors, such as the North Atlantic Oscillation, or internal factors, such as deep convection, or gyre strength, on intraseasonal to interannual and multidecadal time scales (e.g., Bjerknes 1964; Kushnir 1994; Marshall et al. 2001a). The sensitivity of the MOC and meridional ocean heat transport to initial conditions and surface forcing has also been investigated in 3D ocean models using the adjoint method (Marotzke et al. 1999; Bugnion et al. 2006; Sévellec et al. 2008; Czeschel et al. 2010; Heimbach et al. 2011).

General circulation models (GCMs) used in climate studies show multidecadal fluctuations in their Atlantic MOC with pronounced periods ranging from 20 to more than 100 yr (e.g., Weaver et al. 1991; Delworth et al. 1993; Dong and Sutton 2001); SSTs and northward heat transport exhibit similar variability. These fluctuations could be either damped oscillatory ocean modes stochastically excited by atmospheric variability (e.g., Griffies and Tziperman 1995), self-sustained oscillatory modes of the ocean itself (e.g., Winton and Sarachik 1993; Weaver et al. 1993; Greatbatch and Zhang 1995; Chen and Ghil 1995), or of the coupled ocean–atmosphere system (e.g., Timmermann et al. 1998). The driving mechanisms of these MOC fluctuations and what sets their “period” remain largely unresolved despite recent attempts (Frankcombe et al. 2010), and this multidecadal intrinsic variability further complicates climate studies aimed at investigating anthropogenic effects.

While atmospheric dynamics are largely nonlinear and therefore typically unpredictable on time scales longer than a few days, the variability arising from the ocean is potentially predictable on much longer time scales. Evidence from models suggests that interdecadal and multidecadal fluctuations of the Atlantic circulation and SST exhibit some predictability on decadal time scales (e.g., Griffies and Bryan 1997; Pohlmann et al. 2004; Sutton and Hodson 2005; Collins et al. 2006). Generally, these experiments provide relatively large upper bounds on the predictability time of the North Atlantic Ocean and climate, since they assume a perfect model and no errors in the ocean initial conditions.

In this paper, we study the optimal perturbations, defined as singular vectors, leading to the linear amplification of MOC anomalies in an ocean GCM. These perturbations typically arise from the nonnormality of the linearized dynamical system operator (e.g., Farrell 1982, 1988, 1989; Trefethen et al. 1993). Most fluid dynamical systems are nonnormal, meaning that their eigenvectors do not form an orthogonal basis and that transient growth can occur even though the linearized dynamics are stable.

The singular vectors have been successfully used, for example, in numerical weather forecast and ENSO prediction (e.g., Buizza and Palmer 1995; Penland and Sardeshmukh 1995; Moore and Kleeman 1997a,b). Additionally, the nonnormal dynamics of the large-scale ocean thermohaline circulation and MOC has begun to receive some attention in simple box models (Lohmann and Schneider 1999; Tziperman and Ioannou 2002) and 2D models of intermediate complexity (Zanna and Tziperman 2005, 2008; Sévellec et al. 2007; Alexander and Monahan 2009). Tziperman et al. (2008) and Hawkins and Sutton (2009) used linear inverse modeling to approximate the transient growth of MOC anomalies in coupled climate models and subsequently explore the predictability of the MOC.

In this study, we explore the mechanisms for transient amplification of MOC anomalies using the 3D Massachusetts Institute of Technology (MIT) General Circulation Model (MITgcm; Marshall et al. 1997a,b) and its tangent linear and adjoint models (Giering and Kaminski 1998; Marotzke et al. 1999; Heimbach et al. 2005), and we attempt to address the following questions:

Where is the MOC most sensitive to ocean temperature and salinity perturbations?

What is the physical mechanism leading to the growth of the MOC anomalies, and how can it be related to Atlantic Ocean circulation and climate variability?

How fast can errors in ocean initial conditions grow and potentially limit the predictability of the MOC?

By solving an optimization problem and calculating the singular vectors, we find optimal perturbations of deep density anomalies in the northern part of the basin that can amplify MOC anomalies to a maximum that is reached after 7.5 yr. The main growth of MOC anomalies is explained by examining the time evolution of deep zonal density gradients related to the MOC via the thermal wind relation. The growth time scale of the MOC anomalies, and therefore an upper bound on the MOC predictability time, are determined by the cyclonic propagation speed of the density anomalies in the northern part of the basin. The velocity of propagation of the anomalies depends on the mean flow velocity and the mean density gradient.

The paper is organized as follows. The MITgcm and the base steady state reached by the model are described in section 2. In section 3, we present the methodology employed to calculate the singular vectors around the model steady state using the tangent linear and adjoint models, and then we describe the spatial pattern of initial temperature and salinity perturbations found to maximize the growth of MOC anomalies. The amplification mechanism of the anomalies is explained in section 4. We conclude and discuss the results in section 5.

## 2. The base steady state of the MITgcm

### a. Model configuration

To investigate the nonnormal excitation of MOC anomalies, the MITgcm (Marshall et al. 1997a,b), which solves the primitive equations on a sphere under the Boussinesq and hydrostatic approximations for an incompressible fluid, is used. The model configuration is similar to the one described in Zanna et al. (2010), with the exception of the higher horizontal and vertical resolutions used in the present study. A brief model description is given below, and we refer the reader to Zanna et al. (2010) for additional details.

The prognostic model variables are the horizontal velocity **v** = (*u*, *υ*), potential temperature *T*, salinity *S*, and sea surface height *η*. Vertical velocity *w*, density *ρ*, and hydrostatic pressure *ϕ*_{hyd} are all diagnostic quantities. The equation of state is a modified United Nations Educational, Scientific and Cultural Organization (UNESCO) formula derived by Jackett and McDougall (1995), and ocean convection is parameterized by an implicit vertical diffusion scheme. The idealized Atlantic-like double-hemispheric rectangular basin is defined between latitudes 67.5°S and 67.5°N and longitudes 65° and 5°W with a horizontal resolution of 3°. The ocean depth is uniformly set to 5000 m everywhere across the basin, with 15 vertical levels with thicknesses ranging between 50 and 690 m. Time-independent wind forcing and mixed boundary conditions of temperature and salinity are used. The temperature is restored with a time scale of 2 months, and a freshwater flux is prescribed for the salinity surface boundary conditions (Zanna et al. 2010). All relevant model parameters are listed in Table 1.

### b. The steady state

The steady state reached by the model under the above mixed boundary conditions is characterized by SST between 28.5°C at the equator and −0.5°C at high latitudes, and sea surface salinity (SSS) ranging from 35.6 ppt at the equator to 34.7 ppt at high latitudes. The near-surface flow is characterized by a single gyre in each hemisphere, with a similar but reversed and weaker gyre at depth, as shown in Fig. 1. Downwelling occurs at high latitudes near the eastern boundary (Fig. 1). Deep convection is triggered poleward of 50° near the boundaries.

Figure 2 shows the steady-state meridional overturning streamfunction with a strong asymmetric cell and a maximum transport of about 22 Sv in the Northern Hemisphere. The zonally and vertically integrated heat transport has a maximum of 1.2 PW near 30°N (not shown). While the model reached a fairly reasonable steady state, the idealized geometry and forcing lead to perhaps some unrealistic aspects of the ocean state. We should keep in mind that a different steady state may quantitatively affect the results. However, we find here that different steady states of meridional overturning circulation configurations lead to fairly similar optimal initial conditions and growth mechanism, as discussed in section 4d.

## 3. Optimal initial conditions maximizing the growth of MOC anomalies

### a. Defining and evaluating the optimal initial conditions

To evaluate the optimal perturbations, we will be using the tangent linear model derived from the MITgcm equations, which may be represented by

and its adjoint (Giering and Kaminski 1998; Marotzke et al. 1999; Heimbach et al. 2005). The time-dependent vector **P**′(*t*) is a small perturbation from the steady-state solution for temperature, salinity, and horizontal velocity, while the time-independent (autonomous) matrix 𝗔 is the linearized model dynamical operator obtained via the automatic differentiation tool Transformation of Algorithms in Fortran (TAF; Giering and Kaminski 1998; Heimbach et al. 2005). The solution to the linearized model (1) is given by

where 𝗕(*t*) ≡ *e*^{𝗔t} is the propagator matrix, and **P**′_{0} ≡ **P**′(*t* = 0) is the initial perturbation. The system is found to be linearly stable, such that all eigenvalues of 𝗔 have negative real parts (not shown) and every perturbation eventually decays as *t* goes to infinity.

To investigate the sensitivity of MOC anomalies to ocean initial conditions, and therefore the error growth limiting the MOC predictability, we proceed by maximizing the sum of squares of the MOC anomalies at time *τ*, evaluated using a norm kernel 𝗫 as follows:

The cross-sectional area is *dA* = *r dϕ dz*, *ϕ* is the latitude, *z* is the depth, *H* is the uniform ocean depth, and *ψ*′ is the meridional overturning streamfunction anomaly given by

where *λ* is the longitude. The quantity to be maximized, , was chosen such that the integration area *A* in (3) covers the entire Northern Hemisphere. We opt for an averaged value over a given region rather than the value of *ψ*′ at a single latitude or depth to avoid artifacts, as discussed by Zanna and Tziperman (2005) in the context of a simple box model. In addition, high-latitude areas in the model are associated with nonlinear processes, such as convection, and it is therefore preferable to use a quantity reflecting an average over an area rather than a single grid cell to avoid mechanisms relying solely on nonlinear processes. Moreover, the sum of squares of the MOC anomalies permits us to investigate simultaneously the effect of perturbations on the MOC variability and the error growth limiting the MOC predictability, since it can viewed as the error variance of the model forecast for MOC fluctuations anomalies (Lorenz 1982). The norm kernel 𝗫, defined by (3), does not span the entire space and may, in principle, lead to a growth that is not due to the nonnormality of the ocean dynamics. This issue will be addressed in section 4d, and the nonnormality will be shown to be the main cause for the growth of anomalies.

The initial perturbations are constrained to have a unit norm reflecting the energy of the perturbations. Given an initial density anomaly field, the velocity and sea surface height perturbations are expected to adjust to the density field within a few days, leading to a horizontal flow in geostrophic balance with the density gradients. We found in several of our numerical experiments that the initial kinetic energy of the perturbations does not contribute to the nonnormal growth for time scales longer than a few weeks. Since our analysis concentrates on the response of the ocean to perturbations on interannual to decadal time scales, it suffices to only consider initial perturbations of temperature and salinity about the steady-state solution, which maximize the MOC anomalies after a given time *τ* (Zanna et al. 2010). The velocity and sea height anomalies are therefore assumed to initially vanish. This reduces the number of unknown to be calculated and makes the problem more manageable computationally.

The initial temperature and salinity perturbations are constrained to have a unit norm such that ‖**P**′_{0}‖_{𝗘}^{2} = **P**′_{0}^{T}𝗘**P**′_{0} = 1. The energy norm kernel 𝗘 therefore ensures that the contributions of both salinity and temperature to the density, as well as the volume of the grid boxes, are properly accounted for, such that

where *α* = *α*(*T*, *S*) and *β* = *β*(*T*, *S*) are the space-dependent thermal expansion and saline contraction coefficients, respectively; and *T*(*λ*, *ϕ*, *z*) and *S*(*λ*, *ϕ*, *z*) are the steady-state temperature and salinity, respectively.

Maximizing the growth of the MOC anomalies at *t* = *τ* measured by (3) with respect to the initial perturbations of temperature and salinity anomalies using (5) as the initial constraint is equivalent to solving the generalized eigenvalue problem (e.g., Farrell 1988; Zanna and Tziperman 2005) given by,

The optimal initial conditions **P**′_{0} are defined as the fastest-growing perturbations corresponding to the generalized eigenvector of (6) with the largest eigenvalue *γ* (Farrell and Ioannou 1996). In this case, the optimal perturbations correspond to the rescaled first singular vector of 𝗫^{1/2}𝗕(*τ*)𝗘^{−1/2} at time *τ*. Note that an alternative method to define (distinct from the singular vectors) and calculate optimal perturbations of an ocean GCM was used by Sévellec et al. (2008). By avoiding the use of a norm, Sévellec et al. (2008) obtained an explicit solution for the surface salinity perturbations, which depends on the adjoint modes only, leading to the maximum change of MOC anomalies at a single latitude.

In the present study, the optimal perturbations **P**′_{0} are computed using the Lanczos algorithm (Golub and Van Loan 1989) and the routines for symmetric eigenvalue problems of the Arnoldi Package (ARPACK) software (Lehoucq et al. 1998), which requires only the input vector 𝗕^{T}𝗫𝗕**P**′, where the superscript T denotes the matrix transpose (equivalent to the adjoint with respect to a L2-norm). A complete description of the numerical calculation of the optimals can be found in Moore et al. (2004) and additionally in Zanna et al. (2010) for the MITgcm. Zanna et al. (2010) successfully implemented and used this machinery to explore the role of nonnormal ocean dynamics in exciting tropical Atlantic SST variability. The algorithm was found to be robust to the model resolution and to other assumptions.

### b. Spatial structure of the singular vectors

The goal of this study is to find where are the perturbations to which the MOC is most sensitive in order to explore how MOC variability can be optimally excited, and where the predictability of the MOC may be limited by error growth due to uncertain initial conditions and dynamical processes. The optimal initial perturbations of temperature and salinity anomalies **P**′_{0} solutions of the generalized eigenvalue problem (6) and leading to the maximum MOC anomalies at *t* = *τ*, are allowed anywhere in the ocean basin (any latitude, longitude, and depth; Zanna et al. 2010). The maximally amplified MOC anomalies, measured by and defined in (3) under the constraint defined by (5), occur for *τ* ≈ 7.5 yr.

The optimal initial perturbations **P**′_{0}, corresponding to *t* = 7.5 yr, that excite the optimal MOC growth and correspond to the largest *γ* are concentrated at high latitudes in the northern part of the basin below a depth of 1 km with a baroclinic structure. A vertically uniform and weak signal is seen in the upper 1 km (Fig. 3). At depth, the ratio between |*αT* ′| and |*βS*′| is roughly 0.44. While the optimal perturbations show some degree of compensation of the temperature and salinity in the area where is defined, their contributions to the density field do not completely vanish, unlike the findings using 2D models by Tziperman and Ioannou (2002) and Zanna and Tziperman (2005). An estimate of the optimal MOC amplification can be obtained by initializing the linearized model with the optimal perturbations found. Using an initial amplitude of 0.02 kg m^{−3} for the deep density anomalies (approximately equivalent to temperature and salinity perturbations of 0.1°C and 0.05 ppt, respectively, thus keeping their relative amplitude found within the calculated singular vector) results in a 2.4-Sv MOC anomaly after 7.5 yr.

As noted above, the optimal perturbations correspond to the leading singular vector associated with the largest singular value and therefore leading to the largest possible MOC anomaly at 7.5 yrs. The spatial pattern corresponding to the second singular value at *t* = 7.5 yrs is found to have a Southern Hemisphere signal in addition to a Northern Hemisphere high-latitude pattern at both the surface and depth (Fig. 4). The second singular vector, similar to the optimal perturbations, is depth intensified. This spatial pattern excites a much smaller amplification of MOC anomalies, as expected. We find that an initial 0.02 kg m^{−3} deep density anomaly associated with the second singular vector leads to only a 1.1-Sv MOC anomaly. By initializing the model with the parts of the second singular vector corresponding to each hemisphere, we find that about 60% of the growth of MOC anomalies is due to the Northern Hemisphere signal, while 40% is due to the Southern Hemisphere anomalies, such that a Southern Hemisphere anomaly of roughly 0.02 kg m^{−3} leads to a 0.45-Sv MOC anomaly after 7.5 yr. Therefore, despite the significant and major effect of the Northern Hemisphere signal in the first and second singular vectors, a Southern Hemisphere anomaly can still influence the MOC, as found, for example, in previous studies on the effect of Agulhas leakage on the MOC strength (e.g., Weijer et al. 1999; Biastoch et al. 2008). Further investigation is necessary, as the effect on the MOC by the Agulhas Current is primarily found above the thermocline and not at depth (e.g., Biastoch et al. 2008; Heimbach et al. 2011).

## 4. MOC amplification mechanism

### a. Overview of the evolution of MOC anomalies excited by the optimal perturbations

Consider the first singular vector or optimal perturbations (Fig. 3) leading to the maximum MOC amplification at *t* = 7.5 yr. Deep temperature and salinity anomalies appear to evolve similarly in time. Therefore, instead of explaining the mechanism separately for temperature and salinity anomalies, the density is treated below as the independent variable. When the linearized model is initialized with the optimal density perturbations, the time series of the function (3) is characterized by a growing oscillatory-like behavior (Fig. 5a). Some of the eigenvalues and corresponding eigenvectors of the tangent linear model and its adjoint were computed using ARPACK. The dynamics of the ocean model are found to be linearly stable, and the behavior in Fig. 5 can be mainly explained by linear interference of nonorthogonal decaying eigenmodes of the propagator 𝗕. Projections of the optimal perturbations onto the eigenmodes of the propagator 𝗕 reveal that a large number of nonorthogonal eigenmodes participate in the growth of MOC anomalies (not shown). Such an analysis, therefore, does not necessarily shed much light on the physical mechanism responsible for the growth of perturbations, and a global view of the dynamics seems more appropriate. Overall, by analyzing the modal composition of the optimal perturbations at *t* = 0 and *t* = 7.5 yr, we find that the relative phase of the eigenmodes changes dramatically over time and that the spatial structure of the initial excitation and of the response at maximum amplification time are very different, as is typical of nonnormal systems.

Before getting into the details of the growth mechanism and time scale, it is worth mentioning that the optimal initial conditions are the solution to a linear problem and therefore a multiplication by any (positive or negative) constant will also be a solution to the eigenproblem. We arbitrarily pick one sign to simplify the following discussion and to distinguish between the different stages of the growth. The characteristics of the MOC excitation can be explained as follows. The optimal initial perturbations of deep density anomalies lead to a rapid two-week geostrophic adjustment of the velocity and sea surface height anomalies. This is followed by an increase in MOC anomalies between 41° and 60°N and a decrease at other latitudes north of the equator until the first peak in is reached at about *t* = 2 months (Fig. 5). From *t* = 2 months, the MOC anomalies decrease over the entire North Atlantic until reaching a minimum that corresponds to the maximum amplification of the sum of squares of the MOC anomalies, , at about *t* = 7.5 yr. Beyond this time, the MOC anomalies continue to oscillate and decay in amplitude. We now proceed to analyze each of these stages in the evolution of the MOC anomalies.

### b. MOC amplification time scale: Cyclonic propagation of deep density anomalies

The initial density perturbations were projected onto the barotropic and baroclinic modes to understand the initial geostrophic adjustment of the velocities and sea surface height to the optimals. The initial density perturbations in the interior and near the northern boundary project mostly onto the first baroclinic mode, and the geostrophic adjustment of all fields occurs within about two weeks. In the upper ocean, the resulting geostrophic flow leads to a net downwelling near the northern and western boundaries.

Following the period of geostrophic adjustment, despite the contribution of the upper-ocean perturbations to the amplitude of the MOC growth, the evolution of MOC anomalies beyond the first few weeks is dominated by the time evolution of the deep density perturbations. More specifically, from the thermal wind relation

a negative zonal density gradient (∂*ρ*′/∂*x* < 0) in the Northern Hemisphere ( *f* > 0) results in a positive meridional velocity shear (∂*υ*′/∂*z* > 0) and therefore a positive anomaly of the meridional overturning streamfunction *ψ*′. Figure 6 shows that, between 42° and 51°N, the time evolution of the MOC anomalies can be mainly accounted for by the evolution of the 3010-m-depth zonal density gradient, the latter measured between the eastern and western boundaries and averaged over the latitudes 42°–51°N. We observe similar relationships at other latitudes as expected from the thermal wind in Eq. (7).

Our task is therefore reduced to explaining the evolution of the density gradient between the eastern and western boundaries. The first peak in the time series for the MOC anomaly occurs after about 2 months (Fig. 5) and is characterized by positive MOC anomalies between 40° and 60°N and negative MOC anomalies elsewhere in the Northern Hemisphere (Fig. 7, top panels). Near the western boundary, the initial density perturbations between 45° and 60°N are positive (“1” in Fig. 3a) and negative north of 60°N (“2” in Fig. 3a). Near the eastern boundary, the density anomalies are positive north of 23°N (“2” in Fig. 3a) and vanish south of this region. The active convection near the eastern boundary is acting to reduce the density anomalies there. Over the initial 2 months, the negative density gradient between the eastern and western boundaries between 42° and 60°N is thus further decreasing, leading to a positive meridional velocity shear (*υ*′* _{z}* > 0) that increases as the zonal density gradient decreases and becomes negative. The increasing meridional velocity shear corresponds to a positive MOC anomaly (Fig. 7, top panels). Similarly at other latitudes, the positive zonal density gradient explains the negative MOC anomaly. This explains the first peak at about 2 months in the MOC anomaly time series when initialized with the optimals (Fig. 5).

Next, consider the mechanism by which the amplification reaches its main peak at about 7.5 years. Starting at 2 months, the overall overturning decreases until it reaches a minimum at roughly 7.5 years (maximum negative anomaly shown in middle panels of Fig. 7). The negative density anomaly in the deep ocean, initially located at 60°N, 45°W (“2” in Fig. 3a), travels counterclockwise near the northern and western boundaries, reaching the latitude 54°N and longitude 60°W after 7–8 years (“2” in left middle of Fig. 7). If the density perturbation was simply advected by the steady-state velocity, then the initial perturbation should have reached this location after 12 yr. To explain this difference in propagation time scale, we examine the dominant terms in the linearized equation for the density perturbations at 3-km depth. The dominant balance is given by

On the time scales of interest, mixing is required to explain the decay of the perturbations; however, it does not actively participate in the propagation of the anomaly and is therefore neglected for simplicity.

Using the thermal wind equations and assuming that the density anomaly spreads over a vertical scale *h*′, we can approximate the velocity anomalies as

The propagation velocity of the density anomalies thus depends on the velocity of the mean flow and the horizontal component of the mean density gradient (Colin De Verdière and Huck 1999; Te Raa and Dijkstra 2002). Using the (positive) meridional component of the mean density gradient and the (negative) zonal component of the mean density gradient, the direction of propagation is cyclonic, explaining the above optimal perturbation results. The magnitudes of the mean flow and the mean density gradient terms are such that the deep density anomalies should reach the western boundary after roughly 7–8 years. A similar westward propagation of temperature anomaly at the ocean surface was previously analyzed by Colin De Verdière and Huck (1999) and categorized as a special case of nondivergent potential vorticity waves, and later by Te Raa and Dijkstra (2002) and Frankcombe et al. (2008), who named them “SST modes” or “thermal” Rossby waves.

The westward propagation at depth can be explained in terms of the negative meridional density gradient anomaly that develops due to the initial negative density at 45°W near the northern boundary (“2” in Fig. 3a). On the other hand, the meridional component of the mean density gradient is positive. Therefore, after some initial adjustment, the anticyclonic geostrophic anomalous flow around the negative density anomaly causes southward advection of dense water to the east of the anomaly and northward advection of light water to the west of the anomaly. This results in the westward propagation of the light anomaly, as described by Eq. (10). Similarly, the anomalies propagate southward (“1” and “2” in Fig. 7), with a typical velocity given by the southward mean advection velocity and the zonal component of the mean density gradient, *υ* + (*gh*′/*ρ*_{0}*f* )(∂*ρ*/∂*x*). At the maximum amplification time (*t* = 7.5 yr), the positive anomaly previously near the western boundary (“1” in Fig. 7) is now traveling equatorward with a weaker amplitude and is replaced by the westward-propagating negative density perturbation (“2” in Fig. 7). Hence, a positive east–west density gradient at 3010-m depth is created and induces, via the thermal wind relation, a negative meridional velocity shear and thus a negative MOC anomaly at the maximum amplification time (*t* = 7.5 yr) over the entire Northern Hemisphere.

A growing phase of a stable oscillation could be easily confused with nonnormal growth. In a stable linear normal system, the energy (defined as a positive definite quantity) of the initial perturbations decays; however, this is not necessarily the case for a stable nonnormal system in which energy is extracted from the mean flow leading to a transient growth of the perturbations. To confirm that nonnormality is the primary mechanism responsible for the MOC growth, and to rule out a possible growth due to a phase change of the dominant mode, we have diagnosed the kinetic energy anomaly and the density-squared anomaly over the North Atlantic domain. The kinetic energy anomaly reaches a maximum at 7.3 yr with a growth factor of roughly 3.7 (not shown), and the density squared is amplified by a factor of 14 when it reaches its maximum after 6.9 yr (Fig. 5b). We conclude that despite the participation of an oscillatory mode, the growth of all variables shows that the excitation of MOC anomalies is primarily due to the nonnormal dynamics.

Recall that the interaction of multiple linear eigenmodes is necessary for the amplification to occur. After the relatively fast decay of a few eigenmodes, an interdecadal oscillatory mode emerges. In addition to the westward propagation described earlier, this oscillatory eigenmode exhibits similarities to the mode described in Colin De Verdière and Huck (1999) and Te Raa and Dijkstra (2002), such as a phase lag between the north and south temperature gradient and the overturning strength, although some differences exist as well. For example, in the models used by Colin De Verdière and Huck (1999) and Te Raa and Dijkstra (2002), salinity is entirely neglected and the only prognostic variable is temperature, such that the velocities are diagnosed from the temperature. The mode excited here by the optimal perturbations is seen to depend on the mean density gradient (not only the temperature gradient), such that the mean salinity plays a role in the time evolution of the anomalies by setting the background vertical shear, the latter being an active participant in the growth mechanism (section 4c). Therefore, salinity, which is absent from previous studies, plays an important role in both the oscillatory mode that emerges from the optimal excitation and in setting the amplitude and time scale of the growth.

Finally, consider the decay of the optimal perturbations beyond 10 yr. A weaker and positive MOC anomaly peak occurs later around 20 years (Fig. 7a). This peak is not accompanied by growth in either the density anomalies (Fig. 7b) or kinetic energy anomaly, and it is mainly due to the oscillatory behavior of the interdecadal mode described above and not to the nonnormal dynamics. The positive density anomalies initially centered at 15°W and 57°N (“3” in Fig. 3a, which is therefore farther east than anomaly “2” in Fig. 3a) propagate to the western boundary as described above. The relatively slow propagation of the anomalies (“3” in Fig. 7) leaves enough time for the signal to dissipate with time, explaining the weaker peak amplitude of the anomalies at 20 yr (Fig. 5). At this time, the zonal density gradient is negative and therefore the MOC anomaly is positive. The propagation of the negative anomaly (“2” in Fig. 3a) southward leads a positive density gradient and therefore a persisting negative MOC anomaly between the equator and 40°N at depth.

### c. Energy budget and meridional heat transport

The growth of MOC perturbations can also be understood by analyzing the different terms in the energy budget. Figure 5b shows the growth of the squared densities anomalies, related to the available potential energy of the perturbations. Assuming that the anomalous flow is in geostrophic balance, density perturbations alter the available potential energy, some of which is transformed into kinetic energy, as suggested by the strong baroclinic structure of the optimal perturbations leaning against the background velocity shear. From the available potential energy equation, we find that the growth of perturbations is mainly due to the source term *ρ*′**u**′ · **∇***ρ* acting against the damping of perturbations by dissipation and damping of surface temperature by the restoring boundary conditions. This source term is interpreted as the change of available potential energy due to the interaction of the buoyancy perturbation and the anomalous buoyancy advection (Huang 2002) and is transferred to the kinetic energy to increase **u**′_{2}/2 and therefore possibly the MOC as well. While the propagation of the anomalies is due to both the advection of the mean density gradient by the perturbation flow and by advection of the density perturbations by the mean flow, by eliminating different terms in the temperature and salinity equations, we find that the source of nonnormality for the MOC growth mostly comes from advection of the mean density gradient (strongly affected by the salinity) by the perturbed flow (not shown).

The above results suggest that the nonnormal growth mechanism can excite large MOC anomalies. Assuming that the dynamics are linear, a small initial deep density perturbation of 0.02 kg m^{−3} can result in an MOC anomaly of 2.4 Sv after about 7.5 years. Because of the growth of velocity and temperature anomalies, changes in the zonally and vertically integrated meridional heat transport carried by the ocean,

are expected as well, where *c _{p}* is the heat capacity of water. The meridional heat transport anomaly is influenced by the optimal perturbations and the MOC growth over the entire Northern Hemisphere, and its evolution at different latitudes will differ. The heat transport anomalies at 51° and 57°N as function of time, due to the optimal initial conditions leading to maximum MOC growth, are shown in Fig. 8. We find that a 0.02 kg m

^{−3}initial density anomaly results in a 0.08-PW anomaly in the heat transport at the time of maximum amplification of the MOC. The growth of the heat transport is about 7% of the steady-state value, which is again quite large for such a small initial perturbation.

### d. Sensitivity and robustness of the results

Several numerical experiments reveal that the spatial structure of the optimal perturbations is fairly robust to changes in the steady-state solution or the norm kernel used. Using a steady state with a meridional overturning streamfunction that is symmetric across the equator (obtained under restoring boundary conditions rather than mixed boundary conditions) leads to optimals with a strong signal below a depth of 1 km at high latitudes and similar to the ones described above and shown in Fig. 3. Similarly, a steady overturning circulation with two deep and asymmetric cells results again in very similar optimals. Since our optimal perturbations are dominated by a Northern Hemisphere signal at large depth, alternative forms of the norm kernel 𝗫 are explored to ensure that our results are not entirely dictated by its definition in Eq. (3). For a norm kernel reflecting the sum of squares of the MOC anomalies over a different area than the one used in Eq. (3)—for example, between latitudes 21° and 27°N and depths of 200 and 1500 m (where the heat transport is found to be maximal)—again leads to optimals that are similar to the ones shown in Fig. 3. Qualitatively similar results are obtained if the area defining our cost is between 51° and 60°N and 2- and 4-km depth. Overall, we find that the structure of the mean isopycnals sets most of the amplification and its time scale.

We should bear in mind that our results are obtained using a linear analysis where the response is independent of the sign of the initial perturbations, which is not necessarily the case in the fully nonlinear system. There are alternative approaches that allow for the study of nonlinear optimal perturbation growth (Mu et al. 2004); however, this approach is rather expensive to be used in full 3D GCM. The optimal initial perturbations and associated MOC growth mechanism are also relevant as long as the perturbations are small and the linearity assumption is not violated. To ensure that nonlinearities do not become predominant over the 7-yr growth time scale of the MOC anomaly, we initialize several fully nonlinear model experiments with the spatial pattern of the optimal perturbations found. The amplitude of the perturbations for the different numerical experiments is between |0.05|° and |0.4|°C for the temperature and 0.02 and 0.15 ppt for the salinity. These perturbation amplitudes compare with estimates of deep-ocean variability (Forget and Wunsch 2007). For the largest amplitudes of temperature and/or salinity, the nonlinearities lead to deviations on the order of 13% in the MOC response relative to the linear case after 7.5 yr. Overall, the amplitude of the MOC anomalies at time of maximum amplification is slightly reduced by the nonlinearities. While the first peak in the MOC time series after about two months involved convection near the eastern boundary, the main growth mechanism over the following 7 yr did not involve convection. If convection is turned off after the initial 5 months, the main MOC growth is still observed after roughly 7–8 years.

## 5. Conclusions

In this study, the amplification of Atlantic meridional overturning circulation (MOC) anomalies due to the nonnormal linear ocean dynamics is considered. While previous explanations of Atlantic meridional overturning circulation and associated climate variability on interannual to multidecadal time scales relied on the excitation of a damped oscillatory mode (e.g., Griffies and Tziperman 1995; Saravanan and Mcwilliams 1997; Eden and Greatbatch 2003) by stochastic forcing or an unstable oscillatory mode of the ocean (e.g., Quon and Ghil 1995; Te Raa and Dijkstra 2002), the present work suggests the excitation of several nonorthogonal damped modes as an alternative or additional mechanism. More specifically, even though the flow is linearly stable and perturbations eventually decay, we find optimal initial perturbations of temperature and salinity capable of generating significant growth of the MOC anomalies with a time scale of 7.5 yr. These optimal density perturbations leading to the largest possible growth of MOC anomalies, defined as the leading singular vector, are found to be concentrated at high latitudes in the northern part of the basin below a depth of 1 km with a baroclinic structure. The density perturbations are “leaning” against the background velocity shear.

The MOC anomalies’ growth coincides with a growth in both kinetic energy and density anomalies (related to the available potential energy of the perturbations). The growth, reminiscent of baroclinic instability converting potential energy into kinetic energy, involves a cyclonic propagation of deep temperature and salinity anomalies. The propagation speed depends on both the mean flow advection velocity (*u*, *υ*) and the horizontal component of the mean density gradient (*ρ*_{x}, *ρ*_{y}) (similar to Doppler-shifted Rossby waves). The propagation of the anomalies leads to large zonal density gradients between the western and eastern boundaries in the northern part of the basin, thus resulting in a meridional overturning anomaly (via the thermal wind balance). Such a mechanism cannot be found in zonally averaged 2D models (Zanna and Tziperman 2005; Sévellec et al. 2007; Zanna and Tziperman 2008) because of the absence of geostrophic and thermal wind balance. While the propagation of the anomalies is due to both the advection of the mean density gradient by the perturbed flow and by the advection of the density perturbations by the mean flow, the source of nonnormality and therefore amplification mechanism for the MOC comes mostly from advection of mean density gradient (strongly influenced by the salinity) by the perturbed flow. We found that an initial perturbation of 0.02 kg m^{−3} amplitude can result in MOC anomalies of 2.4 Sv (about 11% of the mean circulation) after 7.5 yr and a heat transport anomaly of 0.05 PW (about 7% of the mean value) because of the nonnormal amplification mechanism in both the linearized and full nonlinear models.

If deep temperature and salinity perturbations can be excited by eddies (e.g., Treguier et al. 2006), overflows (e.g., Eldevik et al. 2009), and/or deep convection (e.g., Pickart and Spall 2007), especially at high latitudes near the boundaries, our mechanism may account for some of the MOC variability on interannual time scales. It would therefore be interesting to force the temperature and salinity equations with stochastic forcing at depth to explore what percentage of the MOC variance can be explained by the eddy field, overflows, or deep convection. Moreover, our findings suggest that errors in initial conditions at depth (where observations are very sparse) would strongly limit our ability to predict the climate in the North Atlantic region. Therefore, further observations at high latitudes (e.g., Kohl 2005), but especially of the deep ocean, are crucial to improve our predictions in the coming years and decades.

The optimal perturbations defined in the current study (i.e., singular vectors) allow for both temperature and salinity anomalies at all depths. We find that deep perturbations (rather than surface ones) are far more efficient in creating large growth of MOC anomalies on a time scale of 7.5 yr.

Adjoint-only methods in ocean GCMs have been extremely useful to explore the sensitivity of the MOC to temperature (Heimbach et al. 2011) and salinity (Sévellec et al. 2008) perturbations. Heimbach et al. found, for example, sensitivity of the MOC to high-latitude subsurface temperatures on interannual time scales. Kelvin, Rossby, and continental shelf waves are believed to trigger the MOC response. An alternative method to define and calculate the optimal perturbations, distinct from the singular vectors calculated here, was used by Sévellec et al. (2008). They calculated initial surface salinity perturbations leading to a MOC change. Their mechanism relied on the surface zonal density gradient anomaly produced by the initial salinity anomaly and involved the eastward advection of the mean temperature by the anomalous zonal velocity. They found that a 2-Sv MOC anomaly can be created after about 10 years, if a salinity perturbation of 12.5 ppt was initially triggered over the upper 10 m of their model. Their results suggest that a very large surface salinity anomaly may be necessary to modify the MOC by 11% in their ocean model. These efficient methods based solely on the adjoint of the GCM can be suboptimal compared to the singular vectors, depending on the clustering of the singular values.

In Zanna et al. (2010, manuscript submitted to *Quart. J. Roy. Meteor. Soc.*), restricting the optimal temperature and salinity perturbations to the upper ocean is found to lead to a maximum growth of the MOC after 19 yr, while we find in the present study that deep density perturbations can lead to a faster maximum MOC growth within only 7.5 yr. The excitation of MOC anomalies by surface temperature and salinity in Zanna et al. (2010, manuscript submitted to *Quart. J. Roy. Meteor. Soc.*) is less efficient than exciting perturbations at depth and is found to be due to the interaction of only three nonorthogonal eigenmodes. All three modes are surface trapped, and the temperature and salinity anomalies have rather different behavior from each other over the course of the growth of MOC anomalies. One of these surface-trapped modes does exhibit some similarities with the interdecadal mode excited in the present study. However, the two additional surface salinity modes found by Zanna et al. (2010, manuscript submitted to *Quart. J. Roy. Meteor. Soc.*) are necessary to create the amplification of MOC anomalies.

This higher efficiency of deep perturbations in exciting MOC anomalies is consistent with the results of much simpler 2D models (Zanna and Tziperman 2005, 2008). In light of these results, one may expect that predictability experiments perturbing only the atmospheric state (equivalent in some sense to perturbing the upper ocean only) and leaving the ocean unchanged may lead to an overestimate of the ocean predictability time. Therefore, errors in the deep ocean—either in initial conditions or physical processes—will grow dramatically fast, affecting the potential predictability of MOC. The optimal perturbations studied here need to be taken into careful consideration when initializing models to evaluate ocean predictability on interannual and multidecadal time scales.

This study is a first step toward evaluating singular vectors for the large-scale overturning circulation in a general circulation model arising from the nonnormality of the ocean. A significant strength as well as weakness of the current study is the use of an idealized geometry and model configuration (e.g., Chhak et al. 2006). The simple geometry allowed a better understanding of the results, yet further work is needed to examine the effects of a more realistic representation of the Atlantic Ocean and its mean flow, including potential atmospheric feedbacks. We hope to address these issues in a future study.

## Acknowledgments

We wish to thank Brian Farrell, Jake Gebbie, and Carl Wunsch for very helpful discussions. We are indebted to Chris Walker for his assistance with the numerical simulations on *Swell* and *Odyssey*, and to the MITgcm group for making the model publicly available. We also thank Wilbert Weijer and an anonymous reviewer for their comments, which helped improve the manuscript. ET thanks the Weizmann Institute for its hospitality. This work was supported by the NSF’s Paleoclimate Program (Grant ATM-0455470), and by the NASA ECCO-II and the NOPP ECCO-GODAE projects.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

* Current affiliation: University of Oxford, Oxford, United Kingdom

*Corresponding author address:* Laure Zanna, Clarendon Laboratory, Atmospheric, Oceanic and Planetary Physics, Department of Physics, University of Oxford, Parks Road, Oxford OX1 3PU, United Kingdom. Email: laure.zanna@earth.ox.ac.uk