## Abstract

The role of ocean dynamics in optimally exciting interannual variability of tropical sea surface temperature (SST) anomalies is investigated using an idealized-geometry ocean general circulation model. Initial temperature and salinity perturbations leading to an optimal growth of tropical SST anomalies, typically arising from the nonnormal dynamics, are evaluated. The structure of the optimal perturbations is characterized by relatively strong deep salinity anomalies near the western boundary generating a transient amplification of equatorial SST anomalies in less than four years.

The associated growth mechanism is linked to the excitation of coastal and equatorial Kelvin waves near the western boundary following a rapid geostrophic adjustment owing to the optimal initial temperature and salinity perturbations. The results suggest that the nonnormality of the ocean dynamics may efficiently create large tropical SST variability on interannual time scales in the Atlantic without the participation of air–sea processes or the meridional overturning circulation. An optimal deep initial salinity perturbation of 0.1 ppt located near the western boundary can result in a tropical SST anomaly of approximately 0.45°C after nearly four years, assuming the dynamics are linear. Possible mechanisms for exciting such deep perturbations are discussed. While this study is motivated by tropical Atlantic SST variability, its relevance to other basins is not excluded.

The optimal initial conditions leading to the tropical SST anomalies’ growth are obtained by solving a generalized eigenvalue problem. The evaluation of the optimals is achieved by using the Massachusetts Institute of Technology general circulation model (MITgcm) tangent linear and adjoint models as well the the Arnoldi Package (ARPACK) software for solving large-scale eigenvalue problems.

## 1. Introduction

The variability of sea surface temperature (SST) in the tropical Indian, Pacific, and Atlantic Oceans is commonly believed to be dominated by coupled ocean–atmosphere feedbacks. The coupled El Niño–Southern Oscillation (ENSO) clearly dominates the tropical Pacific variability on interannual time scales. The Indian Ocean variability includes a coupled ocean–atmosphere dipole mode, apparently independent from ENSO (Saji et al. 1999). In the tropical Atlantic, two dominant statistical modes are often used to describe the variability: the Atlantic dipole and the Atlantic Niño mode (e.g., Chang et al. 2000; Xie and Carton 2004; Chang et al. 1997). A positive thermodynamic air–sea feedback through wind-induced surface evaporation is believed to be responsible for the presence of the cross-equatorial Atlantic SST dipole (Carton et al. 1996; Chang et al. 1997). Wind forcing is the main source of SST variability associated with the Atlantic Niño mode, the latter defined as the anomalous SST and heat content in the eastern equatorial basin. In addition to the local forcing, the tropical SST variability in the different basins appears to be strongly influenced by remote forcing, particularly teleconnection mechanisms involving ENSO (e.g., Enfield and Mayer 1997; Penland and Matrosova 1998) and the North Atlantic Oscillation (Chang et al. 2001) for the tropical Atlantic. Unlike in the Pacific, ocean wave dynamics in the tropical Atlantic are usually regarded of secondary importance.

Even though atmospheric excitation of the observed tropical SST variability may be significant, coupled ocean–atmosphere dynamics are not enough to fully explain this variability (e.g., Zebiak 1993). The precise role of ocean-only dynamics in this variability still remains unknown. A few studies invoked ocean dynamics to explain the tropical Atlantic interannual variability. For example, Jochum et al. (2004) suggested that ocean mesoscale eddies due to instabilities of equatorial currents could be a significant source of interannual variability in SST, especially during the Northern Hemisphere spring when the ITCZ is highly sensitive to SST anomalies. It has also been suggested that the Atlantic meridional overturning circulation might affect interannual to decadal tropical ocean variability (e.g., Marshall et al. 2001). The question of whether the ocean passively responds to atmospheric forcing or actively participates in tropical Atlantic variability is yet unanswered. Similarly, understanding the role of the ocean dynamics, even if not dominant, in the tropical Indian and Pacific Oceans, may improve our ability to predict their variability on interannual time scales. Because it is difficult to dissociate the effects of the atmosphere and the ocean on the tropical ocean variability from observations only, we turn to ocean models to address this issue.

Over the past 50 years, the tropical oceans seem mostly governed by linear dynamics (e.g., Penland 1996; Penland and Matrosova 1998). Based on this evidence, can some of the observed tropical SST variability be explained as small-amplitude damped linear ocean dynamics excited by atmospheric or other stochastic forcing? In such stable or damped linear systems, small perturbations can be amplified before their eventual decay, owing to the interaction of multiple nonorthogonal eigenmodes if the linearized dynamical operator is nonnormal (Farrell 1988; Trefethen et al. 1993; Farrell and Ioannou 1996). The initial conditions leading to such transient growth are referred to as optimal initial conditions and are obtained by solving an eigenvalue problem based on the linearized model equations. Several studies applied this concept to climate phenomena such as atmospheric cyclogenesis (e.g., Farrell 1988, 1989), atmospheric predictability (e.g., Buizza 1995; Buizza and Palmer 1995), the oceanic wind-driven circulation (e.g., Moore 1999), the thermohaline circulation (e.g., Lohmann and Schneider 1999; Tziperman and Ioannou 2002; Zanna and Tziperman 2005), and ENSO variability and predictability (e.g., Penland and Sardeshmukh 1995; Moore and Kleeman 1997).

The objective of the present study is to investigate the role of ocean-only dynamics in exciting tropical SST anomalies. To achieve this goal, we are looking for the optimal initial perturbations of temperature and salinity leading to the maximum growth of the sum of squares of the upper-ocean temperature anomalies in the tropics on interannual time scales. More precisely, the current study aims to explore the role of nonnormal ocean dynamics in generating the growth of such perturbations for a better understanding of their role in climate variability and to improve our ability to predict such variability. This study is motivated by the tropical Atlantic SST variability; however, the idealized geometry, model configuration, and results allow for the proposed mechanism to be potentially relevant to other basins as well.

To the best of our knowledge, this work is the first to calculate optimal initial conditions exciting equatorial SST anomalies for an ocean general circulation model (OGCM) that extends beyond the tropics and considers time scales beyond a few months. Note that the goal of the current work differs from the related exploration of optimal SST anomalies in the context of ENSO, even if the tools are similar. Moore et al. (2003), for example, used a hierarchy of hybrid coupled models of the tropics and studied the optimal temperature perturbations (assuming vanishing initial salinity anomalies) leading to an ENSO state on a time scale of 6 months. In their study, the authors found that the optimal temperature perturbations are determined by coupled and atmospheric nonnormal dynamics.

