Abstract

Seasonal hindcasts were made using the COLA model for 16 winter seasons (mid-December through March for 1981/82 - 1996/97). For each season, a nine member ensemble was generated using observed initial conditions in mid-December and observed global sea surface temperature (SST). It is found that in the presence of large tropical SST anomalies, the model is quite successful in simulating seasonal mean height anomalies over the Pacific North America (PacNA) region.

A local spatial pattern correlation field is computed for the ensemble, seasonal mean of 500 hPa height from the GCM and the seasonal means from the reanalyses of the National Centers for Environmental Prediction, hereafter observations. This field exceeds 0.6 over the eastern tropical Pacific and western North America; the maximum values are greatly enhanced during ENSO years. Similar results are obtained for 200 hPa u-winds. The observed enhancement of intraseasonal low pass (10-90 day) variability of 500 hPa height during cold events is simulated, as is the shift in the stormtracks (2-10 day variability of 850 hPa meridional heat flux) during warm events.

Empirical orthogonal function (EOF) analysis is applied in the PacNA region to the ensemble seasonal means of GCM 500 hPa height and the corresponding observed seasonal means; the global combined Arctic / North Atlantic Oscillation is removed from both GCM and observations. The leading EOF mode explains about 50% of the variance for both GCM and observations; the two patterns are nearly identical. Singular Value Decomposition (SVD) analysis between the tropical Pacific SST and 500 hPa height in the PacNA region shows that the nature of the coupling between SST and height is nearly the same in the GCM and observations. The great similarity between the height patterns in the first SVD mode and the corresponding leading mode EOF patterns indicates that the leading height variations are forced by SST.

The SST-forced variance of height was also estimated by regression analysis of (ensemble) seasonal means for the GCM and observations for the 16 years onto an index of tropical Pacific SSTs derived from SVD analysis of a long (30 year) observed record of heights and SSTs. The pattern of percentage variance explained in the GCM and observations are very similar to each other (and to the EOF described above). The higher absolute values in the GCM case reflect the effectiveness of the ensemble in filtering out variability unrelated to SST forcing.

SVD analysis was applied to100 GCM "samples" coupled with the observations; a sample is defined as a single 16-year record obtained by picking one ensemble member randomly for each of the 16 years. Probability distribution functions (pdf) of the pattern correlation for the leading SVD patterns, the percentage explained squared covariance, and the time series correlations all indicated sharp peaks at values of (0.87, 87%, and 0.82) respectively.

The pdf of the projection of individual 5-day means onto the leading EOF described above is quite dramatically shifted during strong warm and cold tropical SST events; the warm [cold] event pdf has almost all its weight in the negative [positive] EOF region. GCM and observations agree well.

The intra-ensemble spread was estimated by computing the PacNA anomaly correlation coefficient (ACC) for each of 36 possible intra-ensemble pairs. For nine years in which histograms of the ACC indicate predominantly positive values, the ACC of the ensemble mean with the observed seasonal mean is also relatively high.

Brier skill scores and reliability diagrams were computed for the "event" of the 500 hPa height being one standard deviation above (or below) the normal, with all such events pooled over the entire Northern Hemisphere or North America only. All skill scores are positive and statistically significant at the 99% level. The North America scores are higher than the whole hemisphere scores; the ENSO year scores are higher than those considering all years.

1. INTRODUCTION

Nearly 20 years ago, it was suggested that the year-to-year variations of the seasonal mean atmospheric circulation and rainfall in the tropical regions are largely determined by slowly varying boundary conditions at the Earth's surface, notably sea surface temperature (SST) and land surface conditions (Charney and Shukla, 1981; Shukla, 1981a). This suggestion was based on sensitivity and predictability calculations using atmospheric general circulation models (AGCM) in which the model simulated variances of circulation and rainfall with interannually varying SST (hereafter referred to as the signal) were compared to the variances generated by the internal dynamics of the atmosphere alone (to be referred to as noise). The noise was estimated by integrating an AGCM with climatological SST. Such calculations were repeated using several AGCMs, and, without exception all models showed that the signal in the tropical regions was far larger than the noise at a very high level of statistical significance. The evidence for the existence of extratropical predictability had been difficult to find. Lau (1981) showed that model simulations with climatological SST could capture 'the most prevalent anomaly patterns observed in the atmosphere.' Estimates of predictability based on model simulated signal and noise gave a generally pessimistic view about the prospects for dynamical seasonal prediction in mid-latitudes (Chervin, 1986). The pioneering monthly prediction experiments by Miyakoda et al. (1983, 1986) had shown the feasibility of dynamical monthly prediction, and studies on the dynamical predictability of monthly means had shown a scientific basis for predictability (Shukla, 1981b). However, it was clearly understood that it was the deterministic predictability of the first ten days that contributed to the dynamical predictability of the monthly mean. It was shown unambiguously that, in the absence of the possible influence of boundary conditions, there was no prospect for predicting monthly means even for the second month. These results further contributed to the skepticism about the prospects for dynamical seasonal prediction in the extratropical regions.

This impression was not consistent with the fact that the observed correlation coefficient (CC) between tropical SST anomalies and mid-latitude circulation anomalies, especially over the Pacific-North America (PacNA) region, was 0.6 or higher (Horel and Wallace, 1981). It was also found that the amplitude of the extratropical height anomalies simulated by some of the early models were too small by a factor of two or three compared to the observed anomalies (see Lau, 1997, for review). For some time it remained an unresolved question as to whether the predictability of seasonal averages in the mid-latitudes is limited by the chaotic nature of the internal dynamics, or whether the deficiency of the early AGCMs produced such a weak boundary forced response that this response was overwhelmed by the noise associated with internal dynamics.

