Upper-Ocean Singular Vectors of the North Atlantic Climate with Implications for Linear Predictability and Variability

The limits of predictability of the meridional overturning circulation (MOC) and upper-ocean temperatures due to errors in ocean initial conditions and model parametrizations are investigated in an idealized conﬁguration of an ocean general circulation model (GCM). Singular vectors (optimal perturbations) are calculated using the GCM, its tangent linear and adjoint models to determine an upper bound on the predictability of North Atlantic climate. The maximum growth time-scales of MOC and upper-ocean temperature anomalies, excited by the singular vectors, are 18.5 and 13 years respectively and in part explained by the westward propagation of upper-ocean anomalies against the mean ﬂow. As a result of the linear interference of non-orthogonal eigenmodes of the non-normal dynamics, the ocean dynamics are found to actively participate in the signiﬁcant growth of the anomalies. An initial density perturbation of merely 0.02 kg m − 3 is found to lead to a 1.7 Sv MOC anomaly after 18.5 years. In addition, Northern Hemisphere upper-ocean temperature perturbations can be ampliﬁed by a factor of 2 after 13 years. The growth of upper-ocean temperature and MOC anomalies is slower and weaker when excited by the upper-ocean singular vectors than when the deep ocean is perturbed. This leads to the conclusion that predictability experiments perturbing only the atmospheric initial state may overestimate the predictability time. Interestingly, optimal MOC and upper-ocean temperature excitations are only weakly correlated, thus limiting the utility of SST observations to infer MOC variability. The excitation of anomalies in this model might have a crucial impact on the variability and predictability of Atlantic climate. The limit of predictability of the MOC is found to be different from that of the upper-ocean heat content, emphasizing that errors in ocean initial conditions will affect various measures differently and such uncertainties should be carefully considered in decadal prediction experiments.