In the present study, we find that perturbations of deep salinity near the western boundary can optimally excite equatorial SST anomalies on a time scale of nearly 4 yr. The optimal deep salinity perturbations project onto the first two vertical baroclinic modes near the western boundary and a geostrophic adjustment of sea surface height (SSH) and velocity anomalies occurs within one week. This results in downwelling at latitudes 40°N/S and upwelling at 20°N/S in the upper levels of the basin near the western boundary. Subsequently, these anomalies excite coastal Kelvin waves traveling equatorward along the western boundary. When the anomalies reach the equator, they excite equatorial Kelvin waves, producing a transient amplification of the equatorial SST anomalies. Although wind anomalies may excite tropical Atlantic SST anomalies, an analysis of the related atmospheric dynamics is beyond the scope of this ocean-only study. The amplification mechanism is shown to rely explicitly on the nonnormality of the ocean dynamics. Overall, we find that ocean dynamics, particularly ocean waves, can play an interesting and dominant role in exciting tropical ocean variability on interannual time scales without the participation of air–sea processes or the Atlantic meridional overturning circulation.

The paper is structured as follows. In section 2, we briefly describe the Massachusetts Institute of Technology general circulation model (MITgcm) used in this study and the steady state reached by the ocean model. In section 3, we explain the methodology for computing the optimal initial conditions with the MITgcm. A description of the optimal initial temperature and salinity perturbations found to maximize the tropical SST anomalies and the associated growth mechanism, including the role of the nonnormal ocean dynamics, are provided in section 4. In the same section, we discuss the sensitivity of the mechanism to the convective parameterization, the linearity assumption, and the model resolution. We then conclude in section 5. Additional details regarding the model configuration and the methodology used to solve for the optimals are provided in the appendixes.

## 2. Idealized base steady state of the MITgcm

### a. Model configuration

To investigate the variability of tropical Atlantic SST anomalies, the MITgcm, which is described in Marshall et al. (1997a,b), is used. The model solves the primitive equations on a sphere under the Boussinesq and hydrostatic approximations for an incompressible fluid.

The configuration is an idealized Atlantic-like double-hemispheric rectangular basin defined between 65°S and 70°N, 70° and 10°W with a horizontal resolution of 5° and solid walls placed at all boundaries. The ocean depth is uniformly set to 4400 m everywhere across the basin. There are seven vertical levels with thicknesses between 40 and 1500 m. The coarse resolution is necessary to reduce the computational cost for the calculation of the optimal modes. Even though certain properties of the mechanisms may not be properly represented owing to the coarse resolution of the model, the results appear to not be overly sensitive to the model resolution (discussed in section 4d). To avoid artificial deep equatorial meridional overturning cells due to the low vertical resolution (Weaver and Sarachik 1990), high values of vertical viscosity and diffusion coefficients were used (Table 1).

The ocean surface is forced by a time-mean wind stress, shown in Fig. 1a and given by (A1), with westerly winds poleward of 20° and easterlies between 20°S and 20°N, with a maximum amplitude of 0.12 N m^{−2} located off the equator to avoid strong Ekman upwelling there (see appendix A for more details). A reasonable steady state was reached by the model using restoring boundary conditions for both temperature and salinity after 8000 years of integration, with a symmetric circulation about the equator. The temperature field was restored to a time-independent latitudinal profile shown in Fig. 1b with an associated time scale of 2 months. At this point, we switched to mixed boundary conditions in which the temperature is still restored to the same field but the salinity restoring term is replaced by a fixed flux equivalent to a freshwater flux defined in appendix A [Eq. (A5)] and shown in Fig. 1c. After another 2000 years of integration and symmetric surface forcing, the model reached a new equilibrium state with asymmetric meridional overturning cells and temperature and salinity fields, consistent with the symmetric overturning being unstable under mixed boundary conditions (Marotzke 1990). Additionally, since the goal of the study is to analyze tropical Atlantic-like variability generated by the ocean-only on time scales ranging from years to decades and discarding the influence of the atmosphere, seasonality was not included in the wind forcing or in the boundary conditions for temperature and salinity. We are aware that the idealized geometry and forcing lead to some unrealistic aspects of the base steady state, such as the lack of a salinity minimum at the equator; however, this allows for a more thorough understanding of the amplification mechanisms explored here, without the full complexity of the ocean, which can mask some of the related ocean feedbacks. We should thus keep in mind that a different steady state may affect the results and that the relatively strong seasonal cycle in the tropical Atlantic may significantly modulate the amplitude of the growth found in this study. Additional details regarding the model configuration, numerical schemes, and model parameters specific to the present study can be found in appendix A.

### b. Steady state

The steady state reached by the model under mixed boundary conditions is now described. The SST, shown in Fig. 2a, varies between 25°C at the equator and −1°C at high latitudes. In Fig. 2b, the sea surface salinity (SSS) decreases from 36 psu at the equator to 34.5 psu at high latitudes. The deep temperature and salinity are fairly uniform with average values of about 3°C and 35 psu. The sea surface height (SSH) exhibits a difference of 80 cm across the western boundary current (not shown).

Because the optimal perturbations (described in section 4) turn out to be concentrated in the western part of the basin, the steady-state temperature and salinity as functions of latitude and depth near the western boundary are shown in Fig. 3. The flow near the surface, shown in Fig. 4a, is characterized by one gyre in each hemisphere with a similar, but reversed and weaker, cell at depth (Fig. 4b), reminiscent of the Stommel and Arons (1960a,b) abyssal circulation. We note that downwelling occurs in the northern and southern part of the basin near the eastern boundary with vertical velocity on the order of 5 to 10 × 10^{−5} m s^{−1}. Moreover, relatively strong upwelling is present near the western boundary with a typical velocity of 10^{−5} m s^{−1} (Fig. 4).