Recently, several modeling studies have suggested that the deficiencies of AGCMs could be primarily responsible for unrealistically low values of simulated interannual variability and the correspondingly low estimates of mid-latitude predictability. A paper by Kumar et al. (1996) showed that a rather small change in the treatment of the physical parameterizations made a significant difference in the amplitude of model-simulated, SST-forced mid-latitude height anomalies. In particular, the amplitude of the SST-forced height anomalies in the revised model nearly doubled over the PacNA region. Bengtsson et al. (1993) and Brankovic et al. (1994) showed that the height anomalies over the PacNA region using two different models were quite well simulated. DeWitt (1996) showed that by replacing the Kuo convection scheme with the relaxed Arakawa-Schubert, or RAS scheme (Moorthi and Suarez, 1992) in the COLA AGCM, the amplitude of simulated mid-latitude height anomalies became comparable to the observed height anomalies, especially for cases with large tropical SST anomalies.

It is in this context that we proposed to conduct a well coordinated study to estimate the predictability of seasonal mean anomalies over the PacNA region, which is also the region of largest interannual variability of height anomalies over the globe. We proposed a multi-institution joint research project called Dynamical Seasonal Prediction (DSP) in the U.S. (Shukla et al, 1999). The DSP, and a similar European project called PROVOST (Prediction of Climate Variations on Seasonal-to-Interannual Time Scales), also contribute to an international project called SMIP (Seasonal Prediction Model Intercomparison Project) coordinated by the CLIVAR (Climate Variability) committee of the World Climate Research Program (WCRP).

There have been several simulations made with observed SST, conducted in a manner similar to that of the Atmospheric Model Intercomparison Project (AMIP; Gates et al. 1999). Some of these studies have carried out detailed and insightful analysis of seasonal predictability (Bengtsson et al. 1996; Kumar and Hoerling, 1997). The motivation for the DSP project is different in one important respect. It has been proposed that, in addition to the identical global SST and sea ice, all models will utilize the same initial conditions of the global atmosphere, including the land surface conditions. The idea is to calculate the theoretical upper limit of seasonal predictability with the current generation of models. These results could then be compared with the AMIP-type simulations, for which the initial conditions of the atmosphere and the land surface can be far different from observations. A systematic comparison of DSP integrations and AMIP integrations will help to estimate the importance of using the most accurate initial conditions of atmosphere and land surface properties for seasonal prediction. The DSP study described herein was unable to use the observed initial conditions of soil wetness and other land surface properties, because such a data set does not yet exist. However, all participating models used observed global atmospheric initial conditions for seasonal integrations.

This paper describes the results of seasonal hindcasts made using the COLA model for 16 winter seasons (1981/82 through 1996/97). For each winter season, the COLA atmospheric model was integrated from the observed global atmospheric initial conditions produced by the National Centers for Environmental Prediction (NCEP) in mid-December, and the integration was carried through 0000 UTC of 1 April of the following year. Section 2 describes the COLA atmospheric model and the initial and boundary data sets. The mean model error is given in Section 3. The ensemble mean anomalies and low and band pass variance patterns simulated by the model are compared to those of the reanalyses of NCEP in Section 4. In addition, Section 4 describes empirical orthogonal function analysis of the 500 hPa height, and singular value decomposition analysis linking the height field with the tropical Pacific diabatic heating. The statistical reliability and reproducibility of the results is described using a variety of methods is described in Section 5, while a summary and discussion is given in Section 6.

2. COLA MODEL, INITIAL AND BOUNDARY DATA SETS

2.1 COLA model

The COLA atmospheric general circulation model (AGCM) is a global spectral model which is truncated rhomboidally at zonal wave number 40 (Kinter et al, 1997). The associated Gaussian grid on which the physical parameterizations are calculated has 128 points in the longitudinal direction and 102 points in the latitudinal direction. The vertical structure of the model is represented by 18 unevenly spaced levels using sigma as the vertical coordinate. The dynamical core of the model is based on a former version of the operational medium range forecast model at NCEP (NMC, 1988). The subgrid scale physical parameterizations have been modified from those which were implemented in the original NCEP code.

This version of the model includes parameterizations of solar radiative heating, terrestrial radiative heating, cloud-radiation interaction, deep convection (relaxed Arakawa-Schubert or RAS scheme), large scale condensation, shallow convection, a turbulence closure scheme for subgrid exchange of heat, momentum and moisture, and biophysically controlled interaction between the vegetated land surface and the atmosphere (the simplified SiB model). The simulated time mean and anomalous wind and geopotential height fields were found to be more realistic when using the RAS scheme compared with the previously used Kuo scheme (DeWitt, 1996).

2.2 Experiments, initial conditions, and boundary conditions

This paper describes the results of seasonal hindcasts made using the COLA model for 16 winter seasons (1981/82 through 1996/97). For each winter season, the COLA atmospheric model was integrated from the observed global atmospheric initial conditions produced by the National Centers for Environmental Prediction (NCEP) in mid-December, and the integration was carried through 0000 UTC of 1 April of the following year. The period January, February, March (JFM) is considered for most of the results presented here, and the winter is labeled by the year corresponding to JFM. More specifically, for each winter season, an ensemble of size nine was generated using observed initial conditions for the nine consecutive 12-hour periods beginning on 0000 UTC, 13 December.

