## Abstract

The climate interannual variability is examined using the general circulation model (GCM) developed at the Laboratoire de Météorologie Dynamique. The model is forced by the observed sea surface temperature for the period 1979–94. An ensemble of eight simulations is realized with different initial conditions. The variability of the Southern Oscillation is studied. The simulated sea level pressure anomalies at both Tahiti and Darwin are realistic compared to observations. It is revealed, however, that the simulated convection activity response to the warm episode of El Niño is too weak over the eastern part of the tropical Pacific. This explains why the simulated Pacific–North American pattern is shifted westward. A global El Niño pattern index is defined and calculated for both the simulation and the National Centers for Environmental Prediction (NCEP) reanalysis data. This serves as a quantitative measure of El Niño’s global impact. A singular value decomposition analysis performed with the tropical Pacific sea surface temperature and the Northern Hemisphere 500-hPa geopotential height shows that the model’s teleconnection between the Tropics and high latitudes is similar to that of the NCEP reanalysis data.

In an exploratory manner, the model’s internal variability versus the external forced variability is studied. It is shown that, except for the equatorial strip, the internal model variability is larger than the external variability. An ensemble mean is thus necessary in order to focus on the model’s response to external sea surface temperature anomalies. An attempt is also made to evaluate statistically the influence of the ensemble’s size on the model’s reproducibility. It is shown that, with this particular GCM, at least five realizations are necessary to correctly assess the teleconnection between the Tropics and the Northern Hemisphere extratropics. This dependency on the number of realizations is less for the tropical circulation.

## 1. Introduction

One of the main objectives of large-scale climate models is to understand the mechanisms that control the present-day climate and future climate change. To get good confidence in a model, it is not sufficient to validate the model against the time-mean climatology, but it is important to examine also the model’s performance in reproducing interannual variability, such as the El Niño–Southern Oscillation (ENSO), a major tropical variability with global impacts.

Numerous studies have reported the use of atmospheric general circulation models with time varying observed global sea surface temperature (SST). Among others, Lau (1985) studied the Geophysical Fluid Dynamics Laboratory (GFDL) model response to the observed SST in 1962–76 and showed that the model was capable of reproducing the main aspects of the ENSO, including its teleconnection patterns at higher latitudes. Latif et al. (1990) reported an experiment using the model of the European Centre for Medium-Range Weather Forecasts (ECMWF) and observed SSTs in mid- and low latitudes for the period 1970–85. The simulated Southern Oscillation index is in good agreement with the tropical Pacific SST variation. Kitoh (1991) studied the interannual variations of the Japan Meteorological Research Institute GCM forced by 1970–89 SST. The simulated variation of the Southern Oscillation is also in good agreement with observations.

Recently in the scope of the Atmospheric Model Intercomparison Project (AMIP; Gates 1992), which is an international effort to evaluate the performance of different climate models and their systematic errors, the study of interannual variability has become routine for an atmospheric GCM. Stephenson and Royer (1995) reported the interannual variability in the Météo-France model with different resolutions. They showed that the model is capable of reproducing the dipole structure of sea level pressure (SLP) over the Pacific, which is, however, shifted eastward, a systematic default in many models. The Tahiti SLP variation is then poorly correlated with the SST anomaly, and Easter Island, which is situated to the east of Tahiti, should give better results. Kleeman et al. (1994) reported the results for the model of the Australian Bureau of Meteorology Research Centre. They emphasized the feedbacks related to the wind stress anomaly and energy budget anomalies during El Niño events. Their result shows a local negative feedback caused mainly by changes in evaporation. However, they found that the radiative variation due essentially to the low-level cloud changes has a small magnitude.

Due to the fact that the internal variability of models can be important, many studies use an ensemble approach to assess atmospheric responses to SST anomaly. The basic idea is that climate simulations are mainly a boundary-forced problem and model internal variability can be considered as noise. Although the cause of model internal variability is not fully determined, it seems that ensemble means can overcome the difficulties of the model internal variability. Harzallah and Sadourny (1995) conducted such an ensemble experiment by using the Laboratoire de Météorologie Dynamique (LMD) GCM with prescribed SST for 1970–88. They performed seven parallel simulations with different initial states and split the total variability into two parts: an internal part due to atmospheric dynamics and an external (or forced) part due to the variability of SST forcing. From the point of view of seasonal forecasts, Stern and Miyakoda (1995) analyzed the ensemble experiment conducted with the GFDL GCM for the time period 1979–88. They define a reproducibility parameter to measure the spread among ensemble members. Reproducibility is, of course, a necessary condition for any seasonal-scale climate prediction that is mainly a boundary condition problem. It was shown that the best reproducibility occurs in the equatorial regions, especially over the Pacific, and that they become poorer at higher latitudes.