Figure 5 shows the meridional overturning streamfunction with two deep cells and a maximum transport of 25 Sv (1 Sv ≡ 10^{6} m^{3} s^{−1}) in the Northern Hemisphere. The barotropic streamfunction (not shown) clearly shows two gyres in each hemisphere. The zonally and vertically integrated heat transport has a maximum of 0.8 PW (1 PW = 10^{15} W) around 30°N (not shown). Deep convection is triggered poleward of 50° near the boundaries.

## 3. Solving for optimal initial conditions resulting in maximum equatorial SST anomalies

Consider now the optimal temperature and salinity perturbations about the steady-state solution, leading to the maximum growth of tropical SST anomalies after a given time *τ*. We start with a brief explanation of the general methodology used to solve such a problem using an OGCM (section 3a) followed by a detailed description specific to the objectives of this study (section 3b).

### a. Finding the optimal perturbations using an OGCM

Let **P**(*t*) be the state vector specifying the time-dependent prognostic variables at all grid points: that is, zonal and meridional velocities *u* and *υ*, respectively; potential temperature *T*; salinity *S*; and sea surface height *η*. The state vector can be written at any given time *t* as the sum of the steady-state solution **P** and a small perturbation **P**′ to the steady-state solution such that **P**(*t*) = **P** + **P**′(*t*). The linearized model for small perturbations (or tangent linear model) is therefore given by

where 𝗔 is the linearized model operator. Because we linearize about a steady state, the linearized operator 𝗔 is time independent and the solution to Eq. (1) is given by

where 𝗕(*t*) ≡ *e*^{𝗔t} is the propagator matrix and **P**′_{0} ≡ **P**′(*t* = 0) is the perturbation at *t* = 0. The system is linearly stable such that all eigenvalues of 𝗔 have negative real parts (Fig. 6). The adjoint model for the small perturbation **P**′ can be expressed as

where the superscript † denotes the Hermitian transpose or adjoint. The minus sign before the time derivative of **P**′^{†} implies a backward integration in time of the adjoint equations. In the following sections, we will drop the primes for convenience. The MITgcm tangent linear and adjoint models are generated with the automatic differentiation tool: TAF (Giering and Kaminski 1998; Heimbach et al. 2005). The adjoint of the MITgcm was used in several studies to investigate for example the sensitivity of the meridional heat and volume transports to surface forcing and model parameters (Marotzke et al. 1999; Bugnion et al. 2006).

Given the state vector **P**(*t*), its norm is defined as **P**^{T}𝗫**P**, where the norm kernel 𝗫 is chosen to reflect the physical quantity to be maximized. The optimal initial conditions **P**_{0} are then defined as the fastest growing perturbations over a given time *τ* and correspond to the eigenmode of the matrix 𝗕(*τ*)^{T}𝗫𝗕(*τ*) with the largest eigenvalue (Farrell and Ioannou 1996). To compute the optimal initial conditions (or optimal perturbations) for a given time *τ*, we proceed as follows: integrating the tangent linear model [Eq. (1)] with initial conditions *δ***P**(*t* = 0) over time *τ*, then multiplying the output by the norm kernel matrix 𝗫, and, finally, using the obtained result as initial conditions to integrate the adjoint model [Eq. (3)] over the interval of time *τ*. This is equivalent to the matrix–vector operation 𝗕(*τ*)^{T}𝗫𝗕(*τ*)*δ***P**. Eigenvectors and eigenvalues of 𝗕(*τ*)^{T}𝗫𝗕(*τ*) are computed using the Lanczos algorithm (Golub and Van Loan 1989) and the routines for symmetric eigenvalue problem of the Arnoldi Package (ARPACK) software (Lehoucq et al. 1998), which requires only the output vector 𝗕^{T}𝗫𝗕*δ***P**. Several iterations (or cycles) of running the combined tangent linear and adjoint models described above are required for a satisfactory convergence of the solution. This machinery was successfully implemented using the automatically differentiated MITgcm tangent linear model and its adjoint combined with ARPACK so as to find the optimal perturbations. Additional details regarding the method are given in appendix B. This method is similar to the one developed by Moore et al. (2004) with the Regional Ocean Model System (ROMS). Another method to define and calculate the optimal perturbations in an OGCM was used by Sevellec et al. (2007) based on sensitivity calculation involving only the adjoint of the model without solving an eigenvalue problem.

### b. Optimal initial conditions of temperature and salinity exciting equatorial SST variability

We are looking for the optimal initial temperature and salinity anomalies **Q** exciting equatorial SSTs at time *τ*, (*τ*). Mathematically, we wish to maximize the sum of the squares of the equatorial SST anomalies between 15°S and 15°N at time *τ*,

under some initial constraint. The norm used for SST anomalies can pick out anomalies of equal or opposite sign around the equator and thus allows the model to choose the optimal SST spatial structure that can develop, dipole or other, without us dictating it in advance. Additionally, it can be viewed as the error variance of the model forecast for surface temperatures near the equator (Lorenz 1982).

Our interest lies in the ocean’s response on interannual and longer time scales, and we expect the initial velocity perturbations to adjust rapidly, leading to a flow in geostrophic balance and therefore entirely determined by the density gradients. It is thus sufficient to consider only initial perturbations of temperature and salinity and assume that the influence of the initial velocity perturbations quickly vanishes. The sea surface height, being evaluated from the horizontal velocity field, was assumed to vanish as well. Note that geostrophy holds for any distance greater than 100–200 km away from the equator (e.g., Lagerloef et al. 1999). Although some of the optimal perturbations (described in details later in section 4a) are located near the equator, the most significant ones are situated sufficiently away from the equator where the flow is nearly in geostrophic balance. In Eq. (4), the vector **Q** is therefore composed of the temperature and salinity anomalies at all grid points. The norm kernel 𝗫 is a diagonal matrix where the entries corresponding to the equatorial SST anomalies are equal to the volume of the corresponding grid cell while all other entries are zero.

For our initial constraint, we require that the optimal initial conditions of temperature and salinity at all model grid points and levels, **Q**(*t* = 0) = **Q**_{0}, have a unit norm such that ‖**Q**_{0}‖_{𝗘}^{2} = **Q**_{0}^{T}𝗘**Q**_{0} = 1. The optimal perturbations are then found from the generalized eigenvalue problem (Farrell 1988; Zanna and Tziperman 2005),

where *γ* are the generalized eigenvalues. The norm kernel 𝗘 reflects the domain-integrated sum of squares of the temperature and salinity anomalies and is given by