The SST and sea ice data used in this study are based on Climate Prediction Center (CPC) weekly 1x1 Optimum Interpolation Sea Surface Temperature (OISST) analyses (Reynolds and Smith, 1994). The data were spatially interpolated by area-weighted averaging to the model Gaussian grid of (x,y) dimension (128,102) and adjusted by the gridded surface topography resulting from the transform of the surface spectral topography used in the model. The adjustment used was -6.5C per kilometer. This gives a closer representation of the corresponding atmospheric temperature in the surface boundary layer. The sea ice locations are set based on the lowest OISST value. The temperature over sea ice points undergoes no height correction. Due to bias corrections added to the OISST, the data were subsequently filtered in time using a 1-2-1 filter (ignoring sea ice points). This time filtered data set is used by the model and is linearly interpolated in time to twice-daily values in all the integrations. All points above -2C are considered open water points, and the value of SST is used for heat flux calculations. Points below -2C are treated as sea ice points with sea ice assumed to extend over the entire Gaussian grid box. The sea ice temperature is computed based on an energy balance calculation and is carried forward in time until SST rises above -2C.

The initial soil wetness was obtained from a climatology calculated from the soil wetness analysis of the European Center for Medium Range Weather Forecasts (ECMWF) analysis-forecast system for 1987-1993 (see Fennessy and Shukla, 1999). After initialization, the COLA model determines the evolution of the soil wetness fields. The operational version of the ECMWF model for that period used three soil layers: the top two layers were prognostic and the bottom layer was prescribed according to time-varying climatological values. Both ECMWF prognostic soil moisture layers, the surface layer with a maximum capacity of 20mm of liquid water and the deep soil layer with a maximum capacity of 120mm of liquid water, were used for initialization. The bottom ECMWF soil moisture layer was not used. Since the land surface parameterization used in the ECMWF model is different from the Simple Biosphere model (SiB) used in the COLA model, the ECMWF soil moisture cannot be used directly. Sato et al. (1989) developed a method to transform soil moisture calculated by other land models to be consistent with SiB. Essentially, the method calculates the time-integral of evaporative demand to which a model grid area would have to be exposed in order to dry down from saturation to a specified level. This same time-integral is then applied to a greatly reduced set of SiB energy balance equations to calculate an equivalent SiB soil moisture level. Soil wetness is the soil moisture content expressed as a fraction of the maximum liquid water capacity for a given grid area. Since the SiB model carries three prognostic soil wetness variables (the surface layer, the root zone and the gravitational drainage zone), the procedure outlined by Fennessy et al. (1994) is used to transform the two combined ECMWF prognostic soil moisture layers into a soil wetness which is used to initialize all three SiB prognostic soil wetness layers.

The definition of the initial snow mass on land points is based on the values of albedo. For albedo values in the range 69%-75%, 49%-69% and 40%-49% the initial snow mass is prescribed to be 20, 10 and 5 kg m-2, respectively. The initial snow depth produced by this procedure is far less than the observed snow depth values. For land points with albedo greater than 75% and land ice points, the initial snow mass is prescribed to be 3x103 kg m-2.



3. MEAN MODEL ERROR

The seasonal means for the January, February and March (JFM) period for all 16 years (1982-1997) were constructed from the twice daily (0000 UTC and 1200 UTC) values from all 144 integrations and compared with observations from NCEP reanalysis for the corresponding 16-year period.

Figure 1a shows the JFM mean model error (generally referred to as the systematic error)

for geopotential height at 500 hPa. The largest errors are found in the polar regions of the

Southern Hemisphere where model tends to produce an anomalous zonally symmetric gradient of geopotential between 40S and the south pole. The largest errors over the Northern Hemisphere are found over the European region. The zonally averaged mean error for height and temperature (not shown) retains the same sign for the entire depth of the troposphere except that, for levels above 500 hPa, the amplitude of error increases with height. The model tends to be much too cold over the poles in the upper troposphere, and, in general, most of the model atmosphere is colder than the NCEP reanalysis. Figure 1b shows the JFM mean error for zonal wind at 200 hPa. The most conspicuous feature is a large mean error of 15ms-1 in the eastern equatorial Pacific area and the northwest part of South America. The magnitude of this error is so large that the observed weak easterlies over the equatorial South America are replaced by spurious westerlies. The main features of the zonal wind error over the northern mid-latitudes reflect a northward displacement of the jet stream in the model compared to the observations. The model simulated core of the jet is also displaced eastward.

Figures 2a and 2b show the model mean JFM error for precipitation and surface (~ 2 m) temperature over land. The observed fields are taken from Xie and Arkin (1996) for precipitation, and Ropelewski et al. (1995) for temperature. The largest precipitation errors (2-5 mm day-1) are found over the regions of largest climatological mean precipitation (equatorial Pacific, South America and southern Africa). Significant errors are also found in the Pacific and the Atlantic storm track regions where the model tends to overestimate precipitation. The model has a cold bias over most of the land areas with the exception of northwest North America and western Eurasia between 40N-55N, where the model can be more than 2-5C warmer than the observations.



4.1 Height and temperature anomaly

Figure 3 shows the composite ensemble mean 500 hPa difference between the averages of three warm ENSO events (1983, 1987, and 1992) and the average of two cold ENSO events (1985, 1989) which occurred during the 16 winter period for the model (top) and the observations (bottom). The model simulates the observed anomalies quite well. There are several other similar model results presented in this volume which clearly show that models have generally improved to such an extent that they can now simulate both the magnitude and the pattern of observed height anomalies realistically. This is in sharp contrast to the earlier models (see Figs. 2 and 4 of Lau, 1997) for which model-simulated anomalies were weaker than the observed anomalies by a factor of two or three. Figure 4 shows the composite 850 hPa temperature anomaly difference for the same two warm and cold events as shown in Figure 3. Both the structure and the amplitude of the temperature anomalies are simulated very well; however, the GCM has a tendency to over-predict the warm anomalies over northeastern North America, and the warmest anomalies in the model are north of their observed location.

