Characterizing Enso with Nonlinear Dynamics

This thesis analyzes three observational data sets related to the El Nino/Southern Oscillation (ENSO) phenomenon to determine if E SO can be considered as a low-dimensional chaotic system. In order to test this hypothesis, I apply Barahona and Poon's [1996] method for detecting nonlinear determinism in short , noisecontaminated time series, and calculate the correlation dimension, De. A slightly modified version of the Vallis [1986] model is used to provide a context for interpreting the results, and to validate the computations done on the observational data. When applied to observational ENSO data, the Barahona and Poon algorithm indicates that low-order nonlinear mapping functions have predictive power which is not significantly different from linear models. In contrast, when the algorithm is applied to correspondingly sampled data from the Vallis model , the algorithm shows the presence of nonlinear determinism, even when these data are strongly contaminated with noise. There is a weak convergence of the correlation dimension in the observational data to a De value between 8 and 10. This indicates that the phase-space dimension of E SO is at least 8. These results suggest that either (1) ENSO is not governed by low-dimensional nonlinear dynamics, or (2) noise related to local physical processes overwhelms the E SO chaotic signal in the observational data.

Vlll El iiio/Southern Oscillation (E SO) is an interannual climatological disturbance centered on the tropical Pacific, and it has global effects and relevance. In this thesis , I investigate the possibility that variability associated with E SO can be attributed to low-dimensional, dissipative, chaotic dynamics.

ENSO
The term "Southern Oscillation" describes a fluctuation in the Southern Hemisphere's atmospheric pressure over the tropical Pacific and Indian Oceans [ Philander, 1990]. When surface pressure is high over the eastern Pacific, it tends to be low over the eastern Indian Ocean and vice versa. The two regions for which this negative correlation is largest are referred to as the "centers of action." "El Nino" refers to an appearance of anomalously warm surface waters off the coasts of Peru and Ecuador. El iiio events, also referred to as E SO events, are coincident with the Southern Oscillation state being such that the atmospheric pressure is low in the east and high in the west. The time interval between successive ENSO events averages about 4 years , but these intervals range from one to eight years in length. E SO and the Sou them Oscillation are considered to be part of a regional cycle with global effects.
Poor prediction of ENSO events has an important effect on the world economy [Philander, 1990]. Conventional models and prediction schemes did not anticipate the recent large E SO event [ CPC, 1997]. If E SO can be shown to behave in a manner consistent with low-dimensional chaos, then it is possible that chaotic prediction techniques can be fruitfully employed to forecast future ENSO events.

Chaos
Chaos refers to the seemingly erratic time-dependent behavior that can occur in otherwise simple, deterministic systems. To be termed chaotic, a bounded time series must not be asymtotically periodic, and it must have at least one positive Lyapunov exponent. 1 For this work, I focus on low-dimensional dissipative chaos, which is operationally defined as a bounded chaotic system having a dimension significantly less than ten. Nonlinearity is a necessary -but not sufficientcondition for chaotic dynamics.
If E SO is chaotic as described above, what would that mean? In reality, E SO is governed by a vast number of equations describing the Newtonian physics which all of the molecules in the atmosphere and the ocean must obey. Even with the 1 The Lyapunov exponents quantify the exponential divergence (in time) of initially close trajectories. If the vector distance between two nearby trajectories goes as s = s-OeLt , then the eigenvalues of the matrix L are the Lyapunov exponents for that system. 2 continuum hypothesis , the effective scale of the smallest fluid motions in the ocean is of order 1 cm, and there are 10 24 such centimeter-size parcels in the ocean. If E SO is governed by low-dimensional nonlinear dynamics , only a much smaller number of equations is important in interpreting and predicting the dynamics. Thus, there would exist a simple model capturing all of the important physics of the system. evertheless, knowing this would not in itself tell us anything about the feasibility of finding that model , or about the uniqueness of that model. Two methods were used to examine the degree of nonlinearity in ENSO data sets. Originally, I also planned to calculate the Lyapunov exponents for the E SO attractor, however , that calculation relies on knowledge of an upper bound for the fractal dimension of the attractor. The results from the two methods cited in this study (Sections 4 and 5) indicate that Lyapunov exponents are essentially incalculable from the available data sets.
The first method used is taken from Barahona and Poon [1996]. This method finds appropriate linear and nonlinear iterated map models for the given data set. It then calculates a "goodness of fit" metric for those linear and nonlinear models. If a nonlinear model has significantly greater predictive power than the linear model , then the data set is said to have underlying nonlinear deterministic dynamics. This method is described in greater detail in Section 4.
The second computation done for this work is the calculation of the correlation dimension [ Grassberger and Procaccia, 1983] for the climate time series. If time series data are taken from all of the phase variables of a chaotic system, they will lie on a chaotic orbit, a "strange attractor." The fractal dimension of this attractor 3 can be calculated on a subset of those variables (or a linear combination of that subset) to reconstruct an analog to the original phase space. The theory of timedelay embedding says that even a single time series can be used in this manner, and the correlation dimension of the reconstructed orbit constitutes a lower bound for the actual number of independent variables controlling the time evolution of that system. For high enough embedding dimensions , the correlation integrals should be approximately linear in log-log space. The slope of that line is the correlation dimension. The details of this method are described in Section 5.
Section 6 summarizes the results and concludes the main part of this thesis. 4

