## Abstract

The vertical structure of the variability in the equatorial Pacific in a high-resolution ocean general circulation model (OGCM) simulation for 1985–94 is investigated. Near the equator the linear vertical modes are estimated at each grid point and time step of the OGCM simulation. The characteristics of the vertical modes are found to vary more in space than in time. The contribution of baroclinic modes to surface zonal current and sea level anomalies is analyzed. The first two modes contribute with comparable amplitude but with different spatial distribution in the equatorial waveguide. The third and fourth modes exhibit peaks in variability in the east and in the westernmost part of the basin where the largest zonal gradients in the density field and in the vertical mode characteristics are found. Higher-order mode (sum of third to the eighth mode) variability is the largest near the date line close to the maximum in zonal wind stress variability. Kelvin and first-meridional Rossby components are derived for each of the first three baroclinic contributions by projection onto the associated meridional structures. They are compared to equivalent ones in multimode linear simulations done with the projection coefficients and phase speeds derived from the OGCM simulation. This suggests that in addition to the first-mode-forced equatorial Kelvin and Rossby waves earlier found in the data, forced waves of higher vertical modes should also be observable. For the first two vertical modes, the anomalies in the linear and the OGCM simulations have a similar magnitude and usually present similar propagation characteristics. Phase speed characteristics are however different in the eastern Pacific with larger values for the OGCM. The effect of zonal changes in the stratification is tested in the linear model for a stratification change located either in the eastern or in the western Pacific. This results in a significant redistribution of energy to higher modes via modal dispersion. In particular the third mode increases to a magnitude closer to the one in the OGCM simulation. Gravest modes are also affected. This suggests that modal dispersion plays an important role and should be considered when interpreting data as a combination of linear long equatorial waves.

## 1. Introduction

Simple linear models have been successful in describing the low-frequency variability of the tropical Pacific (Cane 1984; Busalacchi and Cane 1985). Long equatorial waves as well as their reflection at the meridional boundaries have often been investigated in these models to understand their role in ENSO variability (Battisti 1988; Schopf and Suarez 1988). These ocean models with only one or two baroclinic modes are capable of sustaining ENSO-like oscillations when coupled to simple atmospheric models (Zebiak and Cane 1987; Battisti 1988; Chen et al. 1995; Picaut et al. 1997). They have also demonstrated skill for forecasting ENSO (Cane et al. 1986; Barnett et al. 1988; Kleeman et al. 1995). However, there are significant discrepancies between these models and observations or simulations by ocean general circulation models (OGCMs) raising the issue of the actual role of linear long equatorial waves. Among the most dramatic discrepancies is the tendency of linear models to produce extra-equatorial Rossby waves much larger [by a factor of 4 or 5 in the Graham and White (1988) model] than those observed, which exagerates the importance of reflection on continental boundaries (Périgaud and Dewitte 1996). Note also that some of these simple coupled models failed to predict the recent El Niño events (see CAC bulletins). Furthermore there is a large sensitivity of these models to parameter selection (Zebiak and Cane 1991; Périgaud and Dewitte 1995). Thus caution is required when drawing conclusions on the real world ENSO cycle from those coupled models. For these reasons investigating the role of long equatorial waves and the processes associated with their vertical structure remains of great interest. These can be handled by OGCM simulations for which we have the information on the vertical structure. These models with sufficient vertical and horizontal resolution deal with more complete dynamics than just described by linear equatorial waves (nonlinearities, diffusion, . . .). One needs also to address the relevance of vertical modes to describe GCM dynamics at low frequencies.

The approach proposed here consists in confronting results from a GCM simulation to linear multimode simulations. Giese (1989) also used a combination of linear theory and OGCM to examine the equatorial response to forcing on timescales from days to months. His study sheds light on many processes involving equatorial free waves and their interaction with the local density structure and mean circulation. In contrast with his study, we choose to concentrate on interannual timescales and we will not investigate wave interaction with background flows. We focus on the spatial structure of the variability for the different vertical modes and how this could be influenced by spatial variability in the density structures. Boulanger et al. (1997) analyzed the same OGCM simulation in terms of long equatorial waves, assuming a single baroclinic mode representation of the sea level variability. They detected first-baroclinic-mode Kelvin and first-meridional-mode Rossby propagations but noticed significant discrepancies between simple wave model simulations and wave coefficients derived from the OGCM simulation.

In the present study, instead of assuming an arbitrary vertical structure to decompose sea surface variability in waves for specified baroclinic modes, we perform the decomposition into vertical modes of the model simulation at each grid point and time step. We first investigate the relative importance of baroclinic modes in the OGCM. In particular, the results confirm that more than one baroclinic mode is necessary to accurately hindcast sea level or current anomalies as was already guessed from linear model simulations (Blumenthal and Cane 1989). Propagation characteristics in the OGCM are then analyzed using a multimode linear simulation as a benchmark. The multimode linear model is also used to test the effect of zonally varying stratification on the equatorial Pacific variability. The question of the importance of coupling between vertical modes via modal dispersion is thus addressed.

The paper is organized as follows. Section 2 presents the OGCM and the simulation that will be used. Section 3 is devoted to the analysis of the vertical structure: The approach adopted is presented and the contribution of the different baroclinic modes to surface zonal current and pressure is examined. Parameters necessary for the linear model simulations are derived. Section 4 is devoted to the intercomparison between the OGCM simulation with the standard multimode linear model simulation. The impact of zonally varying stratification is tested. Section 5 discusses the results in light of theory and previous studies. Finally, in section 6, we summarize and provide a perspective for this work.

## 2. Description of the OGCM and the simulation

### a. The model

The simulation was performed with Océan Parallélisé (OPA), a rigid-lid hydrostatic ocean general circulation model [version 7 in Delecluse et al. (1993)] developed at the Laboratoire d’Océanographie Dynamique et de Climatologie. The basic assumptions and the set of discretized equations using tensorial formalism are described in Delecluse et al. (1993) and Marti et al. (1992). Prognostic variables are the horizontal currents, temperature, salinity, and turbulent kinetic energy. There is a turbulent closure scheme with the assumption that small-scale horizontal and vertical transports can be estimated as a function of diffusion coefficients and derivatives of the large-scale flow. The state equation is computed from Millero and Poisson (1981). The equations are discretized with a second-order finite difference scheme on an Arakawa (1972) C grid. The temporal scheme is a leapfrog scheme with an Asselin (1972) filter, except for vertical diffusion that uses an implicit scheme and the horizontal diffusion that uses a forward scheme. The horizontal viscositivity and diffusivity are small positive constants equal to 10^{3} m^{2} s^{−1}. These values were adopted following sensitivity experiments (Maes et al. 1997).

### b. The simulation

The version of the model used here is presented in Maes (1996) and Boulanger et al. (1997). The model domain includes the three ocean basins between 50°N and 50°S (the Indian and Pacific Oceans communicate both south of Australia and in the Indonesian region). The zonal mesh is determined by analytical functions with increased zonal resolution each time that a coastline crosses the equator. The zonal resolution is equal to 0.33° near coastlines increasing to 0.75° in the middle of the basins. The meridional resolution is 0.33° at the equator increasing to 1.5° at the zonal boundaries. The vertical mesh has a variable resolution with a maximum of 30 levels. The vertical resolution is 10 m near the surface, 20 m at 200 m, increasing to 500 m for the deepest level reaching 5000 m (see Table 1). Bottom topography is included based on the 5′ × 5′ ETOPO5 dataset (Marine Geology and Geophysics Division, National Geophysical Data Center, Boulder, Colorado). The vertical eddy coefficients are computed from a 1.5 turbulent closure model in which the evolution of the turbulent kinetic energy is given by a prognostic equation (Blanke and Delecluse 1993). The model is run with a time step of 5400 s and forcing conditions (wind stress and heat and freshwater fluxes) taken from the climate version of ARPEGE (Déqué et al. 1994). The ARPEGE run corresponds to the Tropical Ocean Global Atmosphere period (1985–94) and was forced by prescribed SST (Reynolds and Smith, 1994). The surface heat flux is parameterized as