4.2 Anomaly Correlation Coefficients

Anomaly correlation coefficients (ACCs) between ensemble mean model anomalies and observed anomalies were computed for every year to verify the model's ability to simulate seasonal atmospheric circulation anomalies. The ACCs are the area-weighted correlations between observed and simulated JFM mean anomaly fields over a specific area. The ACCs for 500 hPa geopotential height are presented in Table 1 for several different regions, some of which overlap. The model shows skill in the simulation of the circulation over the Pacific-North America (PacNA) region, which is defined as 15N-70N, 180W-60W. For the five ENSO events during this period (1983, 1985, 1987, 1989, 1992), the ACC values for the PacNA region are 0.91, 0.49, 0.74, 0.84, 0.64 respectively. There are also some non-ENSO years for which the ACC is greater than 0.5.

In order to determine the geographic dependence of the ability of the ensemble mean, seasonal mean simulations to capture the observed anomaly pattern, we introduce the pattern correlation maps (Straus and Shukla, 1997) shown in Figs. 5 and 6. At each point, the pattern correlation between the anomaly of the ensemble mean JFM simulation and the observed winter average anomaly is calculated over a circular region, centered at the base point, with a radius of 2000 km. The resulting map is averaged for the 16 winters using the Fischer Z-transform, and the result shown by the shading for 500 hPa height (Figure 5a) and 200 hPa zonal wind (Figure 5b). The contours in these figures indicate the16 winter average covariance between simulated and observed anomalies. Significant "skill" (pattern correlation values > 0.6) is seen only in the tropical Pacific, and for the upper level winds in the subtropics off the west coast of North America. That this skill is related to ENSO is confirmed in Figures 6a and 6b , which shows the same correlation maps averaged only over the strong ENSO winters during the DSP period: 1983, 1987, 1992 (warm) and 1985, 1989 (cold). Here correlation values above 0.6 (and for height above 0.8) are seen not only in the eastern tropical Pacific but over much of the PacNA region, with clearly enhanced values of the covariances.

4.3 Low pass and band pass variance

The ability of the ensemble mean to reproduce the observed anomaly pattern in the PacNA region suggests that the GCM is also able to reproduce shifts in the transient fluctuations associated with ENSO. Figure 7 compares the low pass variance of 500 hPa height averaged over

the three warm events (1983, 1987, 1992) from the high resolution operational ECMWF analyses (Fig. 7a) to the same quantity averaged over all ensemble members for the same warm winters (Fig. 7c). The cold event low pass variance (averaged over the winters of 1985 and 1987) given in the same manner in Figs. 7b and 7d. (The low pass variance captures fluctuations with periods of about 10-90 days independent of the annual cycle, which is estimated as in Straus, 1983). For purposes of this calculation, the winter season is defined as the 100-day period starting on 20 Dec. The clear enhancement of the variance seen in the ECMWF analyses in the northeastern Pacific during cold events (as in Palmer, 1988) is also noted in the GCM results, although it is not as dramatic. However, since the GCM variance has been averaged over all ensemble members for each year, it is expected that the individual maxima would be smoothed. Figure 8 shows the corresponding results for the transient meridional heat flux on time scales of about 2 - 10 days (band-pass), except that here it is the anomaly from the 16 winter climatology that is shown. During warm events, there is an equatorward shift of the Pacific storm tracks in both ECMWF analyses and in the GCM (as in Trenberth and Hurrell, 1994 and Straus and Shukla, 1997). The opposite shift is consistently seen in the Atlantic. A poleward shift of the Pacific storm track is seen to occur during cold events in the analyses, and again this is captured by the GCM.

4.4 EOF analysis

A pattern which is characteristic of the 16-winter period is obtained by performing an empirical orthogonal function (EOF) analysis (Bretherton, et al., 1992) of winter mean 500 hPa height for observations ( NCEP reanalyses) and the ensemble mean from the GCM integrations for the 16 winters of 1982 - 1997. The EOF analysis was carried out for the region shown in Fig. 9 (150oE - 30oW; 20oN - 90oN) after the combined hemispheric North Atlantic/ Arctic oscillation had been removed (see the companion paper, Straus and Shukla, 1999). The leading EOF mode explains 47% and 50% of the variance, for model and observations which is far greater than the variance explained by any other mode. The patterns, shown in Figs. 9a and 9b are in generally close agreement, as are the associated principal components (time series) given in Fig. 9c. Minor exceptions to this agreement include the small region off the coast of Newfoundland which is involved in the first mode from the observations but not in the GCM's first mode. Also the winter of 1982 is characterized by a large value of the coefficient in the observations, but not in the GCM.

4.5 SVD for tropical SST and mid-latitude height