where *H* is the uniform ocean depth; *λ* longitude, *ϕ* latitude, *z* depth, and *α* and *β* are the thermal expansion and saline contraction coefficients. The coefficients *α* and *β* are calculated at each grid point and depend on the steady-state temperature and salinity values. Although there are a number of possible initial constraints that could be used, the norm chosen here is positive definite (and nonsingular) and ensures that the contributions of salinity and temperature to the density, and the volume of the grid boxes, are taken into account.

## 4. Optimal growth of tropical SSTs

### a. Optimal perturbations leading to the tropical SST growth

Using the implementation described in the previous section for different values of *τ* varying from 1 to 60 yr, the maximum amplification of the tropical SST anomalies, measured by and defined in Eq. (4), occurs for *τ* = 3.8 yr (Fig. 7). The optimal initial conditions **Q**_{0} of temperature and salinity leading to this amplification are dominated by deep salinity anomalies near the western boundary at depths greater than 2 km (Fig. 9). As a result of our problem formulation, the initial sea surface height and horizontal velocity anomalies vanish. Figure 8 shows the zonally averaged contributions of temperature (top panel) and salinity (bottom panel) to the density field (i.e., −*αT* and *βS*). The field of initial density perturbation is dominated by the salinity anomaly contribution with a contribution from temperature anomalies being two orders of magnitude smaller. The amplitude of the initial perturbations is found to be large at depth and is nearly negligible in the near-surface levels. The quantity to be maximized is defined for the surface layers only. It is therefore not surprising that even a relatively modest increase in surface temperatures is expressed as a large amplification factor (final over initial amplitudes) for the surface anomalies (not shown).

To examine the respective contributions of initial temperature and salinity anomalies to the growth of the cost function and thus to the growth of equatorial SST anomalies, two model experiments are performed. First, initializing the linearized model [Eq. (1)] with the optimal initial salinity perturbations only (setting the initial temperature perturbations to zero) results in a tropical SST amplification with roughly the same time scale of 3.8 yr and a cost function attaining 96% of the value found when using both temperature and salinity contributions from the optimals at the time of maximum amplification (as in Fig. 11). In a second experiment, the model is initialized with the optimal initial conditions of temperature only and a maximum growth of the cost function is produced with a time scale of nearly 6 yr. In this second experiment, the cost function reaches less than 1% of its total value found when using both initial temperature and salinity optimals. The contribution of the initial salinity anomalies clearly dominates the excitation of the nonnormal eigenmodes leading to the maximum amplification of the tropical SST anomalies.

To precisely determine the regions in which the optimals of salinity are able to most efficiently excite the growth of , we initialize the linearized model with the initial salinity perturbations limited to each individual vertical layer. The largest amplification can be explained by the salinity anomalies at 2300-m depth, which are shown in Fig. 9. Exciting the model with only the optimal salinity perturbations at this level can reproduce the correct growth time scale and explain over 40% of the cost function value. Adding the contribution of the salinity at a depth of 3650 m (level 7) to the contribution of the salinity at 2300-m depth can explain over 60% of the cost function. The mean flow at these depths (Fig. 4) corresponds to the deep Stommel–Arons-like circulation described earlier.

Our results differ from those of Moore et al. (2003), who found the optimal perturbations to be dominated by coupled processes and therefore the optimal anomalies to be surface trapped. Note also that our study considers interannual time scales, which are an order of magnitude longer than considered in the work by Moore et al.

In summary, our results indicate that the optimal scenario to excite equatorial SST anomalies [defined by in Eq. (4)] on a time scale of nearly 4 yr is by perturbing the ocean with deep salinity anomalies near the western boundary. A crude estimate of the efficiency of nonnormal excitation of tropical SST by deep salinity anomalies in our model was obtained, and an initial salinity perturbation of 0.1 ppt near the deep western boundary results in a temperature anomaly of 0.45°C in the tropical Atlantic after 3.8 yr, assuming that the dynamics are linear.

### b. Transient amplification: A wave mechanism

A schematic of the transient growth mechanism is shown in Fig. 10. The time series of the cost function (4), reflecting the sum of squares of tropical upper-ocean temperature anomalies weighted by their gridbox volume, when the model is initialized with the optimal initial conditions is shown in Fig. 11. This curve exhibits two peaks: the first local maximum occurs after 8–9 months and is followed by a second (and larger) peak after 3.8 yr, when the maximum growth of tropical SST anomalies is reached. It is interesting that the dynamics leads to two peaks of the squared SST anomalies, whereas we have attempted to maximize it only at *τ* = 3.8 yr. Although it is not entirely clear why this occurs, it seems that the nonnormal dynamics excite several oscillatory eigenmodes, leading to the multiple peak behavior. The maxima are due to equatorial Kelvin waves excited near the western boundary, and the corresponding mechanism is now described in detail. It is important to remember that we are solving a linear problem and therefore a multiplication by any (positive or negative) factor will also be a solution to the eigenproblem, thus only the relative growth of perturbations is relevant to the discussion. In the following sections, we arbitrarily pick one sign so as to simplify the discussion of the growth mechanism.

The initial salinity anomalies near the western boundary are concentrated at latitudes of 30° and 50°N/S (Fig. 9). To understand the initial evolution of these anomalies over the course of a few weeks, we have decomposed the vertical structure of the anomalies near the western boundary between latitudes 10° and 60°N/S into normal vertical modes (Pedlosky 1979). The barotropic field is found to be small and the anomalies are initially composed of the first two vertical baroclinic modes near the western boundary. These two first baroclinic modes explain about 80% of the vertical structure of the initial density anomaly and dominate the evolution of the anomalies during the first few days.

A geostrophic adjustment of SSH (Fig. 12a) and velocity anomalies to the initial density occurs within about one week, resulting in downwelling at latitudes 40°N/S and upwelling at latitudes 20°N/S in the upper levels near the western boundary, as shown in Fig. 12b. Subsequently, these anomalies excite coastal Kelvin waves traveling equatorward along the western boundary.