where *Q*_{Arpege} is given by the atmosphere model, SST_{obs} is the observed sea surface temperature of Reynolds and Smith (1994), and SST is the temperature of the upper level of the ocean model. The negative feedback term ∂*Q*/∂*T* is set to −40 W m^{−2} K^{−1}, a typical value for the equatorial Pacific Ocean (Oberhuber 1988). The solar flux is allowed to penetrate below the top model layer with a vertical profile of absorption according to Paulson and Simpson (1977) with Jerlov (1968) optical water Type I. About one-half of the solar flux is absorbed into the first top decimeters while the remaining part is absorbed within the upper 50 m. In order to avoid a slow drift at depth due to an incorrect representation of the deep circulation and of its water mass sources, temperature and salinity are restored toward the Levitus (1982) data (monthly mean for temperature, seasonal mean for salinity) below the mixed layer poleward of 20°N–S.

The model was initialized in January 1984 with Levitus salinity and temperature and was spun up from rest by repeating 1984 once and the year 1985 five times. This simulation reproduces correctly many features of the seasonal cycle (Maes 1996). In particular, simulated climatological surface currents compare well to the climatology established by Reverdin et al. (1994) from buoy drifts and current meter records between January 1987 and April 1992. The simulated interannual variability was not extensively validated with data, but the model reproduces the El Niño and La Niña events in 1987 and 1988. The warming in 1992 is also simulated. Comparison with a few XBT transects in the tropical Pacific (Maes 1996) indicates however that the thermocline is somewhat too diffusive, entailing a small drift at subsurface. This, in particular, results in a trend in the oceanic fields, which was removed before analysis.

## 3. Vertical mode decomposition

### a. Method

A decomposition in vertical modes of the model variables is sought along the equator. The structure functions *F*_{n} are the linear vertical modes solutions of

(Fjeldstad 1933) with the boundary conditions:

where *c*^{2}_{n} is the separation constant and *N* is the Brunt–Väisälä frequency.

We could have discretized the equations on the model grid. Instead, we prefer to treat the profiles as continuous ones. This results in very little difference for the low-order modes. Salinity and temperature profiles are interpolated vertically to a 5-m resolution using 1D cubic splines. The potential density *ρ*(*x, y* = 0, *z, t*) was determined using the state equation from Millero and Poisson (1981). An approximate Brunt–Väisälä frequency *N*(*x, z, t*) is then computed. We define a stretched grid of *M* levels as in Eriksen (1985) such as

*z*′ being the depth of the level *k* on the new grid (*k* = 1, *M*). Here *M* is chosen such as to retain a 5-m resolution in strongly stratified regions. The resolution of the stretched grid is lower in weakly stratified regions below the thermocline (up to 100 m in the deep ocean). We retained *M* = 180, which is sufficient to estimate accurately the vertical modes. Vertical derivatives are estimated on the stretched grid. A minimum value of 10^{−8} s^{−2} for the square root of the Brunt–Väisälä frequency is set to get rid of unstable solutions. Equation (1) is discretized by finite forward differences and the system obtained is solved by a singular value decomposition (SVD).

The most important dynamical area for El Niño is close to the equator. We did the decomposition in vertical modes right at the equator. Results at 1°N or 1°S are very similar, but differences are larger at 5°N or 5°S in the western Pacific, to a large extent because of meridional changes in ocean depth. At 5°N, the presence of a countercurrent trough also contributes to larger differences. Results presented below hold however between 5°S and 5°N.

In this computation, *H* as a function of *x* and *N* depends also on time, thus we obtain a set of *F*_{n} as a function of *x, z,* and *t.* The separation between the vertical and the horizontal structure is based on the assumption that the vertical stratification is steady and meridionally uniform. The zonal and temporal variation in *N*^{2} will couple dynamically horizontal and vertical structures. This expresses the fact that stratification and the meridional deformation scaling change in time and longitude, entailing changes in the horizontal propagating signals. We tested the sensitivity to time and zonal dependency of *N*^{2} by performing the decomposition with a constant *N*^{2} typical of the central Pacific (based on the 1985–94 mean of *N*^{2} at 160°W: 〈*N*^{2}(*x* = 160°W)〉_{t}). We also used *N*^{2} a function of longitude only (〈*N*^{2}(*x*)〉_{t}). Results presented in appendix A show that the major sensitivity originates from the zonal dependency of *N*^{2}.

Monthly means are considered instead of the 5-day output, which has no major impact on the results. For low-frequency variability an alternative approach was used by Boulanger et al. (1997), who applied a 55-day Hanning filter onto the 5-day original fields. We checked that both methods provide similar results. Since we are interested in interannual variability of zonal current {*u*(*x, y, z, t*)} and pressure {*p*(*x, y, z, t*)}, *u* and *p* anomalies are estimated relative to the climatological seasonal cycle of the 1985–94 simulation (the seasonal averages and trends are removed). These fields are then projected on the vertical structure functions *F*_{n}(*x,* 0, *z, t*) (neglecting the meridional structure in the functions). This implies computing the following quantity:

with *q* being *u* or *p.*

Normalization is such that *F*_{n}(*x, z* = 0, *t*) = 1. At latitudes where the ocean depth is different than at the equator, the vertical integration is performed down to the shallower depth common to the equator and to that latitude. Small differences in bottom depth hardly affect the result, as was checked by making other assumptions on the deep portion of the profiles.

### b. Results

#### 1) Zonal current anomalies

In this section, we discuss the contribution of the different vertical modes to the zonal current anomalies (hereafter ZCA). Because the vertical functions for the gravest modes have no structure in the top 40 m, the projections of the current anomalies onto these modes also have no shear in the upper layer. Therefore the surface zonal currents are estimated as the currents averaged over the four uppermost levels of the model (depth 5 m to 35 m). This averaging removes part of the shear within the weakly stratified surface layer associated with incomplete mixing, which cannot be represented in the vertical mode decomposition. The rms variability of surface zonal current anomalies over the 10-yr period is presented in Fig. 1a. Large values are confined within 5°N–5°S, with maximum variability centered near the date line.

The vertical profiles of the simulated zonal current anomalies are then projected onto each baroclinic mode. As a large variability is trapped above the thermocline, it projects preferentially on high-order modes. This variability will be referred to as the upper ocean shear.

##### (i) The two gravest baroclinic modes

Figures 1b and 1c illustrate the contribution of the first two gravest modes to surface currents. Their rms variability near the equator is much weaker than the total rms (Fig. 1a). For mode 1 large values extend meridionally to 5°N–5°S, whereas variability of mode 2 is more confined near the equator. For each of the first 2 modes a meridional scale is estimated by fitting a Gaussian curve onto the rms variability. We find zonally averaged meridional scales of 850 km for mode 1 and 600 km for mode 2. Mode 1 has a peak (10 cm s^{−1}) in the western Pacific and mode 2 presents maxima of a similar magnitude both in the western and eastern Pacific. The contribution of the second baroclinic mode is more important in the eastern Pacific. In the western Pacific mode 1 has higher rms variability than mode 2. In the central Pacific the contributions of the two modes have similar magnitudes. The two gravest modes explain 56% (correlation coefficient *c* = 0.79) of the total surface current variance over NINO4eq (0°, 150°W–160°E) and 68% (*c* = 0.63) over NINO3eq (0°, 90°–150°W).

##### (ii) High-order baroclinic modes

We now concentrate on higher-order modes. Figures 2a and 2b present the rms variability for the contribution of the third and fourth modes and Fig. 2c the rms variability of the sum of modes 3–8. First, we note that the contribution of individual high-order modes is much weaker than for mode 1 or mode 2. However, when summing modes 3 to 8, rms variability reaches 18 cm s^{−1} at the equator just east of 180°, which is larger than for mode 1 or mode 2 (Fig. 2c). This sum presents two equatorial maxima located on both sides of the central Pacific location where the zonal wind variability peaks (Fig. 2d), suggesting local forcing of Kelvin and Rossby waves. The correlation coefficient between zonal wind stress anomalies and the contribution of the sum of the third to the eighth mode reaches 0.7 in the central Pacific [zonal wind stress anomalies have been shifted by 10° to the east; the correlation is significant at the 99% level—the number of degrees of freedom was estimated following Davis (1977)]. Correlation drops to 0.5 (significant at the 95% level) without the zonal shift of the wind stress. Correlation between the sum of modes 1 and 2 and the wind never exceeds 0.5. This pattern fits with the expectation that higher-order modes are more closely related locally to the wind because of their smaller vertical scales and slower horizontal wave speeds. The residual contribution of the modes higher than 8, estimated as the difference between total zonal current anomalies and the contribution of the eight first modes to the zonal current anomalies, present similar characteristics (the correlation between this contribution and zonal wind stress anomalies reaches a significant value of 0.55 near the date line).