To establish that this characteristic pattern is related to the ENSO tropical Pacific SST forcing, a singular value decomposition (SVD) analysis (Bretherton, et al., 1992) was performed linking the 500 hPa height field over the mid-latitude region (150E-30W; 20N-90N) with the tropical Pacific SST field in the region (120oE - 80oW; 20oS - 20oN). Again seasonal (seasonal-ensemble) means were used for observations (the GCM). Figure10 shows the height field patterns and time series associated with the first mode in the same form as in Figure 9. Figure 11 shows the corresponding SST field patterns and time series. The first mode explains 92% of the squared covariance between the height and SST fields for both observations and the GCM. The fact that the two SVD SST patterns are nearly identical, suggests that the nature of the coupling between SST and height is nearly the same in observations and model simulations. The strong resemblance between the first mode pattern for SST (Figs. 11a and 11b) and the first EOF of tropical Pacific SST (not shown), as well as the identity of the extrema in the time series (Fig. 11c) which match the designated SST strong warm and cold events, establish that the SST variation captured by the first SVD mode is that of ENSO. Further, the similarity between the height patterns of the first SVD mode (Figs. 10a and 10b) and that of the first EOF mode (Figs. 9a and 9b) for both observations and GCM indicate that this dominant pattern is linked to the ENSO SST forcing. The time series in Figs. 9 and 10 are also quite similar. The time series of the SST component of the first SVD mode from the GCM and the NCEP reanalyses are nearly identical (Figure 11c).

4.6 Explained variance

The amount of height variance which can be explained by the ENSO signal can be assessed from the heterogeneous correlation pattern for the first SVD mode for the observations and (separately) the GCM. An alternative approach is to calculate the explained variance on the basis of a time series of SST which is based on an observed record longer than 16 years. Figure 12 shows the percentage of seasonal mean variance explained for the 16-year record of observations and the ensemble mean GCM record, from a regression on a time series of SST which emerges from the SVD analysis of a 30-year record (winters of 1958-1997). (Clearly, only the appropriate part of the latter time series was used in the regression.) Locally up to 60% of the variance is explained for the observations (80% for the ensemble mean GCM), and the patterns are not only very similar to each other but also to those in Figs. 9 and 10.

While these results give a strong indication that the ensemble mean of the GCM integrations captures the ENSO SST-forced signal in both pattern and magnitude, it must be recognized that the ensemble mean filters out variations not forced by SST which are due to the sensitivity to initial conditions. That is why the percentage variance explained is higher in the GCM than in observations (Fig. 12). These chaotic fluctuations are, of course, still present in the observations. It may be that individual realizations of the GCM have a larger level of chaotic fluctuations compared to the ENSO SST-forced signal than is realistic, and that only by ensemble averaging can we so clearly distinguish this signal.



5.1 Comparison between model ensembles and a single observed realization

In order to assess the strength of the ENSO signal in single records of 16 seasonal GCM integrations, comparable to the single record available from observations, we have prepared "samples" of GCM integrations. Each sample consists of one integration per calender year, chosen randomly out of the ensemble for each year. (Note that for an ensemble size of N over 16 years there are N16 possible samples.) We can then repeat the 16-year SVD analysis of 500 hPa height and tropical Pacific SST on each such sample. Figure 13a shows the distribution of pattern correlations for the first SVD mode height field between each of a set of 100 samples and the observations. The distribution, presented as a probability distribution function using a Gaussian kernel estimation method (Silverman, 1986), is strongly peaked at about 0.88, and indicates little probability of the GCM/observed pattern correlation falling below 0.80. Similarly, Fig. 13b shows the probability distribution functions for the squared covariance (height / SST) explained by the first SVD mode obtained from the 100 GCM samples (solid line), and for the correlation between the first mode time series of SST and height. The peak in squared covariance explained is at 87%, only slightly lower than the 92% explained by GCM ensemble means, with a very small probability of a sample indicating less than 80% explained covariance. Although the time series correlations are somewhat lower, very few samples give a correlation of less than 0.75.

While the occurrence of ENSO clearly affects intra-seasonal fluctuations (as in Figs. 7 and 8), it is not clear whether the characteristic seasonal mean ENSO pattern defined above leaves a consistent imprint on intra-seasonal fluctuations. Put another way, is the projection of individual 5-day means (pentads) on the SST forced pattern of Fig. 9 predominantly of one sign during a strong ENSO winter, or is that pattern a statistical residue of large episodes of both signs? Figure 14 shows the probability distribution functions (pdfs) for these projections for all GCM pentads (14a) and all observed pentads (14b). In each case the pdf for all pentads is given by the thin line, while those for the winter of 1983 (1989) are given by the thick solid (dotted) line. Clearly, the projections for the extreme warm (cold) ENSO event are predominantly positive (negative). The positive (negative) seasonal mean coefficient for 1983 (1989) in Fig. 9c is thus consistently seen in nearly all pentads. The shifts in the pdf for warm and cold events is less in calculations by Renshaw et al. (1998), perhaps because their EOF was calculated using low frequency filtered daily data.

5.2 Spread within an ensemble

In order to make a simple estimate of spread among each of the nine members of the ensemble for each year, the ACC over the PacNA region for 500 hPa height anomalies was calculated among all 36 possible pairs of ensemble members. Figure 15 shows the frequency distribution of these 36 values of ACC for each year separately. Each panel also contains (at the top) the ACC between the ensemble mean height anomalies and the corresponding observations. It is seen that for nine cases whose intra-ensemble ACC is predominately positive (1983, 1985, 1987, 1988, 1989, 1990, 1992, 1995, 1996), the ensemble mean has a high value of ACC with respect to observations. This implies that the spread among the ensemble members is a good predictor of the skill of ensemble mean. However, there are few notable exceptions. For example, for 1991, most of the members of the ensemble have high ACC values among themselves, but the ensemble mean has a low ACC with respect to the observations. A comparison of individual ensemble members with observations reveals that even when the ensemble mean has a low ACC, some members of the ensemble have high ACC . What needs to be further explored is the question of whether increasing the ensemble size (either by making a larger number of integrations with the same model or by combining integrations from several models), would lead to a better estimate of the probability distribution of predicted states.

5.3 Brier score and reliability diagram