Barnett (1995) also performed an ensemble of climate forecasting experiments by using the University of Hamburg climate model 3 (ECHAM3) GCM (Max Planck Institute, Hamburg). The standard deviation from the different individual experiments was considered as a measure of the internal model variability. This quantity was compared to the model external variability obtained from the forced interannual changes of the boundary conditions (observed SST). He then applied an empirical orthogonal function (EOF) analysis to the internal model variability in order to find the spatial coherent structures. This “spectral” signature of the internal model variability can constitute a good quantitative first guess of the expected model-induced uncertainty in a climate forecast. Barnett found that the internal model variability can be very large, and he concluded that a single model simulation of an interannual climate event or a climate forecast is totally inadequate for judging the model’s abilities or for providing a reliable forecast. Wang and Rui (1996) proposed another methodology for assessing ensemble experiments. They performed the analyses on the relationship between spread and ensemble mean in EOF space. By analyzing an ensemble experiment conducted by the National Centers for Environmental Prediction (NCEP) they concluded that zonal mean responses of ENSO are less predictable, while wavy responses are much more robust in the ensemble.

In this paper, we will present an ensemble experiment performed with a new version of the LMD GCM (LMD-Z model, version 1). The experimental design is very close to that used in Harzallah and Sadourny (1995), but the prescribed SST covers a more recent period (1979–94) in which more observational data are available for comparison. This study is mainly an exploratory one in which we try to explore the ensemble characteristics of our experiment. As in other similar experiments, we will study the model’s spread around the ensemble mean in order to isolate the SST-forced variability (climate signal from the point of view of climate prediction) from the internal model variability (climate noise). We present also a Monte Carlo–based technique in determining the ensemble size necessary to correctly assess a given climate signal. Throughout the paper, we also pay attention to the realism of our model, since it gives the final confidence (or uncertainty) to whatever is deduced from numerical models. We will use mainly the NCEP reanalysis data of the same period and some other observational data to validate the model’s behavior. The validation is performed in terms of mean climate, interannual variability, and global teleconnection structures.

In the next section, we will first describe the experimental design. The main results concerning tropical Pacific variability and global teleconnection patterns are presented in sections 3 and 4. Conclusions follow in section 5.

## 2. Experimental design

The model used in this study (LMD-Z, version 1) is derived from the LMD standard GCM (Sadourny and Laval 1984). The code has been extensively rewritten, has more physical consistency, and is now organized in a more flexible manner for ease of use. In the standard model used in LMD, the meridional discretization was regular in the sine of the latitude. This meant that the model had a poor resolution for high latitudes and the simulation of atmospheric circulation for these regions was less accurate. That is why the model’s grid points in the present study are chosen to be regularly distributed in the latitude–longitude coordinates with resolution of 4° × 5°. A systematic comparison is now under way to study the influence of changing latitudinal distribution and will be reported in a future paper. The model has 11 vertical layers, unevenly spaced to allow for a better resolution in the boundary layer. The time step is 6 min, however, the contributions of physical parameterizations (see the appendix for a brief description) are evaluated only every 30 min.

The model was integrated from January 1979 to December 1994. The SST and sea-ice distributions are the Center for Ocean–Land–Atmosphere (COLA)–Climate Prediction Center (CPC) analyzed data (Reynolds 1988). The original data are interpolated into daily values by a spline scheme and into our model grid by an area-weighted scheme. The ensemble of simulations contains eight members. The initial conditions of the first run are from a 6-month simulation forced by the climatological SST, while the initial states of subsequent realizations were assigned values that were the same as those at the last time step of the preceding integration. The present study hence entails altogether 16 × 8 = 128 yr of model integration.

Before going into the climate variability, we first examine briefly the model’s climatology in relationship with the variability studies presented in the rest of the paper. Figure 1 compares the model’s 200-hPa zonal wind [December–February (DJF) and June–August (JJA)] to the NCEP reanalysis data for the same period. The simulation is successful in capturing most of the principal features of the observed distribution. The model can simulate correctly the position and orientation of the jet stream over the east coasts of Asia and North America but fails to give sufficient amplitude for DJF. Furthermore, the tropical easterlies are slightly weaker. Figure 2 compares the meridional streamfunction for both simulation and observations. The model is capable of correctly simulating the meridional circulation: the direct Hadley cells and the indirect Ferrel cells. However, the dominant winter hemisphere’s Hadley cell for DJF is weaker in the model than in the observations. This is probably the reason why the jet stream over 200 hPa is weaker in the model. The Southern Hemisphere’s winter Hadley cell seems too strong. Figure 3 shows the Northern Hemisphere winter season (DJF) 500-hPa geopotential height for simulation and observation, respectively. The planetary waves are reasonably well simulated given the low resolution of the model, but the ridge over the Rocky Mountains is too weak.