The mechanism is now described more specifically. After one month, the cold anomalies associated with upwelling initially at latitudes 20°N/S and propagating equatorward as coastal Kelvin waves (mostly as first baroclinic modes) begin to arrive in the equatorial region, exciting equatorial Kelvin waves. These waves lead to the first maximum seen in the cost function of tropical upper-ocean temperature anomalies after 8–9 months (Figs. 11 and 13). Similarly, the warm anomalies associated with the downwelling initially at 40°N/S propagate equatorward as warm (downwelling) coastal Kelvin waves (first and second baroclinic modes). Starting at about 10 months, the anomalies continuously excite warm equatorial Kelvin waves leading to a growth of the equatorial SST (Fig. 14) until the maximum amplification time around 3.8 yr. The propagation of the cold (upwelling) and warm (downwelling) anomalies as equatorial Kelvin waves can be seen in the Hovmöller diagrams in Fig. 15. Note that the warm SST peak develops over a longer period of time for several reasons. First, it involves second baroclinic coastal and then equatorial waves that are slower than first baroclinic modes. Second, the downwelling waves slowly arrive at the equator and therefore force a longer and more gradual pulse of warming SST there compared to the upwelling waves arriving earlier in time. Finally, note that the time scale of 3.8 yr does not reflect the Kelvin wave travel time but the evolution of the equatorial SST that is affected by the many Kelvin waves gradually arriving and propagating along the equator.

The most significant terms in the temperature and salinity equations during the 3.8 yr of the growth in the regions governed by wave dynamics, namely the western boundary and near the equator, are the vertical advection of the mean temperature and salinity by the vertical velocity anomaly, *w*′∂*T*/∂*z* and *w*′∂*S*/∂*z*. These terms reflect the changes in thermocline depth due to the propagating waves.

Even though only the deep initial salinity anomalies were necessary to excite the Kelvin waves leading to the growth of the tropical SST anomalies, after several months both temperature and salinity participate equally in the growth mechanism and are necessary to achieve the growth time scale and amplitude shown in Fig. 11. Starting the model with only the temperature field obtained after 1 year of model integration or with only the salinity field after 1 year (with other initial conditions set to zero) shows that both are needed to achieve the optimal growth at 3.8 yr.

At the time of maximum amplification, the density signal is dominated by the contribution of the temperature with an almost negligible effect of the salinity, in contrast to the initial conditions (Fig. 16). We should point out that the maximum values of temperature near the equator are not at the surface but just beneath it because of the damping of temperature by the restoring boundary conditions, as seen in Fig. 16a. In contrast to the rapid appearance of the temperature anomalies and their relatively fast growth, the deep salinity anomalies, which have initiated the amplification, barely decay over the growth time scale of almost 4 yr. On the other hand, the surface salinity, which is initially two orders of magnitude smaller than the deep signal, grows by a factor of a few hundred.

The growth of the SST anomalies is followed by a gradual decay (Fig. 11) due to the asymptotic stability of the model. The warm Kelvin waves responsible for the maximum amplification of the tropical SST anomalies are reflected as westward-propagating Rossby waves at the eastern boundary and largely dissipate because of friction. Only a small fraction of the original warming signal reaches the western boundary. Simultaneously, convective mixing of the perturbation field poleward of 55° in each hemisphere triggers additional cold anomalies that propagate equatorward as Kelvin waves along the western boundary. As the warm westward-propagating Rossby waves associated with the maximum amplification reach the western boundary, they interact with these cold coastal Kelvin waves, and a large degree of cancellation occurs. This process eventually damps the SST anomalies.

### c. Role of nonnormal effects

We performed several tests to verify that the nonnormality of the ocean dynamics indeed contributes to the growth of the tropical SST anomalies found in the present work. We should point out first that, in a stable linear normal system, the energy (if defined as a positive definite quantity) of the initial perturbations decays, whereas 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. We clearly see a strong amplification of the surface temperature anomalies (Fig. 11), but our use of a norm kernel 𝗫 that does not span the entire space may, in principle, lead to a growth that is not due to nonnormality of the ocean dynamics.

Let us define a quantity reflecting the domain-integrated energy-like anomaly of the system such that

where *N* is a constant Brunt–Väisälä frequency and *T*, *S*, *u*, *υ*, *α*, and *β* are functions of space (*λ*, *ϕ*, and *z*). This norm is positive definite, spans the entire space, and no compensation of temperature and salinity can occur. Figure 17 (solid curve) shows a time series of basin-integrated energy-like perturbation when the linearized model is initialized with the optimal perturbations found in the present study (using a final norm reflecting the squared tropical SST anomalies). The growth of is evident from this figure. Moreover, at the maximum amplification time *t* = 3.8 yr, temperature and salinity anomalies are amplified at nearly all levels and their amplification factors for each layer (from the surface to the deepest level) using an L2 norm are as follows:

and

This thus rules out the possibility that the SST anomalies grows simply due an advection of the initial deep perturbations into the surface tropical ocean with no actual nonnormal amplification involved.

To further explore if the nonnormality of the dynamics is the primary cause of the amplification of equatorial SST anomalies, numerical simulations were performed using a norm kernel maximizing the squared temperature and salinity anomalies rather than only the tropical surface temperature anomalies. Therefore, in this numerical experiment, the final norm kernel 𝗫 is equal to 𝗘 defined by Eq. (6) and the initial constraint is unchanged such that the initial and final norm kernels are equal to 𝗘. The norm 𝗘 is not singular, and its amplification can only occur owing to nonnormal dynamics. A maximum growth was found in this case for *τ* ≈ 5 yr, and the growth of the energy-like perturbation is shown in Fig. 17 (dashed curve). The optimal perturbations under the 𝗘 norm (not shown) closely resemble the optimal perturbations found when using the norm kernel reflecting the sum of squares of tropical SST anomalies. While the initial anomalies are not confined to the deep ocean and are amplified both at high latitudes and in the equatorial region, an excitation of coastal and equatorial Kelvin waves leading to a growth of the sum of squares of tropical SST anomalies is obtained. The results under the 𝗘 norm kernel are similar to the ones found by maximizing the tropical SST anomalies, and such similarity points out again the importance of the nonnormality in the amplification seen in the rest of this paper.