SOI
The Southern Oscillation Index (SOI) is a standardized, monthly-averaged sealevel atmospheric pressure (SLP) difference between Tahiti and Darwin, Australia.
Strongly negative SOI values are associated with ENSO events. Tahiti and Darwin are chosen because they are near the centers of action for the Southern Oscillation [Philander, 1990] and because there is a large enough quantity of high-qualitiy data available to make a consistent time series. The version of the SOI used in this study is the Climate Prediction Center (CPC) standardized SOI. The data set starts September, 1932 and it continues through April , 1998. The CPC archives the SOI from 1880 through the present, but there are several large data gaps during 1880-1932, so the set used here starts at the end of the last such gap. Thus, 788 monthly-averaged data points were used for this study. This data set is referred to as the "SOI". For more details, see Appendix A.

Daily T-D
The short length of the SOI limits the accuracy of the computational procedures used in this study. However, lengthy historical measurements are hard to come by, and in any event would not extend significantly farther back into the past than the SOI does. The only other option, then, is to obtain data with a higher time sampling rate than the SOI's monthly rate. To address this need , a daily analog to the SOI was computed for these analyses from a spatially gridded data set.
Although I obtained the base data set through tape archives [Walters, 1997] at the ational Center for Atmospheric Research , it was originally produced at the Australian National Meteorological Research Centre [Seaman, et. al., 1995]. It employs time and spatial averaging to fit all available SLP data (from 4/24/1972 through 12/31/1992) to a 47 x 47 grid extending across the southern Pacific region .
It lists SLP values daily for the grid points at Greenwich Mean Time (GMT) 23:00 until 4/8/73, and twice a day at 11:00 and 23:00 GMT thereafter.
The spatial bins used for this set are about 350 km across (at 17°8) . Because the size of the area represented by a grid point is so large, only one grid point was used to acquire the SLP at Tahiti and one at Darwin. The grid points used are: (9 ,33) centered at 16.4472°8, 149.036°W for Tahiti, and (12 ,10) centered at 13.5319°8, 130.601°E for Darwin.
After extracting the appropriate data from the archive , a daily composite data set was calculated. For days with two available data points, an average of the two points was used. For days with only one data point available, that point was used .
There were 69 missing days with no individual gap being more than 9 days long, so a cubic spline interpolation was calculated to fill in for those days with zero available data points. Finally, the daily SLP value for Darwin was subtracted from that for Tahiti. Altogether, 7557 such "pseudo-SOI" values were computed. These data will be referred to as the "T-D" (Tahiti minus Darwin) data set.

6
To compliment the other two data sets, a daily SLP time series at Darwin was constructed. The Darwin data set is used separately because it is a directly measured data set, and because it has the largest number of available daily data points (16895). However, the Darwin SLP series has a large seasonal component, which could interfere with the calculation of ENSO statistical and chaotic properties.
Eight-times-daily SLP data for Darwin were obtained directly from the Australian Bureau of Meteorology. This data set spans the time period from July 1, 1951 , 2:00 local time, through September 30, 1997, 21:00 local time. To reduce the impact of the diurnal cycle, a daily average was taken by directly averaging all of the available data points between O:OOZ and 23:59Z. Except for one day, every day had at least 3 measurements. That one day had no data, and the SLP for that one day was taken to be the average of the SLP measurements for the previous and following days. Finally, 1000 mbar was subtracted from each SLP value. This processing resulted in a (residual) SLP series 16895 days long. The data set will be referred to as the "Darwin SLP" data.