In the context of ensemble predictions, it is of interest to consider skill scores which measure the occurrence (or non-occurrence) of an "event", e.g., seasonal mean 500 hPa geopotential height Z being one standard deviation greater (or less) than normal. The probability forecast for this happening at each gridpoint in the GCM can be simply given by the fraction of ensemble integrations, while for observations such an event either happens or does not happen (probability of 1 or 0). The observed frequency of events over the Northern Hemisphere (NH) is plotted against the forecast probability in Fig. 16a (height below minus one standard deviation) and Fig. 16b (height above one standard deviation). For these calculations, 17 years (1982-1998) were used. Similar diagrams for events over North America are plotted in Fig. 16c and Fig. 16d. Such "reliability diagrams" (see Stanski et al., 1989 for details) are used to give a graphical depiction of probability forecast performance. All of the corresponding Brier skill scores (Stanski et al., 1989; Palmer et al., 1999) are positive and statistically significant at the 99% level (based on a Monte Carlo approach), suggesting that the skill of ensemble forecast has higher skill than a climatological probability forecast (see the figure caption). When all gridpoints in the NH are counted, the below-normal events are significantly over-forecast (Fig. 16a), leading to a lower Brier skill score than that of the above-normal events (Fig. 16b). Over North America, the reliability curves lie closer to the "perfect reliability" 45° line, with correspondingly higher skill scores than their NH counterparts. It thus suggests that the model performs better over North America than it does over the whole NH domain. By comparing the Brier skill scores for the ENSO years (i.e., 1983, 1987, 1992, 1998 for warm events and 1985, 1989 for cold events) with the normal years, it is seen that the model performs better during ENSO. The high skill scores over the North America, therefore, are likely to be attributed to the ENSO SST forcing.

6. SUMMARY AND DISCUSSION

The COLA model is successful in simulating seasonal mean height anomalies over the PacNA region in the presence of large tropical SST anomalies. This is not a unique property of this model, but this is a property of the atmosphere because as shown in Shukla et al, (1999) and several other papers in this volume, several models have produced similar simulations. However, in the absence of significant SST anomalies, there is a large spread in the values of ACC among the members of an ensemble and there is little, if any, predictability of extratropical seasonal mean anomalies. The SVD analysis linking 500 hPa height field over the PacNA region and the tropical SST clearly suggests the predominant role of tropical forcing in producing mid-latitude circulation anomalies.

In the presence of large SST anomalies, the spread in the values of ACC among the members of each ensemble is small, suggesting that at least in this limited number of cases, if SST were predicted accurately, the winter seasonal mean over PacNA region could be predicted with a high degree of confidence. The prospects for prediction should not be considered entirely hopeless even for the years when SST anomalies are not large, because models still have quite large systematic errors and the probability distributions of predicted states have not been examined in ensembles of large size (50-100).

A preliminary analysis of summer season hindcasts (not shown) with this model has revealed that the mean model errors in rainfall are far larger than the observed anomalies of rainfall especially over the land areas. In fact, even with prescribed global observed SST, the model could not simulate the rainfall difference over the U.S. even for two years of opposite extreme rainfall anomalies: the drought in 1988 and the floods in 1993. There is a clear need to improve models to simulate mean summer rainfall over land areas and only then it can be determined whether summer rainfall anomalies have any predictability.

Acknowledgements

This research was supported by NSF grant ATM-93-21354, NOAA grant NA-76-GP0258 and NASA grant NAGW-5213. We gratefully acknowledge the computing resources provided by the Climate Simulation Laboratory.



7. REFERENCES

Bengtsson, L., K. Arpe, E. Roeckner, and U. Schulzweida, 1996: Climate predictability experiments with a general circulation model. Clim. Dyn., 12, 261-278.

Bengtsson, L., U. Schlese, E. Roeckner, M. Latif, T.P. Barnett, and N. Graham, 1993: A two-tiered approach to long range climate forecasting. Science, 261, 1026-1029.

Brankovic, C., T.N. Palmer, and L. Ferranti, 1994: Predictability of seasonal atmospheric variations. J. Climate, 7, 217-237

Bretherton, C. S., C. Smith, and J.M. Wallace, 1992: An intercomparison of methods for finding coupled patterns in climate data. J. Climate, 5, 541-560.

Brier, G.W., 1950: Verification of forecasts expressed in terms of probabilities. Mon. Wea. Rev., 78, 1-3.

Charney, J.G., and J. Shukla, 1981: Predictability of monsoons. Symposium on Monsoon Dynamics, New Delhi, December 1977. Monsoon Dynamics, Cambridge University Press, Sir James Lighthill and R.P. Pearce (eds.), 99-109.

Chervin, R.M., 1986: Interannual variability and seasonal climate variability. J. Atmos. Sci., 43, 233-251.

DeWitt, D.G., 1996: The effect of cumulus convection on the climate of COLA general circulation model. COLA Tech. Report #27. COLA, 4041 Powder Mill Road, Suite 302, Calverton, MD, 20705.

Fennessy, M. J., and J. Shukla, 1999: Impact of initial soil wetness on seasonal atmospheric prediction. J. Climate, to be published.

Fennessy, M. J., J. L. Kinter III, L. Marx, E. Schneider, P. J. Sellers, and J. Shukla, 1994: GCM simulations of the life cycles of the 1988 drought and heat wave. COLA Tech. Report No. 6, 53 pp. COLA, 4041 Powder Mill Road, Suite 302, Calverton, MD, 20705.

Gates, L., and Coauthors, 1999: An overview of the results of the atmospheric model intercomparison project (AMIP I). Bull. Amer. Meteor. Soc., 80, 29-55.