Finally, one can project the optimal perturbations onto the eigenmodes of the propagator 𝗕. Some of the eigenvalues (shown in Fig. 6) and corresponding eigenvectors of the tangent linear model were computed using ARPACK. The optimal initial perturbations can be described as a linear superposition of the eigenmodes of 𝗕. As the optimal perturbations evolve in time and some of the eigenmodes decay, the spatial pattern of the response at the maximum amplification time emerges. We find that the modal composition of the initial excitation and of the response at maximum amplification time are very different, as is typical of nonnormal systems. By projecting the optimal perturbations onto the eigenmodes obtained using ARPACK, about 25 eigenvectors appear to significantly contribute to the initial perturbations and their evolution (not shown). Their degree of nonnormality as defined by Farrell and Ioannou [1999, Eq. (38)] is found to be relatively high (between 10 and 99). These are again clear indications that the nonnormality does participate in the tropical SST amplification mechanism found in this work.

The source of nonnormality in this study is shown to come mostly from advection of tracers (temperature and salinity) similarly to results obtained with models of intermediate complexity (Zanna and Tziperman 2005, 2008). In addition, surface boundary conditions appear to contribute to the nonnormality of the system. We find that the response of the system is larger when switching from restoring to mixed boundary conditions in both the MITgcm and models of intermediate complexity. Moreover, box models showed that interaction with an evolving atmosphere amplifies to nonnormal response of the ocean (Zanna 2009) and similarly for GCM studies related to ENSO dynamics (Moore et al. 2003).

### d. Sensitivity of the results

#### 1) Convective parameterization

The parameterization of convection in ocean models can lead to numerical artifacts when calculating optimal initial conditions (Zanna and Tziperman 2005). To validate our results and ensure that convection does not play a crucial role in the growth mechanism, especially in the propagation of the initial deep salinity anomalies up to the surface, we initialized the tangent linear model with the optimals with the convective parameterization switched off. The growth mechanism for the SST anomalies and its associated time scale are found unchanged, whereas some properties of the decay mechanism are altered when convection is turned off.

#### 2) Linear versus nonlinear dynamics

The solutions of the eigenvalue problem defined by Eq. (5) are valid under the assumption that the perturbations are small enough such that the nonlinear terms do not affect the physical mechanism. Furthermore, the linear response is independent of the sign of the initial perturbations, which is not necessarily the case in the fully nonlinear scenario. Bearing this in mind, we examine the response of the fully nonlinear model to the optimals obtained above. We find that introducing the optimal perturbations characterized with deep salinity anomalies with amplitude ranging from 0.05 to 0.14 ppt into the nonlinear model leads to a similar excitation of Kelvin modes as discussed in section 4b. These perturbation amplitudes were chosen to compare with estimates of variability in the deep ocean (Forget and Wunsch 2007). For these amplitudes, the error introduced by using the linear instead of fully nonlinear model are *O*(8%) only. The growth time scale changes by up to several months depending on the strength of the initial perturbations introduced. In general, we find the decay of the perturbations to be somewhat faster due to nonlinear processes.

#### 3) Adjoint sensitivity versus optimals

The sensitivity to model parameters or initial conditions can be calculated by using the adjoint model only (Marotzke et al. 1999; Bugnion et al. 2006; Sevellec et al. 2007). The adjoint equations, obtained by differentiating the original model, allow the calculation of the linear sensitivity (about the full nonlinear model state) of the cost function to initial conditions in one single calculation. The output of the adjoint model calculation provides the sensitivity of the same squared equatorial upper-ocean temperature measure, which was used in the calculation of the optimals, to the initial conditions of temperature and salinity over the entire ocean domain. The sensitivity is simply expressed as the gradient of the cost function at final time *τ* with respect to the initial conditions. It is useful to compare our results to the sensitivity of the squared equatorial upper-ocean temperature anomalies to initial conditions of temperature and salinity as calculated by the adjoint only, without solving the full eigenvalue problem. Similarly to the optimal initial conditions found, the adjoint shows a larger sensitivity of squared tropical SST anomalies to the salinity relative to the temperature. However, the adjoint sensitivity solution is found to be closer to the surface (with a maximum at 750 m), compared to the optimals which have a maximum at a depth of 2300 m. The anomalies near the western boundary are centered around the same latitudes as those of the optimals. Additionally, the adjoint sensitivities show large perturbations along the equator where the cost function is defined. We estimate that a salinity perturbation of 0.1 ppt at a depth of 750 m results in an equatorial SST anomaly of 0.33°C after 3.8 yr. Therefore the adjoint solution leads to about 75% of the amplification found by exciting the model with the optimal initial conditions. It seems that the adjoint sensitivity solution was able to obtain much of the nonnormal effects—but not be completely optimal as the solution of the full eigenvalue problem. The optimality of the adjoint solution in general depends on the clustering of the eigenvalues of 𝗕^{T}𝗫𝗕 for the optimal initial condition. There are still significant advantages to using the adjoint for exploring the growth mechanism of initial anomalies, which may involve nonnormal dynamics (Sevellec et al. 2007). Among the advantages are computational efficiency and the freedom to choose a quantity to be maximized that is not necessarily positive definite and may be more appropriate for exploring certain aspects of the dynamics.

#### 4) Model resolution and growth time scale

At such coarse resolution, wave dynamics is not appropriately resolved since the gridbox size is larger than the Rossby radius of deformation, and we find Kelvin waves with reduced speed. Nevertheless, several studies showed that the phase speed of coastal and equatorial Kelvin waves on a C grid is independent of model resolution (Hsieh et al. 1983; Ng and Hsieh 1994). However, it was shown that the phase speed of Kelvin waves in coarse resolution models is reduced owing to high viscosity and diffusion coefficients used in such configurations (Wajsowicz and Gill 1986). Furthermore, the dispersion properties of the effective C–D grid have been studied in detail by Adcroft (1995) and are shown to resemble those of the B grid. The latter grid, in turn, is known to exhibit reduced wave propagation speeds. In addition to their reduced speed, the waves are found to decay much faster because of the strong effects of friction and diffusion. To get a better sense of the effect of resolution as well as diffusion and viscosity on the results, we have first increased the number of vertical levels to 15 (instead of 7). This allowed a dramatic decrease in the vertical viscosity coefficient to a more reasonable value of 1 × 10^{−3} m^{2} s^{−1}. We find that the vertical and horizontal structure of the optimal initial conditions are identical to those described in the previous section, and so is the mechanism (geostrophic adjustment leading to the excitation of coastal Kelvin waves followed by equatorial Kelvin waves). The only difference is the time scale, which we find to be slightly faster, with a maximum amplification occurring after 3.4 yr. Using horizontal resolutions of 4° and 3° (instead of 5°) with 15 vertical levels, we found similar results with amplification time scales of 3.3 and 3.1 yr, respectively. Therefore, as the model resolution increases and the viscosity and diffusivity decreases, the time scale for the growth slightly decreases. Moreover, by increasing the model resolution and decreasing its eddy viscosity and diffusivity, we may expect that the number of degrees of freedom increases, and this may lead to a larger number of nonorthogonal eigenvectors and therefore a greater transient amplification. This is merely a speculation at this point, as such higher-resolution experiments could not be carried out owing to the computational cost involved.