Random numbers
For the analysis techniques used in this work , it is helpful to compare the results with an equivalent random data set. Therefore, for each data set, a surrogate data set was constructed following the method described in Tsonis and Elsner [1993]. The Fourier transform of the data set was computed, and then the frequency amplitudes 7 were multiplied by a random phase factor , ei<f>, where <P is uniformly distributed in the interval 0 to 2n. The inverse transform was then used as a surrogate series.
This process produces a random series with a mean, variance, and autocorrelation function similar to the original series. Although in theory the inverse transform of the newly randomized data should be real if </;(!) = <P( -!) , in practice, numerical errors induced a significant imaginary part to the inverse. Thus, for this analysis two additional computational steps were performed. First, the imaginary part of the inverse transform was dropped (the real part was used). Secondly, the standard deviation, a , was manually adjusted to match the a of the original time series.
Inspection of the random time series obtained in this way indicated that their autocorrelation functions match those of the original data sets well. Figure 1 shows the autocorrelation functions for the SOI and its corresponding random surrogate data set.   Vallis [1986] proposed a simple model for ENSO which captures many of its qualitative features. The Vallis equations are

.1 Model Physics
The free variables are Te (upper-ocean temperature in the East) , Tw (upper-ocean temperature in the West) , and u (eastward velocity of the upper-ocean). The other parameters, A , B , C , T, T* , and u* are predetermined. A is the Newtonian ternperature damping constant, B is a coupling constant, C is the frictional damping constant, T is the mean deep-ocean temperature, T* is the relaxation temperature of the ocean, and u* is the mean trade wind velocity. In this research , A , B , C , and f' are taken to be the values suggested by Vallis [1986] (see table below) which are chosen to correspond to the real ocean dynamics, and to make one time unit correspond to about one month , but T* and u* are adjusted to match the observed E SO dynamics more closely. Figure 2 shows 2000 months of Vallis model data. This figure , as with all Vallis model data sets used in this thesis , was calculated using a fourth-order Runge-Kutta integration with 6t = 0.05 month. The model output with this time interval is refered to as the Vallis "daily" data set.

Seasonal Forcing
Because the tropical Pacific Ocean is strongly influenced by seasonal forcing , for this research a seasonal forcing is included. In his original paper, Vallis [1986] suggested that a seasonal forcing could be added to the model by replacing the constant parameter, u* , with u*(t) = u*(l + 3sinwt) , where w = 1~; ar· However, this level of seasonal variation implies that the trade winds reverse for almost half of the year. According to Philander [1990], the real trade winds only reverse for about two months in the spring. In order to make the seasonal forcing more realistic , this research uses u* ( t) = u* (1 +a sin wt) and a is a constant amplitude chosen to make the trade wind reversal occur for only two months in the spring (see Figure 3). Using Vallis's value u* = -0.45 ms-1 and taking t = 0 to be the beginning of the year, this requires that a = -1. l.

11
Both the computational methods used here assume that only one time series is available, but the Vallis model has three independent phase variables. To make comparisons between the model and the data, one needs to choose which of the phase variables to use as a proxy to the data. I decided to use the phase variable which showed the clearest difference between the model "ENSO" and "normal" states. Figure 4 shows histogram plots of the variable values when they were at relative maxima. These were calculated for T* values which ranged from 2° to 20°C (with a step size of .01°C). At each T* value, the model was iterated to produce 120,000 months of data and all of the model runs were combined to produce There are two distinct groups in the histogram for u corresponding to the ENSO and non-E SO conditions. Hence, it is easy to define an ENSO criterion in terms of the relative maximum of u . I have chosen the relative maxima of u being greater than 1 to be my criterion for defining an E SO condition, but any value between .9825 and 1.825 would have given identical results.