## 3. Tropical variability

The main tropical circulation is related to the Southern Oscillation (SO), which was discovered by Walker (Walker and Bliss 1932). Its relationship to the Pacific El Niño events was suggested by Bjerknes (1969). Figure 4 shows the SST anomalies of the Niño-3 region (5°S–5°N and 92°–150°W), which is a good proxy (Niño-3 SST index) for El Niño events. Since the main focus of this paper is on the interannual variability, we apply a low-pass 11-term filter (Trenberth 1984) to the anomaly field in Fig. 4 in order to eliminate signals with periods less than 8 months. Several warm and cold episodes occurred over the experiment period. We can note the El Niño events of 1982–83, 1987, 1991–92, and the most important La Niña event of 1988.

Figure 5a and 5b show, respectively, the observed (dashed line) and simulated (solid line) SLP anomalies at Tahiti (eastern Pacific) and Darwin (western Pacific). The simulation is the *ensemble mean* (see section 4d for a more rigorous definition) calculated for the two model points that are, respectively, closest to Tahiti and Darwin. The associated error bars are from the standard deviation among the ensemble members. For both Tahiti and Darwin, the standard deviation is about 0.4 hPa with a small seasonal variation. However, the simulated interannual variation is between 1 and −1 hPa. It is clear that the internal variability of a model is not negligible, even in the Tropics.

The SLP anomaly at Darwin is positively correlated with the Niño-3 SST index, and that at Tahiti is negatively correlated with the Niño-3 SST index. For both Tahiti and Darwin the simulated variations are very close to the observed ones. However, the simulated amplitude is smaller than the observed one for the 1982 El Niño event and the 1988 La Niña event. The mean correlation coefficients are, respectively, 0.71 for Tahiti and 0.75 for Darwin.