Horel, J.D., and J.M. Wallace, 1981: Planetary-scale phenomena associated with the Southern Oscillation. Mon. Wea. Rev., 109, 813-829.

Kinter III, J.L., D.G. DeWitt, P.A. Dirmeyer, M.J. Fennessy, B.P. Kirtman, L. Marx, E.K. Schneider, J. Shukla, and D.M. Straus, 1997: The COLA atmosphere-biosphere general circulation model. Volume 1: formulation. COLA Tech. Report, 51, 46pp. COLA, 4041 Powder Mill Road, Suite 302, Calverton, MD, 20705.

Kumar, A., and M.P. Hoerling, 1997: Interpretation and implications of the observed inter-El

Nino variability. J. Climate, 10, 83-91.

Kumar, A., M. P. Hoerling, M. Ji, A. Leetma, and P. D. Sardeshmukh, 1996: Assessing a GCM's suitability for making seasonal prediction. J. Climate, 9, 115-129.



Lau, N. C., 1997: Interactions between global SST anomalies and the mid-latitude atmospheric circulation. Bull. Amer. Meteor. Soc., 78, 21-33.

Lau, N.C., 1981: A diagnostic study of recurrent meteorological anomalies appearing in a 15 year simulation with a GFDL general circulation model. Mon. Wea. Rev., 109, 2287-2311.

Miyakoda, K., J. Sirutis, and J. Ploshay, 1986: One month forecast experiments-without boundary forcings. Mon. Wea. Rev., 114, 2363-2401.

Miyakoda, K., T. Gordon, R. Caverly, W. Stern, J. Sirutis, and W. Bourke, 1983: Simulation of a blocking event in January 1977. Mon. Wea. Rev., 111, 846-869

Moorthi, S. and M.J. Suarez, 1992: Relaxed Arakawa-Schubert: A parameterization of moist convection for general circulation models. Mon. Wea. Rev., 120, 978-1002.

NMC Development Division Staff, 1988: Research version of the medium range forecast model. NMC Documentation Series #1. Development Division, NMC, Washington, D.C., 20233.

Palmer, T. N., 1988: Medium and extended-range predictability and the stability of the Pacific- North America mode. Quart. J. Roy. Meteor. Soc., 114, 691 713.

Palmer, T.N., C. Brankovic, and D.S. Richardson, 1999: A probability and decision-model analysis of PROVOST seasonal multi-model ensemble integrations. Quart.J.Roy. Meteor. Soc. In preparation.

Renshaw, A.C., D.P. Rowell, and C.K. Folland, 1998: Wintertime low-frequency weather variability in the North Pacific-American sector 1949-93. J. Climate, 11, 1073-1093.

Reynolds, R.W., and T.M. Smith, 1994: Improved global sea surface temperature analyses using optimum interpolation. J. Climate, 7, 929-948.

Ropelewski,C.F., J.E. Janowiak, and M.F. Halpert, 1995: The analysis and display of real time surface climate data. Mon. Wea. Rec., 113, 1101-1107.

Sato, N., P. J. Sellers, D. A. Randall, E. K. Schneider, J. Shukla, J. L. Kinter III, Y.-T. Hou, and E. Albertazzi, 1989: Effects of implementing the Simple Biosphere Model in a general circulation model. J. Atmos. Sci., 46, 2757-2782.

Shukla, J., and Coauthors, 1999: Dynamical seasonal prediction. Bull. Amer. Meteor. Soc.

Shukla, J., 1981a: Predictability of the tropical atmosphere. NASA Tech. Memo. #83829, NASA Goddard Space Flight Center, Greenbelt, MD, USA, 51pp.

Shukla, J., 1981b: Dynamical predictability of monthly means. J. Atmos. Sci., 38, 2547-2572.

Silverman, B. W., 1986: Density estimation for statistics and data analysis. Chapman and Hall, 175pp.

Stanski, H.R., L.J. Wilson and W.R. Burrows, 1989: Survey of common verification methods in meteorology. World Weather Watch Tech. Report #8. (WMO/TD, #358), WMO, Geneva, 114pp.

Straus, D. M., 1983: On the role of the seasonal cycle. J. Atmos. Sci., 40, 303-313.

Straus, D. M., and J. Shukla. 1997: Varations of midlatitude transient dynamics associated with ESNO. J. Atmos. Sci., 54, 777-790.

Straus, D. M., and J. Shukla, 1999: Distinguishing the mid-latitude response to ENSO from internal variability: Analysis of observations and GCM simulations. Quart. J. Roy. Meteor. Soc.

Trenberth, K. E., and J. W. Hurrell, 1994: Decadal atmosphere-ocean variations in the Pacific. Climate Dyn., 9, 303-319.

Xie, P. and P. Arkin, 1996: Analysis of global monthly precipitation using gauge observations, satellite estimates, and numerical model predictions. J. Climate, 9, 840-858.

8. FIGURE CAPTIONS

Figure 1. JFM mean model error for 144 seasonal integrations for (a) 500 hPa height (meters) and (b) 200 hPa zonal wind (ms-1).

Figure 2. JFM mean model error for 144 seasonal integrations for (a) precipitation (mm day -1) and (b) surface (2 m) temperature (C).

Figure 3. Ensemble mean JFM 500 hPa height difference between the average of 1983 and 1992, and the average of 1985 and 1989. (a) GCM difference ( b) Observed difference

Figure 4. Ensemble mean JFM 850 hPa temperature difference between the average of 1983 and 1992, and the average of 1985 and 1989. (a) GCM difference (b) Observed difference