Justification for T* change
In order to make the frequency of ENSO event occurrences in the model match that of the SOI, the model T * parameter was adjusted.
First, the SOI E SO event rate was calculated as follows. The SOI data set was smoothed using a 12-month-period low-pass Butterworth filter (applied in both the forward and reverse directions). Then as with Figure 4 , a histogram of the relative 12 SOI minima was plotted (see Figure 5) to find a suitable threshold SOI value for determining E SO events. The pronounced minimum at SOI= -0.9 suggests this as an appropriate threshold value. Finally, the number of SOI relative minima below this threshold were considered to be ENSO events (with a threshold value of -0.9, there were 15 ENSO events) , and the ENSO rate was calculated by dividing the number of SOI relative minima by the total number of years in the data set (788 months= 65.67 years). This resulted in an SOI ENSO rate of 0.2284 year-1 . Second, the model E SO rate was calculated as a function of T*. This is shown in Figure 6. Again , a model E SO event occurs (by definition) when a relative maximum of u is greater than 1. The number of relative maxima which were ENSO events is divided by the total number of years (12 model months) to calculate the model ENSO rate for each T* parameter value between 2° and 20°C in increments of 0.01°C. The model ENSO rate matches that of the SOI when T* is 7.32°C.   18 Barahona and Poon [1996], describe a method for detecting nonlinear determinism in short, noisy time series. That algorithm is used here and referred to as the "B/P algorithm" or "B/P analysis" . Barahona and Poon [1996] use the discrete Volterra-Weiner-Korenberg (VWK) series as an iterated mapping function , where Yn = Yn(Yn -1, Yn-2, ... , Yn-1t) · So,

Description of B/P algorithm
(1) K is the memory of the system (this corresponds to an embedding dimension) , d is the degree of the polynomial, and M = (K + d)!/(d!K!) is the total number of terms in the complete expansion. The coefficients, am , are then calculated to find a best-fit VWK mapping function to the data. The parameters K and d are chosen to minimize Akaike 's information criterion, where r :::; M is the number of mapping function terms used to make short-term predictions on the data, and E is the standard one-step-ahead prediction error of (1 , 1) -complete 3 ao + a1Yn-1 + a2Yn-2 (2 , 1) -complete 4 ao + a1Yn-1 + a2Yn-2 + a3y;_l (2 , 2) -truncated 5 ao + a1Yn-1 + a2Yn-2 + a3y;_1 + a3Yn-1Yn-2 (2 , 2) -truncated 6 ao + a1Yn-1 + a2Yn-2 + a3y;_ 1 + a3Yn-1Yn-2 + a3y; _2 (2 , 2) -complete  As is apparent from Figure 9, the linear and nonlinear mapping functions , and the mapping function based on a random surrogate to the SOI, have similar predictive capabilities. This indicates that E SO is not driven by a nonlinear, deterministic dynamical system with a small number of degrees of freedom, or that small-scale phenomena (random noise) overwhelm ENSO, or that monthly averaging of the climate data (from compiling the SOI) obscures the results. To help distinguish between these possible explanations I use the algorithm on Vallis model data.

Vallis Model
The Vallis model (as described in Section 3) was used to generate a monthly-sampled data set. The length of the data set is chosen to correspond to the length of the

Noise
Surrogate noise with standard deviations , a , up to 100% of the signal standard deviation was added to the model data before the B/P algorithm was applied (Figure 11). When the noise strength was 60% or lower, then the system still displayed identifiable nonlinear determinism. The surrogate C ( r) functions in these plots are 22 surrogates of the noisy Vallis model data. Adding the noise to the Vallis daily data set before forming monthly averages produced essentially the same results. So, if the proposed E SO climate attractor is similar in essence to the Vallis model , the fact that nonlinear determinism was not detected in the observational data puts a lower limit of 60% on the noise level in the SOI data set. In other words , either the SOI is not characterized by low-dimensional nonlinear dynamics, or the noise in the data set is greater than 60%.

Smoothing
The monthly-sampled Vallis model data set described above with and without added noise (a= 100% of the signal a) was smoothed with a low-pass Butterworth filter (run forwards and backwards) before applying the B/P algorithm. The B/P analysis failed to detect nonlinearity for all possible period cutoffs greater than the yquist period (2 months) , thus smoothing cannot be reliably used to recover nonlinear determinism. This was observed in both smoothed pure data sets ( Figure 12) and in smoothed noisy data sets ( Figure 13). One possible explanation for this is that a smoothing algorithm increases interdependence among neighboring data points, and this interdependence is linear in nature. Thus smoothing enhances the predictive capabilities of linear mapping functions.

7 Sub-sampled Vallis model
A Vallis model data set was constructed to correspond in sampling interval and length with the T-D data set subsampled to 1 data point per 7 days. The data set was subsampled because the B/P algorithm does not work well with oversampled data. By trying various sampling intervals, I found that 1 point per 7 "days" seemed to be the optimal interval in resolving the nonlinear determinism of the Vallis model data when limited to 7557 daily points. 1079 such data points were analyzed ( Figure 14) .