Let us now examine the statistical significance of these correlation coefficients. In general, to test the statistical significance of the linear correlation coefficient *r*_{xy} obtained from two time series (*x* and *y*) of length *N,* one can calculate (Press et al. 1992, chap. 14) the statistic *t* = *r*_{xy}[*μ*/(1 − *r*^{2}_{xy})]^{1/2}, where *μ* is the number of degrees of freedom. If the distribution of the two time series forms a two-dimensional Gaussian distribution around their mean values, the statistic *t* is distributed in the null hypothesis (no correlation) like Student’s t-distribution with *μ* degrees of freedom, whose two-sided significance level is given by 1 − *A*(*t*/*μ*), where *A* is the Student’s distribution. It is evident that the number of degrees of freedom may not be the same as the number of data points that enter into the calculation of correlation coefficient. We use here the method of Davis (1976) to estimate *μ*:

where

and *r*_{xx} and *r*_{yy} are the autocorrelation coefficients of *x* and *y.* Here *τ* indicates the average number of data required to gain a new degree of freedom in the estimation of *r*_{xy}. The results show that the correlation coefficients, for both Tahiti and Darwin, are statistically significant with a confidence level greater than 99%.

Figure 6 shows the correlation map of the SLP anomaly over the equatorial strip with the Niño-3 SST index (as shown in Fig. 4) for both simulation and observation. In the case of the model, all individual members of the ensemble are used in calculating the correlation. The dipole structure is dominant in the Tropics with a negative center over the eastern Pacific and a positive center over the western Pacific and Indian Oceans. There is, in general, a good agreement between the model and the reanalysis. The model’s structure is, however, more zonally elongated, compared with NCEP data, which are more meridional and have larger extension to higher latitudes.

The 850-hPa trade wind over the equatorial Pacific is also a good index of the ENSO strength. In general, El Niño events imply a weakening of the trade wind. For three selected regions (western Pacific: 5°N–5°S and 135°E–180°; central Pacific: 5°N–5°S and 175°–140°W;eastern Pacific: 5°N–5°S and 135°–120°W), both the simulated (solid line) and observed (dashed line) 850-hPa trade wind anomalies (zonal component) are drawn in Fig. 7. The correlation coefficients between simulation and observation for the western and central Pacific reach very high values (0.84 and 0.86, respectively). However, it is less for the eastern Pacific (0.70). The amplitude for the western Pacific is slightly larger in the model than in the observations. As one moves to the east, the model’s amplitude becomes too small, especially for the eastern Pacific. This could be related to our convection parameterization, which may not respond correctly to the SST anomaly for the eastern Pacific. The error bars plotted in Fig. 7 indicate the standard deviation among the 8 members. We can see that it is larger over the western Pacific than over the central and eastern Pacific. For the western Pacific, the convective activity is strong and its variation is also large.

Outgoing longwave radiation (OLR) is an efficient measure of the tropical convection and precipitation rate. Negative OLR anomalies are a manifestation of the increase of cloudiness and/or enhancement of cloud tops. Positive anomalies are usually indicative of suppressed convection. Figure 8a shows the OLR variation in a longitude–time diagram, averaged over the zonal belt of 5°S–5°N. Figure 8b shows the National Oceanic and Atmospheric Administration (NOAA) satellite observations. In general, the agreement between simulation and observation is very good, especially for the El Niño–La Niña cycle of 1987–89. However, the model fails to reproduce the amplitude of the strong El Niño event of 1982–83, where the SST anomalies occurred mainly over the colder oceans in the eastern part of Pacific. In the next section, we will show that this lack of atmospheric convective response is mainly responsible for the westward shift of the Pacific–North American pattern.

## 4. Global teleconnection patterns

El Niño is essentially a tropical phenomenon but its impacts are global. Many studies in the literature discuss teleconnection patterns—how anomalies propagate to other regions of the globe. By using long-term observation data, Wallace and Gutzler (1981) and Barnston and Livezey (1987) documented the main teleconnection patterns. They used, respectively, the regression technique and the rotated EOF analysis. Hoskins and Karoly (1981) provided a theoretical explanation of the global teleconnection pattern by using the principle of Rossby wave dissipation over a rotating sphere. Recent studies based on modeling (Lau and Nath 1994, 1996) or observations (Deser and Blackmon 1995; Zhang et al. 1996) emphasize more and more the role of an “atmospheric bridge” in relating climate variabilities over different regions.

In the scope of the AMIP experiment, most of the GCMs available in the international community have been run over several El Niño episodes. It is important to verify the GCMs’ performance in reproducing the teleconnection pattern. This kind of diagnostic study should become standard in the future. In this paper, we use two statistical approaches to validate the model’s performance against observations (NCEP reanalysis data). The first approach consists of examining the correlation between the Niño-3 SST anomaly and the anomaly field of different climate variables over the globe. Since the Niño-3 SST anomaly is a proxy for El Niño events, the correlation maps give a good representation of structural differences between a composite El Niño event and a composite La Niña event. The second approach is based on the empirical orthogonal function analysis, which searches for the main spatially coherent structures of a single field and the singular value decomposition analysis, which is a powerful tool in determining the correlation structures between two fields. The basic idea is to first identify the spectral information of the atmospheric circulation and then compare model simulation against observations.

### a. Niño-3 SST indexed global correlations

Figures 9a and 9b (contour) show, for both the simulation and observations, the correlation map for 200-hPa geopotential height. In the case of the model, the individual members are used in calculating the correlation. As in Chen (1982), we can also apply the above-mentioned statistical significance test (section 3) to the correlation map. The shaded regions indicate confidence levels larger than 90%. For the NCEP data, the correlation for the entire tropical belt (20°N–20°S) is statistically significant. The whole North Pacific region exceeds the 90% level with an extension to Siberia and China. Two other regions are the south and east of the United States. They correspond to the active centers of the PNA pattern. Several significant regions are also observed in the Southern Hemisphere around 30°–60°S. The model’s 90% significant level covers larger areas. This is the direct effect of ensemble simulation in which the number of degrees of freedom is largely increased.

Both the observations and simulation give positive correlations for the whole tropical zone. The structure is of a weak wavenumber-two pattern, with one maximum over the eastern Pacific and the other over the western Indian Ocean. The Pacific maximum has a specific structure, two maxima south and north of the equator, which suggests the solution given by Gill (1980) for a heating symmetric about the equator. The maximum over the western Indian Ocean is centered over the equator. Both of the structures are tightly related to a shift in convective activity, a decrease of diabatic heating over Indonesia, and its migration to the central Pacific. This modification affects the equatorial Kelvin wave and the disturbance propagates poleward and eastward through a Rossby wave dissipation mechanism. We can easily recognize two wave train structures for both of the hemispheres. These results are in agreement with other published observations by Horel and Wallace (1981) and Karoly (1989). This gives us confidence regarding the performance of the model in the Tropics and subtropics. However, for high latitudes, especially over the PNA region, the disagreement between simulation and observations is large: the model seems to give a shift of correlation structure to the west, and the correlation coefficient is also smaller than in the observations.

We can now examine the regression map (Figs. 9c and 9d) between the 200-hPa geopotential height and the Niño-3 SST index. The map is obtained by multiplying the geopotential height standard deviation with the correlation map. The map indicates the geopotential height variation corresponding to a typical amplitude of the Niño-3 SST index. We can observe that the signals are essentially over the tropical regions, the whole Pacific, and North America. For the latter region, it is clear that the model’s structure is shifted westward; in particular, the positive anomaly center in the eastern part of North America is now completely shifted to Alaska in the model. This westward shift is probably related to the fact that the model’s tropical response to the El Niño events, as illustrated by the OLR anomaly (Fig. 8) or 850-hPa trade wind anomaly (Fig. 7c), is weaker over the eastern Pacific. So the local Hadley circulation does not move close enough to the North American coast for the El Niño events. We hope to overcome this model discrepancy in the future by using other more recent convection schemes (e.g., Tiedtke 1989; Hack 1994).

If we consider the regression map (or the correlation map, which is biased to the tropical regions) as the global structure of the El Niño impact, we can introduce now the notion of “El Niño global pattern index.” It is defined as the projection of the anomaly field to the observed regression map (Fig. 9d):

where *υ*_{m} is the regression map loading (normalized to unity), *Z*_{m}(*t*) is the anomaly field at time *t,* and *M* is the number of spatial points. We plotted *q*(*t*) in Fig. 10 for both observation and simulation (in the case of the simulations, the ensemble mean and associated error bars are drawn). The correlation coefficient between observation and simulation reaches 0.67 for the whole period with a statistical confidence level at 93%. This illustrates that the model, in the global scale, can simulate well the El Niño–type interannual response, although the main contribution to the good agreement comes from the Tropics. Notice, however, that since it is not weighted by area there is an underestimation for low latitudes. The weighted index (not shown here) increases further the agreement between observation and simulation. Careful inspection of this figure seems to reveal two different regimes. The first one is from 1982 to 1991, where ENSO episodes are very strong and the statistical dispersion is much smaller than the climate signal. The agreement between observation and simulation is very good, although the model’s response for the 1988/89 La Niña event is too weak. The second regime covers the beginning and the end of the figure for which the ENSO signal is weak, the model’s spread is large, and the discrepancy between observation and simulation is more important. This behavior of our model is coherent with recent results obtained by K. Miyakoda, A. Navarra, and N. Ward (1997, personal communication) showing that the tropical predictability increases when the SST anomaly is large and decreases when the SST anomaly is small.

### b. Northern Hemisphere winter circulation and its connection to the tropical SST

We now study the Northern Hemisphere winter circulation by using the 500-hPa geopotential height. We concentrate on the relationship between the Northern Hemisphere winter circulation and the central Pacific SST. The statistical method employed is the singular value decomposition (SVD) analysis, which is a powerful technique for identifying the dominant coupling modes between two data fields. The formulation of this statistical method was given in Bretherton et al. (1992), to which the interested reader is referred for further details. The application of this tool to the extratropical air–sea coupling problem has been demonstrated in Wallace et al. (1992), Lau and Nath (1994), and Zhang et al. (1996), among others.

Note that the SVD technique is applied here to the cross-correlation matrix as in Lau and Nath (1994); that is, the anomaly fields are normalized by their standard deviation. Furthermore, in order to take into account the spatial inhomogeneity of the grid, the anomaly fields are also normalized by cos*ϕ,* where *ϕ* is the latitude of the grid. We should emphasize that the NCEP data have been objectively interpolated to the model’s grid before the statistics are applied, so that the analysis was done in the same manner for both model and NCEP reanalysis. Note also that the model’s results are obtained by using the *ensemble mean* (in the next section, we will study the influence of internal model variability on the results).

Figure 11 shows, respectively, the leading SVD mode from the model and from the reanalysis data. The bottom panels give the expanded time coefficients, which are calculated by projecting each data field (SST or 500-hPa geopotential height) on its corresponding singular vector at each of the time points. The bottom panels represent the temporal evolution of the amplitude and polarity of the characterstic patterns. The top (middle) contour panels give the heterogeneous correlation maps, which are the temporal correlation coefficients between the data field of 500-hPa geopotential height (SST) and the above-mentioned expanded time coefficient for SST (500-hPa geopotential height). They highlight the spatial features in one field that exhibit strong temporal relationships with the companion field. For the tropical Pacific SST, it is a typical El Niño structure for both observation and simulation. The resemblance between the model and the reanalysis is also good for the 500-hPa geopotential height. Both show a coherent structure in the Tropics and several active centers for the high latitudes. The dominant structure is the PNA pattern, but other teleconnection patterns are also merged onto this leading structure. Unlike EOF technique, we cannot apply any rotation to SVD in order to better isolate different structures. On the other hand, all these patterns can be interpreted to be tightly correlated to the tropical Pacific SST. Several discrepancies should be noted between Figs. 11a and 11b: the negative centers over the north of Japan and the North Atlantic are absent in simulation. The pair of centers over eastern and central Asia is eastward shifted to China. The Aleutian negative center of the PNA structure is, however, westward shifted.

To be sure that the obtained SVD structures are not spurious results (see Newman and Sardeshmukh 1995 and Cherry 1996 for some cautions regarding SVD analysis), we also performed EOF analysis separately for every field: the tropical Pacific SST and Northern Hemisphere 500-hPa geopotential height, the leading EOF structures (not shown here), are very close to the leading SVD structures.

Some statistical quantities are summarized in Table 1 to characterize the strength of coupling [fraction of the squared covariance (SCF); temporal correlation coefficient between the two expanded time coefficients, *γ*] and the significance of the heterogeneous correlation patterns [fraction of domain-integrated variance (VARF)]. We can see that all the indicators are similar for the model and reanalysis data. So we could conclude that our model is capable of reproducing the linkage between the Tropics and the extratropics.

### c. Internal model variability versus external variability

Previous authors (Harzallah and Sadourny 1995; Stern and Miyakoda 1995; Barnett 1995) have shown that the internal model variability can be very large for the extratropical circulation. This variability can make a large impact on the statistically based results of model output. We now perform, once again, the SVD as in section 4b. But instead of using the ensemble-mean data, we use the eight simulations (16 yr each) as a “grand” simulation of 128 yr. So the total signal is now analyzed. The last line of Table 1 gives the main statistical indicators for the leading mode. They are in general smaller than those obtained by the external signal, except for the value of SCF. The leading SVD structure is shown in Fig. 12. The comparison with Fig. 11a reveals that the spatial patterns are roughly the same, but the actual values are smaller. The comparison with NCEP reanalysis (Fig. 11b) is less satisfactory in Fig. 12 than in Fig. 11a. It is clear that the presence of internal model variability can modify quantitatively the results. But the main results hold qualitatively and remain robust.

One important question in the approach of climate ensemble simulations is the following: what is the minimum number of runs necessary to overcome the internal model variability and under study, with confidence, the model’s response to external forcing? It is clear that the answer depends upon the signal to study and the method used to assess the signal. To give a quantitative answer to the question, we use a Monte Carlo method, already used by Cheng et al. (1995) and Richman (1986) for a similar purpose, and repeat the SVD analysis for subsets of the eight simulations. We can construct eight classes of subset by taking, respectively, 1, 2, . . . , or 8 runs. For the subset *n,* the number of possible combinations is *C*^{n}_{8}. For convenience, the structure of the subset 8 (the total ensemble) will be referred to as master structure (Fig. 12). The structures obtained from other subsets will be called perturbed structures. We can now define formally a congruence coefficient to measure the resemblance between the master structure and the perturbed ones:

where *m*_{i} and *p*_{i} denote elements of the respective (master or perturbed) structure vectors. The form of this coefficient is similar to a correlation coefficient, as the possible values range from 1 for total agreement through 0 for no relationship to −1 for total disagreement. As stated in Richman (1986), the congruence coefficient is a better measure than the pattern correlation coefficient in the context of quantitatively assessing the goodness of match of loading vectors since the former measures the similarity between two patterns without removing their respective means. However, the congruence coefficient is biased toward a higher value than the correlation coefficient. According to Richman (1986) and the references therein, it is reasonable to attach the following goodness of match names to specific ranges of absolute congruence coefficients: 0.98–1 (excellent match), 0.92–0.98 (good match), 0.82–0.92 (borderline match), and less than 0.82 (poor match). Table 2 shows the results of our calculation. We can see that, for a single realization, five cases out of the eight are under the borderline. However, for a combination of five simulations, all the congruence coefficients are in the good or excellent range. We conclude that at least five or six simulations are necessary to obtain a correct assessment of the teleconnection between the Tropics and the 500-hPa circulation of the whole Northern Hemisphere. The time period of our study (1979–94) includes roughly three ENSO cycles, so our conclusion can be then extrapolated to another form; that is, at least 15 ENSO cycles are necessary to get a good response for the Northern Hemisphere teleconnection patterns. Here we used implicitly the assumption that the climate system is an ergodic one.

To demonstrate that the minimum number of realizations is signal dependent, we replace the field of whole Northern Hemisphere 500-hPa geopotential height by that of the 200-hPa velocity potential over the region of 45°S–45°N and repeat the same analysis procedure. The results are given in Table 3. All coefficients pass the borderline level. With a minimum number of 3 (even 2), the assessment of the signal is thus robust. This is coherent with our general knowledge that the 200-hPa velocity potential field, especially in the Tropics, is very strongly related to the ENSO cycle, and the associated internal model variability is small.

### d. Analysis of variance

To complete the studies on the model’s external and internal variability, we now employ another widely used technique, analysis of variance. It has already been used by Harzallah and Sadourny (1995), Zwiers (1996), and Wang and Rui (1996), among others, in the context of ensemble simulations. Assume that *x* is a time series of a monthly mean (or seasonal mean) variable, *x*(*n, y*). Here *x*(*n, y*) represents the *n*th member among an ensemble of size *N* (8 in the present study) for the *y*th year of a simulation of length *Y* (16 in the present study). We can define two averages: ensemble mean and climatological mean. The ensemble mean is taken with respect to the simulation counter *N*:

And the climatological mean is taken with respect to years *Y* and simulations *N*:

For climate forecasting purposes, the most interesting signal is the interannual variability and is represented by the year-to-year variation of the ensemble mean *x*_{e}(*y*). We can use the standard deviation *σ*^{2}_{e} to measure the external, or forced, signal. However, the dispersion of the simulations (measured by the standard deviation *σ*^{2}_{i}) indicates the internal model variability that is, in some sense, noise for climate forecasting. From the definition of the two averages, we have

and

It is easy to show that the total variability, defined by the total standard deviation,

is the sum of the external and internal variabilities:

To illustrate the above definitions, we take as an example the geopotential height at 500 hPa. We examine only the Northern Hemisphere winter season DJF. Figure 13 shows, respectively, the distributions of *σ*_{e}, *σ*_{i}, and the ratio *λ* (=*σ*_{e}/*σ*_{i} − 1). For the external variability, we find maximum values over the North Pacific with extensions to North America and to Asia. Another high-value region is over the North Atlantic, near the European coast. For the internal variability, the largest values occur in the vicinity of the stationary troughs over the northern oceans, where we also observe large external signals. The overall pattern shows increasing scatter as one progresses toward the pole. We can see from Fig. 13c that the ratio *λ* is positive for the tropical regions and negative poleward of 20°. One exception is in Mexico where the external signal is strong and the internal variability is relatively small.

As in Barnett (1995), we can now assess the spatially coherent structures contained in the internal variability by performing an EOF analysis. This gives the spectral signature of the model internal variability. Note that the EOFs here are obtained by resolving the eigenvalue problem of the correlation matrix in order to include the structures at low latitudes, since the variance at low latitudes is small. Figure 14 shows the first two spatial EOF structures of the internal variability. EOF-1 explains 25% of the total variance and the corresponding large-scale coherent structure is mainly confined to the Tropics. EOF-2 accounts for 13% of the variance and shows a clear wave train structure originating in the North Pacific. It bifurcates into two parts, one which propagates to the North Pacific and then returns equatorward and eastward, another which crosses the North Pole and reaches Europe. EOF-3 (not shown) explains 7.3% of the total variance and shows coherent structures related to the subtropical jets. EOF-4 and the remaining EOFs show less coherent structures and mainly represent the model’s random noise.

From Fig. 14b, it is clear that the internal variability presents a coherent spatial structure over the PNA sector. So it is plausible that the PNA pattern is not totally a response to the tropical SST anomalies, but an intrinsic characteristic of the high-latitude circulation (Simmons et al. 1983). The number of our ensemble simulations is not large enough to confirm this hypothesis, since the above-calculated internal variability could be contaminated by the SST-forced variability.

## 5. Summary and conclusions

By using an atmospheric GCM forced through the observed SST from January 1979 to December 1994, we have completed an ensemble experiment with eight members. Our results confirm the conclusion of Barnett (1995); that is, a single model simulation is totally inadequate and misleading for studying the interannual climate variability, since at latitudes higher than 20°N, the internal model variability is larger than the forced variability. Furthermore, by using a Monte Carlo approach, we evaluated the influence of the internal model variability on the external climate signal and showed that five or six is the minimum number of ensemble simulations to assess the teleconnection between the tropical Pacific SST and the 500-hPa circulation for the whole Northern Hemisphere. However, this minimum number can be reduced if we are only interested in tropical signals.

We analyzed the 500-hPa geopotential height for the Northern Hemisphere winter; the internal variability did have coherent spatial structures, which shows that it is not totally a random noise. Our results suggest that the PNA pattern could be a response to the tropical SST anomalies, but it could be also an intrinsic characteristic of the Northern Hemisphere circulation.

Since the results presented may well be model dependent, we increase confidence in our evaluation by also presenting the performance of our GCM in reproducing forced interannual variability. The variation of the Southern Oscillation was extensively studied in this paper. A well-developed dipole structure was obtained in the correlation map for sea level pressure: a negative center over the eastern Pacific, and a positive center over the western Pacific and Indian Oceans. The signal was strong for both Tahiti and Darwin. The variation of the 850-hPa trade wind was also well correlated to observations, especially for the western and central Pacific.

By using the correlation map between the Niño-3 SST index and the climate variation at 200 hPa, we have evaluated the upper-tropospheric response to El Niño and the global teleconnection patterns. The most significant response was over the Pacific. It is related to the displacement, from the Indonesian region to the central Pacific, of the intensive convective activity and diabatic heating center. This agrees with the response of the atmospheric circulation to tropical heating predicted theoretically by Gill (1980). It was also consistent with observation-based analyses given by Horel and Wallace (1981) and Karoly (1989). A wave train pattern propagated poleward and eastward for both hemispheres. An induced response over the tropical Indian Ocean and Africa was also evident in our model. It was different from the Pacific response, since the maximum correlation center was over the equator. The western Indian Ocean perturbation also propagated to high latitudes in a wave train structure. This behavior of the model was in agreement with Rasmusson and Mo (1993) using 1986–89 National Meteorological Center (now the National Centers for Environmental Prediction) operational analysis.

By combining the EOF and SVD methods, we studied the coupled structures of tropical Pacific SST and the Northern Hemisphere winter 500-hPa geopotential height. The same method was also applied to the reanalysis data. The model was capable of reproducing the teleconnection between the tropical Pacific SST and the Northern Hemisphere circulation.

The present study also revealed several problems with our model. The main one is that the atmospheric response to the warm episode of the ENSO oscillation is too weak over the eastern part of tropical Pacific. This probably induces a westward shift of the PNA pattern through the local Hadley circulation (Rasmusson and Mo 1993; Tyrrell et al. 1996). It will be interesting to study the influence of convection schemes in future versions of our model.

## Acknowledgments

This work is partly supported by the Environment Program of the European Commission and the French national program on climate dynamics (PNEDC). The author would like to thank David Stephenson, Tim Barnett, Bryan Weare, and Bryant McAvaney, whose suggestions greatly improved the manuscript. Useful discussions with Robert Sadourny and Hervé Le Treut are acknowledged. Comments from referees helped to clarify many points of the text. Technical assistance from the colleagues of LMD and IDRIS was indispensable for the realization of this work. Contour graphics were drawn by using GrADS developed at COLA. Observation data in Figs. 5, 7, and 8 were obtained through the ftp server of the NOAA/CAC. The NCEP reanalysis data were supplied by NOAA/CDC and were available through the LMD Climate Data Base Project.

## REFERENCES

### APPENDIX

#### Brief Description of Model’s Physical Parameterization

The convection parameterization is a combination of the moist adjustment scheme (Manabe et al. 1965) and the Kuo (1965) scheme. The first is used to eliminate the conditional instability when saturation occurs, and the second to take into account large-scale water vapor convergence and evaporation from the surface. These quantities are then considered as the convection forcing to enable the vertical mixing of energy and water.

Cloud liquid water is a prognostic variable of the model, which is calculated by considering the source and sink terms. At present, the clouds created by convection are not considered. Only those associated with stratiform precipitation are taken into account. Partial condensation is allowed through a statistical approach (Le Treut and Li 1991).

The radiative transfer code has been designed by Fouquart and Bonnel (1980) for the solar radiation and by Morcrette (1991) for the terrestrial radiation. In the current experiments, the diurnal cycle of insolation is not taken into account. The radiation code is called at every physical time step to calculate the influence of clouds. However, the clear-sky transmission coefficients are evaluated only every 4 h.

The surface model is a bucket model for which we consider only a homogeneous layer of 150 mm. The calculation of surface temperature is incorporated in the boundary layer and based on the surface energy balance equation. For the surface moisture, a holding capacity is fixed at 150 mm of water, and all the water above this value is lost as runoff.

## Footnotes

*Corresponding author address:* Dr. Z.-X. Li, Laboratoire de Météorologie Dynamique du CNRS, Ecole Normale Supérieure, 24 rue Lhomond, 75231 Paris Cedex 05, France.

Email: li@lmd.ens.fr