Although high-order modes can be explained by local forcing, it is worth noticing that the third and fourth modes present variability peaks in the far eastern and western Pacific reaching 8 cm s^{−1} where there is little local forcing (see Figs. 2a and 2b). For mode 3 and to a lesser extent for mode 4, there is strong variability in the eastern Pacific from 120°W to the coast of Peru and in the western Pacific from 170°E to the Indonesian archipelago, regions where we find large zonal changes in the stratification inside the thermocline (see Fig. A1). These features are no longer present when the projections are done on modes estimated using a spatially uniform density structure (see Fig. A3). Rms variability for the modes 3 and 4 is then maximum in the west-central Pacific and reaches peak values of only 4.5 and 3.3 cm s^{−1}, respectively.

##### (iii) Spectrum

Figure 3 presents the power spectra of ZCA, of the first three mode contributions to ZCA and of the zonal wind anomalies, all quantities being averaged over NINO4eq. No peak significant at the 90% confidence level emerges. Note however that all the modes present low-frequency variability and that mode 1 has relatively more variance at higher frequencies (between 5 and 10 months) than at lower frequencies compared to the other modes. Although the mode 1 spectrum does not exhibit a specific high-frequency peak, this is consistent with linear simulations suggesting that the first baroclinic mode exhibits energetic high frequencies basin modes whose characteristics are determined by phase speed and basin geometry (Périgaud and Dewitte 1996).

##### (iv) Correlation and explained variance

The correlation between the individual modal contributions and the surface current are presented along the equator in Fig. 4a. The number of degrees of freedom involved in estimating significance of correlation ranges from 20 to 80 depending on the longitude. Zonally averaged it is larger than 45 for all the modes so that correlation larger than 0.28 could be considered to be significant at the 95% level. Mode 2 is highly correlated with the total signal (0.7) at the equator, whereas correlation is relatively low (of the order of 0.45) for mode 1, except in the far western Pacific. Correlation between mode 3 and the total signal is relatively high in the central Pacific where Ekman current anomalies are important, and low in the eastern Pacific and the far western Pacific. In the central Pacific mode 3 can be explained to a large extent by local forcing similarly to what we find for higher-order modes (cf. earlier discussion). What happens in the western and eastern ends of the basin is less clear. Although mode 3 has its maximum rms variability in these regions, correlation with zonal wind stress is very low there, indicative of high-frequency fluctuations not directly related to local forcing. An analysis of the influence of zonal change in density profiles on linear modes suggest that there is probably significant modal dispersion in these regions (discussion in section 4c). With the summed eight gravest modes correlation is very large everywhere along the equator.

Figure 4b presents the explained variance of the various modal contributions relative to the total ZCA. Whereas the modes 1 and 2 explain almost the same variance in the central Pacific, the variance is higher for mode 1 in the west and for mode 2 in the east where it reaches 70% at 100°W. Mode 3 explains up to 20% of the variance in the central Pacific, which is identical to mode 1. The percentage of variance explained by the sum of the three gravest modes is rather high (∼70%) except around 160°W where it drops to 50%. In this strongly wind-forced region the contribution of higher-order modes is required to reach 90% of explained variance.

#### 2) Sea level anomalies

We now discuss the vertical decomposition of pressure anomalies. Figure 5 presents time–longitude plots along the equator of the first, second, third, and fourth baroclinic mode contributions to the simulated sea level anomalies (hereafter SLA). The joint contributions of the first two modes and of modes 3 to 8 are also presented. Whereas variability for the first mode is largest in the western Pacific, the second mode is more active in the central and eastern Pacific. The sum of the first two modes captures the spatial structure of the simulated sea level anomalies although with lesser amplitude. Note that modes 2 and 3 contribute significantly to the sea level decrease during La Niña 1988. This is also true for the contribution of higher-order modes (sum of third to the eighth mode), which in mid-1988 contributes as much as mode 2. During this period high-order modes with finer vertical scales tend to be more active because the thermocline is shallower after an upwelling Kelvin wave has crossed the basin (Périgaud and Dewitte 1996).

##### (i) Correlation and explained variance

Figure 6 is similar to Fig. 4 but for SLA. Correlations larger than 0.3 are significant at the 95% level. Correlation is high between mode 1 and the total signal as well as between mode 2 and the total signal. Correlation between mode 1 and mode 2 is therefore high (averaging 0.7 along the equator) confirming results from previous studies (Busalacchi and Cane 1985; Giese 1989; Reverdin et al. 1996). Mode 3 is also highly correlated with the total signal in the central and eastern Pacific. For mode 3 both the correlation and the explained variance are very low west of 165°E. Note that correlation for these modes decreases by 0.2 near the date line. This is not the case for the contribution of higher-order modes (correlation between the sum of the fourth to the eighth mode and sea level is 0.8 at 180°) suggesting that near the date line part of the sea level variability is due to changes in the density profile with rather small vertical scales. This is consistent with variability associated with surface current convergence above the barrier layer near the salinity front often located close to the date line (Vialard and Delecluse 1998).

The percentage of explained variance is the same for modes 1 and 2 in the eastern basin. It is twice as large for mode 1 than for mode 2 west of 160°E. Mode 3 explained variance is weak in the western and central Pacific but increases eastward to the same level of explained variance as for mode 1 or 2 (∼20%). With 3 modes, the explained variance reaches 92% in the western Pacific and 80% in the eastern region but is as low as 35% around the date line. Even with eight modes, the percentage of explained variance does not exceed 60% there. This differs from ZCA and indicates that to recover the sea level variability in the central Pacific up to a 80% level, a large number of high-order modes is required.

#### 3) Phase speeds and projection coefficients

##### (i) Phase speed

Low-frequency linear waves have horizontal propagating phase speeds proportional to the separation velocity *c*_{n}. In the OGCM simulation the speed is a function of longitude and time (see section 2). The mode-1 phase speed *c*_{1}(*x, t*) along the equator is presented in Fig. 7a. It shows that temporal variability is not very large. The phase speed is usually increasing (decreasing) only by up to 8% from its average value during a warm (cold) event as expected from the temporal changes in stratification of the upper thermocline (see Fig. A2). This is also true for higher-order modes (not shown). For each mode, there is a sharp decrease in phase speed from the central toward the eastern Pacific (and also to the west) that is highly correlated with bottom topography. This is in agreement with results from Picaut and Sombardier (1993), who performed a similar computation using Levitus climatology (Levitus 1982) and a realistic bottom topography. Figure 7b presents zonal equatorial plots of the phase speed averaged over the 10-yr period for the first five modes. The pattern is in close agreement with the Picaut and Sombardier (1993) values, but the OGCM values are higher. This is to some extent due to the spatial smoothing involved in deriving the Levitus dataset and less so to differences between model and observations (the model thermocline is however slightly more diffuse in the simulation than in observations). The first-vertical-mode phase speed estimated from CTD profiles by Picaut and Sombardier is 2.84 m s^{−1} at 180° and 2.2 m s^{−1} at 110°W (see their Fig. 3a), whereas model mean values for these locations are, respectively, 2.78 m s^{−1} and 2.18 m s^{−1} (differences are less than 3%). The zonal gradient in phase speed around 120°W is less pronounced in their case than in the model simulation. Here the model ocean depth is used, which presents a large change near that longitude corresponding to one model vertical grid point (see Table 1). However, the difference in phase speed between the central Pacific and the eastern Pacific in the simulation is to a large extent due to zonal changes in stratification, as was demonstrated by performing the calculation maintaining a uniform 3000-m ocean depth (Fig. 7b). A similar conclusion is reached for the western Pacific. These results confirm Picaut and Sombardier’s conclusion that changes in phase speed due to density structure dominate those caused by bottom topography in the tropical Pacific.

##### (ii) Projection coefficient

Another parameter needed for linear model simulations that we can derive from this decomposition is the projection coefficient of the wind stress: It is defined as in Lighthill (1969) assuming that the wind stress acts like a body force in a mixed layer of depth *d* and not below:

where H_{1} is a dimensionalizing depth equal to 150 m in this study, and *d* is the depth of the near-surface layer of weak stratification in which momentum and buoyancy are quickly transported vertically by turbulent processes. We used a value for *d* of 50 m. For the gravest modes, the vertical functions do not vary much between 0 and *d* and are assumed to be equal to *F*_{n}(*x,* 0, *t*) = 1.0 so that