#### 5) Norm kernel

We have calculated the optimal initial conditions using different norm kernels 𝗫 such as the domain-integrated energy-like or the sum of squared temperature and salinity over the entire domain (see section 4c for additional explanations). The optimal initial conditions were again mostly confined to the western boundary, similarly to the equatorial SST norm kernel described above. The velocities adjust to the density gradients, whereas temperature and salinity field barely change. This confirms our hypothesis that calculating optimal initial conditions for temperature and salinity alone suffices for studying the excitation of tropical SST anomalies. In addition, the mechanism does not involve only waves propagating the temperature and salinity signal but also horizontal advection and convection processes.

## 5. Conclusions

We investigated the role of the nonnormal ocean dynamics in exciting equatorial SST anomalies in an Atlantic-like ocean general circulation model. Even though the flow is stable and any perturbation eventually decays, we found that optimal initial conditions of deep salinity near the western boundary can lead to a significant growth of the sum of squares of equatorial upper-ocean temperature anomalies after about four years.

The optimal perturbations of deep salinity anomalies near the western boundary initially project onto the first two vertical baroclinic modes. This is followed by a rapid geostrophic adjustment of sea surface height and surface velocity anomalies after about one week, resulting in downwelling at latitudes 40°N/S and upwelling at latitudes 20°N/S in the upper layers of the model near the western boundary. Subsequently, these anomalies excite equatorward traveling coastal Kelvin waves near the western boundary. Finally, once the anomalies reach the equator, they excite equatorial Kelvin waves, leading to a maximal amplification of equatorial upper-ocean temperature anomalies. The maximum amplification time of about 3.8 yr is reached when upper-ocean temperature anomalies are spread in the equatorial region at most longitudes. The main contribution to the growth of temperature and salinity anomalies is the vertical advection of mean temperature and salinity gradients by the anomalous flow. We also find that the optimal excitation of the tropical upper-ocean temperature anomalies involves several nonnormal modes of the linearized model operator. The optimal initial conditions and the associated mechanism leading to maximum amplification of the tropical SST anomalies are found to be insensitive to the model resolution. However, the growth time scale slightly decreases as the model resolution increases and viscosity and diffusion decrease.

The amplification process described in this study is excited by salinity anomalies in the deep western boundary area. The relevance of this process to variability in the real ocean can be addressed by projecting observed eddy noise onto the spatial structure of the modeled optimal initial conditions. Deep ocean mesoscale eddies are observed to be strong in the western boundary areas (Dengler et al. 2004) where the optimal initial conditions are concentrated, and Seo et al. (2006) showed the important role of eddies in controlling the mean tropics and its variability, suggesting that our calculated optimals may be relevant to observed Atlantic Ocean variability.

Furthermore, the amplitude of observed deep salinity variability compares with observed tropical Atlantic SST variability in a manner fairly consistent with our results. Using hydrographic data, Forget and Wunsch (2007) estimated the variability of temperature and salinity with the seasonal cycle removed. They found the range of salinity variability near the deep western boundary to lie between 0.06 and 0.1 ppt, and the Atlantic surface equatorial temperature variability to be between 0.44° and 1°C. In the linear mechanism described here, a deep salinity anomaly of 0.1 ppt is associated with a tropical SST excitation of 0.45°C at the time of optimal growth. Therefore, a nonnegligible part of the interannual variability of the equatorial Atlantic SST could perhaps be initiated near the deep western boundary via nonnormal growth without participation of the meridional overturning circulation or air–sea fluxes.

An issue that the current study did not address quantitatively is how important is the role of the large-scale coherency of the optimals in exciting a significant SST signal in the tropics. One does not expect the stochastic salinity anomalies due to eddies to strongly project on the large-scale symmetric structure of the optimals obtained. Therefore, even if the amplitude of deep salinity anomalies is appropriate for exciting tropical SST anomalies based on our finding, the suboptimal spatial structure may significantly reduce the efficiency of the proposed mechanism. It would be interesting to quantitatively examine the degree of deterioration of the excitation of the SST due to perturbations that are limited, for example, to either the Northern or Southern Hemisphere.

The nonnormal mechanism leading to the amplification of equatorial SST relies mostly on wave dynamics and could potentially be relevant to other ocean basins as well. However, as mentioned in the introduction, the strong tropical Pacific ENSO variability may dominate the nonnormal amplification discussed in this paper, and the mechanism proposed here seems therefore more relevant to tropical Atlantic variability.

The highly idealized configuration and steady state used in this study can alter the nonnormal growth mechanism. For example, bathymetry can provide an additional source of nonnormality (e.g., Chhak et al. 2006) and amplify the excitation of anomalies. In addition, wave reflection from complex boundaries such as the ones present in the Atlantic can affect the excitation and dissipation of the SST anomalies. The model used in this study has specified mixed surface boundary conditions, which may be expected to damp surface perturbations more strongly than a fuller atmospheric model. Finally, stochastic atmospheric forcing by winds, heat, and freshwater fluxes are likely to force the surface ocean at a larger amplitude than the deep stochastic forcing due to ocean eddies. It would be useful, therefore, to study the nonnormal amplification of stochastic surface forcing in a three-dimensional model, and we hope to address this in a future study.

## Acknowledgments

