## Abstract

Singular vector analysis and Floquet analysis are carried out on a linearized variant of the Zebiak–Cane atmosphere–ocean model of El Niño–Southern Oscillation (ENSO), hereinafter called the nominal model. The Floquet analysis shows that the system has a single unstable mode. This mode has a shape and frequency similar to ENSO and is well described by delayed oscillator physics. Singular vector analysis shows two interesting features. (i) For any starting month and time period of optimization the singular vector is shaped like one of two nearly orthogonal patterns. These two patterns correspond approximately to the real and imaginary parts of the adjoint of the ENSO mode for the time-invariant basic-state version of the system that was calculated in previous work. (ii) Contour plots of the singular values as a function of starting month and period of optimization show a ridge along end times around December. This result along with a study of the time evolution of the associated singular vectors shows that the growth of the singular vectors has a strong tendency to peak in the boreal winter. For the case of a stochastically perturbed ENSO model, this result indicates that the annual cycle in the basic state of the ocean is sufficient to produce strong phase locking of ENSO to the annual cycle; it is not necessary to invoke either nonlinearity or an annual cycle in the structure of the noise.

The structures of the ENSO mode, of the optimal vectors, and of the phase locking to the annual cycle are robust to a wide range of values for the following parameters: the coupling strength, the ocean mechanical damping, and the reflection efficiency of Rossby waves that are incident on the western boundary. Four variant models were formed from the nominal coupled model by changing the aforementioned parameters in such a way as to (i) make the model linearly stable and (ii) affect the ratio of optimal transient growth to the amplitude of the first Floquet multiplier (i.e., the decay rate of the ENSO mode). Each of these four models is linearly stable to perturbations but is shown to support realistic ENSO variability via transient growth for plausible values of stochastic forcing. For values of these parameters that are supported by observations and theory, these results show the coupled system to be linearly stable and that ENSO is the result of transient growth. Supporting evidence is found in a companion paper.

## 1. Introduction

Since Zebiak and Cane (1987), the most credible hypothesis to explain El Niño–Southern Oscillation (ENSO) is that it is a large-scale deterministic system, the seemingly random behavior of which can be explained by nonlinear dynamics. That is, ENSO can be characterized by a set of unstable equations where the growth of the ENSO mode is limited by nonlinearity and the erratic behavior is caused by chaos. Much literature has been devoted to studying the chaotic behavior of the Zebiak–Cane model (ZC model). See, for example, Jin et al. (1994) or Tziperman et al. (1994).

However, an alternate hypothesis is now gaining favor: that the ENSO dynamical system is approximately linear, and that its erratic behavior is caused by external stochastic forcing. For example, Penland and Sardeshmukh (1995) built a statistical model of the Indian and Pacific basins and drove it stochastically as a simulation of ENSO. Also, Chang et al. (1996) run several computer models, some of which are unstable and chaotic, and others that are stable and stochastically driven. Similarly, Moore and Kleeman (1999) took the intermediate model of Kleeman (1993), subjected it to stochastic forcing (in both stable and unstable configurations), and showed that it exhibited ENSO-like statistics.

In the stochastic hypothesis, the main source of variability in the ENSO system would be due to the occurrence of “optimal initial conditions” in the external forcing that will cause transient growth. In a series of papers, Farrell (1988a,b and others) shows that for a non-self-adjoint dynamical system, there sometimes exist initial conditions that temporarily induce a large amount of growth—even in a linear damped system. Several studies have looked at the optimal initial conditions (also called the singular vectors) in either the ZC model or the Battisti (1988) version of the ZC model (hereinafter called the B88 model). Specifically, Chen et al. (1997), Xue et al. (1997a,b), and Thompson (1998a, hereinafter T98a) explored the singular vectors in the context of predicting the error growth due to uncertainties in initialization. In this paper, the characteristics of the optimals are explored in the context of constructing a plausible stochastic model of ENSO.

The four possible scenarios for ENSO dynamics are (i) ENSO is an unstable chaotic mode; (ii) ENSO is an unstable, saturated mode, where external noise causes transient growth that either cancels or reinforces the warm and cold events so that they take on their irregular appearance; (iii) ENSO is a neutral or slightly damped mode where noise alters the phase and amplitude of the peaks so that they take on their irregular appearance; and (iv) ENSO is a highly damped mode where noise and transient growth create all warm and cold events.

What is meant by the terms “highly” damped versus “slightly” damped in scenarios iii and iv? Loosely speaking, a highly (slightly) damped system will have a short (long) decorrelation time relative to the time between events. Hence, a slightly damped system would be one where an ENSO mode that is initiated will last for at least a few cycles, unless modified by subsequent noise, while in a highly damped system, transient growth leading to an event would not last beyond that event. Under scenario iv, if the time spectrum of the noise is white, ENSO events will occur independently of each other.