(using the varying vertical function within the mixed layer results in less than 5% deviation in the *P*_{n} for the gravest five modes).

Figure 8 presents average values of *P*_{n} along the equator for the first five modes. Not surprisingly, *P*_{n} is low in the central Pacific for individual higher modes. A large number of these modes is therefore required to reproduce the correct amplitude of the upper ocean shear. For high-order modes (3–5) *P*_{n} reaches similar values as for modes 1 and 2 near the meridional boundaries. Although zonal wind stress anomalies are usually weak there, mode 3 may be strongly wind forced along the equator near the eastern boundary. As pointed out by Eriksen (1988), when the thermocline is thin and shallow as in the eastern Pacific, forcing projects most efficiently onto the higher baroclinic modes. This is also consistent with results from Giese (1989), which demonstrate that high-frequency forcing projects preferentially on high-order modes. Westerly wind burst activity in the western Pacific (Maes 1996) acts toward increasing *P*_{n} for the high-order modes in this region. The similar amplitude of *P*_{n} for modes 1 and 2 is consistent with the simulation results and with other studies (Gill 1982; Giese and Harrison 1990) indicating that the second baroclinic mode is as important as the first one in contributing to the eastern Pacific variability. The temporal dependency is such in the central Pacific that *P*_{n} increases (decreases) during El Niño (La Niña) event for the first three baroclinic modes. In the east, it presents a seasonal cycle for which amplitude increases with mode number. High-order modes in the east present high-frequency fluctuations that do not show any coherency with local wind forcing.

The modal decomposition exhibits different characteristics for the first three modes. Higher-order modes behave rather similarly to locally trapped modes forced by the wind (with the exception of mode 4, which bears strong resemblance to mode 3). We will therefore concentrate in the following on the behavior of the first three modes considering that their sum accounts for a large part of the explained variance both in terms of ZCA (82.35% over NINO3eq, *c* = 0.91) and SLA (79.64% over NINO3eq, *c* = 0.92) (see also Figs. 4b and 6b).

## 4. Comparison with linear simulations

We are investigating whether the projection of the model fields on individual vertical modes behaves dynamically in a similar way to what is expected from linear theory, that is, with wave dynamics characterized by the phase speed and by the projection coefficients estimated from the OGCM simulation. There are many reasons why this should not be the case: 1) linear theory assumes spatially uniform density profiles for separation in vertical modes, 2) OGCM includes nonlinear terms and diffusive processes that also couple horizontal and vertical equations. The OGCM simulation will be refered to as OPA hereafter.

The approach proposed here is based on an intercomparison of OPA with simulations of a three-mode linear model. The model parameters are tuned to fit the modal phase speeds and projection coefficients *P*_{n} from the OGCM simulation. The linear model is from Cane and Patton (1984), an adiabatic shallow-water model on the equatorial *β* plane. The code is written for a rectangular basin extending from 28.75°S to 28.75°N, 124°E to 80°W with a resolution of 0.5° lat and 2° long. For each baroclinic mode, the model first solves the equatorial Kelvin mode propagating eastward and then derives the long low-frequency Rossby and anti-Kelvin modes propagating westward. The Kelvin and the first meridional Rossby components of the motion contribute to a large share of the equatorial variance in ZCA and SLA (see Table 2). We will therefore focus the comparison on the Kelvin and first meridional Rossby components. The linear model will also be a tool to investigate the effect of a zonal change in the phase speed and the density profile.

The wind stress forcing is the same as for OPA: wind stress anomalies from ARPEGE computed relative to the climatology over the 10-yr period. The linear model is forced by wind stress, but not by density fluxes. We found that density fluxes contribute little to SLA and ZCA variability in this simulation. A five-year spinup from rest is obtained by forcing the linear model with the year 1985 repeated five times.