We thank Brian Farrell and Jake Gebbie for very helpful discussion. We are also grateful to Chris Walker for his assistance with the numerical simulations and to the MITgcm group for sharing the source code. Comments from two anonymous reviewers significantly improved the manuscript. ET thanks the Weizmann Institute for its hospitality. This work was supported by the NSF Paleoclimate Program Grant ATM-0455470.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**(

**,**(

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**(

**,**

**,**

**,**

**,**(

**,**(

**,**(

**,**

**,**(

**,**

**,**

**,**(

**,**(

**,**(

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX A

#### Model Description

The MIT general circulation model (MITgcm) was described by Marshall et al. (1997a,b). In this appendix, we only briefly summarize the relevant variables and modifications made for the present study. In addition, all model parameters used are summarized in Table 1.

The MITgcm solves the primitive equations such that, at any time *t*, the prognostic variables for the flow are the horizontal velocity **v** = (*u*, *υ*); potential temperature *T*, salinity *S*, and the sea surface elevation *η*, which is solved using a linearized implicit scheme (Dukowicz and Smith 1994). The model solves diagnostically for the vertical velocity *w*, the density *ρ*, and the hydrostatic pressure *ϕ*_{hyd}. This study uses a relatively coarse resolution and the C–D scheme is therefore used for the computation of the Coriolis terms (Adcroft et al. 1999). The equation of state is a modified UNESCO formula derived by Jackett and Mcdougall (1995) that uses a horizontally and temporally constant pressure *p*_{0} = −*gρ*_{0}*z*. Ocean convection is parameterized by an implicit diffusion scheme.

No–normal flow and no-slip conditions are prescribed at the bottom and sides boundaries, and the vertical boundary conditions are given by *ż* = *w* = 0 at the ocean bottom *z* = −*H* and *w* = *dη*/*dt*.

The time-independent wind stress forcing applied at the ocean surface, shown in Fig. 1a, is given by

and

where 1 ≤ *x* ≤ *N _{x}*; 1 ≤

*y*≤

*N*;

_{y}*N*and

_{x}*N*are the number of grid points in the zonal and meridional direction, respectively; and the other model parameters are given in Table 1. This expression describes westerlies north (south) of 20°N (20°S) and easterlies between 20°S and 20°N. The presence of the term

_{y}reduces the magnitude of the wind stress near the equator and therefore diminished a too strong Ekman upwelling that, otherwise, occurs right at the equator.

The model was first run to steady state using restoring boundary conditions for both temperature and salinity. The fields are restored to time-independent profiles symmetric across the equator, given by

in which *ϕ* is latitude and *ϕ*_{0} is the latitude of the northernmost grid point. The restoring time scales for temperature *τ _{T}* and salinity

*τ*are 2 and 6 months, respectively. We then switched to mixed boundary conditions, with the temperature still restored to (A3) and with the salinity restoring term replaced by a fixed flux equivalent to a freshwater flux shown in Fig. 1c and defined by

_{S}where Δ*H*_{1} is the depth of the surface layer, SSS* is the restoring salinity field, and SSS is the steady-state sea surface salinity field reached under restoring boundary conditions.

Several assumptions were found necessary so as to obtain smooth sensitivity fields from the adjoint model, which in turn is needed to avoid numerical artifacts in the computation of the optimal initial conditions: flat ocean bottom, idealized basin geometry, and diffusion of tracers along grid coordinate directions rather than along isopycnals (Gent and McWilliams 1990; Redi 1982).

The model is run to steady a state with an asynchronous time stepping (see Table 1). However, the calculation of the optimals and the entire analysis of the growth mechanism is done with a synchronous time stepping.

### APPENDIX B

#### Solving for the Optimal Initial Conditions Using the Tangent Linear and Adjoint Models and ARPACK

We now describe the methodology used to solve for the optimals including technical solutions to the numerical challenges encountered. We wish to evaluate the initial conditions leading to the maximum amplification of the positive definite cost function at time *τ*,

where the vector **P** includes all prognostic variable anomalies at all grid points (potential temperature *T*, salinity *S*, zonal and meridional velocities *u* and *υ*, and sea surface height *η*), 𝗕 is the propagator matrix, and **P**_{0} is the initial state vector. We wish to maximize (*t* = *τ*) subject to the constraint that **P**_{0} has a unit norm under some possibly different norm kernel 𝗘 such that ‖**P**_{0}‖_{𝗘}^{2} = **P**_{0}^{T}𝗘**P**_{0} = 1. This optimization problem is solved using Lagrange multipliers, denoted here by *γ* (Farrell 1988), leading to the generalized eigenvalue problem,

The initial conditions **P**_{0} leading to an optimal growth of at a time *t* = *τ* are therefore the generalized eigenvectors of 𝗕(*τ*)^{T}𝗫𝗕(*τ*) (Farrell 1988; Zanna and Tziperman 2005), and the generalized eigenvalue *γ* is equal to the amplification of (*t* = *τ*) if the two kernels 𝗫 and 𝗘 are equal.

Assuming that the norm 𝗘 is not singular and therefore invertible, it is convenient to solve this problem by transforming the generalized eigenproblem (B2) into a standard eigenvalue problem by multiplying both sides of (B2) by the inverse of 𝗘, denoted by 𝗘^{−1}, such that

To efficiently solve this problem with ARPACK, it is best to transform it to a symmetric problem by defining (Moore et al. 2004)

allowing the eigenvalue problem to be written with a standard symmetric matrix on the left-hand side as

The eigenvalues *γ* are equal to the generalized eigenvalues of the original problem, and the solution for the eigenvectors **s**_{0} can be transformed back to the generalized eigenvectors **P**_{0} using **P**_{0} = 𝗘^{−1/2}**s**_{0}.

As explained above, we limit the optimal initial conditions to temperature and salinity only, setting the other prognostic variables [horizontal velocity (*u*, *υ*) and sea surface height *η*] to zero. We are interested in interannual time scales and, therefore, expect the horizontal velocity and sea surface height to quickly adjust to the initial density to achieve a geostrophic balance (except right at the equator where the Coriolis parameter is identically zero). For this reason, initial values of horizontal velocity and sea surface height are not expected to affect nonnormal growth on the longer time scales of interest to this study. Moreover, we note that the pressure and therefore the sea surface elevation are solved using an implicit solver in the MITgcm. The velocities are then reevaluated to satisfy the sea surface height equation, leading to further complications when trying to include the initial conditions of horizontal velocity and sea surface height in the optimal initial equations.

## Footnotes

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