This study is limited to linear models, and so scenarios i and ii, which require nonlinearities, will not be explored. Scenarios iii and iv can be addressed using a linear model, and the linear ocean–atmosphere model developed in T98a will be used. However, this model is unstable when using the background states and parameterizations of either the ZC model or the B88 model and so will need to be stabilized if it is to be run as a stochastic model.

Most previous intermediate ENSO models use parameterizations that put them in an unstable regime. This class of models is distinguished by the use of a steady-state atmosphere, either empirically derived (e.g., Barnett et al. 1993) or based on a particular, reasonable dynamical scenario [e.g., Gill (1980), and Lindzen and Nigam (1987) models]. As a result these models are devoid of atmospheric noise and therefore produce the ENSO cycle by normal mode growth and not by transient growth. Certain key parameter values were prescribed in these early coupled atmosphere–ocean models, but without the benefit of guidance from observations and theory. However, recent observational and theoretical results can now be used to better constrain these parameter values. By including the more realistic values in the models developed in this study, the coupled system becomes linearly stable, which suggests that noise in the system is important for maintaining the variance in the model ENSO mode.

This study consists of two papers. In this, the first paper, the linear model from T98a will be extended to include an annual cycle in the background states. This model will be analyzed, and then four linear stochastic models will be designed using the results of the analysis. The main characteristics to be controlled in the model design will be the decay rate of the principal eigenmode (the ENSO mode) and the amount of transient growth due to the singular vectors. In the companion paper these four “candidate” models will be run and their output statistics compared to observations. Also in the companion paper the implications of these stochastic models for the predictability of the ENSO system will be analyzed.

The structure of this paper is as follows. Section 2 contains a brief review of the model physics. In section 3, Floquet analysis (analysis of the modes of a system with a cyclic background) and singular vector analysis of the model (configured to be a linear version of the B88 model) are performed. It is shown that the shape of the principal singular vector is approximately one of two orthogonal shapes for all start months and optimization periods, and that these two shapes correspond to the real and imaginary parts of the adjoint of the annually averaged background case of T98a. It is also shown in section 3 that independent of start months and optimization period, the singular vectors evolve into an ENSO mode that tends to peak during the boreal winter. This implies that in a stochastic model, the phase locking of ENSO to the annual cycle would occur even under forcing that is white in time.

In section 4 sensitivity studies are performed to determine the effect of several key parameters on the model stability and transient growth. Specifically, the coefficient of the western boundary reflection is lowered, and the ocean mechanical damping is increased to values that are in agreement with recent studies. Four candidate models are then designed to have stabilities that range from slightly damped to highly damped. The leading mode in each of the four models is still ENSO-like and evolves in a manner consistent with the delayed oscillator physics. However, it is shown that through these parameter changes, the modal growth is substantially reduced with only a moderate reduction in the transient growth. It is also shown that the tendency of the system to phase lock to the annual cycle is robust to these parameter changes. Last, the Niño-3 index for a simulation of each candidate model is compared with the same index for observational data. However, the bulk of the model comparisons are left to the companion paper. Section 5 is the summary.

## 2. Model description

The model used in this study is a linearization of Battisti’s version of the Zebiak and Cane (1987) model of the atmosphere–ocean system in the tropical Pacific. The atmosphere is the Gill (1980) two-layer model of the tropical atmosphere, with heating a function of SST (evaporation) and surface wind convergence, following Zebiak (1986). The ocean is composed of a 1½-layer linear dynamic model of the Pacific basin (30°N–30°S, 124°E–80°W); the interface between the layers is equated with the depth of the 20°C isotherm. The equation for the sea surface temperature includes the effects of horizontal advection, surface fluxes, and a parameterization for the effects of upwelling. The thermodynamic changes in the ocean are constrained to the mixed layer (50-m depth), where modified Ekman dynamics take place. See Zebiak and Cane (1987) or Battisti (1988) for further details.

The original model is a nonlinear anomaly model: solutions are calculated as perturbations from a prescribed basic state, which includes the climatological annual cycle of SST, surface currents, upwelling, and surface winds. For this study, these background fields are taken to be the same as those used by Battisti (1988). Besides linearity, the major difference between the model used in our study and that used in Battisti (or Chen et al. 1997) is in how the model equations are solved. Here, in a manner similar to Hirst (1988) the equations were cast in matrix form (see section 3), with discretization of the equations in the zonal direction and expansion onto the equatorial modes in the meridional direction. The resolution in the *x* direction is 11° longitude; 14 meridional modes are kept. The basic results are insensitive to increasing the *x* resolution or the number of meridional modes that are retained. A full discussion of the equations and the numerical techniques is found in Thompson (1998b).

## 3. The ENSO mode and the spatial and temporal structure of optimal growth

### a. Eigenvalue analysis (Floquet analysis)