Different linear simulations were done with modal phase speeds and projection coefficients corresponding to different areas. The characteristics of the simulations are summarized in Table 3. The standard simulation (refered to as runREF has constant phase speeds and projection coefficients characteristic of the central Pacific (150°W; Figs. 7b and 8, Table 3). RunPX is the linear simulation where we chose zonally varying projection coefficients (from Fig. 8) but retained the phase speeds of runREF. This run tests the impact of zonally varying projection coefficients on linear simulation. However runPX is very close to runREF and will not be discussed thoroughly (see Tables 4 and 5 for comparison with runREF). The linear formalism is also used to test the effects of zonally varying phase speed, assuming a zonal discontinuity of the vertical density profile at a prescribed longitude. Details on the analytical treatment of a zonally varying stratification are given in appendix B. Two cases will illustrate the effects of a density front: runCXE with the discontinuity at 120°W, runCXW with the discontinuity at 162°E (longitude independent projection coefficients as in runREF). The linear simulations are also sensitive to prescribed Newtonian damping friction formulation and to reflectivity properties of the zonal boundaries, which will not be tested thoroughly. Linear models with different formulations of friction allow vertical separation where the modal decay time (1/*r*_{n}) is of the form: *r*_{n} ∼ *c*^{−q}_{n} with *q* ranging from 0.5 to 2.5 (Gent et al. 1983). These choices result in more damping along the characteristics for the higher modes. We arbitrarily retained *q* = 0.5, which favors propagation of high-order modes. Studies of the sensitivity of the simulations on *q* indicate that the conclusions presented in this section also hold for other choices of *q*. The spindown time for the first vertical mode is chosen as 30 months. Alternatively for runCXW, we retain a stronger friction in the western part of the basin (see Table 3 and discussion in section 4c).

### a. Baroclinic mode contributions

#### 1) Zonal current anomalies

In this section, we present runREF ZCA rms variability maps, as well as the correlation and rms differences between the OGCM simulation and the linear simulation runREF along the equator. Variability maps for ZCA are displayed in Fig. 9, which can be compared to Fig. 1. Total surface current anomalies are derived from the contribution of the first three baroclinic modes added to Ekman current anomalies estimated as in Blumenthal and Cane (1989) assuming a 50-m mixed layer. The variability peaks in the western Pacific for modes 1 and 2 are very close to the OGCM simulation. However modes 2 and 3 do not present the variability extrema of OPA in the eastern Pacific. Mode 3 variability in runREF is also much lower (by a factor of 2) in the central and western Pacific leading to weaker total surface current anomaly variability in runREF than in OPA. High-order modes also contribute to the differences in total ZCA between the two simulations. For each of the first three modes the meridional scales of runREF and OPA rms current variability are very similar in the western and central Pacific. They increase east of 120°W for the three modes in runREF. This was not the case in OPA for modes 2 and 3; their meridional scales do not vary much over most of the basin. The meridional scales of runREF ZCA at 100°W estimated by fitting a Gaussian curve are 1000 km for mode 1, 700 km for mode 2, and 500 km for mode 3 (to be compared respectively to 800 km, 400 km, and 500 km for OPA). Another difference between OPA and runREF is that runREF rms spatial variability distribution for mode 2 is similar to the one for mode 1. It presents however larger values for mode 2, which has a shorter meridional scale and a similar input of wind stress forcing than mode 1 (*P*_{1} = *P*_{2}). Mode 3 does not peak at the eastern and western ends of the ocean as in OPA. RunREF mode 3 clearly is a locally forced mode with dynamics very close to the one of a frictional layer.

Figures 10 and 11 are equatorial sections of the correlation and rms differences between runREF and OPA for the zonal current and sea level anomalies for each baroclinic mode and for the total fields. For mode 2 ZCA correlation presents a maximum of 0.65 in the central Pacific whereas mode 1 correlation never exceeds 0.4. The ratio of high-frequency variability versus low-frequency variability is larger in runREF than in runOPA explaining the low correlation for mode 1. For mode 3 ZCA, there is no correlation in the far western Pacific. In this region, amplitude is small in runREF and relatively large in OPA. In runREF the relatively strong friction and small phase speed damped the solution, which is largely wind forced, whereas in OPA the response to wind forcing is not well identified in the western Pacific. Total surface current anomalies from runREF and OPA are well correlated in the central Pacific. Correlation is low near the meridional boundaries, which can be caused by the following reasons. First, the linear model does not have realistic coastlines near the western boundary, which can be detrimental to the solution in the western Pacific region (Du Penhoat et al. 1983) and might overestimate wave reflection. Furthermore, these boundary regions are characterized by a different density structure than the one in the central Pacific, which was used for runREF. Note that runPX does not perform better than runREF in these regions, suggesting that this low correlation is not due to inadequate wind forcing. The rms difference for the total ZCAs along the equator is 20 cm s^{−1} on average (Fig. 10). It is larger in the west than in the east where variability is less.

#### 2) Sea level anomalies

SLA correlation between OPA and runREF mode 1 contributions is highest in the west (Fig. 11). It reaches 0.85, which is as large as the correlation between OPA and the combined contribution of the first three modes of OPA, illustrating that sea level anomalies are dominated by the first baroclinic mode in this region as shown by other studies (Gill 1982, 1983; Delcroix et al. 1991). As expected mode 2 correlation is also high because the first two baroclinic modes are highly correlated in multimode linear simulation (Reverdin et al. 1996). Correlation between OPA and runREF mode 3 remains relatively high in the central and eastern Pacific but is insignificant in the far western Pacific. Although mode 3 amplitudes are small both for OPA and runREF in this region, OPA mode 3 exhibits more high-frequency variability relative to the low frequency compared to runREF mode 3. Along the equator correlation usually increases by 0.1 everywhere compared to mode 1 when summing the contributions of the three gravest modes. The correlation and rms difference zonally averaged along the equator are 0.7 and 2.7 cm, respectively (Fig. 11).

Considering the simplicity of the linear model compared to the OGCM the agreement between the two simulations for the first two baroclinic modes is encouraging. Discrepancies are larger for mode 3.

### b. Propagating characteristics

The previous comparison has shown overall good agreement between a multimode linear simulation and the OGCM simulation for the magnitudes of the anomalies of the two gravest modes. We now investigate the propagation characteristics in both models. In the linear model, the variability for a given vertical mode can be described as the sum of forced equatorial waves. At the equator, the Kelvin and first meridional Rossby (hereafter R1) waves provide the main part of the signal (Boulanger and Menkes 1995; see also Table 2). Their meridional characteristics depend on the vertical mode. For each vertical mode component, we estimate in the OGCM equivalent linear Kelvin and R1 mode, which are compared to the linear simulation. For a given vertical mode Kelvin and R1 coefficients are obtained by projecting the corresponding contribution of sea level and zonal current anomalies onto the complete set of meridional modes. For OPA the scaling coefficients (*c*_{n} for ZCA and *c*^{2}_{n}/*g* for SLA) and the determination of the horizontal modes meridional structures are estimated from the zonally varying mean phase speed given in Fig. 7b. The projection is done onto the (10°S–10°N) data. Land areas are filled by bilinear interpolation. Similar results are usually obtained if zero values are put instead on land areas. However, for Papua New Guinea, which is close to the equator, differences between the two methods are large and the projection coefficients will be uncertain. The issue of determining wave coefficients from sea level data has been addressed by Boulanger and Menkes (1995). The difference with their work is that both *u* and *p* surface fields are used instead of the *p* field alone, so this does not require additional assumptions. We also ignore the far western Pacific, which is a complex area of islands and passages to the South China Sea and the Indian Ocean. In the linear model simulations, the meridional boundaries are ideal walls.

#### 1) Lag-correlation for the R1 and Kelvin time series

##### (i) In the central Pacific

The propagation characteristics can be estimated from a lag-correlation analysis of the R1 or Kelvin components at different longitude. Figure 12 presents the results for the R1 component with a reference longitude at 180° for the first three baroclinic modes of OPA and runREF. The plots for runREF exhibit a clear low-frequency first meridional Rossby wave signal with high coherency (correlation >0.6) extending longitudinally over 4000 km. The pattern associated with positive correlation is rather broad in time for modes 2 and 3 compared to mode 1, expressing that the latter modes have more energy at low frequency than mode 1. The R1“beam” is particularly clear for mode 1 with high correlation lasting less than 5 months. OPA also presents evidence of propagation at low frequency for the first three vertical modes. However the high correlation zone does not extend as far longitudinally. The slope of the linear best fit to the points with maximun correlation at each longitude gives an estimate of *c*_{n}/3 (theoretical speed for unforced *n*th baroclinic first-meridional Rossby wave) as indicated for runREF. We find for runREF values for *c*_{n} of 2.73, 1.64, 1.04 m s^{−1}, respectively, for vertical modes 1, 2, and 3. This method applied to OPA yields 2.71, 1.45, and 1.11 m s^{−1} respectively for mode 1, 2, and 3. These values are all in close agreement with the free linear wave speeds obtained from the decomposition in vertical modes (see Fig. 7b, values at 180° are 2.75, 1.72, and 1.05 m s^{−1}).

There are other differences between the two simulations particularly for modes 1 and 3. In OPA mode 3, there are negative coherency zones on each side of the R1 beam that are not present in runREF at the same time lags, and suggest high-frequency fluctuations. A similar comment can be made for runREF mode 1 and to a lesser extent to OPA mode 1. Note also the sharp change in the propagating signal west of 160°E for OPA mode 3 where we find slower propagation (0.6 m s^{−1}), a characteristic not shared by runREF.

##### (ii) In the eastern Pacific

Similar plots are presented for a reference location at 110°W in Fig. 13. Propagation characteristics for the linear model are even clearer than in Fig. 12. This region presents less local wind forcing compared to the central Pacific so that we are mostly in the presence of free propagating waves. Very high coherency extends over 4000 km for runREF. Negative correlation beams for runREF mode 1 are due to multiple wave reflections onto the meridional boundaries. The secondary maxima at ±9 months result from the resonant basin mode Kelvin+R1 (Périgaud and Dewitte 1996). Mode 1 in OPA is similar to mode 1 in runREF except at lags exceeding 7–8 months (there is no energetic resonant mode for OPA first baroclinic component). For mode 2, the high correlation lasts longer than for mode 1 both for OPA and runREF. Zonal coherency is however less in OPA than in runREF. Differences between OPA and runREF are larger for mode 3. Phase speed values derived from a linear best fit to the correlation maxima gives 2.95, 2.25, and 1.32 m s^{−1} respectively for mode 1, 2, and 3, which is 10%–30% larger than the values derived from the vertical mode decomposition (values for runREF are respectively 2.71, 1.65, and 1.03 m s^{−1}, closer to the ones expected from vertical mode decomposition).

We draw similar conclusions from the Kelvin component analysis, but they are less certain because Kelvin waves propagate faster and the monthly resolution is too coarse for estimating with confidence phase speed values, in particular for mode 1. The phase speed values averaged for the two longitudes 180° and 110°W derived from the Kelvin waves are 2.70, 1.69, and 1.55 m s^{−1} for the OPA mode 1, 2, and 3. We find 2.1, 1.47, and 1.01 m s^{−1} for runREF vertical mode 1, 2, and 3. This also suggests that phase speed of the OGCM Kelvin contribution is about 20% larger than the linear phase speed. Along with the fact that high correlation bands across the Pacific do not extend longitudinally as much for OPA as for runREF, this suggests that processes other than linear long-wave propagation contribute to the OGCM variability. Modeling work by Giese (1989) indicates that nonlinear self-advection makes an important contribution to the phase speed increase of model Kelvin pulses. The increase in phase speed has also been described as Doppler shifting by the equatorial undercurrent (McPhaden et al. 1986). It is also possible that modal dispersion is involved in explaining such increase in phase speed for modes higher than the gravest one.

#### 2) Extended EOF analysis

To investigate further the characteristics of the different modes in OPA and runREF, we decompose the R1 and Kelvin time series in empirical orthogonal functions to extract the most significant features and realize frequency sampling (Kutzbach 1967). Such techniques have been used routinely for meteorological time series. A suitable technique for a case where propagating structures are expected consists in introducing the temporal structure in the basic observation vectors by padding together data taken at adjacent or distant times. This is the extended EOF (E-EOF) technique introduced by Weare and Nasstrom (1982), which is identical to the Multi-Channel Singular-Spectrum Analysis (M-SSA) analysis technique (Vautard and Ghil 1989). It has been used in numerous studies (Lau and Lau 1990; Jiang et al. 1995; Robertson et al. 1995a, 1995b; Vintzileos 1996, among others). The E-EOF decomposition has been performed considering both Kelvin and R1 time series. Because of the length of the time series and the timescales of the phenomena investigated (seasonal to several years), we will use a time window of 2 yr with a 2-month interval. This leads to a 2 × 12 × *N*-dimensional “observation” vector where *N* is the number of longitude grid points (fields have been interpolated on a 5° grid extending from 135°E to 80°W so that *N* = 30). The variance maximization procedure for these extended vectors leads to the diagonalization and eigenvalue problem of the lagged-covariance matrix. This often produces pairs of modes with similar spatial structure and explained variance but with a quadrature phase shift. The similar characteristics imply that these modes essentially describe the same propagating phenomenon. Modes that are not paired often correspond to nonpropagating phenomena.

In the following we will concentrate on the first dominant E-EOF mode pair for each vertical mode because this pair is the one most representative of the El Niño–La Niña oscillation. All other modes are at higher frequencies. Most of them exhibit characteristics that are more difficult to interpret and correspond to peculiarities of the ARPEGE wind forcing. The only “high frequency” mode we need to consider is the resonant mode for the first vertical mode of the linear simulations. This mode contributes to a significant share of the total variance and is well extracted from the time series by E-EOFs.

Figures 14a–c and 15a present the results of the decomposition for the Kelvin (left-hand side) and R1 (right-hand side) coefficients of the first three baroclinic modes. For clarity only one element of the pair in each statistical mode is displayed (statistics for both elements of the pair are given in Table 4). The associated time series (middle) are shifted in time so that the time-lag correlation is maximum between the OPA and runREF time series. Then an initial time for the spatial structures is selected corresponding to the earliest common date within all the time series in each figure (January 1985 for Fig. 14 and June 1986 for Fig. 15). We find high correlation between the OPA and runREF time series (see Table 4) implying that we are in the presence of the same statistical modes for OPA and runREF. We display the spatial structures for three consecutive time sequences (separated in time by 2 months). Thus for Fig. 14 the spatial structures are for January, March, and May 1985. The overlaying period between the two simulations is limited because of the limited span of the time window used for the decomposition. For the first two baroclinic modes time–space correlation is estimated based on seven consecutive sequences (Fig. 14a and 14b) both for the Kelvin and R1 coefficients. It exceeds 0.80 (see Table 4) indicating good agreement between OPA and runREF both for the Kelvin and R1 coefficients for each of the first two baroclinic modes. This is not the case for the third baroclinic mode with negative and insignificant correlation both for the R1 and Kelvin spatial structures. Another discrepancy between OPA and runREF for vertical mode 3 (Fig. 15a) is the percentage of variance explained by this first pair of E-EOFs: 60% for runREF and only 32% for OPA (see Table 5). This originates in that mode 3 in runREF is mainly a low-frequency mode wheras in OPA it also includes high frequencies.

The leading E-EOF mode for vertical mode 1 in runREF explains less variance than the one in OPA. This is due to the linear simulation near-resonant basin-scale oscillation at 9–10 months, which accounts for a large amount of the total variance. The time series indicate that it contributes to the simulation over the whole period and explains 18% of the total variance (see Table 5). It is extracted as a subharmonic of the 9–10 month mode by the E-EOF decomposition in a specific pair of E-EOFs (see Fig. 14c). Its associated zonal wavenumber is given by the spatial structures (1.5 wavelengths across the basin for R1 and 0.5 for Kelvin, which clearly indicates its resonant character). Note also that the phase speed corresponds to first baroclinic free wave speeds 2.7 m s^{−1} for the Kelvin component and 0.9 m s^{−1} for the Rossby one.

### c. Impact of zonally varying stratification: RunCXE and runCXW

OPA and runREF present fairly good agreement for the two gravest vertical modes, at least at low frequencies. The major difference for the gravest modes between the two simulations is the excessive high-frequency variability present in vertical mode 1 for runREF compared to OPA. This could be an issue for coupled models involving linear multimode model of the ocean. Moreover the previous analysis suggests that linear models will not produce realistic high-order modes, which, albeit of weaker variability than modes 1 and 2, can affect the surface ocean signature. Périgaud and Dewitte (1995) and Périgaud et al. (1997) showed how sensitive a simple coupled model could be to small and localized variations in the oceanic fields. For instance, they showed that parameterizations resulting in changes of the downwelling (upwelling) anomaly in the far eastern Pacific in Zebiak and Cane’s (1987) simple coupled ENSO model modify drastically its behavior within two years of the initialization.

The regions where long-wave propagation is likely to be most altered correspond to the largest zonal gradient in the vertical stratification. Figure A1 (average buoyancy profile as a function of longitude) indicates that this happens around 160°E and 120°W. Thus we want to know how sensitive the discrepancies of runREF with respect to OPA are to the specificied stratification. In particular what is the effect of zonally varying stratification?

The problem of a propagating wave in a zonally varying stratified fluid has been addressed analytically by Busalacchi and Cane (1988, hereafter BC). They considered two cases: 1) when the Väisälä frequency *N* may be written as *N* = *N*_{1}(*x*)*N*_{2}(*z*), preventing any modal dispersion from one region to the next (dispersion coefficients are then given by *γ*_{m,n} = *δ*_{m,n}—see appendix B for a definition of *γ*_{m,n}) and 2) *N* derived from theoretical density profiles. Their results suggest that zonally varying stratification does not produce substantial changes in the energy flux of propagating equatorial waves. However, the second case indicates significant modal dispersion related to relatively large values of nondiagonal *γ*_{m,n} (see their Table 2 and discussion). We adopted their approach using density profiles from the OGCM simulation so that solutions are obtained for realistic forcing. The methodology used to solve the equations for the linear model with a density front is summarized in appendix B. Whether the density front is close to the eastern or to the western boundary results in differences in the dispersion coefficients.

For runCXE, the density front is set at 120°W and density profiles for the western and eastern parts of the basin are chosen from the OGCM simulation, respectively at 135°W and 110°W. This provides the coefficients for the dispersion matrix and the phase speeds for each part of the basin (see Table 3). Other locations were tested on each side of 120°W resulting in changes for the coefficients smaller than 10%. The vertical structure in the eastern basin (east of 120°W) is subject to relatively strong annual variability, which appears most prominently for the high-order vertical modes. Thus the coefficient coupling the third mode to the west with the second mode to the east is (*γ*_{3,2}) = +0.2 in winter and spring and it drops to −0.1 for August and September. We chose to keep the annual cycle for the dispersion matrix in this experiment because the simulation is very sensitive to those coefficients. Friction is the same for each basin. The eastern meridional boundary has been shown to be a good reflector for long equatorial waves (Boulanger and Fu 1996). After reflection onto the eastern meridional boundary free waves for different modes may combine to propagate vertically as shown for the annual cycle by Kessler and McCreary (1993). This tends to reduce the surface amplitude of the anomalies. A way to take that into account in our model would be to add more friction. We prefer not to do that in the absence of appropriate studies for the interannual timescale. Note also that, if we use a strong friction in the eastern part of the basin, it would also damp the impinging Kelvin wave, which would be unwarranted.

For runCXW, the density front is now placed at 162°E and density profiles for the western and eastern parts of the basin are derived from the OGCM simulation at 160°E and 160°W. They provide coefficients for the dispersion matrix and the phase speeds for each basin (see Table 4). The far western Pacific is a complex area where many processes other than long equatorial waves propagation are at work. Due to the presence of numerous islands, reflection of long equatorial waves is also altered (du Penhoat and Cane 1991; Clarke 1991), which could result in the dissipation of a significant part of the incoming Rossby wave energy. We model this by introducing in the linear model a large Rayleigh type friction west of the density front.

Results are presented as before in terms of the dominant E-EOF pairs for runCXE (Fig. 15b) and runCXW (Fig. 15c) baroclinic mode 3. The results for the first two baroclinic modes in runCXE and runCXW are summarized in Table 4. Compared to what was found with runREF, E-EOF dominant modes for vertical modes 1 and 2 still present a high correlation with equivalent E-EOFs of OPA both for Kelvin and R1 coefficients. However rms difference between OPA and the linear simulations runCXE and runCXW is on average a little bit larger except in the far eastern Pacific for runCXE, where the Kelvin coefficient for each baroclinic mode undergoes an increase compared to runREF. This is consistent with a decrease in phase speed east of 120°W and in closer agreement with OPA in the east.

The Kelvin spatial structures for vertical mode 3 in runCXE and runCXW are in much better agreement with OPA compared to runREF. Correlation between runCXE and OPA is now significantly positive. However the R1 spatial structures are negative west of 120°W for runCXE, whereas they are almost zero and positive for OPA so that the space–time rms difference between OPA and runCXE is larger than the one between OPA and runREF. RunCXW performs better in reducing rms differences together with a very high correlation. The high-frequency variability relative to the low-frequency one is increased. The percentage of explained variance of the leading E-EOF for vertical mode 3 is now 41% for runCXE and 36% for runCXW (see Table 5) closer to the OPA value. Nevertheless, for runCXE there is still too much high-frequency variability for the gravest modes, in particular for vertical mode 1. The percentage of explained variance for the resonant mode in runCXE (23%) is even larger than in runREF (18%). It contaminates vertical mode 3 via modal dispersion suggested by the transmission coefficient of 3.25 for a third baroclinic transmitted Kelvin wave resulting from an incident first baroclinic Kelvin wave of unit amplitude (see Table 6a and appendix B). Note also that the associated coefficient for the third baroclinic first meridional reflected Rossby wave is also large in absolute value (−3.19), which may explain the change in the shape of the R1 spatial structures in runCXE compared to runREF. RunCXW is more successful in producing a percentage of variance for the first and third baroclinic modes comparable to the one in OPA together with a diminution in the energy of the resonant mode for vertical mode 1 (see Table 5).

## 5. Discussion

Linear theory predicts that the baroclinic response to forcing can be expressed as a sum of vertical standing modes. In this study we decomposed fields from an OGCM simulation in the linear vertical modes. As the primitive equations in the OGCM do not fullfill all the required assumptions for dynamical separation between the horizontal and the vertical, coupling between vertical modes is expected. The sensitivity of linear simulations to the spatial variability of stratification discussed in section 4 illustrates an important process responsible for coupling vertical modes. In our study, mode 3 is the most affected and presumably higher-order modes will also be affected. Note that the simulations runCXE and runCXW are sensitive to the number of vertical modes we retain. The continuity equations (see appendix B) at the interface were arbitrarily truncated. However retaining more modes acts toward more intense coupling between them so that we expect effects comparable to the ones identified with three modes (simulations with four modes were carried out and present similar characteristics to the ones with three vertical modes). Impact onto modes 1 and 2 is difficult to assess because the two modes are highly correlated for all the simulations, including runREF for which stratification is spatially uniform. Correlation for the temporal functions of the dominant E-EOF modes between vertical mode 1 and 2 exceeds 0.8 in runREF as well as in runCXE or runCXW (see also the time series in Fig. 14). Analyzing other statistical E-EOFs for the gravest vertical modes is not so conclusive because of their sensitivity to parameter selection (friction, phase speed, projection and dispersion coefficients). As already mentioned, they may also correspond to peculiarities of the ARPEGE wind forcing.

Previous theoretical studies investigated the effect of zonally varying stratification. Most of them used energy flux conservation equations in the WKB limit (Gill 1983; Giese and Harrison 1990), which involve solely two integrated quantities, the phase speed *c*_{n} and the vertical structure norm *∫*^{H}_{0} *F*^{2}_{n}(*z*) *dz.* In this limit there is no scattering of energy from one vertical mode to another and there is no reflection of the Kelvin wave into westward propagating Rossby waves. Alternatively Busalacchi and Cane (1988, hereafter BC), using the pressure and zonal current continuity equations, yields in addition the *γ*_{i,j} coefficients, which brings further insight about the effect of zonal changes in the density structure. However, as pointed out by Giese and Harrison (1990), BC used exponential density profiles that were not realistic enough [see Fig. 15 and discussion in Giese and Harrison (1990)] and they also assume a flat bottom ocean. In our study we used density profiles from an OGCM simulation that are more realistic (cf. Fig. A2 with Fig. 8 in Lukas and Firing (1985)] along with spatially varying ocean bottom. The derived nondiagonal dispersion matrix coefficients suggest stronger coupling between vertical modes compared to the BC study (see appendix B for details). Results are presented in Tables 6a and 6b in terms of dimensionalized wave coefficients, which suggest a large modal dispersion. For instance, for runCXW a contribution from modes 2 and 3 to mode 1 and a contribution from modes 1 and 2 to mode 3 is non-negligible, consistent with the results just presented. Note that runCXW would not perform so well if we had not used a strong friction in the western part of the basin (see discussion in section 4c), suggesting that a standard linear model overestimates reflections onto the meridional western boundary. In runCXW, the major part of the energy passing eastward through the density front reaching the central basin either is generated by wind forcing close enough to the western edge of the density front or results from part of the gravest baroclinic modes waves being reflected onto the western boundary that is not completely dissipated. This process can strongly modify the response of standard multimode linear models and may explain that observational evidence for the delayed oscillator theory (Schopf and Suarez 1988; Battisti and Hirst 1989) remains equivocal.

We also expect that vertical mixing and advection terms in the momentum equations will also result in coupling between the vertical modes. In the presence of vertical energy propagation this will occur along WKB ray paths into the deep ocean as discussed in a linear model by McCreary (1984). There has been observational evidence that this latter process takes place on an annual timescale (Kessler and McCreary 1993), but we have not investigated whether vertical energy propagation is a significant process in OGCMs and whether it also happens at other frequencies.

## 6. Conclusions

A primitive equation high-resolution model of the tropical Pacific was forced by the daily surface fields of the ARPEGE AGCM for the 1985–94 decade to investigate the contribution of baroclinic modes to surface zonal current and sea level anomalies on interannual timescales. We find that the two gravest modes contribute with similar magnitudes to surface currents and sea level. This differs from the results of a recent study by Fukumori et al. (1998), which shows from another OGCM simulation that sea level variability in the tropical Pacific is dominated by the first baroclinic mode. Note however that the OGCM used in their study has a much lower resolution (12 vertical levels with 4 levels in the top 325 m and 1° in latitude) suggesting that vertical and meridional resolutions have a strong impact on such calculations.

Although the first two modes appear largely correlated, they present different spatial patterns. Mode 1 variability is the largest in the western and central Pacific, whereas mode 2 is more active in the eastern Pacific. The zonal current and pressure anomaly contributions from each baroclinic mode are decomposed in terms of long equatorial waves and are compared to the results of a multimode linear model forced by the ARPEGE wind stress fields. A similar spatial distribution of variance is found for the first two modes in the OGCM and in the linear model simulations. Kelvin and first-meridional mode Rossby wave propagations are observed in the OGCM for the first three baroclinic modes at phase speeds comparable in the central Pacific to those of unforced linear waves but with speeds larger than those of linear waves in the east. In the linear simulations, similar propagation characteristics are observed but at speeds closer to unforced linear waves. Modal dispersion due to zonal change in the density field is shown to improve the mode 3 Kelvin and first-mode Rossby waves linear simulation. The two gravest modes are also affected.

This study thus suggests that the analysis of sea level height data assuming a single baroclinic mode has strong limitations. This may explain some of the difficulties encountered by Boulanger et al. (1997) in recovering the correct amplitude for the first-mode Rossby wave with a single mode wave model both in the central and eastern Pacific. In our study, with the “standard” multimode linear simulation (runREF) we find close agreement between the linear model and the OPA Kelvin and first meridional mode Rossby coefficients for the first two baroclinic modes both in terms of phase and amplitude. Considering a density front in the linear simulations successively in the east and in the west leads to better agreement between the linear simulations and OPA wave coefficients for vertical mode 3. To pursue further the investigation of Boulanger et al. (1997), the role of wind forcing and the reflection efficiency of both eastern and western boundaries in the OGCM will have to be investigated for each baroclinic mode taking into account the modal dispersion process.

Another result of our investigation is that the contribution of higher-order modes depends significantly on the local density structure. Modal dispersion caused by zonal change in stratification occurs mainly for these modes, but there is also redistribution of energy from these high-order modes to the gravest modes. Salinity structure contributes to the density structure, in particular near the sea surface in the west-central Pacific where a barrier layer is present [see the discussion of its dynamics in Vialard and Delecluse (1998) from another OGCM experiment]. A change in the salinity structure due to wind or heat flux forcing will have an effect on the vertical modes, which could result in complex modal dispersion. This, also, will have to be investigated. As simple to intermediate ocean–atmosphere coupled models are being developed, this study also calls for the investigation of the impact of zonally varying stratification in a coupled context.

Finally, a question raised by this work is whether separate contributions of individual baroclinic modes can be extracted from the time sequence of sea level observations, for instance, from satellite altimetric measurements. Sea level from altimetry is indeed usually interpreted as “one vertical mode” to separate Kelvin and Rossby contributions. This study points out that it can be an inaccurate assumption. We will need to investigate further the OGCM zonal current and sea level baroclinic contributions to determine if this is feasible.

## Acknowledgments

Computing supports for running the OGCM was provided by the Institut du Developpement et des Ressources en Informatique Scientifique. Corinne Michau (CNES) provided invaluable technical assistance. The authors would like to thank Pascale Delecluse and Jean François Minster, who were instrumental in the inital stage of this work. We particulary appreciated stimulating discussions with Rui Ponte during his stay at LEGOS. Helpful comments and suggestions were provided by Claire Périgaud during Boris Dewitte’s visits at JPL. We also thank Vincent Echevin for interesting discussions and Eric Guilyardi for presenting some results of this study at the IAPSO conference in Melbourne (July 1997). Thanks are adressed to Ken Clarke, Denise Kriesel, Eric Houplain, and Markus Durstewitz for their support during the course of this study. Fruitful discussions with Yves du Penhoat, Joël Picaut, and Jean-Philippe Boulanger were also appreciated. Finally, the authors wish to thank Paul S. Schopf, Julian P. McCreary, and an anonymous reviewer for their helpful comments. Funding for B. Dewitte and G. Reverdin was provided by Centre National de la Recherche Scientifique. Funding for C. Maes was provided by AGORA (Project of the environment program of the EEC).

## REFERENCES

### APPENDIX A

#### Influence of *N*^{2} Temporal and Spatial Variability on to the Vertical Mode Decomposition

The linearization of the equations about a basic state along with the assumption that the vertical stratification is horizontally uniform and assuming specific formulations of vertical and horizontal mixing leads to a Sturm–Liouville system. Its eigenfunctions are the vertical normal modes and its eigenvalues the corresponding phase speeds. In standard linear theory the result of this computation is constant in space and time.

As pointed out by Busalacchi and Cane (1988), hydrographic observations indicate, however, a rich structure in the density field at low latitudes. They particularly investigated the zonally varying stratification whose most prominent expression is the zonal slope of the thermocline (Meyers 1979). Zonal variation in the density field is present in the OGCM simulation as illustrated by the mean buoyancy profile along the equator (see Fig. A1). There is also temporal variability of *N*^{2} in the OGCM thermocline. Figure A2 presents profiles in the central Pacific for El Niño (January 1987) and La Niña (September 1988) periods indicating respectively an increase and a decrease in stratification. All these changes have an impact on the vertical mode decomposition, which we discuss in the following.

##### Decomposition with 〈N^{2}(x)〉_{t}: Impact of time dependency

The decomposition is performed with the time mean *N*^{2}(*x*). Figure A3 present maps of rms variability of the vertical mode contributions to the zonal current anomalies (to be compared with Figs. 1b–d and Fig. 2a). It shows little difference with the decomposition for which the temporal variability in *N*^{2} is retained. The rms difference is less than 4 cm s^{−1} for the first three baroclinic modes all over the basin.

##### Decomposition with 〈N^{2}〉_{t}: Impact of zonal dependency

The decomposition is now performed using the profile for *N*^{2} at longitude 160°W: 〈*N*^{2}(*x* = 160°W)〉_{t}. Figure A4 is similar to Fig. A3 but for this new decomposition. Rms variability is less (except at 160°W) for all modes compared to the results using *N*^{2}(*x, t*) (see Figs. 1b–d and Fig. 2a) or 〈*N*^{2}(*x*)〉_{t} (see Fig. A3). Differences are largest for high-order modes. Mode 3 variability is now largest in the central Pacific and reaches only 4.5 cm s^{−1}.

### APPENDIX B

#### The Linear Model and the Inclusion of a Density Front

We start from the linear shallow-water equations on an equatorial *β* plane written for the *n*th vertical mode (see Cane and Sarachik 1981):

with (*τ*^{x}, *τ*^{y}) being the surface stress, *r*_{n} equals *r*(*c*_{n}*β*)^{−q} where *r*_{n} is the Rayleigh friction and *P*_{n}, the projection coefficient defined as in Lighthill (1969) assuming that the stress acts as a body force in a mixed layer of depth *d* and not below:

The equations are nondimensionalized in the usual equatorial way with velocities ≅ *c*_{n}, lengths ≅(*c*_{n}/*β*)^{1/2}, frequency ≅(*c*_{n}*β*)^{1/2}, and pressure ≅*c*^{2}_{n}. Knowing *P*_{n} and *c*_{n}, the set of Eqs. (B1) can be solved numerically with appropriate boundary conditions for the tropical Pacific (see Cane 1984; Cane and Patton 1984).

Although the effect of zonally varying stratification on equatorial wave phenomena is neglected in standard linear model, it is possible to consider a single zonal discontinuity in the density field of a shallow water system without loosing too much of its simplicity. The formalism was developed by Busalacchi and Cane (1988). We will use their notations in the following. Their method is based on Cane and Sarachik (1981) and is similar to the technique used by Cane and du Penhoat (1982) to solve for the flow past an equatorial island. It uses the continuity of zonal velocity and dynamic pressure at the meridional density front, the truncation at some upper limit of the infinite sum of Rossby structures, and neglects the short Rossby waves contribution. The reader is invited to refer to BC for more details.

The method requires the determination of the transmission coefficients of a Kelvin wave and of a Rossby wave across an interface. This is done by solving linear systems [see Eqs. (2.12) and (3.3) in BC] knowing the dispersion coefficient *γ*_{m, n} = ∫^{0}_{−H} *F*_{m}*F*^{′}_{n }*dz,* where *F*_{m} is the vertical structure of the *m*th baroclinic mode in the western ocean (the one west of the meridional density front) and *F*^{′}_{n}, the vertical structure of the *n*th baroclinic mode for the eastern ocean (the one east of the meridional density front). The calculation in BC assumes a flat bottom ocean and dispersion is only due to zonally varying stratification. In the context of the OGCM simulation with realistic topography and possible depth jumps at the density discontinuity we use an artifact to extend the method. We extend the vertical structure profile on the shallow side to the depth of the deeper one as follows: assuming *H*_{W} > *H*_{E} (depth of the western ocean > depth of the eastern ocean), we define *F*^{′}_{cn} as follows:

Thus the new set (*F*^{′}_{cn})_{n=1,M} is still orthogonal and the formalism in BC holds.

We solved the linear systems [Eq. (2.12) and Eq. (3.3) in BC] in the case of an incident Kelvin and first-meridional Rossby wave of unit amplitude impinging onto the density fronts used in runCXE and runCXW (see Table 3 for the values of the *γ*_{i,j}). Tables 6a and 6b display the amplitude of the transmitted and reflected Kelvin, first- and third-meridional Rossby waves for the various baroclinic modes considered (discussion in section 4c).

There is a constraint on the location and the number of density fronts that the model can handle. Indeed the numerical scheme (Cane and Patton 1984) imposes that we have *c*DT < *L,* where *c* is the phase speed, DT is the time step of the model (typically 10 days), and *L* is the longitudinal extension of the portion of the ocean where *c* is constant. For waves propagating at 2.9 m s^{−1}, fronts have to be separated by at least 22°. The westernmost density front could only be at 148°E.

Note that we can also easily include a longitude-dependent projection coefficient *P*_{n}. This coefficient only intervenes in the projection of the wind onto the meridional structures and in the calculation of the Ekman shear. Thus it consists in multiplying the wind by a zonally varying coefficient instead of a constant one. This does not affect the method used for solving the equations.

The algorithm is based on Cane and Patton (1984). Starting from the eastern boundary, we integrate the Rossby component for the eastern ocean. At the interface, transmitted Rossby waves are computed imposing the continuity of pressure and zonal velocity. The reflected Rossby part from the incident Kelvin wave computed at the eastern edge of the western ocean is also derived from the continuity equation (see BC). Both components are added providing the condition at the interface to start the integration of the Rossby and Kelvin components in the western ocean. We can then compute the Kelvin part along characteristics in the eastern basin.

## Footnotes

*Corresponding author address:* Dr. Boris Dewitte, LEGOS/UMR 5566/CNES, 14, Ave. Edouard Belin, 31401 Toulouse, Cedex 4, France.

Email: bxd@pontos.cst.cnes.fr