Figure 5. (a) Maps of spatial correlation and covariance between the JFM mean observed and ensemble mean COLA GCM 500 hPa height for 16 winters (1982-1997). The correlation (covariance) map gives a sliding spatial correlation (covariance) between two fields over a region of radius 2000 km centered at the given point (see text for details). Spatial correlation (covariance) is given by the shading (contours, interval of 102m2. (b) As in (a) but for 200 hPa zonal wind. Contour interval of 3 m2s-2.

Figure 6. (a) Maps of spatial correlation and covariance between the JFM mean observed and ensemble mean COLA GCM 500 hPa height (a) and 200 hPa zonal wind (b) for the strong ENSO winters of 1983, 1985, 1987, 1989, 1992. Otherwise as in Fig. 5.

Figure 7. Low pass variance of 500 hPa height for: (a) warm ENSO winters from ECMWF analyses; (b) cold ENSO winters from ECMWF analyses; (c) Warm ENSO winters from COLA GCM ensemble average; (d) Cold ENSO winters from COLA GCM ensemble average. See text for details. Contour interval of 3 x 103 m2.

Figure 8. Band pass variance anomaly of 850 hPa transient meridional heat flux for: (a) warm ENSO winters from ECMWF analyses; (b) cold ENSO winters from ECMWF analyses; (c) Warm ENSO winters from COLA GCM ensemble average; (d) Cold ENSO winters from COLA GCM Ensemble average. See text for details. Contour interval of 1 oK ms-1.

Figure 9. Leading empirical orthogonal functions of 500 hPa height (in the region shown) for 16 winters for: (a) seasonal mean observations (NCEP reanalyses) and (b) COLA GCM ensemble seasonal mean. Percent variance explained given in parentheses (see text for details). (c) The principal component coefficients associated with the patterns in (a) and (b) as a function of time. For any winter, using the coefficient to scale the pattern gives the height in meters.

Figure 10. Leading mode of singular value decomposition of 500 hPa height in the region shown with tropical Pacific SST for 16 winters for: (a) seasonal mean observations (NCEP reanalyses) and (b) COLA GCM ensemble seasonal mean. Percentage squared covariance explained given in parentheses (see text for details). (c) The first mode coefficients associated with the patterns in (a) and (b) as a function of time. For any winter, using the coefficient to scale the pattern gives the height in meters.

Figure 11. Leading mode of singular value decomposition of tropical Pacific SST in the region shown with 500 hPa height in the Pacific-North America region for 16 winters for: (a) seasonal mean observations (NCEP reanalysis) and (b) COLA GCM ensemble seasonal mean. Percentage squared covariance explained given in parentheses (see text for details). (c) The first mode coefficients associated with the patterns in (a) and (b) as a function of time. For any winter, using the coefficient to scale the pattern gives the SST in oK.

Figure 12. Regression of seasonal mean (Jan-Mar) 500 hPa height for 16 winters (1982 - 1997) on time series of tropical Pacific SST. The time series was derived from singular value decomposition of tropical Pacific SST with 500 hPa height for 30 winters (1968 - 1997). (a) Variance explained for observations (NCEP reanalyses); (b) Variance explained by the GCM ensemble seasonal means. Light (dark) shading for values above 60% (80%).

Figure 13. (a) Probability distribution functions of pattern correlation between leading singular value decomposition (SVD) patterns for observations (NCEP reanalyses) and 100 GCM samples, each sample being obtained by randomly choosing one GCM integration per calendar winter. The SVD analysis used seasonal mean (Jan-Mar) tropical Pacific SST and 500 hPa height in the Pacific / America region for 16 winters (1982-1997). See text for details. (b) Probability distribution functions for the GCM squared covariance explained (solid) and correlation between time series of height and SST for the first SVD mode (dash) for 100 GCM samples.

Figure 14. Probability distribution functions (pdfs) for projection of pentad means on the leading empirical orthogonal function (EOF) of (a) ensemble seasonal means for the COLA GCM and (b) seasonal means for observations (NCEP reanalyses) for 16 winters. The thin line is the pdf of all pentad means for 16 winters, the heavy (dotted) line is the pdf for the winter of 1983 (1989). (see text for details).

Figure 15. ACC along x-axis among the 36 pairs from nine ensemble members for 500 hPa height for each year. Percentage of total pairs along y-axis. ACC between the ensemble mean model and observed 500 hPa height anomalies at the top (left) of each box denoted as skill.

Figure 16. Reliability diagram and frequency distribution for the COLA DSP 17-year ensemble forecast of the seasonal mean (JFM) 500 hPa height Z. The horizontal axis gives the forecast probability and the vertical axis gives the observed frequency. The probability density function (pdf) is plotted as histograms. Two events are forecast: 1) Z is one standard deviation lower than normal, denoted as Z<(-1), and 2) Z is one standard higher than normal, denoted Z>(+1). The upper panel shows results derived from all gridpoints in the NH (NH, 0N-90N, 0E-360E) and the lower panel from gridpoints in North America (20N-70N, 125W-65W). Also shown in the top left hand corner of each diagram is the Brier skill score B (BENSO for ENSO years and Bnormal for normal years; score enclosed in parenthesis does not pass significance test at the 99% level), defined as B= (bcli- b)/bcli, where b is the Brier score (Brier 1950) and bcli is the Brier score for a climatological forecast. We define the climatological probability of an event at a gridpoint for each year as the climatological frequency of occurrence of the event, calculated from the observations of all other years at the same gridpoint, and a climatological forecast of the event as the forecast in which the climatological probability is predicted.