In T98a the eigenvectors and singular vectors of the model were found assuming a time-invariant background. Specifically, in T98a the annual cycle was approximated by its average value. This proved to be an adequate simplification to find the general shape and mechanism of the optimal patterns,^{1} but it tells nothing about the seasonal dependence of the singular vectors. In this study the model is generalized to include the annual cycle in the background states.

The basic equations controlling the dynamics can be written in matrix form as

where *T*(*t*) is the discretized SST, and (*t*) are the discretized components of the equatorial Rossby and Kelvin modes (i.e., the ocean dynamics). In general, Eq. (1) has the solution

where exp( ) is the matrix exponential. Because **M** varies in time, the matrix exponential is not constant, and the system does not have the simple eigenvalue decomposition used in T98a.

To solve Eq. (1) numerically, the problem will be discretized in time by holding the matrix operator constant over each month and calculating the linear propagator for that month. These propagator matrices are then multiplied in order to get the propagator matrix between any two times. In other words,

where **M**_{1} is the matrix operator with the January climatological background, **M**_{2} has the February background, etc., and where *t*_{m} is one month. In this model all months are 365.25/12 days long. The linear propagator for 1 yr from January to December is therefore

The eigenvalues of the matrix **R**_{year} are sometimes called the Floquet multipliers of **M**(*t*) (Iooss and Joseph 1980, 126–133). The Floquet multipliers indicate the frequency and growth rate of the eigenvectors of **R**_{year}, which are the modes of the system described by Eq. (1). Recall the definition of eigenvalues

where _{j} is the *j*th eigenvector and *λ*_{j} is the *j*th eigenvalue of **R**_{year}. From Eqs. (1)–(5) it can be seen that

Keeping in mind that both *λ*_{j} and _{j} are complex quantities, the frequency of the mode is determined by mod(*λ*_{j}), which is the angle through which _{j} rotates in 1 yr.

Because this eigenvector is determined by the 1-yr propagator, it indicates only what the fields look like at 1-yr intervals. From just an eigenvalue analysis of **R**_{year}, nothing can be said about the states that the fields go through in getting from to *λ*_{j}_{j}. In fact, a yearly propagator going from July to June,

will have (in general) eigen*vectors* that are different than those of the January-to-December propagator. However, the eigenvalues are the same for both systems. Hence, the growth rate over one year is the same, no matter the start month.

To illustrate, Fig. 1 shows the Niño-3 index for a pure ENSO mode propagating freely, but with the exponential growth suppressed. In Fig. 1, the solid line indicates Niño-3 as calculated by multiplying the current state by the appropriate monthly linear propagators, and the circles show Niño-3 at yearly intervals (at the beginning of January). Note that the system’s state at the yearly intervals can be calculated either by multiplying it by **R**_{year} or by multiplying it by *λ*_{j} and taking the real part of the resulting complex fields. As seen in Fig. 1, the ENSO mode is more complicated with an annual cycle than the simple sinusoidal ENSO mode from T98a (see Fig. 3 from T98a); the Niño-3 signal contains several frequencies and the states through which it passes cannot be fully described by a single complex eigenvector. However, it does have a dominant frequency that is determined by (any one of) the yearly propagators. The period of this dominant frequency will be considered to be the period of the ENSO mode.

Similar to the annually averaged case, the model now produces a single unstable mode that looks like ENSO. This ENSO mode has a period of 2.74 yr and grows by a factor of 1.82 per year. Figure 2 shows the real and imaginary parts of the sea surface temperature anomalies (SSTA) and thermocline fields of the leading eigenvector for the January-to-December linear propagator of Eq. (4). In Fig. 2 the eigenvector has been rotated so that the real part shows the ENSO at the peak of a warm event, while the imaginary part represents the transition phase. The eigenvectors for the other 11 yearly propagators are not identical to Fig. 2 but are still recognizably versions of ENSO.

### b. Singular value analysis

Before singular value analysis can be performed on the system, a suitable norm must be chosen. For this analysis, the following two properties would be desirable. (i) At the initial time, the norm should reflect the probability distribution of the stochastic forcing on the system, or in the case where singular value decomposition (SVD) is being performed to analyze error growth, the desirable norm is picked based on the assumed probability distribution of the initial error. According to Palmer et al. (1998) the preferred norm is the covariance norm (sometimes called the Mahalanobis metric) that is derived from the probability distribution. In the case of spacially uncorrelated forcing, the covariance norm is simply the *L*_{2} norm. In the stochastic models to be developed later in this study, spacially uncorrelated noise will be assumed and so the *L*_{2} norm is appropriate. (ii) At the optimization time, the norm should weigh heavily the mature phase of the ENSO cycle. The objective of the parameter studies carried out in section 4 will be to maximize the growth of ENSO events. The direct method of achieving a high value for ENSO events would be to weight the area of the Niño-3 index heavily in the norm. However, for the ENSO mode in this model, the *L*_{2} norm peaks at the mature phase of the cycle and so no special weighting is required.