Sub-sampled T-D
The B/P algorithm was applied to a data set consisting of 1079 points from the T-D data subsampled to 1 data point per 7 days. The C(r) function for nonlinear mapping functions was not significantly different than the C ( r) function for the linear mapping functions or for mapping functions predicting the surrogate data set of the weekly sampled T-D ( Figure 15). Thus, the T-D data set is either not characterized by low-order nonlinearity, or there is large localized noise obscuring the low-order nonlinear signal. For this data set, small-scale physical processes have been averaged out by the spatial binning process, so dynamical noise attributed to localized non-ENSO dynamics is not a likely explanation for the null result.
Unless the binning process itself destroys the nonlinear signal, this indicates that the Southern Oscillation is a physical system with many significant variables, which this analysis cannot distinguish from "random. " 24

Bi-weekly-subsampled Darwin SLP
The daily-averaged SLP values from Darwin, Australia were subsampled on a biweekly time scale. The resultant time series has 1207 points, and the parameter values used here were r;, = 7, and d = 3. The B/P analysis was then applied, and the C(r) for the linear (and surrogate data) mapping functions were smaller than the c ( r) for the nonlinear mapping function (Figure 16). Thus nonlinear determinism is not found in this data set. The parameter values used here are "' = 7 and d = 3.

35
One of the earliest methods for characterizing a chaotic time series is the Grass berger and Procaccia [1986] method of calculating the correlation dimension, De, from the slope in log-log space of the correlation integral.

Description of Grassberger-Procaccia method
Consider a dynamical system with state vector, X, which evolves in ad-dimensional phase space. The correlation integral, C(l) , is defined as where N is the number of vector pairs, Xi , Xj , used in the calculation. The distance function can be any vector norm; in this study I use the Euclidian norm unless otherwise specified. If or equivalently logC(l) , . __, vlogl , then vis called the correlation dimension , De.
For instance, consider a line (a one-dimensional object), in an n-dimensional phase space. Equation 4 says that as the radius , l , of an n-dimensional hypersphere grows, the number of points encompassed by that hypersphere will grow as l1.
When working with experimental data, one does not often have access to all of the components of the state vector, X. Therefore, Grassberger and Procaccia [1986] suggest using the method of time-delay embedding to apply their method to time 36 series consisting of a single variable.
First one constructs a set of d-dimensional vectors from the original time series, x(t) as follows: (5) 7 is a fixed time interval, and d needs to be larger than the attractor's dimension.
Furthermore, for a bounded data set of non-infinite length, C(l) -+ 1 as l -+ oo (eventually, all the data are encompassed by a hypersphere of radius l). Also , if the time series is finite , there is some minimum distance between the embedded phase space vectors; thus for all l less than that minimum distance , C(l) = 0. So for real data sets, equation ( 4) only holds for a finite range of l values. This range of appropriate l values is referred to as the "scaling region".
For this work, I calculate correlation integral derivatives , dC(d)(l)/dlogl , for embedding dimensions , d, ranging from 2 to 14. A 0.41 log-unit width box average is applied to smooth the correlation integral derivatives. The scaling region is determined to be the region , in log l space, where the correlation derivatives converge and are horizontal for large enough d. The derivative value in the scaling region is the correlation dimension , but the scaling region can be narrow.

Vallis Model
The Vallis model (as described in Section 3) is used to substantiate the correlation integral calculations on the data sets. Although a similarly long data set is unavailable for comparison, the "true" correlation dimension for the model is cal-37 culated from 50,000 monthly sampled points. From the correlation derivative plot  (Figure 18, right) , so monthly averaging, similar to that done in calculating the SOI, does not significantly affect the algorithm. However, Figure 18 also shows that the De may be underpredicted in these cases. Note that the slopes appear to be about 2-3, rather than 3-4 as in Figure 17.
A twenty-times-monthly, "daily", Vallis model data set was generated to test the plausibility of doing correlation integral calculations on the daily T-D data ( Figure 19). The apparent correlation dimension of the daily set is significantly lower than that of the monthly-sampled model data. Because this probably indicates an oversampling problem , the calculation was repeated using a maximum norm in the integration instead of the Euclidian norm , as suggested by Grassberger [1986], but no improvement over the original calculation was found ( Figure 20).