Introduction
The past century North Atlantic surface climate variability is often described as a combination of an interannual and a multi-decadal mode (e.g. Deser and Blackmon, 1993). The analysis of the limited available observations led to the hypothesis that interannual variability is mainly driven by the atmosphere, while multi-decadal variability is mostly due to ocean dynamics and its large-scale overturning circulation (e.g. Bjerknes, 1964;Kushnir, 1994). The Atlantic meridional overturning circulation (MOC), defined as the zonally averaged meridional volume transport and forced by both wind and buoyancy fluxes, transports large amounts of water and heat from low to high latitudes and therefore plays a crucial role in surface climate variability.
Many studies involving numerical models show potential predictability of the North Atlantic sea surface temperatures (SSTs), upper-ocean heat content and MOC on decadal timescales (Griffies and Bryan, 1997;Boer, 2000;Pohlmann et al., 2004;Sutton and Hodson, 2005;Collins et al., 2006). These studies are performed by constructing model ensembles in which each run corresponds to a slightly perturbed atmospheric state while leaving the ocean state unchanged. The spread of the individual model trajectories gives a measure of predictability. Since the ocean initial conditions are not perturbed, such results can be considered an upper limit on the predictability. One therefore wonders how realistic this predictability limit is and how errors in ocean initial conditions or in model parametrizations can affect the predictability in the North Atlantic region. While the development of initialization of decadal prediction systems has been the focus of many groups over the past several years and some improvement in this area has been made (Smith et al., 2007;Keenlyside et al., 2008), our understanding of the processes governing the error growth of anomalies is still incomplete.
Due to the non-normality of the ocean dynamics, the linear interference of non-orthogonal eigenmodes may result in a rapid transient growth even if the system is linearly stable (e.g. Farrell, 1988Farrell, , 1989. Such transient growth of anomalies can therefore provide a limit on the predictability of the large-scale ocean circulation, upper-ocean heat content and SST fluctuations, before the nonlinearities become important. The purpose of this study is to investigate the variability and predictability limits (in a linear sense) of the North Atlantic Ocean to upper-ocean perturbations and air-sea forcing in an idealized configuration of a three-dimensional ocean general circulation model (GCM). More specifically, we explore the fastest growing (optimal) perturbations, constrained to the upper ocean, leading to the largest amplification of MOC and upper-ocean temperature anomalies. The anomalies leading to the largest excitation are given by the leading singular vector. The analysis of the leading singular vectors or optimal perturbations provides information on (i) the sensitivity of MOC and upper-ocean temperatures to perturbations, (ii) the physical mechanism by which the perturbations grow and impact the North Atlantic climate, and (iii) the growth time-scale relevant for climate variability and predictability.
In other words, the leading singular vectors help identifying regions where observations could be targeted to improve our understanding of climate variability and our ability to predict these fluctuations, provide linear estimates of uncertainties in ocean initial conditions and model parametrizations, and provide time-scales for error growth limiting the predictability in the region.
The related sensitivity of the MOC and meridional heat transport to initial conditions and surface forcing has been investigated in 3D ocean models using adjoint methods (Marotzke et al., 1999;Bugnion et al., 2006;Sevellec et al., 2008;Czeschel et al., 2010;Heimbach et al., 2011). However the singular vector approach provides a stricter bound on error growth than an adjoint sensitivity experiment or excitation by adjoint modes. Singular vector analysis has been successfully used for example in numerical weather and El Niño-Southern Oscillation prediction (e.g. Buizza, 1995;Penland and Sardeshmukh, 1995;Moore and Kleeman, 1997) and applied to various climate variability problems (e.g. Moore, 1999) including the MOC (e.g. Lohmann and Schneider, 1999;Tziperman and Ioannou, 2002;Tziperman et al., 2008;Alexander and Monahan, 2009;Hawkins and Sutton, 2009).
In a previous study, the authors used an idealized configuration of the Massachusetts Institute of Technology ocean GCM (MITgcm) and its combined tangent linear and adjoint to calculate the singular vectors and explore the role of the non-normal ocean dynamics in MOC variability and predictability (Zanna et al., 2011, hereafter ZHMT11). The largest amplification of MOC anomalies resulted from deep density perturbations at high latitudes, with a growth timescale of about 7.5 years. The amplification of anomalies, reminiscent of baroclinic instability converting available potential energy into kinetic energy, involved a cyclonic propagation of deep temperature and salinity anomalies. The density anomalies propagated under the influence of the velocity of the mean flow and the steady state density gradient, determining the growth time-scale of the MOC anomalies and therefore an upper bound on the MOC predictability time (error growth time-scale).
While eddies, deep convection and overflows can excite deep density anomalies, the atmosphere is likely to excite perturbations at the surface of the ocean. In the present study, we investigate the non-normal excitation of the MOC and upper-ocean temperature dynamics due to upperocean forcing alone, using the same model and underlying methodology presented in ZHMT11. The optimal excitation of MOC and upper-ocean temperature variability by atmospheric forcing is of great interest. An important question is whether ocean feedbacks can amplify surface temperature perturbations or whether the ocean can merely dampen the atmospheric thermal forcing. We find that MOC variability can be excited by upper-ocean initial temperature and salinity on a time-scale of 18.5 years due to the interference of three eigenmodes of the linearized dynamics. Similar initial conditions of upper-ocean temperature to those identified when maximizing the MOC are found to amplify the North Atlantic upper-ocean temperature by a factor of 2 on a time-scale of 13 years due to the interaction of two eigenmodes. The implications for predictability of the North Atlantic Ocean are discussed by focusing on the time-scales for error growth.
We briefly describe the model and the ocean steady state in section 2. In section 3, we present the methodology used to evaluate the singular vectors with the MITgcm. The initial perturbations of upper-ocean temperature and  50, 70, 100, 140, 190, 240, 290, 340 , 390, 440, 490, 540, 590, 640 salinity found to maximize the MOC anomalies and the associated growth mechanism are discussed in section 4. This is followed by a similar analysis of upper-ocean temperature perturbations leading to the growth of North Atlantic temperature anomalies in section 5. We conclude and discuss the possible implications of this study for variability and predictability of the North Atlantic Ocean in section 6.

The MITgcm and the ocean steady state
In this work, we use the MITgcm (Marshall et al., 1997b,a) in an idealized Atlantic-like double-hemispheric rectangular basin configuration between latitudes 67.5 • S and 67.5 • N and longitudes 65 • W and 5 • W with a 3 • horizontal resolution, as in ZHMT11. The ocean depth is uniformly set to 5200 m with 15 vertical levels of thicknesses varying between 50 and 690 m. No-normal flow and no-slip conditions are prescribed at the bottom and side boundaries. The prognostic variables are the horizontal velocity v = (u, v), potential temperature T and salinity S, with all other quantities being diagnosed. The equation of state is a modified UNESCO formula (Jackett and McDougall, 1995), and ocean convection is parametrized by an implicit vertical diffusion scheme. The model is forced with timeindependent wind and mixed boundary conditions for temperature and salinity: 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 I. The steady state reached by the model 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 nearsurface horizontal velocity field is shown in Figure 1 with a strong subtropical gyre and a relatively weak and small subpolar gyre north of 60 • N. Figure 2 shows the steady state meridional overturning streamfunction with a strong Northern Hemisphere asymmetric cell. The overturning cell has a maximum transport of about 22 Sverdrup (Sv; 1 Sv= 10 6 m 3 s −1 ). While the steady state reached by the model is fairly reasonable, the idealized geometry and forcing perhaps lead to some unrealistic aspects of the ocean state. We should keep in mind that a different steady state and geometry may quantitatively affect the results. Further details on the model and its steady state are provided in ZHMT11.

Evaluation of upper-ocean singular vectors
The tangent linear model derived from the MITgcm primitive equations, which may be represented by and its adjoint (Giering and Kaminski, 1998;Marotzke et al., 1999;Heimbach et al., 2005), are used to compute the singular vectors under different norm kernels. The matrix A is the time-independent linearized dynamical operator obtained via the automatic differentiation tool TAF (Giering and Kaminski, 1998;Heimbach et al., 2005) and P (t) is a small perturbation from the steady state solution for temperature, salinity and horizontal velocity such that The matrix B(t) ≡ e At is the propagator representing the tangent linear model, and P 0 ≡ P (t = 0) is the initial perturbation.
To calculate the singular vectors leading the maximum growth of the MOC or upper-ocean temperature anomalies, we proceed by maximizing subject to initial unit norm ||P 0 || = 1. The norm kernel X reflects the quantity to be maximized, which in the present study is either the sum of squares of MOC anomalies or the sum of squares of upper-ocean temperature anomalies. Using the sum of squares of anomalies permits us to investigate simultaneously the impact of upper-ocean perturbations on the variability and error growth limiting the predictability, since X can be viewed as the error variance of the model forecast for MOC or upper-ocean temperature fluctuations (Lorenz, 1982). The maximization of Eq. (3) is performed for lead times τ varying from 1 to 40 years, while constraining the initial temperature and salinity perturbations to have a unit norm, reflecting the available potential energy of the upper-ocean perturbations, λ is the longitude, ϕ the latitude, z the depth, α = α 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. The initial perturbations are constrained to the upper ocean by setting h = 360 m. Entries of the norm kernel E corresponding to anomalies below h = 360 m are set to zero. The norm kernel 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. The initial kinetic energy of the perturbations is neglected in all calculations as it does not contribute to the non-normal growth for time-scales longer than a few weeks. The velocity and sea surface height perturbations adjust to the density field within a few days, leading to a horizontal flow in geostrophic balance with the density gradients. Using only temperature and salinity perturbations is therefore sufficient to explore growth on interannual and decadal time-scales (Zanna et al., 2010.
Maximizing anomalies as measured by Eq.
(3) at t = τ with respect to the initial perturbations of temperature and salinity anomalies, using Eq. (4) as the initial constraint, is equivalent to solving the generalized eigenvalue problem given by (Farrell, 1988;Zanna and Tziperman, 2005) The optimal initial conditions P 0 are defined as the fastest growing perturbations corresponding to the generalized eigenvector of Eq. (5) with the largest eigenvalue γ (Farrell and Ioannou, 1996). In this case, the optimal perturbations correspond to the rescaled first singular vector of X 1/2 B(τ )E −1/2 at time τ . For clarification, it is worth mentioning that the term optimal perturbations in this work refers to the leading singular vector, similar to early papers discussing the subject; the Introduction gives references. The initial perturbations P 0 were computed using the Lanczos algorithm (Golub and Van Loan, 1989) and the routines for symmetric eigenvalue problems of the Arnoldi Package software (Lehoucq et al., 1998) which requires only the input vector B T XBP , where the superscript T denotes the matrix transpose (equivalent to the adjoint with respect to a L2-norm). The eigenvalues and eigenvectors of the tangent linear B and adjoint models B T are evaluated through the same procedure except that non-symmetric routines are used. The system is found to be linearly stable such that all eigenvalues of A have negative real parts and every perturbation eventually decays as t goes to infinity. A complete description of the numerical calculation of singular vectors can be found in Moore et al. (2004) and additionally in Zanna et al. (2010) specifically for the MITgcm. In Zanna et al. (2010), the algorithm was found to be robust to the MITgcm model resolution and to other model assumptions.

Maximizing MOC anomalies
In order to determine the effect of atmospheric forcing on the MOC excitation, and estimate the growth time-scale of anomalies limiting the predictability of the MOC due to nonnormal ocean dynamics, we consider initial perturbations of temperature and salinity constrained to the upper 360 m of the ocean. The quantity (Eq. (3)) to be maximized is given by J MOC (τ ) is the sum of squares of MOC anomalies at time τ north of the Equator (as in ZHMT11) and dA = rdϕ dz(z) is the cross-section area element. The meridional overturning streamfunction anomaly ψ is given by where H is the ocean depth. Note that a sensitivity study to the choice of the domain shows that the main results are not overly sensitive to this choice (ZHMT11).  The generalized eigenvalue problem defined by Eq. (5) is solved for lead times τ between 1 and 40 years. The leading singular value γ , reflecting the largest possible amplification of the MOC anomalies squared at a given lead time τ , is shown in Figure 3. The grey curve shows the leading singular value when solving the generalized eigenvalue problem for MOC anomalies with the optimal perturbations constrained to the upper ocean, while the black curve shows the leading singular value when allowing the perturbations to be at any depth (numerical experiments from ZHMT11). The maximum amplification occurs around 18.5 years when the singular vectors are constrained to the upper ocean, compared with 7.5 years when they are allowed over the entire ocean depth (ZHMT11).

Spatial structure of the leading singular vector
Consider the leading singular vector of upper-ocean perturbations, often called (global) optimal perturbations, resulting in the maximum MOC amplification (corresponding to the largest γ ) and occurring for τ =18.5 years. The upper-ocean depth-averaged optimal temperature and salinity perturbations weighted by their respective contributions to density, defined as  Figures 4(a, b). Both temperature and salinity patterns appear to contribute to the initial optimal density anomaly field (Figure 4(c)). The signal, with relatively large amplitude at high latitude, is located in the northern part of the basin (similar to ZHMT11), primarily at the boundary between the subtropical and subpolar gyres and in the subpolar gyre. The sensitivity of the MOC to high-latitude perturbations is consistent with other studies (e.g. Tziperman et al., 2008;Hawkins and Sutton, 2009;Czeschel et al., 2010;Heimbach et al., 2011;Zanna et al., 2011). The perturbations are located near the downwelling branch of the steady state MOC and in the region where the slumping of the isopycnals is the steepest (ZHMT11). Note that the optimal perturbation patterns are to some extent a mathematical construction using the physical properties of the dynamical system to initially minimize the MOC and energy anomalies in order to excite the largest possible growth over 18.5 years.

Eigenmodes of the linearized propagator
The growing singular vectors are evidence of linear interference of non-orthogonal decaying eigenmodes which give rise to the transient growth of MOC anomalies which can be further investigated. Let us denote by s i and r i the eigenvectors of B and of its adjoint (Hermitian transpose) B † , respectively. The initial perturbation P 0 can be written as a linear superposition of all the eigenmodes of B, such that The coefficients a k , which are the projections of P 0 onto the eigenmodes s i , are given by a k = (r T k P 0 )/(r T k s k ). The set of largest coefficients a k corresponds to the main eigenmodes contributing to the initial conditions and participating in the transient amplification. A linear superposition of mainly three eigenmodes of the linearized dynamical model, responsible for the MOC growth during the initial 18.5 years, was identified. Note that the eigenmodes are not the leastdamped modes of B, indicating that mode selection is important and that a simple adjoint mode calculation would not be optimal to pick out the growth of anomalies. Some of the properties of the three eigenmodes are summarized in Table II. One eigenmode (mode 1) is found to be mainly characterized by temperature anomalies, while the two additional eigenmodes (modes 2 and 3) are dominated by salinity anomalies.
Given that the eigenmodes are oscillatory and that the norm kernel X MOC (defined by Eq. (6)) does not span the entire space, small transient growth can occur without the participation of non-normality of the ocean dynamics. A useful measure of the non-normality of a particular eigenmode s k is given by (Farrell and Ioannou, 1999) ν provides a measure of the projection of s k on the remaining eigenmodes of B, where L is such that X = LL T . Therefore the larger the value of ν, the higher is the degree of nonnormality. The values of ν for the three eigenmodes of the linearized dynamical model responsible for the MOC growth are given in Table II. The value of ν for each eigenmode is at least five times larger than the value of ν for the remaining eigenmodes. Therefore the three eigenmodes responsible for the MOC growth are significantly more non-normal than the remaining elements of the eigenspectrum. Using only a linear combination of these three highly non-normal eigenmodes results in an amplification timescale of 17 years for the MOC anomalies and a growth factor similar to the one obtained with the singular vector (within 10% error). Therefore the optimal growth of the upper-ocean singular vectors can be mainly explained in terms of the three non-normal eigenmodes identified. To explain the growth of MOC anomalies, which is mainly determined by the evolution of the zonal density gradient via the thermal wind relation, we must consider separately the propagation of the salinity and of the temperature during the amplification for each eigenmode. Considering the temperature and salinity separately in the current study differs from ZHMT11 in which the evolution of the temperature and salinity anomalies were well correlated and the density could be treated as the controlling variable.
The first eigenmode (mode 1) is characterized by temperature with a structure similar to the upper-ocean optimal temperature anomalies in the Northern Hemisphere shown in Figure 4(a). The salinity component of this eigenmode is negligible and the evolution of the salinity does not contribute to the density and MOC growth (Table II). In the upper ocean, the temperature anomalies propagate westward around 60 • N ( Figure 5(a)) against the eastward mean flow (Figure 1). The westward propagation can be explained by diagnosing the different terms in the linearized temperature equation. The different terms are averaged over time and volume and their magnitudes are found to be: Retaining only the dominant terms in the temperature equation, the time evolution of the temperature anomalies is therefore determined by where the steady state quantities are denoted by the bar and the primes are the perturbations from the steady state. Using the thermal wind relation to rewrite v , we obtain under the assumptions that the anomalies are in geostrophic balance and spreading over a depth h of 360 m (due to our initial constraint for the singular vectors), and neglecting the small effect of salinity compared to that of the temperature on the density anomalies. The zonal velocity of propagation for the temperature is therefore given approximately by For the temperature to propagate westward against the mean flow, the necessary condition is Colin De Verdiere and Huck (1999) and Te Raa and Dijkstra (2002) analyzed an unstable oscillatory eigenmode of surface-trapped temperature anomalies with a similar spatial structure and westward propagation (e.g. compare Figure 4 here with Figure 4 in Te Raa and Dijkstra, 2002) to explain self-sustained MOC variability. In our model regime, the oscillatory eigenmode is damped with a decay time-scale of 68 years and period of about 34 years (Table II) and the interaction with the two additional salinity modes is necessary to create the amplification of MOC and energy anomalies.
The leading singular vector projects onto two additional eigenmodes which are characterized by upper-ocean salinity anomalies (Table II). One eigenmode (mode 2) with anomalies near the western part of the basin ('1' in Figure 4(b)) has a decay time-scale of 16 years. The different terms in the salinity equation for mode 2 are evaluated and found to be . The time evolution of the salinity anomalies of mode 2 is therefore determined by The third and last eigenmode participating in the optimal growth (mode 3) has a decay time of 4 years and oscillation period of 10 years with salinity anomalies east of 45 • W ('2' and '3' in Figure 4(b)). The different terms in the salinity equation of mode 3 are evaluated and found to be: . The time evolution of mode 3 is thus determined by The sum of modes 2 and 3 obviously closely resembles the spatial pattern of optimal initial conditions of salinity shown in Figure 4(b) and their evolution is therefore given by The overall decay of salinity anomalies is relatively fast compared to the decay of temperature anomalies associated with mode 1 previously described. For all the modes considered, mixing and surface damping of temperature anomalies are required to explain the decay of the optimal perturbations, but do not actively participate in the propagation and are therefore neglected for simplicity.

Transient amplification of MOC anomalies
The main growth of MOC anomalies is explained by examining the time evolution of zonal density gradients related to the MOC via the thermal wind relation.
The first peak at about 6 years ( Figure 6, grey curve) is mostly driven by a density gradient dominated by salinity anomalies, yet with a contribution of the temperature gradient (in particular due to the temperature anomaly located near the western boundary). The stationary negative temperature anomalies in the northwestern region decay over the initial 5 to 6 years. The positive temperature anomalies initially situated in the northeastern region of the basin around 20 • W (Figure 4  The evolution of the salinity anomalies is rather different. The salinity anomalies located in the northeastern part of the basin, especially near the boundaries ('3' in Figure 4(b)), are rapidly damped by the advection of the salinity anomalies by the mean vertical velocity, w∂S /∂z. However, the salinity anomalies located in other parts of the basin ('1' and '2' in Figure 4) are slowly advected by the eastward mean flow (see the direction of the mean velocity in Figure 1) toward the eastern boundary ( Figure 5(b)). Decay and propagation lead to a nearly vanishing salinity anomaly in the vicinity of the western boundary, and to a negative salinity anomaly near the eastern boundary. The newly emerged density gradient at around 6 years ( Figure 5(c)), creates a vertical shear of meridional velocity anomaly, therefore explaining the local maximum in the timeseries of the MOC (Figure 6).
Consider the main MOC peak at maximum amplification time (about t = 18.5 years, Figure 6) which is entirely dominated by the temperature signal. As already mentioned, the positive temperature anomalies initially situated around 20 • W at high latitudes (Figure 4(a)) propagate westward. The anomalies start reaching the western boundary after 10 years, yet the large amplitude of temperature anomalies continues to build over the following 8 to 9 years. As a result of the westward propagation of the warm anomaly, the zonal temperature gradient anomaly becomes negative, with a maximum amplitude around 18-19 years. Consider now the evolution of the salinity anomalies. The overall decay of the salinity anomalies in the basin leads to a vanishing salinity anomaly near the eastern boundary after 16 years. Therefore the negative zonal temperature gradient induces a negative meridional overturning flow after 18.5 years ( Figure 5(d)), which is further enhanced by the decay of the salinity near the eastern boundary. Note that the meridional anomalous advection of the mean salinity by the perturbations, v ∂S/∂y, is negligible in the time evolution of salinity anomalies. This explains why the salinity does not show a clear westward propagation against the mean flow ( Figure 5(b)) during the amplification of the MOC anomalies, unlike the temperature anomalies ( Figure 5(a)). The amplification time-scale of the MOC is therefore determined by the travel time of the initial warm anomaly toward the western boundary and the decay time-scale of the salinity modes. Since the model is linearly stable, the eigenmodes decay at longer times and so do the MOC anomalies ( Figure 6).

Energy growth and further evidence of non-normal effects
In an asymptotically stable normal system, the energy of the initial perturbations always decays and no growth of anomalies can be supported. However, due to the presence of asymmetries and shear, the ocean dynamics are non-normal and the energy of perturbations can grow significantly. Given that the eigenmodes participating in the growth are oscillatory and that the norm kernel X MOC , defined by Eq. (6), does not span the entire space, small transient growth can occur without the participation of the non-normality of the ocean dynamics. We have already described the growth as a linear combination of highly non-normal eigenmodes, showing that the fast decay of the salinity anomalies and the advection of the temperatures anomalies lead to a large MOC anomaly after 18.5 years. In order to further investigate the role of the non-normality in the growth of perturbations, we have diagnosed the energy growth in the basin when initializing the model with the upper-ocean leading singular vector. Figure 6 shows the square root of the growth of the MOC anomalies defined by J MOC and the growth of the energy of the anomalies (E, defined as a sum of the kinetic and available potential energy) when the model is initialized with the optimal perturbations. The sum of the kinetic energy and the available potential energy reaches a maximum around 21 years with a growth factor of roughly 17 ( Figure 6, black curve). Therefore the overall anomalies in the model grow as a function of time and roughly correspond to the growth of MOC anomalies, which can only occur due to the non-normality of the dynamical operator. Similarly to ZHMT11, we find that the growth of perturbations is mainly due to the source term ρ u · ∇ρ, the extraction of potential energy from the mean stratification, acting against the damping of perturbations. This source term is interpreted as the change of available potential energy due to the interaction of buoyancy perturbation and the anomalous buoyancy advection (Huang, 2002). We conclude that, despite the participation of oscillatory modes which can cause resonance and lead to small to moderate growth, the bulk of the inferred amplification of all variables shows that the excitation of MOC anomalies is primarily due to the non-normal dynamics.

Implications for predictability and variability of the MOC
The ability of the system to undergo rapid transient amplification has important implications for how the North Atlantic climate will respond to atmospheric or oceanic stochastic noise, but also for its predictability. The model steady state is identical to that of ZHMT11, and we can therefore compare the present results with the results from ZHMT11. In both studies, the leading singular vectors are characterized by perturbations at high latitudes in the Northern Hemisphere. The relatively high sensitivity of the model to high latitude perturbations suggests that the MOC variability can be efficiently excited on interannual to decadal time-scales by high-latitude stochastic noise at the surface (i.e. due to heat, freshwater and wind forcing) but also at depth by overflows, convection or eddies (ZHMT11). The growth of these perturbations suggests that errors in ocean initial conditions or model parametrizations at high latitudes could potentially limit the predictability of the MOC and the North Atlantic climate. Additional observations at high latitudes could be crucial to improve our predictions (e.g. Tziperman et al., 2008;Hawkins and Sutton, 2009) but also improvements in non-explicitly resolved ocean processes such as convection, overflows or eddies.
The growth time-scale of the MOC anomalies, excited by upper-ocean perturbations of temperature and salinity anomalies and analyzed in this study, can be explained by the fast decay of the two salinity modes and relatively slow decay of the temperature mode. This time-scale of 18.5 years is longer than the 7.5 years obtained when the perturbations are allowed over the entire ocean depth (ZHMT11). This result implies that the predictability time-scales of 10 to 20 years obtained when only atmospheric perturbations are used to initialize ensemble experiments (e.g. Griffies and Bryan, 1997;Pohlmann et al., 2004) may be overestimates. In addition to the difference in growth time-scales, the MOC anomaly appears to be less sensitive to upper-ocean perturbations than to deeper ones. We find here that a density perturbation of 0.02 kg m −3 in the upper ocean leads to an MOC anomaly of 1.7 Sv compared to 2.4 Sv in ZHMT11 when the anomalies are mostly located in the deep ocean.
It seems also that interannual ocean variability can be excited by deep density anomalies while multi-decadal variability can be excited by near-surface anomalies due to the non-normal dynamics in this ocean model. A simple (but only partial) explanation of the difference in the MOC anomalies growth time-scale between the present study and ZHMT11 is given by the amplitude of the westward velocity of propagation. At depth, the westward mean flow, u < 0, and the meridional density gradient, ∂ρ/∂y > 0, acted to amplify the magnitude of the westward propagating signal, given by hence the relatively rapid time-scale of 7.5 years found in ZHMT11. In the current study, the perturbations near the surface propagate against the mean flow, such that there is a competing effect of the meridional density gradient and the eastward mean flow slowing down the westward propagation of the anomalies, leading to a relatively longer time-scale of 18.5 years compared to ZHMT11.

Excitation of upper-ocean temperatures
The ocean is a major player in determining the heat content of the climate system due to its large heat capacity. The ocean dynamics can induce large fluctuations in the ocean heat content on a decadal time-scale, as seen in the upper layers of the ocean (Levitus et al., 2004) with spatial and temporal variations (Lozier et al., 2008). Many studies have shown that MOC variability is associated with changes in the meridional transport of heat to high latitudes and therefore leads to SST and upper-ocean heat content variability with implications for the potential predictability of North Atlantic climate on decadal time-scales. Model experiments find that SST fluctuations are lagging the MOC by roughly 5 to 10 years (Pohlmann et al., 2004), suggesting that upperocean temperatures could be used to reconstruct the MOC (without using the salinity). Recent studies have analyzed the benefit of using the temperature record to initialize decadal prediction systems (Keenlyside et al., 2008). This raises several important questions: To what extent are changes in upper-ocean temperatures associated with MOC changes? How useful is the upper-ocean temperature record in determining the predictability of the North Atlantic climate? How do uncertainties in upper-ocean initial conditions of temperature limit the predictability of North Atlantic temperatures? Our numerical analyses can also serve as a preliminary experiment for investigating the optimal perturbations of North Atlantic SSTs using the observed record (Zanna, 2011).

Maximizing upper-ocean temperature anomalies
In this section, we calculate the singular vectors of upperocean temperatures which maximize the sum of squares of temperature anomalies in the upper 360 m in the Northern Hemisphere at time t = τ . The upper-ocean temperature anomalies to be maximized are given by where dV(λ, ϕ, z) is the volume of each grid cell. The initial conditions are constrained to satisfy J T (0) = 1. The upper-ocean temperatures are maximized for different lead times τ and the largest amplification is found to be for τ = 13 years. The leading upper-ocean singular vector resulting in the maximum growth of upper-ocean temperatures is used to initialize the linearized model and a timeseries of the growth of J 0.5 T is shown in Figure 7 (grey curve). The sum of squares of the temperature, J T , is amplified by a factor of about 4 due to the non-normal ocean dynamics, and the temperature itself is amplified by a factor of roughly 2 (Figure 7). The evolution at different times of the upper-ocean temperature as function of latitude and longitude is shown in Figure 8. The leading singular vector of temperature is shown in Figure 8(a). The temperature pattern in Figure 8(a) shows some similarity with the one found to maximize the MOC anomalies in Figure 4(a). The signal is again concentrated at high latitudes, although a relatively cold anomaly is present in the northwestern part of the basin ('1' in Figure 8(a)).

Transient growth of the upper-ocean temperatures
We can describe the amplification mechanism in two complementary ways, first by analyzing the dominant terms leading to the rate of change of J T , and second by discussing the interference of two non-orthogonal eigenmodes in building up the temperature anomalies in the upper ocean north of the Equator in this model. The full linearized equation for temperature anomalies is given by where κ is the diffusion coefficient and τ T the SST restoring time-scale (such that τ T −1 vanishes below the surface level). Multiplying Eq. (11) by T and integrating over the volume defined in Eq. (10), we obtain 1 2 where ϕ=4.5N x z(0) dA , with dA being the cross-section area element and z(0) the upper-level thickness.
During the 13 years of growth of J T , the dominant term near the surface is − T u · ∇T V . All other terms in Eq. (12) are several orders of magnitude smaller than T u · ∇T V . The interaction between the temperature perturbations T v ∂T/∂y V , is found to be particularly large north of 35 • N. In the western part of the basin, the interaction of temperature anomalies and vertical temperature advection, T w ∂T/∂z V , is found to be of a similar magnitude to T v ∂T/∂y V . It indicates that the sum of squares of upper-ocean temperature growth is driven by energy being extracted from the mean meridional temperature gradient and mean temperature stratification. Consider next the perspective based on the interference of eigenmodes resulting in the amplification of J T . Similarly to section 4, the leading singular vector for τ = 13 years is projected onto the eigenmodes of the linearized propagator B. Two oscillatory decaying eigenmodes of temperature anomalies were found to participate in the transient growth of upper-ocean temperature anomalies. The modes are again highly non-normal and their evolution is dominated by the temperature anomalies with a negligible contribution of the salinity anomalies over the time period considered (Table III).
One eigenmode (Mode 1 in Table III) is similar to the oscillatory mode involving the westward propagation of temperature as discussed in section 4, with an e-folding time of 68 years and a period of about 34 years. The time evolution of the eigenmode excited is dominated by u∂T /∂x and v ∂T/∂y, leading to a westward propagation of the warm anomaly ('2' in Figure 8) with a westward velocity The time evolution of mode 1 is therefore determined by Eq. (8).
A second decaying oscillatory mode (Mode 2 in Table III) is characterized by the cold anomaly in the western part of the basin ('1' in Figure 8). The time evolution of mode 2 with an e-folding time of 15 years and a period of about 54 years is dominated by the advection of zonal and vertical temperature gradients of the mean flow by the velocity perturbation, namely u ∂T/∂x and w ∂T/∂z. All other terms diagnosed in the temperature equation are four to five orders of magnitude smaller than u ∂T/∂x and w ∂T/∂z. The time evolution of the temperature anomalies for mode 2 is therefore given by Consider now the evolution of the different anomalies ('1' and '2' in Figure 8) and the growth of temperature anomalies. The evolution of the cold anomaly ('1' in Figure 8) in the northwestern part of the basin is driven by the evolution of mode 2. The cold anomaly, initially in the northwestern region of the basin ('1' in Figure 8(a)), moves southward. In the western part of the basin, a positive meridional density gradient anomaly (∂ρ /∂y > 0) is created due to the presence of this cold anomaly. The thermal wind relation induces a positive meridional density gradient anomaly and a positive vertical zonal shear anomaly (∂u /∂z > 0) in the western part of the basin at high latitudes associated with upwelling near the western boundary. The creation of an anomalous upwelling velocity therefore enhances the cold anomaly ('1' in Figure 8(b)) as it moves southward during the initial 4 years.
Simultaneously, the warm initial anomaly ('2' in Figure 8(a)) creates an anomalous anticyclonic geostrophic flow with northward velocity to the west of the anomaly and southward velocity to the east. As the temperature anomaly propagates westward, southward velocity perturbations induced by the temperature anomaly are replaced by a northward velocity anomaly to the west of the warm anomaly. The southward velocity anomaly to the west of the warm anomaly creates a negative advection of the gradient of mean temperature by the meridional velocity perturbations, such that v ∂T/∂y < 0, which therefore enhances the positive temperature anomaly ('2' in Figure 8(b c)) as it travels westward. As the higher temperature propagates and reaches the vicinity of the western boundary, it creates a negative meridional density gradient in the northwestern part of the basin. This reverses the sign of ∂u /∂z (which was positive around t = 4 years) creating a negative vertical velocity anomaly in the western part of the basin. Therefore w ∂T/∂z is now negative and relatively strong such that the amplitude of anomalies is further enhanced (Figure 8(b, c)).
We can summarize the upper-ocean temperature growth as follows. The initial warm anomaly in the northeastern part of the basin propagates westward as a Dopplershifted Rossby wave (Eq. (8)) with the mean meridional temperature gradient serving as the background meridional vorticity gradient. This anomaly induces an anticyclonic circulation around it, which moves with the anomaly further amplifying the temperature anomaly by advecting the mean temperature gradients (meridional and vertical) when reaching the vicinity of the western boundary. Therefore, the temperature anomaly increases during the first 13 years due to v ∂T/∂y and w ∂T/∂z via the interaction between two non-normal eigenmodes.
The growth of upper-ocean temperatures has an impact on the mid-and high-latitude circulation (horizontal and vertical). However, despite the similarities between the initial conditions of temperature found when maximizing the MOC anomalies and those found when maximizing the North Atlantic upper-ocean temperatures, the excitation of upper-ocean temperature anomalies is not well correlated with the MOC variations at any lag between 1 and 10 years ( Figure 9). The lack of initial salinity anomalies in the upperocean singular vectors could possibly explain the differences between the mechanisms for the growth of MOC and upperocean temperatures. The upper-ocean temperature growth influences to some extent the MOC and vice versa, although not dramatically. Therefore large amplification of upperocean temperatures can occur without the participation of the MOC. This implies that the use of upper-ocean temperature fluctuations for reconstructing the past MOC may perhaps be more difficult than believed by other studies (e.g. Pohlmann et al., 2004), at least in the current configuration of the MITgcm, and possibly in other models, and in the real ocean.
In terms of predictability, the time-scale of upperocean temperature growth is relatively slow (13 years), suggesting that uncertainties in upper-ocean temperatures alone, compared to deep perturbations for example, may not be the dominant barrier to making decadal predictions of upper-ocean temperatures. This is further investigated by evaluating the upper-ocean singular vectors of temperature and salinity maximizing the upper-ocean temperatures (Eq. (10)). The largest amplification is reached after 11 years and the amplification factor is roughly 10. If upper-ocean temperatures are maximized allowing full-depth singular vector perturbations, the maximum amplification occurs after 5 years and is about 15. The perturbations are then located in the upper 2 km of the ocean basin. As for the MOC growth, deep perturbations of temperature and salinity limit the predictability of upper-ocean temperatures in this model.

Discussion and conclusions
We examined the active role of the ocean in amplifying the meridional overturning circulation (MOC) and upperocean temperature anomalies due to excitation of optimal initial temperature and salinity anomalies (given by the leading singular vectors) constrained to the upper ocean in an idealized configuration of the MITgcm. We first computed optimal initial perturbations of upper-ocean temperature and salinity maximizing MOC anomalies. The results showed an amplification of MOC anomalies on a time-scale of 18.5 years. By comparing the current results to ZHMT11, where deep initial perturbations were also allowed and the growth time-scale was found to be only 7.5 years, we concluded that deep density perturbations can yield to a faster and more efficient growth of MOC anomalies than upper-ocean perturbations. This conclusion is supported by 2D model results (Zanna andTziperman, 2005, 2008) and indicates that predictability experiments in which only the atmospheric state is perturbed (equivalent to perturbing the upper ocean only) may strongly overestimate the ocean predictability time. While eddies, deep convection and overflows can excite deep density anomalies, the atmosphere is likely to primarily perturb the surface of the ocean. Therefore, both deep and near-surface ocean perturbations need to be taken into careful consideration when initializing models to evaluate ocean predictability on interannual and multi-decadal time-scales.
We then proceeded to the analysis of optimal excitation of upper-ocean temperature anomalies. We found that North Atlantic upper-ocean temperature anomalies can be amplified by a factor of 2 within 13 years. The optimal perturbations found are very similar to the ones obtained when maximizing MOC anomalies. Despite this similarity, upper-ocean temperature anomaly growth is not associated with MOC growth. Hasselmann (1976) and Frankignoul and Hasselmann (1977) showed that the heat capacity of the upper ocean acts to passively integrate atmospheric stochastic forcing and therefore amplify the low-frequency response of SST. The amplitude of upper-ocean temperature anomalies in this scenario is limited by the atmospheric forcing amplitude. In contrast to their scenario, we showed in this study that ocean dynamics play an active role in amplifying temperature anomalies in a North Atlantic-like ocean on relatively long time-scales.
While we were able to identify distinct growth mechanisms for the upper-ocean temperature and MOC anomalies, the idealized geometry and model configuration used here make the application of our conclusions to the analysis of North Atlantic observations difficult. However, we note that a preliminary analysis of observed optimal SST perturbations using linear inverse modelling (e.g. Penland and Sardeshmukh, 1995) in the North Atlantic finds an amplification by a factor of 1.5 to 2 after about 4 to 6 years, and a westward propagation in the northern part of the basin similar to our findings (Zanna, 2011).
The singular vectors are solutions of a linear problem, therefore a multiplication by any (positive or negative) constant will also be a solution to the eigenproblem. In the above analysis, we arbitrarily picked one sign in order to explain the growth mechanisms. Moreover, the amplitudes used for the initial perturbations to evaluate the response of the non-normal dynamics are comparable to estimates of ocean variability (Forget and Wunsch, 2007). To check that the linear mechanism is valid over such long time-scales (at least in our model), we initialized the full nonlinear model with the singular vectors found in both numerical experiments, namely for the MOC growth and the upper-ocean temperature growth. We varied the perturbation amplitude and found the mechanism to be robust and the growth amplitude to be within 20% of the linear approximation for density anomalies of 0.05 kg m −3 ( Figure 10).
This study and ZHMT11 are the first numerical experiments calculating the singular vectors for the largescale overturning circulation and upper-ocean temperatures using a GCM and its combined tangent linear and adjoint. The use of a relatively coarse resolution and an idealized geometry GCM permitted us to explore the non-normal effects extensively at relatively low computational cost and avoid several numerical artifacts in the computation of the singular vectors (Zanna et al., 2010). A sensitivity analysis to the resolution of the model was performed. Using an horizontal resolution of 1 • and 2 • , the singular vectors and their growth were found to be qualitatively similar to the one described in this study (ZHMT11). ZHMT11 found that the optimal patterns leading to the largest amplification of MOC anomalies are at high latitudes at depth, suggesting that non-resolved processes such as overflows, convection or eddies could lead to large errors in evaluating the predictability time-scales of the North Atlantic Ocean. While the effect of bottom topography was found to damp interannual fluctuations of the MOC in some models (Huck and Vallis, 2001), Chhak et al. (2006) showed that bathymetry can provide an additional source of non-normality and further amplify the growth of MOC and upper-ocean temperatures anomalies. The influence of atmospheric and sea-ice feedbacks, the role of eddies and overflows in addition to a realistic geometry is currently under investigation. This will enable a detailed comparison with observations and determine how these processes affect the variability and predictability. In addition, the impact of stochastic forcing on the non-normal ocean dynamics is being explored to estimate the stochastic optimals. The stochastic optimals are the spatial structure of the forcing leading to the maximum variance of the MOC and upperocean temperatures.