The *L*_{2} norm therefore seems to be a reasonable norm, the only difficulty being how the SST and ocean dynamics term should be weighted relative to each other. Referring to Eq. (1), the *L*_{2} norm for the system can be written:

where *N*_{r} and *N*_{T} are the number of points in the discretization of *r* and *T,* respectively, and *w*_{1} and *w*_{2} are constants that weight the ocean dynamics relative to the SST. For this study, *w*_{1} and *w*_{2} were set so that the magnitude of the *r* and *T* portions of Eq. (8) are approximately the same at the mature phase of the ENSO mode. T98a showed that for time periods *τ* greater than 3 months, the shape of the *r* and *T* portions of the singular vectors are approximately independent of the weights used. However, the magnitude of the *r* and *T* portions of the singular vectors are strongly dependent on the weighting. For convenience, the variables were nondimensionalized so that *w*_{1}/*N*^{1/2}_{r} = *w*_{2}/*N*^{1/2}_{T}, and so a simple unweighted *L*_{2} norm can be used in the calculations below. A more complete discussion of the norm used here can be found in appendix F of Thompson (1998b), hereinafter called T98b.

Using the *L*_{2} norm, the conditions that promote optimal growth over time period *τ* can be computed using the SVD of **R**(*τ*), the propagator. Singular value decomposition breaks **R**(*τ*) into three matrices:

where **U** and **V** are orthornormal matrices and **S** is a diagonal matrix. The columns of **V** and **U** are called the right and left singular vectors, respectively. (Hereinafter, the right singular vectors will be simply be called the“singular vectors.”) The diagonal elements of **S** are the singular values. The column of **V** with the largest corresponding singular value is the optimal initial field. Its magnitude grows by a factor equal to its singular value (in time *τ*), and its final shape after time *τ* is given by the corresponding column of **U**. This can be seen by noting that **V** is orthogonal; it has the property that **V**^{T}**V** = **I**, where **I** is the identity matrix. Equation (9) can be rewritten as **RV** = **US** or **R***υ*_{i} = *s*_{i}_{i}, where _{i} and _{i} are the *i*th columns of **V** and **U**, and *s*_{i} is the *i*th diagonal element of **S**. Because **R**(*τ*) is the propagator, it takes the field of shape _{i}, to a field of shape _{i}, and size *s*_{i} at time *τ.* More details about SVD can be found in Noble and Daniel (1988, 338–345).

Figure 3 shows a contour map of all the singular values for the first singular vectors (SV1s) with *τ* between 3 and 18 months, starting at each month of the year. This kind of graph will be called an SV map. From this figure note that there is a local maximum in the transient growth. This maximum occurs at around *τ* = 9 months, starting at the beginning of May. Hereinafter the notation SV1(*m, τ*) will be used to mean the first singular vector for starting month *m* and optimized over the time period *τ.* Note that the maximum value of the singular value is about 12, the same as the maximum singular value for the annually averaged case from T98a. Figure 3 has an maximum on the edge of the graph at (6, 18), which deads to a second local maximum that is off the graph. Because this system is unstable, the singular values grow increasingly large as *τ* → ∞. However, for *τ* larger than the first local maximum, these large singular values are due to the exponential growth of the ENSO mode rather than the initial transient growth. For this reason, the SVD analysis will be limited to *τ* less than 18 months.

The singular vectors produced are very similar to those of the annually averaged case. Figure 4 shows the pattern for SV1(5, 9).^{2} The SST field shows the east–west dipole along the equator and the same north–south dipole in the eastern Pacific, which were characteristic of the SST optimals from T98a and Chen et al. (1997). The thermocline field consists of a trough extending across the basin along the equator, spreading slightly at the western boundary, and also displays the dipole located near the eastern boundary and south of the equator. This basic structure of the thermocline of SV1(5, 9) is also found in the optimals of T98a. The similarity in the patterns indicates that while the cyclic background effects the magnitude of the growth of the singular vectors, the mechanisms of growth are the same as detailed in T98a and Chen et al. (1997).

Figure 5 shows a contour map of the correlation coefficients between SV1(5, 9) and SV1(*m, τ*) for all *m* and *τ.* Over a large area, the shape of the first singular vectors are not particularly sensitive to the time of year or period of growth. For time periods of less than 12 months, the correlation is greater than 0.8 for almost the entire domain. This is consistent with the results of others studies (e.g., Chen et al. 1997; Xue et al. 1997a,b). For time periods greater than 12 months, the correlation drops precipitously near (*m, τ*) = (8, 15). Figure 6 shows SV1(8, 15), and Fig. 7 shows a contour map of the correlation coefficient between SV1(8, 15) and SV1(*m, τ*) for all *m* and *τ.* Comparing Figs. 5 and 7, note that one or the other of these two patterns (Figs. 4 and 6) represent that optimal growth conditions over the entire range of *τ* from 3 to 18 months. In fact, the first singular vectors for all *τ* greater than 3 months are approximately represented by one of these two patterns;the patterns repeat with a 12-month period for *τ* greater than 18 months (not shown).

There is a relationship between the patterns SV1(5, 9) and SV1(8, 15). They each represent (approximately) the second singular vector of the system at the other’s (*m, τ*). That is, SV2(5, 9) ≈ SV1(8, 15) and SV2(8, 15) ≈ SV1(5, 9). The correlation coefficient between SV1(8, 15) and SV2(5, 9) is 0.84, and the correlation between SV1(5, 9) and SV2(8, 15) is 0.92. This is not a coincidence; these two nearly orthogonal patterns correspond (approximately) to the real and imaginary parts of the adjoint from the annually averaged version of the model. (See Fig. 3.1 of T98b for a picture of the adjoint.)

Figure 8 shows the contour maps of the second and third singular values of the system. Comparing Fig. 8b, the third singular values, to Figs. 8a and 3a, the second and first singular values, respectively, it can be seen that the first two singular vectors, as represented nearly everywhere by the two patterns Figs. 4 and 6, should explain most of the transient growth of the system.

### c. The inspiration for a stochastic model

There is an interesting feature of the SV map for the model; appearing in Fig. 3 is a dashed line that marks the locations where the starting month plus the time period *τ* yields 1 January. Along the lines of the fixed start month, the singular values tend to peak close to this dashed line. In other words, whatever time of year the system is perturbed, it has a tendency to peak during the (northern) winter months.

The real ENSO system of course has a tendency to peak during the winter months (e.g., Rasmusson and Carpenter 1982). For instance, Fig. 9 shows the evolution of the Niño-3 index for the six largest warm ENSO events from the Consolidated Ocean Atmosphere Dataset (COADS) data between 1950 and 1992. In all six of these events, the ENSO matures around December or January. (Not shown are two ENSO warm events that extend over two consecutive winters.) A statistic that quantifies this tendency to peak during winter is the monthly standard deviation, displayed in Fig. 10 for the Niño-3 index from the COADS data.

In the model, the tendency for transient growth (growth from the singular vectors) to peak during the winter is even more pronounced than indicated by Fig. 3. Recall from T98a that when an SV1(*m, τ*) evolves it does not necessarily peak at time *τ.* For example, Fig. 11a shows the growth of SV1(May, *τ*) for various *τ.* The circles indicate the value of the *L*_{2} norm at time *τ* for each curve. As seen in this figure, each singular vector produces a mature ENSO event in December, January, or February, regardless of the value of *τ. *Figure 11b shows the same diagram for SV1(Sep, *τ*). Here the results are similar, except for the larger *τ*’s where the optimals start to peak in the *following* winter. Figure 11c shows the same set of curves (with the circles removed) for SV1(*m, τ*) over *τ* = 3, 6, 9, 12, and 15 months, starting every other month. This graph clearly shows that the evolving singular vectors tend to peak during the winter, regardless of the starting month or optimization period.

Comparing Fig. 11 with Fig. 10, it appears that the model peaks about 2 months late relative to the real system. At least part of this discrepancy is due to Fig. 10 using the Niño-3 index, while Fig. 11 uses the *L*_{2} norm. The Niño-3 index of the model ENSO peaks about 1 month earlier than the *L*_{2} norm. The rest may be due either to deficiencies in the model or to structure in the noise that drives the real system.

## 4. Stochastic model development

### a. Transient growth versus modal growth

Because the class of models being studied here is highly parameterized, there are many coefficients within the basic equations that are known only approximately. The objective of the parameter studies in this section will be to look for those parameters that can be altered (within realistic ranges) to damp the modal growth of the system, while preserving the transient growth. For this purpose the transient growth-to-modal growth ratio (TMR) will be defined as

The modal growth used will be the growth of the complex vector that represents the ENSO mode. The growth of the complex ENSO vector is used because this growth rate is constant in time, unlike the real part of the ENSO mode, which grows at different rates over different phases of the ENSO cycle (see T98a). The TMR is not an attempt to quantify the relative importance of the two types of growth for the system, but rather a convenient benchmark for the parameter studies.

The TMR for the model with the B88 parameters is 6.97, which will be considered the baseline. For all of the parameter studies, the TMR will not be listed directly, but rather the change in the TMR from the baseline will be shown. From the definition of the TMR it would seem that it has been assumed that the maximum growth in the model is always at (*m, τ*) = (5, 9) = (May, 9 months). This assumption has been made to avoid calculating all the SVs for each point of the parameter studies. This assumption has been shown to be reasonable by checking the full SV map at the extremes of the parameter ranges; the maximum optimal remains consistent at (5, 9) to within the resolution of the SV maps.

### b. Parameter studies

The most straightforward way to suppress the growth in the model is to reduce the coupling strength between the ocean and atmospheric components of the model. Specifically, the constant *γ* in Eq. (A.2) of T98a controls the strength of the atmospheric heating (and therefore the amplitude of the atmospheric reaction) to a prescribed SST anomaly. Alternatively, the drag coefficient *C*_{D} from Eq. (A.6) of T98a controls the strength of the ocean response to a given wind anomaly. Changing either of these coefficients has the same effect on the model stability.

Table 1 shows the period and growth of the ENSO mode, the (5, 9) optimals growth, and the TMR as the coupling is reduced from its nominal value. As might be expected, both the strength of the optimal and the growth rate of the ENSO mode are reduced. By interpolation, it is found that the ENSO mode becomes neutral at a coupling multiplier of about 0.76. At this neutral point the TMR is reduced by a factor of about 0.92. So reducing the coupling, while neutralizing the model as desired, does not increase the TMR, but reduces it slightly. This result is not surprising, since reducing the coupling damps the model, but without targeting any particular mechanism.

One characteristic of the transient growth versus modal growth difference that should be exploitable is the timescale difference. The *τ* optimal experiences its growth over a 9-month period, which is too short a period of time for reflection from the western boundary to have much effect. Both the B88 model and the ZC model use a 100% reflective boundary, yet Clarke (1991) and duPenhoat and Cane (1991) have estimated using linear theory that the western boundary reflection (WBR) coefficient is in the neighborhood of 0.8 for interannual Rossby signals. Mantua and Battisti (1994) estimate that the WBR may be as low as 0.7. Since duPenhoat and Cane (1991) show that the eastern boundary reflection coefficient is close to 1.0, only the western boundary coefficient was varied. Table 2 shows the effect on the model of reducing the reflection efficiency of the western boundary. As expected, the ENSO growth rate is reduced, while the SV remains essentially unchanged.

A second modification to the ocean equations that should affect the TMR positively would be an increase in the mechanical damping of the ocean. The damping is set to 2.5 yr in the ZC model and the B88 model, but this damping time is probably too long. Picaut et al. (1993) compares the output of an intermediate ocean model driven by observational winds to the tide gauge, moorings, and Geosat data and finds a damping time of between 6 months and 1 yr for the frequencies of interest to this study. While increasing the damping will affect both types of growth, it should affect modal growth more owing to its longer timescale and since it is made up of equatorial waves that propagate back and forth across the basin.

Table 3 shows the effect of the increased mechanical damping on the model in the range from 2.5 yr to 6 months. As expected, the ENSO mode is damped more quickly than the optimal, with the TMR increasing by 77% for the 6-month damping coefficient. The ENSO mode becomes neutral when *a*^{−1} = 0.703 yr, at which point the transient growth still retains 77% of its strength.

### c. Stochastically driven ENSO models

Using the parameter studies from section 4b, four alternative ENSO models were built. The first model will be the nominal model, with all parameters left at their B88 model values, except that the coupling between the ocean and atmosphere is reduced until the model is slightly damped. This model will be called N.97: “N” for “nominal” and “.97” because 0.97 is the yearly growth of the ENSO mode. Remember that these models are to be driven stochastically and so the B88 model parameters could not be used (since B88 model is unstable). For this nominal model the TMR was kept nearly the same as in the B88 model.

The second model was built to match the ENSO growth rate of N.97, but the transient growth will be kept as large as possible. To do this, the WBR, which preserves the transient growth most effectively, is set to a reasonable value of 0.8. The ocean mechanical damping is then increased until the ENSO growth reaches 0.97 per year; this damping turns out to be 10 months. This model will be called T.97 where “T” stands for“transient growth.”

The last two models are designed to explore the region where the ENSO mode is not nearly neutral, but either moderately or highly damped. Both models are built to preserve the transient growth. The third model, designated T.80, has an ENSO growth rate of 0.8 per year. It was built by reducing the WBR to a value of 0.7, and increasing the ocean mechanical damping to a realistic value of 8.7 months. The fourth model has an ENSO growth of 0.6 per year. This was achieved by setting the WBR to 0.7, and increasing ocean damping until the ENSO mode achieved the desired growth rate. The required ocean damping turned out to be 6.5 months—a high damping, but still plausible.

Table 4 shows the parameter values for each model, and the basic characteristics of each model such as the ENSO growth rate and the model singular values. The shapes of the ENSO modes and singular vectors are very similar to those for the linear B88 model (see Figs. 4 and 6) and will not be shown.

The SV maps for all four models are presented in Fig. 12. All four of the alternate models produce SV maps with the same characteristic ridge along the winter end times as was found for the model in its B88 model configuration (Fig. 3), and therefore the hypothesis that time-invariant stochastic forcing will lead to model variance with an annual cycle is unchanged by these parameter modifications. One important difference between the SV maps of Figs. 12 and 3 is that in the damped models, the local maximum at (5, 9) is an absolute maximum. This is because the models are damped, and therefore the transient growth is truly transient. In the unstable parameter regime, the singular vectors evolve into the ENSO mode and so grow without bound.

It should be noted that all these models were built with the primary design criteria of manipulating the transient and ENSO growth rates. Unfortunately, it is not possible to hold all other characteristics constant. Specifically, the period of the ENSO mode ranges from 2.9 yr in N.97 up to 4.3 yr in T.60. An examination of the parameter studies (Tables 1, 2, and 3) shows that all the parameters that were varied in producing the models above have the characteristic of lengthening the period as they lower the growth rate. In fact, none of the parameters had the property of shortening the period while both lowering ENSO growth and enhancing the TMR. Therefore, it was not possible to control the ENSO period using these parameters.

### d. Stochastic simulations

In this section, the four alternate ENSO models developed in section 4c were integrated while undergoing stochastic forcing. The analysis of the behavior and statistics of the models is deferred to the companion paper. Figure 13 shows the Niño-3 index of the COADS data along with a 42-yr sample from each of the four models. The model samples come from a 500-yr stochastic simulation where the forcing consists of monthly perturbations to the state variables, *r* and *T.* The random perturbations were time invariant and spacially uncorrelated (i.e., white in space and time). Each of the stochastic simulations was made with the same type of random forcing, except the amplitude of the forcing was changed from model to model so that the annual averaged variance of the model output is the same as the annual average variance of the 42 yr of COADS observations. Table 5 shows the amplitude of the (normally distributed) forcing for both the state variables where the equatorial ocean mode amplitudes *r* are given in terms of the thermocline perturbation.

A visual examination of the model output reveals that the nearly neural models, N.97 and T.97, are too regular as compared with the COADS data. However, the moderately damped models, T.80 and T.60, which use realistic parameter values for the western boundary reflection efficiency and ocean mechanical damping, display irregular behavior more like the real world. In the companion paper, a more complete comparison will be made between the statistics of the various models and the statistics of the observed system, including an assessment of how dependent the coupled system may be on the exact form of the stochastic forcing.

## 5. Discussion and conclusions

Floquet analysis is performed on a linearized variant of the Zebiak–Cane model of the coupled atmosphere–ocean system in the tropical Pacific. Consistent with the results of Battisti and Hirst (1989), the leading mode of the coupled system is very similar to the observed ENSO mode, and the physics of the mode is described by the so-called delayed oscillator physics. With the parameters that are traditionally used in the intermediate class of coupled atmosphere–ocean models [e.g., the Zebiak and Cane (1987) and the Battisti (1988) models], the system is linearly unstable (see, e.g., Battisti and Hirst 1989; and T98a).

Singular vector analysis is also performed on the linearized model and shows that for any start month and for any time period of optimization, the singular vector is shaped like one of two nearly orthogonal patterns. These two patterns correspond approximately to the real and imaginary parts of the adjoint of the ENSO mode for the time-invariant basic-state version of the system that was calculated in Thompson (1998a). However, the growth of the first singular vector, applied as an initial condition to the model, is sensitive to the start month and optimization period. This optimal transient growth has a strong annual cycle due to the annual cycle in the basic-state climatology. Maximum growth occurs when perturbations are allowed to evolve over boreal summer, with perturbations reaching their peak at the end of the year, regardless of the month in which the perturbation is applied. These results are consistent with those of Chen et al. (1997). This tendency for perturbations to peak in the boreal winter implies that if stochastic perturbations are applied to the model, the resulting ENSO events will be phase locked to the annual cycle, as per the observations. This phase locking would not require there to be an annual cycle in the stochastic forcing of the model, nor would it require nonlinearities in the system.

In order to convert this study’s linear model into a stochastically forced model, it is necessary that the model be moved into a stable regime. Two model parameters seemed particularly appropriate for modification to this end: the ocean mechanical damping *a*^{−1}, and the western boundary reflection efficiency WBR. Note that at the time the Zebiak–Cane model was constructed, there was insufficient data to constrain the values of the *a*^{−1} or the WBR. The traditional values of these two parameters in the Zebiak and Cane (1987) and Battisti (1988) models are as follows: *a*^{−1} = 2.5 yr, WBR = 1. However, the recent results from Picaut et al. (1993) suggest the value of *a*^{−1} for annual and longer timescales is about 6–14 months. Mantua and Battisti (1994) performed a comparison study between the model and sea level observations. This study suggests that the efficiency of reflection of annual and longer timescale Rossby waves off the Indonesian archipelago is about 0.75, which is consistent with the theoretical results of du Penhoat and Cane (1991) and Clarke (1991).

Parameter studies carried out in section 4 show that moving *a*^{−1} and WBR into these more realistic regimes stabilizes the model while preserving the system sensitivity to perturbations. Changes are made in the values of various parameters to construct four candidate coupled models that are asymptotically stable but fall into different sensitivity regimes: nearly neutral but with nominal transient growth (N.97), nearly neutral with enhanced transient growth (T.97), moderately stable with enhanced transient growth (T.80), and very stable with enhanced transient growth (T.60). The nominal candidate model was produced by reducing the atmosphere-to-ocean coupling until the system stabilized, (reducing the transient growth by approximately the same amount as the modal growth), whereas the enhanced transient growth models were produced by altering *a*^{−1} and the WBR. Analysis of the candidate models shows that the leading mode (the ENSO mode) and the optimal perturbation patterns are very similar to those from the coupled model with the standard parameter set. In addition, the tendency of optimal perturbations to peak during the boreal winter is also preserved. The four candidate models were then run with external stochastic forcing, and it was noted that the two more highly damped models, T.80 and T.60, produced ENSO events irregularly in a manner similar to observations.

These results suggest that several of the basic characteristics of these coupled models are robust to a wide range of values of the parameters and to whether the parameter selection renders the coupled system asymptotically stable or unstable. In particular, the following characteristics of the coupled system are robust: the shape and dominance of the ENSO mode, the tendency of the ENSO events to peak at the end of the calendar year, the shapes of the principal singular vector (optimal initial conditions), and the tendency of these optimal initial conditions to peak the end of the year (regardless of when the perturbation is applied). Therefore, it seems that intermediate coupled models such as the Zebiak and Cane model and its variants are capable of producing ENSO-like events under a wide range of parameter space, even if the system is linearly stable, providing that some external stochastic forcing is present.

Interestingly, many of the coupled atmosphere–ocean models that employ atmosphere models that are constructed to be in equilibrium with the SST (e.g., the Gill-type models or empirical atmosphere models) require unrealistically large values of the coupling strength to sustain ENSO-like variability. Zebiak (1984) argued the large drag coefficient helped to compensate for the lack of synoptic activity in the these equilibrium atmosphere models. In effect, Zebiak (1984) argued that the details of the unresolved atmospheric variability are not important to the physics of the model ENSO events, and the results from the paper demonstrate this to be the case.

The research community is beginning to explore the hypothesis that the irregularity of warm and cold ENSO events may result from stochastic forcing. Models at a variety of levels have been used. Specifically, simple stochastic ENSO models have been built by Burgers (1999) and by Wang et al. (1999); an intermediate, stochastic, Markov model was built by Penland and Sardeshmukh (1995); intermediate, stochastic, physical models are used by Moore and Kleeman (1999) and by Flugel and Chang (1999); a stochastically forced hybrid coupled model is used by Eckert and Latif (1996); and Chang et al. (1996) use a fully coupled GCM (in addition to hybrid and intermediate models). All these studies support the plausibility that ENSO irregularity, and possibly ENSO variability, are due to external stochastic forcing.

In part two of this study, evidence will be provided that the more realistic sets of parameters (those used in the T.80 or T.60 models) also afford model responses that are in good agreement with observations for a number of statistics. It will be shown that the observed level of variance in the coupled system can be produced by a modest (realistic) amplitude of stochastic forcing that is temporally white and spatially uncorrelated, and that the pattern of the first singular vector is responsible for the bulk of the system variability. The limits of predictability (the potential predictability) of the coupled system will be examined, and are more modest than was previously estimated using the traditional deterministic models.

## Acknowledgments

We would like to thank the many denizens of the Climate Palace for their support during the years of research that lead to this paper. We also thank Cecile Penland and two anonymous reviewers for many helpful comments on the original manuscript. This work has been supported by a grant from the NOAA Office of Global Programs to the Hayes Center at the University of Washington.

## REFERENCES

## Footnotes

*Corresponding author address:* David Battisti, JISAO, University of Washington, Box 354235, Seattle, WA 98195-4235.

Email: david@atmos.washington.edu

^{1}

In Thompson (1998a, b) the term “*τ* optimal” is used to mean the “first right singular vector for time period *τ.*” In this paper the“optimal” term will be used only in reference to the results of this previous work.

^{2}

A 1-2-1 meridional filter was used on the SV patterns that suffer from grid noise. The nature of this grid noise is explained in appendix B of T98a. The filter degrades the performance of the singular vectors (SVs) by about 10%; that is, the filtered SVs shown will grow to about 90% of the amplitude of the unfiltered SVs at the optimization time *τ.*