Noisy Vallis Data
To evaluate the effects of noise on the determination of the correlation dimension , various levels of normally distributed noise were added to the 788-point, monthlysampled, Vallis data set. The correlation dimension calculation is inconclusive if the noise has a standard deviation over 103 of the signal standard deviation (Figure 21). Nevertheless, with 5-103 noise , the correlation integral slopes display weak convergence similar to the boundedness observed for the noise-free data ( Figure 22).

SOI
There is a weak indication in Figure 23 that De ~ 10 for the SOI data. This indicates that the supposed climate attractor has a fractal dimension greater than or equal to 10, that the noise in the data is sufficiently large to destroy the convergence in the correlation integral calculation, or that the data set is too short to yield a reliable result. In any case, this result is consistent with the results from Section 4: the Southern Oscillation does not appear to be an example of low-dimensional chaos.

Daily T-D
For embedding dimensions higher than 11 , the calculation indicates that the correlation dimension for the T-D data set is roughly 8.5 (Figure 24). The same De estimate is obtained using the maximum norm ( Figure 25). However , in the light of the analysis on the Vallis model cited above (Figures 17, 19, and 20) this may be an underestimate due to oversampling.   values sampled 20 times monthly. This calculation was done using a maximum norm rather than the standard Euclidian norm.
a. a. a. Surrogate noise is added as shown.

-d=3
::::::..   Thus measurement errors, which are almost certainly less than 60% of the signal standard deviation, are not causing the nonlinear signal to be lost. (3) The B/P algorithm detected nonlinear determinism in a model data set which mimicked the monthly time-averaging of the SOI. Thus the SOI's time averaging is not masking 51 any nonlinear signal in the data. ( 4) Spatial averaging, such as that done in processing the T-D data set, should average out any signals from small-scale physical processes unrelated to the large-scale dynamics. The fact that nonlinear determinism is not detected in the T-D data set suggests that those small-scale processes are also not responsible for the null result.
Thus, the most likely conclusion to be drawn from the B/P analysis in Chapter 4 is that E SO is not an example of low-order chaos. The short lengths of the data sets put computational limits on the size of K, , consequently the B/P analysis can only be used to discern low-order nonlinearity. For these data sets, /' \, = 7 seems to be the largest possible embedding dimension , so this result leaves open the possibility that ENSO is driven by mid-or high-order nonlinear dynamics.
Section 5 shows that the derivatives of the correlation integrals for the SOI and T-D data sets weakly indicate De ~ 8 -10. This calculation is less immune to noise than the B/P calculations in Chapter 4: the model correlation integral derivatives only converge when less than about 10% noise is added. Furthermore, the model correlation integral derivatives apparently converge to De values which are less than the real model De , when the model data set has a length and sampling rate comparable those of the data sets. However, monthly averaging does not affect the model De calculations noticeably.
Assuming that measurement noise is not significant, the correlation integral calculations, nevertheless, place a lower bound on the number of dimensions in which ENSO evolves. Any false results here would be an underestimate of the real De, not an overestimate. Thus, it appears that the correlation dimension of the SOI 52 is at least eight, i.e., not "significantly less than 10" . Taking this as our criterion for "low-order" chaos we conclude from Sections 4 and 5 that the SOI, and hence ENSO , cannot be usefully described as low-order chaos.

Conclusions about ENSO
If one accepts the idea presented here that ENSO is not low-order chaos, what ramifications does this have? Most simply, one can say that currently available chaos prediction techniques will not enhance our ability to predict ENSO events.
The lack of low-order nonlinear determinism implies that the tropical Pacific Ocean and atmosphere processes responsible for ENSO are complicated. Any models for ENSO which exhibit low-order chaos are probably too simple to produce reliable , accurate predictions of E SO episodes.
If daily Tahiti SLP data from the 1950's to the present could be obtained, they could be combined with the Darwin SLP data set to produce a nearly 17,000 term "daily Southern Oscillation Index" . With such a long series it would be possible, in principle, to investigate the hypothesis that E SO is governed by medium-order chaos by increasing"' in the B/P method to values in the range of 10-20.

53
For this research, the Southern Oscillation Index data set was obtained from the OAA Climate Prediction Center in Washington DC through their web site at http : //nic . fb4 .noaa . gov :80/data/cddb/ This data set was downloaded May 6, 1998. Two data files from the site were used in compiling the complete SOI data set; these are "soi" and "soi . his" .