** Section 5 on Potential Predictability and Model Validation was collaborative
research with Dr Xiaogu Zheng (NIWA) **

1. Introduction

2. Statistical Forecasting of Sea Surface Temperatures

    2.1 Datasets

    2.2 Analysis Techniques
        2.2.1 Spectrum Analysis
        2.2.2 AR Processes

    2.3 Results
        2.3.1 M-SSA of Sea Surface Temperature EOFs (5° Grid)
        2.3.2 M-SSA of 2°Sea Surface Temperature EOFs
        2.3.3 M-SSA of Weekly Sea Surface Temperature EOFs
        2.3.4 AR Process Modelling

    2.4 Discussion and Conclusions

3. Dynamical Seasonal Forecasts During the 1997/1998 El Nino

    3.1 Model and Experiments

    3.2 Skill Scores – LEPS

    3.3 SSTA Forecasts

    3.4 Climate Forecasts
        3.4.1 Rainfall
        3.4.2 Surface Air Temperature
        3.4.3 200hPa Geopotential Height
        3.4.4 Mean Sea Level Pressure

    3.5 Discussion and Conclusions

4. AGCM Simulation of Climate 1958-1991

    4.1 Model and Experiments

    4.2 Simulation Skill Scores

    4.3 SOI Relationships

    4.4 Discussion and Conclusions

5. Potential Predictability and Model Validation

    5.1 Results from Monthly Values (1950-1991)

    5.2 Results from Daily Values (1958-1991)

        5.2.1 Potential Predictability
        5.2.2 Correlation and Model Validation

6. Conclusions

REFERENCES

1. Introduction

The impact of global sea surface temperatures (SSTs), and in particular tropical SSTs, on global weather and climate systems has been recognised for some time. Observational studies have identified areas around the world where climate is strongly influenced by the El Nino/Southern Oscillation (ENSO) phenomenon (e.g. Ropelewski and Halpert, 1987; Philander, 1990). Over the Australian region, a number of empirical studies have provided evidence of relationships between the interannual variability of Australian rainfall and ENSO (e.g. Pittock, 1975, 1984; McBride and Nicholls, 1983; Nicholls, 1988; Whetton, 1988; Drosdowsky, 1993b; Ropelewski and Halpert, 1987, 1989) and sea surface temperature anomalies (SSTAs) in the east equatorial Pacific and around Australia (Streten, 1981, 1983; Whetton 1990b; Nicholls, 1984, 1989; Drosdowsky, 1993b). Some studies have also examined the relationship between this temporal variability and large scale circulation anomalies over the Australian region (e.g. Pittock, 1975, 1984; Streten, 1981; Whetton, 1988, 1990a; Wright, 1988a,b; Nicholls, 1989; Drosdowsky, 1993b; Ropelewski and Halpert, 1987, 1989).

Pittock (1975, 1984), Nicholls (1989) and Drosdowsky (1993a,b) showed that the interannual variability of Australian winter rainfall was dominated by two spatial patterns. One was centred over the eastern third of the continent and the other a broad northwest-southeast band with maximum in the south/southeast and associated with extensive cloudbands, with origin off the northwest coast, commonly referred to as Australian northwest cloudbands (NWCBs). The former pattern is related to Pacific SSTs and the latter to the gradient in SST between the Indonesian archipelago and the central Indian Ocean. Such relationships have, with varying success, also been confirmed in

atmospheric general circulation models (AGCMs) forced by idealized SST anomalies (e.g.Voice and Hunt, 1984; Simmonds and Smith, 1986; Simmonds and Trigg, 1988; Simmonds et al., 1989; Simmonds 1989, 1990; Simmonds and Rocha, 1991; Frederiksen and Balgovind, 1994; Frederiksen and Frederiksen, 1996; Frederiksen et al., 1999).

Frederiksen et al. (1999) showed that a version of the Bureau of Meteorology Research Centre (BMRC) AGCM, forced by SSTs for the period 1950-1991, was able to simulate much of the observed Australian rainfall variability, especially over the east and southeast of the continent. In particular, the model was able to capture the observed SOI relationship and the dominant patterns of seasonal and interannual variability, and their relationships with global SSTs. In addition, the model captured correctly the associated global anomalous circulation patterns. These results suggest that there might be some useful forecast skill from this model, provided reasonably good forecasts of the SST forcing can be generated.

Recently, statistical approaches, based upon these relationships between Australian rainfall and the Southern Oscillation Index (SOI) and global SSTAs, have been developed to make seasonal forecasts (Drosdowsky and Chambers, 1998). Similar methods have also been developed for forecasting surface temperatures (Jones, 1999). These techniques have shown significant forecasting skills over some parts of the Australian continent. However, the impact of the SST forcing on the weather variability is dependent on the position, timing and strength of these boundary forcings. For this reason, it is important to investigate the feasibility of employing AGCMs in providing dynamical seasonal forecasts.

Theoretically, dynamically extended seasonal forecasts, beyond the limit of predictability

of individual synoptic weather systems (typically 10 – 15 days), are possible because atmospheric boundary forcings, such as, for example, SST and soil moisture, which evolve on slower time scales than that of the synoptic systems, largely determine the averaged weather features on a seasonal (1 – 3 month averages) time scale. However, even on this time scale and for seasonal averages, the climate system is still inherently chaotic and it is of interest to determine how much of the climate variability, at different geographical locations, is in fact predictable. A number of studies have considered the question of potential long-range predictability of climate variability. These have involved ensembles of multidecadal simulations using atmospheric general simulation models (AGCMs) with observed (SSTs) and different initial conditions (e.g., Dix and Hunt 1995; Harzallah and Sadourny 1995; Kumar et al. 1996; Rowell et al. 1995; Stern and Miyakoda 1995; Rowell 1998; Zwiers 1996; Zheng and Frederiksen, 1999). With such an approach, it has been possible to separate interannual variability into chaotic components (due to the sensitivity to initial conditions) and a potentially predictable component (based on the ensemble average). Thus, it is possible to identify areas where different climate variables are potentially predictable.

In this study, we investigate the possible use of sophisticated atmospheric climate models in forecasting seasonal climate anomalies one to three months ahead. Such models require the sea surface temperature (SST) as a lower boundary forcing. On this time scale, statitistical models for SST prediction are expected to do better than dynamical models, and therefore we consider in section 2 a number of possible statistical approaches. In section 3, we present results from trial seasonal forecasts conducted during the 1997/98 El Nino. Results from multidecadal simulations of the climate during 1950-1991 are presented in section 4, and in section 5 we consider the question of potential predictability and model validation using correlation analysis. Our conclusions are presented in section 6.

2. Statistical Forecasting of Sea Surface Temperatures

In order to develop a dynamical seasonal forecasting scheme using an AGCM, we need first to predict the SSTs which will act as the lower boundary forcing for the atmospheric model. Given the SST climatology, it is only necessary to predict the SST anomaly (SSTA) from the climatological mean value. On a time scale of 1 – 3 months ahead, statistical models for SST prediction are likely to outperform model based SST predictions. Barnston et al. (1999) discuss the predictive skill of a number of dynamically based SST prediction models during the 1997/98 El Nino episode. Here, we investigate two statistical methods for predicting the SSTAs. One is based on a multichannel singular spectrum analysis (M-SSA) and the other on time series analysis using a low order autoregressive technique. The usefulness of these methods is assessed by comparing forecasts with SSTAs obtained using persistence, that is, using the latest observed SSTAs and holding them fixed throughout the forecast period.

2.1 Datasets

Two datasets were considered in this study, comprising a monthly and a weekly set. For the initial analysis SSTA data was derived from the UK Meteorological Offices GISST 1.1 data set (Drosdowsky and Chambers 1998). The monthly data from January 1949 to December 1991 was used. This data was extended, from December 1991 to July 1997 using the National Centre for Environmental Prediction (NCEP) optimum interpolation analysis, which is available from 1982 to the present. Both of these data sets have a 1° by 1° spatial resolution and were reduced to a 5° by 5° grid covering the region from 55° S to 60° N and 30° E to 70° W. Analyses were also performed on a 2° by 2° grid and using the region from 60° S to 60° N and 30° E to 70° W.

However, because the derived forecasts lacked some of the smaller scale features seen in the observed, it was decided to also evaluate these techniques using weekly SSTAs obtained from the weekly NCEP SST optimum interpolation analysis. This data set covers the period from November 1981 to the present. Again this data was converted from a 1° by 1° grid to a 5° by 5° grid covering the same spatial domain 60° S to 60° N and 30° E to 70° W.

2.2 Analysis Techniques

Both the monthly and weekly SST fields were analysed using rotated S-mode principal component analysis as described in Drosdowsky (1993c) using IMSL subroutines. The principal component analysis was primarily performed as a data reduction measure. For a discussion on the use of rotated principal components see Drosdowsky and Chambers (1998). Twelve components were retained, for the monthly SST data, explaining around 46% of the variance. In the case of the weekly SST data 13 components were retained which explain about 46% of the variance. In both cases the first principal component represents the mature phase of an El Niño/Southern Oscillation (ENSO) event.

Two methods of forecasting SSTs were compared. The first involved spectrum analysis and the second low order autoregressive time series analysis.

2.2.1 Spectrum Analysis

Spectrum analysis was developed to evaluate the non-randomness of a time series without assuming the nature of the non-randomness. The analysis is based on the assumption that the time series does not necessarily consist of a finite number of oscillations with discrete wavelengths, but that it is comprised of practically infinite numbers of small oscillations covering a continuous distribution of wavelengths. The resulting spectrum gives a measure of the distribution of variance in the time series for all possible wavelengths. The longest wavelength is the linear trend (infinite wavelength) and the shortest obtainable is equal to twice the interval between successive observations in the time series. Here, because there are a number of input variables, or time series, to be analysed, we consider the use of multichannel singular spectrum analysis (M-SSA) (see Plaut and Vautard, 1994, for further details on the method of multichannel singular spectrum analysis).

Multichannel singular spectrum analysis (M-SSA) is a method of identifying coherent space-time patterns when given a regularly sampled archive of maps. Plaut and Vautard (1994) comment that M-SSA is mathematically equivalent to the extended EOF analysis of Weare and Nasstrom (1982). However, M-SSA differs from an extended EOF analysis as it uses a large number of lags in the state vectors, which are used to obtain the spectral properties.

To obtain the period of each EOF output from the spectrum analyses the maximum-entropy method was used. The maximum-entropy method (MEM) is designed to choose the spectrum that corresponds to the most random, most unpredictable, time series whose autocorrelation function agrees with a set of known values. That is, the autocorrelation function of the available time series is extrapolated by maximising the entropy of the process. The entropy is defined as a measure of the average information content. Hence, unlike spectrum analysis, MEM does not assume a periodic extension of the data or that the data outside the available record length is zero.

An important issue in the maximum-entropy method is the choice of filter order, M. It has generally been observed that, for a given record length N, small values of M can give spectral estimates which have insufficient resolutions, whereas large values of M yield estimates that are statistically unstable with spurious details (Haykin 1979). The optimal value of the filter order used in this study was obtained using the ‘final prediction-error (FPE) criterion’. According to Haykin (1979) this method generally gives slightly better results than the Information Theoretic Criterion (AIC) and Autoregressive Transfer Function Criterion (CAT) (see Haykin 1979 p42 for definitions of these criteria). M_opt is then defined as the value of M corresponding to the minimum value of FPE(M) where

and S_M² is the residual sum of squares (see Ulrych and Bishop 1975 for computational method).

2.2.2 AR Processes

Time series analysis was only performed on the weekly SSTA data.

Each of the PC time series of the weekly sea surface temperature anomalies were forecast at 1, 2 and 3 week lead times using the IMSL procedure ARMME. The PC time series are stationary since they were originally constructed from anomalies with zero means. The model is given as:

Where p is the order of the model, determined by the MEM described in section 2.2.1. Method of moments parameter estimates of b₁ to b_p were obtained from the subroutine ARMME, which uses the extended Yule-Walker equations, and the equation above was used to obtain forecast values for each of the SSTa PCs.

2.3 Results

2.3.1 Multichannel Singular Spectrum Analysis (M-SSA) of Sea Surface Temperature (SST) EOFs (5° Grid)

The SST data set, consisting of the monthly time series of the first twelve rotated SST eigenvectors, had dimensions 12 by 516 (the number of months of data) (see Figure 4 of Drosdowsky and Chambers, 1998). As we are particularly interested in shorter term trends, the number of lags in the M-SSA was set at 50. This gave a new variable of dimensions 466 by 600 (516-50 by 12x50). The M-SSA consisted of essentially performing a principal component analysis on this new data set.

A plot of the resulting percentage of variance explained by the eigenvectors indicated that only 10 of the eigenvectors needed to be retained (Figure 2.1).

From here on we will refer to the original 12 principal component EOFs are S-PCs (spatial principal components) and the EOFs obtained from the M-SSA as ST-PCs (space-time principal components).

Reconstructing the S-PCs using the EOFs obtained from the M-SSA, ST-PCs, for the corresponding channel results in smoothed versions, containing less noise. In most cases the reconstructed data from the channels represents the S-PCs fairly well. This is particularly true for the reconstructions of the first channel.

By plotting pairs of eigenvectors, over time, for each channel it is possible to look for periodicity within the eigenvectors. For most ST-PC pairs the amplitude in channel 1 is larger than in other channels; hence the associated oscillation is dominated by fluctuations of the first S-PC, which is mostly the ENSO pattern. Channels 2 and 6 appear to represent trends in the data.

Table 2.1 gives the percentage of variance that each of the pairs of ST-PCs explains. Also given in Table 2.1 are the periods of the oscillations which were obtained using Maximum Entropy Spectral Analysis. Again, we used the Burg method.

For ST-PC pairs 1 & 2 and 5 & 6 the peak frequency occurred at zero. For these pairs it was not possible to resolve the period over this time scale (50 years).

ST-PCs 3, 4, 7, 8, 9 and 10 were used to predict SST time series patterns five months ahead. Correlations between the forecast values and the actual time series for the same period were computed for each of the channels. For most channels the correlations were quite low and it was decided that more ST-PCs probably need to be retained. In all further analyses the number of ST-PCs used was set to 30. This resulted in the following pairs: 3 & 4, 7 & 8, 10 & 11, 29 & 30, and the triple 25, 26, 27. Reconstructions of the time series were based on ST-PCs 1-8, 10, 11, 25-27, 29 and 30 since ST-PC 1 represented the mean, ST-PCs 2, 5, and 6 represented trends and the other ST-PCs were pairs or triples. The forecast time series now fitted the original time series fairly well, particularly for the first two channels. However, some of the smaller scale patterns in the SSTs were not picked up very well when the time series for the channels were recombined to give spatial patterns. It was thought that more detail could be obtained by using 14 rotated global S-PCs initially.

The new spatial region covers all areas between 60° North and 60° South. After re-performing the M-SSA we were left with 19 ST-PCs for use in the reconstructions of the original SST patterns. Forecasts of 1, 2, 3, 4 and 5 months ahead were compared with persistence using RMS errors. It was found that, if the forecast period was only out to 2 months ahead that there was no really advantage in using M-SSA over persistence. However, once the forecast lead time is greater than 3 months M-SSA tends to out perform persistence. Figure 2.2 shows the bias and root mean square error (rmse) at 3 months lead, when both have comparable rmse. At shorter lead times the M-SSA has larger bias and rmse.

The above method of forecasting SSTs appeared to perform worst in regions where the original SST EOF patterns were quite large and extensive, for example in the ENSO region of EOF1. It was thought that not enough fine scale patterns were being picked up when using only 12 EOFs and 5° gridded data. To attempt to retain the fine scale patterning a 2° analysis was performed.

2.3.2 M-SSA of 2°Sea Surface Temperature (SST) EOFs

Due to the large number of 2 (spatial principal components) gridded data points the principal component analysis of the 2° SST temperature data needed to be performed on subsets of the full data set. Three such regions were selected; the Atlantic Ocean, the Pacific Ocean and the Indian Ocean. The PCA analysis resulted in retaining 6 EOFs for the Atlantic Ocean and 8 each for the Pacific and Indian Oceans, and are shown in Figure 2.3, Figure 2.4 and Figure 2.5. Comparing Figures 2.3-2.5 and Figure 4 of Drosdowsky and Chambers (1998) we see that a number of the patterns are contained in both analyses. For example, the first Pacific EOF is present in the first Global EOF (see Figure 4 of Drosdowsky and Chambers, 1998), as is the first Indian EOF. Note that we have excluded a region in the Pacific Ocean where the data is considered to be unreliable before the use of satellite data in the 1970's.

A separate M-SSA was conducted for each of the ocean regions. The number of ST-PCs used by the Atlantic, Indian and Pacific regions were 9, 12 and 14 respectively. Forecasts three months ahead were made and the data recombined to form a global picture. Tested over the same period as used in previous analyses the M-SSA generated forecasts had lower rmse over most of the globe than a persistence forecast. This was particularly true for the Indian Ocean region. Any improvement over the previous global analysis was not so clear cut.

The use of the 2° SST data set for the prediction of Australian rainfall is limited at present due to the large number of resulting EOFs (6 for the Atlantic, and 8 each for the Indian and Pacific Ocean regions) and the lack of reliable past data. There are just too many predictors and not enough data.

Another method of increasing the amount of small scale patterns but not the resolution is to look at weekly data rather than monthly data.

2.3.3 M-SSA of Weekly Sea Surface Temperature EOFs

1° Weekly SSTA patterns were converted to a 5° by 5° grid before performing a PCA. The first 13 rotated components were retained for use in further analysis. Despite using a number of different lags there was not really enough years of data to isolate the oscillatory patterns properly, indicating that the use of M-SSA for prediction of SSTA may not be the best method in this case.

The use of a low order AR process may be more beneficial for such a short data set.

2.3.4 AR Process Modelling

The order of the AR process determined by the MEM was stable for each of the SST PC time series. The first four SST PCs required a higher order, or greater number of parameters, than the remaining SST PCs (Table 2.2) However, no PC required more than 6 parameters with most of the order of 2. The parameter estimates for each of the AR models are given in Table 2.2. The parameter estimates were insensitive to the number of lags used.

Forty weeks of independent data, from the week beginning 29/12/96 to the week beginning 04/10/97, were available for forecasting and verification.

When forecasting up to 3 weeks ahead the bias in both the forecasts using the AR process and using persistence averaged close to zero over the entire domain. However, the bias in the persistence forecasts were more consistently close to zero over the spatial domain than the time series forecasts. The time series forecasts also had generally higher rmse, around twice that of the persistence forecasts (see Figure 2.6 for an example of the errors associated with the forecasts).

2.4 Discussion and Conclusions

The statistical prediction of SSTs was examined using two methods, spectral analysis and time series analysis. For the weekly SST data the spectral analysis predictions were more accurate than those obtained using time series (AR) methods. However, when predicting SSTs on a weekly time scale, particularly for short lead times, the use of persistence as the forecast method leads to an improvement in forecast skill over the statistical methods.

For the monthly SST data, the statistical model (M-SSA), performs better than persistence when the lead time is over 3 months. When the lead time is less than this then persistence is a better forecast method.

The results from this study indicate that little is to be lost by using persistence in models requiring forecasted SST values, particularly when short lead times are involved.

3. Dynamical Seasonal Forecasts During the 1997/1998 El Nino

In this section, we present some results from an AGCM seasonal forecasting experiment which was conducted in the BMRC during the 1997/1998 El Nino event. The aim is to study to what extent the model can reproduce some the observed features of seasonal rainfall and surface temperature anomalies over the Australian region and globally by persisting earlier observed SST anomalies into the forecast period. We also investigate the capability of the model to capture the observed atmospheric circulation anomalies. Two different convection schemes are applied in the model to assess the sensitivity of the forecasts to the convection parameterization. NCEP reanalysis data are used to evaluate, or verify, the model forecasts.

3.1 Model and Experiments

The global AGCM used in this study is similar to the one used in Frederiksen et al. (1999) except that the model vertical resolution is increased. This is a BMRC climate model (version 3.7). It is a spectral model with rhomboidal R31 horizontal resolution, together with a 17 sigma level vertical representation. A single layer Manabe bucket land-surface model is applied to represent the soil moisture storage. In this study, we have concentrated on assessing the success of using persisted observed global SST anomalies in the dynamically extended seasonal forecasts, thus the impacts of soil moisture anomalies are excluded by initializing the model with observed soil moisture climatology. To investigate the sensitivity of the model forecasts to the convective scheme used in the model itself, two different convective parameterizations are employed which include a version of the Kuo convective scheme (Kuo, 1974) and a version of the Tiedtke mass-flux convective scheme (Tiedtke, 1989). Further details of the model dynamics and physics in the BMRC AGCM are given in McAvaney and Colman (1993) and references therein.

In this section, we summarize the model forecasts made from March 1997 to March 1998, covering the developing and mature phases of the 1997/98 El Nino event. These forecasts were conducted in a quasi-operational mode, with the results presented each month at the National Climate Centre’s monthly outlook meeting. Typically, a series of forecasts would commence in the second week of each month (designated month zero) and forecasts would be made for the next three succeeding months. Thus, following on from the conclusions of the section 2, we used a simple SSTA forecasting approach employing simply the latest observed weekly SSTAs, for the first week of month zero, which would be persisted into the forecast period. These SSTAs are then superimposed onto a twenty-week climatology extending into the forecast period. This twenty-week climatology is formed from the monthly climatology of Reynolds and Smith (1995) using the harmonic analysis technique of Epstein (1991) to derive daily SST values. The skill of using such SST forecasts will be discussed below in section 3.3.

In evaluating the model forecasts, we will concentrate only on the skill of the model in predicting climate anomalies, that is, with respect to either a model climatology or a control run using only an SST climatology. Two sets of runs were conducted for the model with the Kuo convection scheme (hereafter referred to as the K-forecasts); one set uses the climatological weekly SSTs and the other uses the aforementioned method. Thus, the first set of runs constitutes a set of control runs against which the second set can be compared, and differences, or anomalies, can be determined. Only one set of runs was performed for the model with the mass-flux scheme (hereafter called the M-forecasts) which only uses the forecast SST. Anomalies for this latter case are determined with respect to a ten year Atmospheric Model Intercomparison Project (AMIP) run (using the same model configuration and the observed global SSTs from 1979 to 1988).

An ensemble approach was employed to overcome the model’s sensitivity to initial conditions (see also section 5, and the discussion on potential predictability). In each set of experiments, an ensemble of six runs was generated with different initial conditions from the 6-hourly BMRC analysis system near the time when the weekly SSTAs were obtained. We were restricted to six runs per monthly forecast because of the quasi-operational nature of the project, with forecasts needing to be prepared for the monthly NCC outlook meeting. For the K-forecasts, the monthly averaged anomalies were calculated as the differences between ensemble averages of the six runs using forecast SSTs against the ensemble averages with climatological SSTs. As mentioned previously, the M-forecast anomalies were derived as differences between the ensemble averages and the model climatology from the AMIP run.

These anomalies were compared with anomalies derived for the forecast period using observed data for the period and observed climatologies. For rainfall and surface temperature, two observational datasets were used; one from the observational rainfall and temperature data collected for the Australian region by the NCC, and the other from the daily averaged National Center for Environment Prediction (NCEP) reanalysis dataset. As both showed very similar features in most of the months examined, and because the NCEP dataset is global, we shall concentrate here on the comparison of the model forecasts with the NCEP reanalysis. For the comparison of model and observed large-scale circulation anomalies, we also will use the NCEP datasets. Monthly and seasonal anomalies of the NCEP analysis, during the forecast period, were calculated with respect to a 40-year climatology derived from NCEP daily values for the period 1958-1997.

Before describing the results of this study, we consider, in the next section, appropriate measures of forecast skill.

3.2 Skill Scores – LEPS

The most commonly used measures of skill, when assessing climate model simulations or forecasts, are the rmse and the anomaly correlation. Here, in addition to these, we used as our principal skill measure the linear error in probability space (LEPS) score developed by Ward and Folland (1991) and revised by Potts et al. (1996). This score is related to the difference between the position of the forecast and the observation in the cumulative probability distribution of the particular climate variable under consideration.

Thus, for individual forecasts, if the position of the forecast in the cumulative distribution is  (ranging from 0.0 – 1.0) and that of the observation is , then the LEPS skill score is defined to be

From Eq.(3.1), it follows that, for a given observation , the largest skill score occurs when . In addition, most skill is attached to correctly forecasting at the extremes of the cumulative distribution (giving a score near 2). Correctly forecasting near the middle of the cumulative distribution (), gives the smallest correct skill score (approximately 0.5). The worst possible score (-1.0) is assigned to forecasting at the exact opposite extreme for an extreme observation (e.g.  and ).

It is desirable to have a measure of overall skill, given a set, or ensemble, of forecasts. To achieve a skill score range from 100% to –100%, average skill (SK) may be defined (Potts et al., 1996) as

Here the summation is over all pairs of forecasts and observations,  is the individual score for each forecast, and  depends on whether the numerator is positive or negative. For a positive numerator,  is the sum of the maximum possible scores given the observations (obtained by setting  in Eq.(3.1)). If the numerator is negative,  is the sum of the moduli of the worst possible scores given the observations, obtained by setting setting either  or  in Eq.(3.1) and taking the negative value with the largest modulus.

Here, we shall be mainly concerned with testing the skill of the forecast anomalies. This has the advantage of removing some of the biases seen in the model simulation of climate variables. We shall also use the cumulative probability distribution derived from the 40 year NCEP dataset 1958-1997 for each climate anomaly we consider. The question of the significance of the model skill score also arises. To determine a 95% significant level for the overall skill score SK we have used a Monte Carlo approach in which we have generated 1000 random forecasts covering the entire period March 1997 – March 1998. These forecasts have been generated by a random choice from the 40-year NCEP dataset of anomalies. The corresponding skill scores have been ranked with the 50^th top score providing the 95% significant level.

3.3 SSTA Forecasts

As mentioned in section 2, the use of persistent SSTAs, whether monthly or weekly, tends to give improvements in forecast skill compared to the M-SSA and AR methods, for short lead times extending out to three months. It is of interest therefore to calculate the LEPS scores for an SSTA forecasting scheme with persistence out to about three months ahead. In Figure 3.1 we show the overall skill score SK for forecasts 4, 8 and 12 weeks ahead. Here, the score is based on forecasts using each week of the period March 1997 – March 1998 to forecast for week 4, 8 and 12 ahead. The SSTA cumulative probability, used in the skill score, is based on the Reynolds 40 year monthly dataset extending from 1950-1989, for which anomalies are given only for the region 45S-65N.

From Figure 3.1, we see that there is high skill over much of the world’s oceans during the forecast period. In particular, there is very high skill in the eastern equatorial Pacific, the eastern north Pacific, the western Indian Ocean and parts of the Atlantic. As the forecast period is extended, skill decreases in most areas with regions of negative skill appearing, especially at higher latitudes and in the western Pacific. However, there is still appreciable skill in the eastern Pacific and along the Pacific-rim, the western Indian Ocean and around Australia. The statistical significance of these results is shown in Figure 3.2 where we show only those areas that are significance at 95%. The forecasts still show significant skill in those areas thought to be important regions of SST-forcing during the period.

3.4 Climate Forecasts

In this study, we have restricted our evaluation of the model forecasts to an in-depth analysis of four climate variables: rainfall, surface temperature, 200hPa geopotential height and mean sea level pressure (MSLP). We also concentrate on three-month seasonal mean anomalies, although individual month forecasts were also evaluated. For rainfall and surface temperature we highlight the Australian region as well as discuss the global forecasts. The 200hPa geopotential height and MSLP are used to evaluate forecasts of the large-scale circulation anomalies.

3.4.1 Rainfall

Overall, model forecasts of rainfall show very limited skill over the Australian region during this period. Figure 3.3 and 3.4 show the observed and K-forecast rainfall anomalies (positive shaded, negative unshaded) for the first twelve of the thirteen seasonal forecasts. While there are regions where both show anomalies of the same sign, the differences are quite obvious, especially during the latter half of the forecast period over the north and east of the continent, when the model gives substantial anomalies of opposite sign. The skill is quantified in Figure 3.5, where we have applied the LEPS score of Eq.(3.1) to each forecast (positive skill shaded). While there are areas of skill, for each forecast, it is clear that there is no consistent systematic pattern of skill throughout the thirteen forecasts. This is also true for monthly forecasts at 1-3 month lead times.

However, globally there are regions where our forecasting strategy does show some consistency in skill. This can be clearly seen in Figure 3.6 which shows the global LEPS skill score (only areas of positive skill are plotted). Throughout the forecast period, there is appreciable skill in the tropical eastern Pacific, tropical western Indian Ocean and over Indonesia. This diagram also showing increasing skill as the El Nino develops with strengthening SSTAs in the eastern Pacific. By May 1997, there were clear indications in the surface and sub-surface thermal structure of the tropical eastern Pacific that an El Nino was imminent. Prior to the May forecast, skill scores are globally relatively moderate.

The overall skill score, for all the forecasts combined, is shown in Figure 3.7, with only the skill at the 95% level plotted. As in Figure 3.6, significant skill is confined largely to the tropics, with maxima in skill in the east tropical Pacific, over Indonesia and southeast Asia, the western Indian Ocean and central east African coast, and northeast Brazil. Over Australia, there are some isolated areas of significant skill. Skill over the middle and high latitudes does not generally exhibit a coherent pattern and probably reflects the more stochastic nature of the rainfall process in the extratropics.

It should be noted that we are dealing here with only a small number of forecasts and consequently it is to be expected that randomly generated forecasts may prove to be fairly skillful, and will therefore raise the 95% significance level. The same procedure applied to a larger number of forecasts would be expected to yield more skill.

In Figure 3.8, we show three other measures of the model skill, being the root mean square error (rmse), the mean anomaly error (me, positive shaded), or anomaly bias, and the anomaly correlation (ancor, positive shaded) for all thirteen anomaly forecasts. As expected, the largest rmse occurs in those regions where the rainfall anomalies are largest (that is, the tropics). Even for the anomalies, there are some quite large mean errors. In particular, the model over-estimates the rainfall anomalies throughout much of the equatorial Pacific and Indonesian area. Over Australia, there is generally a negative bias and some slight positive bias over the southern parts. The anomaly correlation between the model forecast and the observed is, by and large, consistent with the pattern of overall skill shown in Figure 3.7. That is, regions of significant positive skill correspond with regions of positive correlation.

Individual LEPS scores for each M-forecast over Australia is shown in Figure 3.9. While there is some consistency between these and those shown in Figure 3.5 for the K-forecasts, it is clear that there are also some marked differences in some regions. Globally (Figure 3.10), the M-forecast skill scores are qualitatively similar to Figure 3.7, but with a more coherent structure in the region of the South Pacific Convergence Zone (SPCZ), where the model response appears to be improved with the mass-flux convective scheme. Over Australia, there is some significant skill over the northwest and southeast corners of the continent.

The corresponding rmse, me and ancor are shown in Figure 3.11. It is clear that model errors are generally larger when the mass-flux convective scheme is used. Otherwise, there are many qualitative similarities with the K-forecast results.

3.4.2 Surface Air Temperature

In the case of surface temperatures, both the K-forecasts and M-forecasts show some forecasting skill over Australia, especially over the eastern part. Figure 3.12 and 3.13 show the observed and K-forecast surface temperature anomalies (positive shaded). In this case there appears to be much closer qualitative and coherent similarities, with the exception of the southern parts of the continent. This is quantified in Figure 3.14 which shows the individual forecast LEPS scores. Throughout much of the forecast period there is consistent skill over much of the eastern coast.

Globally (Figure 3.15), the major skill is over the tropical oceans, particularly over the eastern Pacific Ocean, the Indian Ocean and the Atlantic Ocean. As with rainfall, the skill shows an evolving pattern, becoming progressively larger and more coherent as the El Nino strengthens. In the extratropics, there is skill over the northeastern Pacific Ocean, the southeastern Pacific Ocean and the southern Indian Ocean. There is also appreciable skill over tropical/subtropical landmasses, notably central/southern Africa, South America and Southeast Asia. The K-forecasts also demonstrate skill over the North American continent, with a pattern reminiscent of a PNA-like pattern. These same features can also be seen in the overall skill score (Figure 3.16).

Between 60S-60N, the root mean square temperature errors and biases are relatively small (Figure 3.17), with the largest errors occurring at the higher latitudes. As with the rainfall, the pattern of anomaly correlation is basically consistent with the overall LEPS score, with areas of positive and negative correlation corresponding, respectively, to areas of significant skill or lack of skill.

Figure 3.18 shows the individual LEPS score for the M-forecasts over Australia. Especially in the latter half of the forecast period, these are qualitatively similar to the corresponding K-forecasts (Figure 3.14). This is also true of the overall global LEPS scores (compare Figure 3.16 and Figure 3.19). The most noticeable differences can be seen in the rmse and me error plots (Figures 3.17 and 3.20). The M-forecasts have much larger errors and biases over the continents, and especially over mountainous regions. Over Australia, the rmse is about twice that for the K-forecasts. The global anomaly correlation pattern is, however, similar in both cases.

3.4.3 200hPa Geopotential Height

Results for the model forecasts of the 200hPa geopotential height suggests that both sets of forecasts have very good skill in the low latitudes (roughly between 30S-30N) and moderate skill extending into the extratropics. This can be seen, to some extent, in Figures 3.21 and 3.22 which show the seasonal anomalies for the NCEP analysis and the K-forecasts. It is clear from this sequence of diagrams that both observed and model show similar features, especially in the latter half of the forecasts period when the El Nino is well established. Thus, for example, the evolution of features such as the bi-polar positive anomaly in the tropical eastern Pacific, the Pacific North American (PNA) anomaly pattern and Southern Hemisphere wave structures, especially south of South America, are seen in both figures. In the first half of the forecast period, there is also evidence, in both figures, of the characteristic Pacific South American (PSA) anomaly pattern with high centred south of South America.

As with rainfall and surface air temperature, the skill of the forecasts evolves throughout the forecast period, as the Pacific and Indian Ocean SSTAs establish themselves and grow in magnitude. Thus, in Figure 3.23, we see that it is not until JJA that the K-forecasts show coherent and quite substantial skill throughout the subtropical band 30S-30N. In the latter half of the period, the model shows very high skill in the tropics/subtropics, over the North Pacific and North America (PNA pattern) and in the southern Pacific Ocean south of South America.

The overall skill score (Figure 3.24) suggests that the model has significant skill over much of the globe with the most noticeable exception of the high latitudes of the Southern Hemisphere. Figure 3.25 shows that this is also the region where the K-forecasts have the largest root mean square errors and biases. Consistent with Figure 3.24, the positive anomaly correlations also occur where the model has maxima in skill.

The individual skill scores for the M-forecasts are shown in Figure 3.26. It shows essentially the same features as seen for the K-forecasts in Figure 3.23. However, there appears to less skill in the extratropics (including the Northern Hemisphere) and over Australia, in this case. This is reflected also in the overall skill score (Figure 3.27), and in the larger root mean square errors and biases seen in the M-forecasts (Figure 3.28) over these regions.

3.4.4 Mean Sea Level Pressure

With significant positive SSTAs over the eastern Pacific Ocean during the El Nino event, the model successfully simulated the shift in the tropical Walker circulation. This feature is reflected in the increase of MSLP over the western Pacific Ocean and the decrease of MSLP over the eastern Pacific Ocean (see Figures 3.29 and 3.30, for the observed, and K-forecasts of, MSLP anomalies). Associated with the eastward shifting of the ascending branch of the Walker circulation, convection activities also move towards the central Pacific Ocean. In both model results, there is a significant increase in convective cloud formation over the eastern Pacific Ocean, and a significant decrease over the western Pacific Ocean (not shown).

Figure 3.31 shows the individual skill scores for the K-forecasts. Again, it is not until the latter half of the forecast period that large regions of coherent and substantial skill appear on these maps. The overall skill score appears in Figure 3.32. These results indicate that the Southern Oscillation during the course of the 1997/98 El Nino event is well forecast by the model. Good forecasting skill is seen over large areas of the eastern and western Pacific Ocean and over the SPCZ. There is also some skill over the tropical Atlantic Ocean, parts of the North and South Pacific Ocean, and parts of the Indian Ocean. Except for the northern half of Australia, Indonesia, parts of central Europe, the Middle East, Alaska and eastern Russia, there is little significant skill over land and the extratropical oceans, where the root mean square errors and biases (Figure 3.33) tend to be largest.

The individual and overall skill scores for the M-forecasts are given in Figures 3.34 and 3.35, respectively. Compared with Figures 3.31 and 3.32, for the K-forecasts, there are some noticeable differences. Firstly, there appears to be additional skill, in this case, over much of the Indian Ocean, central and southern Africa, the Middle East, India and Southeast Asia, even in the early forecasts. Compared with the K-forecasts, the overall skill in the eastern Pacific Ocean is confined to the far east. Also, as was true for each of the other three climate variables, the root mean square errors and biases, are generally much larger than those for the K-forecasts, in regions of low skill.

Both the K-forecasts and M-forecasts show a failure to forecast MSLP changes in the Southern Hemisphere storm track region. Zheng and Frederiksen (1999) have shown that the model does not simulate the storm tracks well, especially during austral summer and this may lead to enhanced weather noise in the model. This could be responsible for the very limited skill of the model rainfall forecasts over the Australian region and the temperature forecasts in the southern part of the continent.

3.5 Discussion and Conclusions

Using persistent SSTA forecasts, the skill of model forecasts of rainfall, surface temperature and the large-scale atmospheric circulations have been assessed against NCEP reanalysis data. Two convection schemes have been employed to investigate the sensitivity of the model forecasts to the convective parameterization. Overall, the model results show skill over large areas of the tropics. The model does not show systematic patterns of rainfall forecasting skill in the extratropics. Skillful rainfall forecasts are mainly over regions where significant SSTAs are located. Surface temperature forecasts are mainly skillful over the tropical oceans. In the extratropics, skill is seen over the northern and southern Pacific Ocean, and the southern Indian Ocean. Over land there is skill over parts of Australia, central and southern Africa, India, Southeast Asia, South America and over part of the North American continent.

Both model forecasts capture the Southern Oscillation and changes in the tropical convection associated with an eastward shift of the ascending branch of the Walker circulation. This is reflected in high MSLP forecast skill in the west and east Pacific. The 200hPa geopotential height is perhaps the most skillfully forecast, with significant skill over much of the globe, excepting mainly the high latitudes of the Southern Hemisphere.

Model forecasts indeed show sensitivities to the convection schemes used, although there are some broad similarities in the patterns of overall skill. Differences are most obvious when the focus is on a specific region. Generally, the M-forecasts of all four climate variables, considered here, have larger root mean square errors and biases than for the corresponding K-forecasts in those areas of low skill. This is primarily associated with the different mechanisms of convection parameterized in each of the schemes.

4. AGCM Simulation of Climate 1958-1991

In this section, we examine the ability of the BMRC AGCM to simulate inter-annual variations of rainfall, surface temperature, 200hPa geopotential height and the MSLP when the model is forced by observed SSTs over many decades. We also investigate the extent to which the model is able to reproduce the observed SOI relationships.

4.1 Model and Experiments

The model used here is that of Frederiksen et al. (1999), and is described briefly in section 3.1 above. It has only nine vertical sigma levels and only the Kuo convective scheme is implemented in this set of experiments. An ensemble of five runs were conducted with different initial conditions and were forced by observed SSTs (taken from the Hadley Centre’s Global Ice and Sea Surface Temperature dataset version 1.1) for the period December 1949 to November 1991. The NCEP daily reanalysis data for the period 1958-1991 is used to verify the model simulations and calculate skill scores for the five run ensemble average of each climate field. As in section 3, we concentrate here only on climate variable anomalies.

4.2 Simulation Skill Scores

Here we present plots of the overall skill score SK (Eq.(3.2)), significant at the 95% level, for each calendar month and calculated on the basis of 34 years of simulation (1958-1991). Figure 4.1 shows the simulation skill for rainfall anomalies over this period. As in Figure 3.7, significant and coherent skill is confined largely to the tropical oceans and has a maximum in the eastern Pacific Ocean during most of the year. In austral winter and spring, there are secondary maxima in skill in the region of the Indonesian archepelago and eastern South Africa. In the extratropics there appears to be no systematic pattern of skill. Over Australia, there is some significant skill over northern Queensland and parts of the central east of the continent during austral winter and spring. In autumn and winter, there are some months with significant skill in a band stretching from the northwest across the continent.

In Figure 4.2, we present the 95% skill score based on all months taken together. Again, significant skill is confined to the subtropical band 30S-30N. However, there is significant skill over much of Australia, South Africa, along the eastern coast of Africa, the Middle East, Southeast Asia and northeast Brazil.

The simulation of surface air temperature anomalies is skillful over most of the global oceans (Figure 4.3). This is due mainly to the response of the near surface temperatures to the precribed SSTs. However, as in Figure 3.16, maxima in skill occur in the eastern Pacific Ocean, the North Pacific, the Indian Ocean and the Atlantic Ocean. The location of the skill maxima also seems to have some dependence on the annual cycle of solar radiation. This is especially true for the Atlantic Ocean where maximum in skill occur in the summer hemisphere, and in the subtropics during spring and autumn.

Except for the tropics, there is no significant skill over land. As in Figure 3.16, skill is confined to the central African continent, Southeast Asia, Central America and northern South America. Over Australia, there is some skill on the eastern coast during some months of the year.

For the 200hPa geopotential height field, significant skill is very much confined to the 30S-30N band during much of the year (Figure 4.4). This contrasts with Figure 3.24, where, during the 1997/98 El Nino, there was significant skill in the extratropics. This can be partly explained by the fact that the simulations involve both ENSO and non-ENSO years, and, as we saw in Figure 3.23, we might expect to have the best skill in peak ENSO years.

Figure 4.5 show the corresponding plots for the simulation skill of MSLP anomalies. Again, as in Figure 3.32, significant skill occurs predominantly in the subtropical band, with maxima generally in the eastern and western Pacific Ocean, and Atlantic Ocean. Thus the model clearly simulates the SOI during 34 years considered here. In fact, the correlation between the SOI from each of the 5 runs conducted and the observed, is better than 0.8. Given this fact, it is of interest to see to what extent observed SOI relationships of climate variables are reproduced in these simulations.

4.3 SOI Relationships

A number of studies (e.g. Ropelewski and Halpert, 1987; Philander, 1990) have identified regions around the globe where climate variability is strongly influenced by ENSO. Here, we investigate the observed and modelled relationships between rainfall, surface air temperature, 200hPa geopotential height and MSLP with the corresponding observed and modelled SOI. We have calculated monthly averaged anomalies, with the annual cycle removed, of each of the simulated and observed (NCEP) climate variables, for the period 1958-1991. To these, and the corresponding monthly SOI time series, we have applied a five month running filter, and have correlated the observed and simulated time series with, respectively, the observed and simulated SOI time series. Thus, we are essentially correlating seasonal mean time series with the SOI.

Figure 4.6 shows global plots of these correlations for the rainfall. The similarity between the observed and modelled correlation patterns is quite remarkable, and shows similar features to those seen in Ropelwski and Halpert (1987). The magnitude of the correlations are also comparable between the observed and modelled. With the exception of western and central Africa, parts of Europe and the very high latitudes of the Northern Hemisphere, the SOI relationship is correctly simulated over most continents. Thus, for example, a positive SOI (La Nina), is associated with enhanced rainfall over Australia, Southeast Asia, northeast Brazil, Canada, southern and eastern Africa, and India, and reduced rainfall over the United States and southern South America. However, over many of the continental areas, the linear SOI relationship explains a relatively small fraction of the rainfall variation. Consequently, as we show in the next section, internal, or chaotic, variations can play a large role in many of these locations. This may explain, for example, why over Australia, even though the SOI relationship is strong, the model forecasts during 1997-98 were not very skillful.

In Figure 4.7, we show the corresponding diagram for the surface air temperature, and again, with the exception of the high latitudes of the Northern Hemisphere, the similarities between the observed and modelled is remarkable. Both show a correlation pattern, over the oceans, which has features of the typical ENSO SSTA pattern (La Nina phase), with negative correlations in the Indian Ocean, Atlantic Ocean, throughout most of the equatorial Pacific Ocean and the eastern Pacific rim. Positive correlations are seen in the western, northern and southern Pacific Ocean. For a negative SOI, as occurred during 1997/98, positive surface temperature anomalies would be expected over, for example, the eastern and western coasts of Australia, parts of central and southern Africa, subtropical South America and Central America, and Canada. From Figure 3.16, we see that these are the areas where the model forecasts showed significant skill during this latest event. Interestingly, however, over eastern North America, the cooler than normal surface temperatures were not skillfully forecast.

The relationship between variations in the general circulation and the SOI is also very well reproduced in the model simulations. Figure 4.8 show the SOI/200hPa geopotential height relationships. They are virtually indistinguishable and show the PNA and PSA teleconnections, and the bi-polar structure of the subtropical Pacific associated with ENSO. For negative SOI, these features would have opposite sign. All of these features were seen in the 1997/98 model forecasts with significant skill (Figure 3.24). Similarly, changes in the Walker Circulation associated with the SOI (Figure 4.9), were also largely captured by the model forecasts, although restricted mainly to the tropics/subtropics where the correlations are largest.

4.4 Discussion and Conclusions

The results of this section show that, in simulations with observed SSTs, the model has some significant skill in simulating climate anomalies, especially in the 30S-30N latitude-band, not unlike the model forecasts for 1997-98. The model is also able to reproduce the observed SOI relationships over the same 1958-1991 period. These relationships go some way to explaining the skill seen in the model forecasts of section 3, especially in those regions with very high SOI correlation. However, it is clear that even when these correlations are high, the model may not exhibit significant forecast skill (as happened for example with rainfall over Australia). This may partly be due to the role played by internal, or chaotic, processes in climate variability. This is an issue that we investigate in the next section.

5. Potential Predictability and Model Validation

The results of sections 3 and 4 beg the question what limits there are to predictability of climate variables. In particular, how much of the interannual variation in seasonal means can be attributtable to the external SST forcing. In an endeavour to determine this we have applied the technique of Rowell et al.(1995) to monthly values of each variable from the AGCM simulations of section 4. However, because this technique can only be applied to ensembles of model runs, we have developed in collaboration with Dr Xiaogu Zheng (National Institute of Water and Atmospheric Research) a new technique, using daily values. This technique can also be applied to observed data (Zheng and Frederiksen, 1999). With this approach, it is possible to estimate the correlation between the forced (or predictable) components of the simulated and observed seasonal means, thus providing a method of validating the model results.

5.1 Results from Monthly Values (1950-1991)

The technique of Rowell et al. (1995) uses monthly data from model ensemble runs to determine SST-forced and internal, or chaotic, components of variation of climate variables. Thus, seasonal anomalies of a given climate variable in simulation s and year y are represented as

Here, s = 1,…, S is the simulation number, y = 1,…, Y is the year;  is the overall mean,  is the deviation from induced by external forcings (e.g. SST) where this deviation is identified with year y;  represents the interannual variability in simulation s induced by internal or chaotic sources, where this variability is assumed to be purely random and independent from one simulation to another. Following Zwiers (1996), an average over a given subscript will be represented by an ``open circle’’. Thus, represents the average over all years and the average over all years and all simulations.

If we let  denote the estimated variance of a variable, then, for example, represents the variance of . Hence, the following estimates are possible (Rowell et al., 1995),

Thus, /  and /  are estimates of the proportions of the total seasonal variance which is attributable to chaotic and forced components. For the simulations of section 4, this latter ratio represents the fraction due to SST-forcing, i.e. the potentially predictable part of the total variance.

Eqs.(5.2)-(5.4), are equally valid if y represents the number of months in each simulation and monthly anomalies with the annual cycle removed. Figure 5.1 shows, for the five simulations of section 4, the fraction of the total variance of the monthly mean of global precipitation, surface air temperature, 200hPa geopotential height and MSLP, which is attributable to SST-forcing over the period January 1950-November 1991.

From Figure 5.1, it can be seen that there is appreciable SST-forcing in precipitation throughout the tropics with local maxima located in the eastern Pacific, to the north of Australia, throughout the Indonesian archipelago, in the south Atlantic, and throughout the Indian Ocean. Over land, better than 20% of the variation in rainfall over eastern and western central Africa, northwest Australia, Indonesia, Northeast Brazil and central America is SST-forced. These results are very similar to the patterns of forecast and simulation skill shown in Figures 3.7, 3.10, 4.1 and 4.2. That is, the model appears to have skill in those areas where SST-forcing plays a large, if not dominant role. For rainfall, this tends to be predominantly over the oceans. Over Australia, less than 20% of the variation is due to SST-forcing or, alternatively, more than 80% is due to internal variability.

In general, the same conclusions hold for the other three climate variables shown in Figure 5.1. However, for surface air temperature, the SST-forcing, and associated model skill (Figures 3.16, 3.19, and 4.3), extends into the mid-latitude oceans of both hemispheres. There are also substantial land areas throughout the sub-tropics (central Africa, northwest Australia, Southeast Asia, Central America and Brazil) with appreciable SST-forced variation.

The largest model skill for the 200hPa height (Figures 3.24, 3.27 and 4.4) occurs predominantly throughout the latitude band 30S-30N, and this is reflected also in Figure 5.1. For MSLP, maxima in SST-forcing occur in the tropical eastern Pacific, Atlantic Ocean, Indonesian archipelago, and Indian Ocean. Again, these are regions where the model shows skill (Figures 3.32, 3.35 and 4.5), although, in this case the convective scheme used in the model seems to affect the results (compare Figure 3.32 and 3.35).

5.2 Results from Daily Values (1958-1991)

The techniques discussed in the previous section can only be applied to ensembles of model simulations. In this section, we propose a new technique (Zheng and Frederiksen, 1999) which can be applied to both model and observed data for which daily values are available. Thus, it is possible to compare results from both model and observations. To quantify this comparison, we also define estimates of the correlations between model and observed components of variability.

5.2.1 Potential Predictability

In deriving this new technique, it is useful to consider a seasonal mean, derived from daily meteorological values, as comprising three parts: the forced component, the internal source component and the weather noise component. The forced component is the response to slowly varying external boundary forcing (e.g., SST, ice coverage and vegetation coverage) and radiative forcing (e.g., greenhouse gas concentration and sulphate aerosol loading). The internal source component is the atmospheric variability induced by internal dynamics, land surface processes and soil moisture. The weather noise component is the response to rapidly varying day-to-day weather events. While this may be an artificial distinction, it allows useful progress to be made on the analysis of the long-term variability.

The forced component is potentially predictable at the long-range, because the forcings themselves are potentially predictable. The internal source component may be potentially predictable up to three months in advance (Madden 1983; Shukla 1983), but appears to be unpredictable beyond that. Historically there are two definitions of potential predictability: the potential predictability defined by Madden (1976), i.e. the variability of the forced and internal source components, and the potential predictability which is defined by Rowell (1998), i.e. the variability of the forced component only. The weather noise components are unpredictable beyond a couple of weeks, because weather events are unpredictable beyond about 10 days (Lorenz 1973). Therefore, it is of critical importance, in the context of AGCM-based long-range forecasting, to skillfully simulate the variability of the forced component.

A traditional problem of model validation is whether the interannual variability of the components of seasonal means is being correctly simulated. Some variability, such as the weather noise component and the sum of the forced and internal source components, can be estimated for both observational and simulated data. For a skillful AGCM, the simulated variability should correspond closely with the observed (Zwiers 1996). Some variability, such as, that due to the forced or internal source component, can be estimated from simulated data but not from observations.

We represent a daily value from an ensemble of AGCM simulations forced by observed SST as

Here, s = 1,…, S is the simulation number, y = 1,…, Y is the year and t = 1,…, T is the day in a specified season with length T. The terms in the model have the following interpretation:  is the overall mean,  is the deviation from induced by forcings where this deviation is identified with year y;  represents the low frequency interannual variability in simulation s induced by internal sources, where this variability is assumed to be purely random and independent from one simulation to another.  represents the weather noise at time t of the specified period in year y of simulation s. Theare assumed to comprise a stationary normal stochastic process with mean zero and are independent and identically distributed with respect to s and y (Zwiers 1996). A daily observed value is represented as

Here, the definitions of , ,  and  are similar to the definitions of , ,  and  respectively (Madden 1981; Zwiers 1987; Zheng 1996). Furthermore, , ,  and  are assumed to be statistically independent of each other.

Following Zwiers (1996), a "circle" represents a subscript when an average is taken over that index. For example,  indicates the seasonal mean of  in year y and simulation s and  is the estimated overall mean computed by averaging  across all days and simulations, etc. With this notation, a simulated seasonal mean of daily values can be expressed as

where ,  and  are referred to as the simulated forced component, internal source component and weather noise component (of the seasonal mean), respectively. Similarly, an observed seasonal mean of daily values can be expressed as

where ,  and  are referred as the observed forced component, internal source component and weather noise component (of the seasonal mean) respectively. Again, the symbol  denotes the estimated variance of a variable. Formulae for the estimated variability of simulated data are listed in Table 5.1.

Name	Notation	Formula
Variability of forced component
Variability of seasonal means
Variability of weather noise component
Variability of internal source component
Variability of forced and internal source components

The top five formulae are adopted from Rowell et al (1995). The second formula () is the key to deriving the following three, and it only requires seasonal means, not daily data. However, the procedure used by Rowell et al. (1995) to derive

is based on the assumption that the simulated forced components are statistically independent of each other and of the other two components and. No such assumption is made here. Since we have assumed that the forced components are fixed factors, the Kolmogorov strong low of large numbers can be used to derive (5.9) (Zheng and Frederiksen, 1999). The sixth formula () is a key to deriving the last four. The last three formulae (,  and ) are related to both key formulae. A precondition for the sixth formula is that the true spectral density function of weather events {} is smooth at the zero frequency (Zwiers 1987). An approach to check this precondition is given by Zheng and Renwick (1998).

Formulae for the estimated variability of observed data are listed in Table 5.2.

Name	Notation	Formula
Variability of seasonal means
Variability of weather noise component
Variability of forced and internal source components
Variability of forced component		if is assumed.

The second formula () is similar to the key formula for  in Table 5.1. However, no formula in Table 5.2 is similar to another key formula for  in Table 5.1, because observations can only have one sample. The last formula () relies on the assumption that the ratio of variability of forced component over variability of internal source component is correctly simulated. If it is not, the last formula is a biased estimate. However, if the ratio is incorrectly simulated, the model is already wrong in this aspect.

We have applied this method to the 200hPa geopotential height field for JJA and DJF, using daily values from the AGCM simulations of section 4 and the NCEP reanalysis for the period 1958-1991. In both seasons, the model tends generally to underestimate the total inter-annual variance, with maxima typically reduced by a factor of a half. Figures 5.2 and 5.3 show, respectively, for the AGCM, the components of the JJA inter-annual variance and the proportion of the total variance due to each. The corresponding results for the NCEP data is shown in Figures 5.4 and 5.5.

The spatial patterns of variance are very similar and this is reflected in the weather noise and the combined low frequency and SST components, which can be derived for both model and observed data. Interestingly, the weather noise maximum in the Southern Ocean has similar magnitude in both cases. In the mid- to high-latitudes of the Northern Hemisphere, the model underestimates the weather noise variance. In the sub-tropics and tropics, both model and observed have the smallest weather noise variance.

The deficit in model variance in the Southern Hemisphere appears to arise from an underestimation of the combined low frequency and SST component, and this varies with geographical location. Thus, for example, while the model variance, in this component, in the Southern Ocean/South Pacific, east of the dateline, is about a half the observed, it is severely underestimated in the region to the southwest of Australia and in the south Atlantic. In the Southern Hemisphere, the SST-forced component tends to be larger than the low frequency component for both the model and observed (if the assumptions of Table 5.2 are assumed). In the Northern Hemisphere, where the variance during JJA is generally smaller, both components are, by and large, of comparable magnitude. Throughout the tropics, the low frequency component is very small for both model and observed.

Figures 5.3 and 5.5 show the proportions of the total variance attributable to each component. As in Figure 5.1, the chaotic components dominate the middle and high latitudes, with the weather noise component generally the larger of the two. The SST-forced component dominates the 30S-30N latitude-band and this is the region where the potential for predictability is highest.

During DJF (Figures 5.6 ,5.7,5.8 and 5.9), the deficit in model variance arises from the combined low frequency and SST-forced variance being generally considerably underestimated in middle and high latitudes, except importantly in the north Pacific. Consequently, the weather noise component in the model is relatively much more dominant at these latitudes (see Figures 5.7 and 5.9). The low frequency component appears also to be larger in the tropics during DJF than JJA. As for JJA, the SST-forced component is largest in the 30S-30N latitude-band.

5.2.2 Correlation and Model Validation

In this study, an estimation procedure for the correlation between the forced components for simulated and observed seasonal means is proposed. It represents the correlation between the simulated and real-world response to the observed SST variations. If the correlation coefficient is high, it suggests that the SST-forced components for simulated seasonal means are realistic. If the correlation is low, it may indicate some deficiency in the AGCM’s climatology.

No forced component can be calculated exactly. In practice, the seasonal mean is used as an unbiased estimate, with the seasonal mean of the internal source and weather noise components being considered as the random error term of the estimate. The procedures proposed in this paper enable us to estimate the temporal relationship between the forced components for the simulated and observed seasonal means without knowing the numerical values of these forced components.

Historically, for an ensemble of runs, the correlation between simulated ensemble seasonal means and observed seasonal means is often used as an index of AGCM skill in simulating the SST-forced variability. In fact, this represents the real-world forecast skill for an ensemble of finite size (or simulation number) given a perfect SST forecast. A low correlation could be due to the strongly chaotic nature of the climate system, so it is not the optimal index of AGCM skill in this context. The correlation between the forced component of simulated seasonal means and the observed seasonal means represents the real-world forecast skill for infinite simulations, because the ensemble mean of simulated seasonal means converges towards the simulated forced component as the simulation number increases. The difference between the two correlation coefficients discussed in this paragraph is a measure of the possible improvement of the real-world forecast skill by adding more simulations.

Here, the symbol  will also be used to denote the estimated covariance between two variables. For example,  represents the estimated covariance between  and . Furthermore, the symbol  will be used to denote the estimated correlation coefficient between two variables.

The formulae for the correlation between the simulated and observed forced components () and the correlation between the simulated forced components and the observed seasonal means () are listed in Table 5.3.
 
 

Name	Notation	Formula
Correlation between simulated ensemble seasonal means and observed seasonal mean
Correlation between simulated forced components and observed seasonal means
Correlation between simulated forced components and forced and observed internal source components
Correlation between simulated and observed forced components

depends on the assumption imposed on the last formula in Table 5.2 which may not be true. Therefore the correlation between the simulated forced components and the observed forced and internal source components
() is introduced as its estimated lower bound. If the variability of the observed internal source component
() is significant, the lower bound estimate is always less than the true estimate. Otherwise, the two estimates are equal. The variability of the observed internal source component cannot be estimated as observations afford only a single sample, and therefore its significance cannot be tested. Zwiers (1996) proposed a procedure for testing the significance of the simulated variability of the internal source () which can be used as a reference.

The correlation between the simulated ensemble seasonal means and the observed seasonal means () is also listed in Table 5.3. It is an ordinary Pearson correlation coefficient, while the last three correlation coefficients in Table 5.3 are inflated Pearson correlation coefficients, because the variances of related variables are not estimated as the sample variances. So it is possible that an estimated correlation is greater than one. In this case, the correlation is truncated to one. If the variability of the simulated forced component , or the variability of observed forced and internal source components , is too small, then the three estimated inflated correlation coefficients could be severely distorted. Analysis of variance can be used to test whether these variances are erroneously inferred as significant (Zwiers 1996). These tests may be too conservative for the current purpose, because the variability can be erroneously inferred as insignificant by only using these tests. In practice, we use the following rule: when the ratio of variability of simulated forced component over the variability of simulated seasonal means (/), or the ratio of variability of observed forced and internal source components over variability of observed seasonal means (/) is less that 0.05, the three inflated correlation coefficients are not further calculated. If the variability of the forced component is significant and well simulated, the correlation between the simulated ensemble seasonal means and the observed seasonal means should be positive. Otherwise, there is no reason for further estimating the three inflated correlation coefficients. The key point in proving the formulae in Table 3 is that the boundary forcing is the only possible source of the covariability between the simulated and observed seasonal mean temperature, that is,

       ».                            (5.10)

Its detailed proof is given in Zheng and Frederiksen (1999).

All the estimates in Tables 5.1-5.3 are strongly consistent, that is, they converge to the true values almost everywhere as the number of data points increases to infinity. Most formulae in Tables 5.1-5.2 only require large Y´ T , but the formulae in Table 5.3 require large Y. Because data are often limited in practice, it would introduce error. Another error source is that the statistical model may not exactly fit real meteorological data. Therefore, all the estimates should be treated with caution. Particularly we suggest that the inflated correlations in Table 5.3 be treated qualitatively and not quantitatively. Where possible checks of their consistency with other climatological information should be made.

Figure 5.10 shows for JJA plots of the correlation of the model total field and observed total field and the correlation of the model SST-forced component with, respectively, the observed total field, SST-forced + low frequency and SST-forced components. The corresponding plots for DJF appear in Figure 5.11. As might be expected, the correlations increase as we progressively remove the chaotic components of variation. This is especially true in the mid-latitudes. The largest correlation between the total fields occurs in the 30S-30N latitude-band in both seasons. Comparing the plot for the correlation between the total fields and that for the correlation between the model SST-forced component and the observed total field (equivalent to taking an infinite ensemble), we see that, for JJA, little extra skill would be gained by increasing the ensemble size. This is not the case for DJF where further improvement is possible.

Figures 5.10 and 5.11 show that the model has considerable skill in simulating the SST-forced component of the interannual variation of 200hPa geopotential height. This skill extends from the tropics into the extratropics. In the extratropics, especially, it is the chaotic components that prevent a more skillful simulation of the total field. In the tropics, during JJA, there is little difference between the plot of the correlation of the total fields and that for the SST-forced components. This is due to the fact that the chaotic components in the tropics are relatively small, and differences arise because of the weak weather noise component. During DJF, the situation is quite different because there is also a non-negligible low frequency component of chaos.

6. Conclusions

The main aim of this study has been to assess the use of a stand-alone AGCM, forced by predicted sea surface temperatures, in forecasting seasonal climate anomalies one to three months ahead with particular reference to Australian rainfall. Because statistical models for SST prediction are likely to outperform dynamical model based SST schemes in this situation, we have investigated two possible statistical schemes: one based on a multichannel singular spectrum analysis and the other based on a low order autoregressive technique. Neither of these techniques was able to perform better than simple persistence when using either weekly or monthly time scales with lead times less than 3 months. Thus, we have employed a simple SSTA forecast scheme where the latest observed SSTA persist into the forecast period. This approach was shown to have significant skill (based on LEPS scores) throughout the forecast period during the 1997/98 El Nino. This was especially true in the eastern equatorial Pacific Ocean, the eastern north Pacific Ocean, the western Indian Ocean and parts of the Atlantic.

A series of quasi-operational forecasts were conducted using the BMRC AGCM throughout the 1997/98 El Nino period, covering mainly the development and mature phases of this event. Two different convective schemes ( Kuo, 1974; Tiedtke, 1989) were used to investigate the sensitivity of the model forecasts to convective parameterization. In general, the model forecasts showed significant skill in forecasting seasonal anomalies of rainfall, surface air temperature, 200hPa geopotential height and MSLP over large areas of the tropics and subtropics, but predominantly over oceans. There was little systematic skill in the rainfall and MSLP forecasts in the extratropics. In contrast, skill in surface temperature and 200hPa geopotential height forecasts extended into the extratropics of both hemispheres.

Model forecasts were able to capture the Southern Oscillation and changes in tropical convection associated with the observed eastward shift of the ascending branch of the Walker circulation. Important teleconnections in the global circulation, such as the PNA and PSA anomaly patterns were also skillfully forecast. Interestingly, model skill for all climate variables evolved throughout the forecast period, showing increasingly more widespread and larger skill with the maturing of the El Nino. Thus, most skill occurred in the latter half of the forecast period.

The model forecasts showed sensitivities to the convective schemes used. However, except for MSLP, there was not much different in the overall patterns of skill as determined by combined LEPS scores. Differences were most obvious when specific locations were targeted. There was, however, a tendency for the forecast using the mass flux scheme to consistently have larger root mean square errors and biases. In particular, the response to the SST-forcing appears much more concentrated and localized over areas of largest SSTA.

Over Australia, there was little significant and systematic skill in the rainfall forecasts. This may be due to excessive internal variations (or sensitivity to initial conditions) in our forecasts. Zheng and Frederiksen (1999) show that this model does not simulate the storm tracks well, tending to be too baroclinically unstable, in the Southern Hemisphere, especially in summer. This may be a source of additional weather noise in the model. Alternatively, as we discuss below, this may be due to the relatively low SST-forced component of rainfall variation over Australia. The forecasts did, however, show some significant skill in the surface temperature, 200hPa geopotential height and MSLP forecasts over Australia. Here, and generally globally, significant skill tended to coincide with areas that show high correlation with the SOI.

An analysis of the model’s ability to simulate climate anomalies, over the period 1958-1991, when forced with observed SSTs, showed patterns of skill consistent with those seen in the model forecasts discussed above. The model was also able to reproduce the observed SOI relationships during this period. These relationships were shown to explain, to a large extent, the skill seen in the forecasts, especially in those regions highly correlated with the SOI.

Finally, in this study, we have looked at the roles played by SST-forcing and internal dynamics, or chaos, in forcing climate variations, as a means to understanding what might be the limits of predictability. That proportion of the climate variation that is SST-forced is potentially predictable because the SST-forcing is itself potentially predictable. Using the technique of Rowell et al. (1995) with monthly data for the period 1950-1991, we have shown that SST-forcing predominates mainly in the 30S-30N latitude-band for all four climatic variables considered here. That is, climate variations are potentially predictable mainly in this region. Surface temperature also showed high SST-forcing in the extratropics. Again, as with the SOI relationships, these results go some way to explaining the skill of model forecasts during 1997/98.

For this project, we have also developed a new technique for determining potential predictability and for validating model simulations. This technique, using daily data, can be applied to both model and observations. With this method, model and observed climate variations are split into components consisting firstly of a weather noise component, due to weather events, and a combined low frequency and SST-forced component. For the model, it is then possible to further split the latter component into its two constituent parts. Using assumptions about the model, these two components can also be estimated for the observed.

We have applied this method to the 200hPa geopotential height field for JJA and DJF taken from AGCM simulations and NCEP reanalyses for the period 1958-1991. During JJA, we found that the weather noise and low frequency internal components dominated the extratropics and gave little contribution in the tropics, where the SST-forced component prevailed. This was also generally the situation during DJF, with the exception that the low frequency component gave a significant contribution in the tropics as well. In both seasons, the model tended to underestimate the total inter-annual variance. We found that the deficit in the model variance, in the Southern Hemisphere during JJA, appears to arise from an underestimation of the combined low frequency and SST-forced component. During DJF, this component is generally considerably underestimated in middle and high latitudes, with the exception of the north Pacific. As a result, the weather noise is proportionately much higher than in the observed, at these latitudes. The analysis suggests that some improvement is possible in the model simulation if these deficits can be corrected.

An estimation procedure for the correlation between the forced components of the simulated and observed seasonal means has also been proposed. This is a fairer measure of model skill than that currently used, that is, the correlation between the total fields. Applying this to the JJA and DJF model and observed data we found that, as expected, the correlation between model and observed increase as we progressively remove the chaotic components of variation. This was particularly true for the mid-latitudes. Largest correlations tended to occur in the 30S-30N latitude-band where the model forecasts and simulations showed highest skill scores. This technique also allows us to see the impact of increasing the ensemble size of simulations. We found that for JJA, little skill would be gained by increasing the ensemble size from five members. This was not the case during DJF where further improvement was possible.

Our results show that the model has considerable skill in simulating the SST-forced component of the interannual variation of the 200hPa geopotential height. This skill extended from the tropics to the extratropics, where it was the chaotic components that prevented a more skillful simulation of the total field. During JJA, in the tropics, we found that there was little difference between the plot of the correlation of the total fields compared with the SST-forced components, due principally to the small contribution of the chaotic components. In DJF, this is not the case, and correlations between the SST-forced components were considerably higher because of the non-negligible contribution of the low frequency component to the total field.

Continued efforts are being taken to study the model forecasts in more detail to see if improvements can be made in future. One area we hope to investigate is the impact of soil moisture anomalies on the seasonal forecasts. Here, the initialization of soil moisture was with climatology. This may not be sufficient, and some longer spin-up strategy may be necessary to initialize the soil moisture to more realistic values. We also hope to investigate moisture fluxes during the 1997/98 El Nino to try to understand how these varied from earlier episodes, and what influence this may have had on Australian rainfall during this period. The analysis here has also pointed out model deficiencies that need to be investigated in more detail. Another avenue of future research is to investigate global patterns of variation of the weather noise, low frequency and SST-forced components of the total climate fields. Techniques for doing this are currently being developed within the Climate Group.
 
 

REFERENCES

Barnston, A.G., M.H. Glantz, and Y. He, 1999: Predictive skill of statistical and dynamical climate models in SST forecasts during the 1997-98 El Nino episde anf the 1998 La Nina onset, Bull. AMS, 80, 217-243.

Dix, M. R. and H. G. Hunt, 1995: Chaotic influences and the problem of deterministic seasonal predictions. Int. J. Climatol., 15, 729-752.

Drosdowsky, W., 1993a : An analysis of Australian seasonal rainfall anomalies : 1950-

1987. I : Spatial patterns. Int. J. Climatol., 13, 1-30.

Drosdowsky, W., 1993b : An analysis of Australian seasonal rainfall anomalies : 1950-

1987. II : Temporal variability and teleconnection patterns. Int. J. Climatol., 13, 111-

149.

Drosdowsky, W. 1993c Potential Predictability of Winter Rainfall Over Southern and

Eastern Australia using Indian Ocean Sea-surface Temperature Anomalies. Australian

Meteorological Magazine 42:1-6.

Drosdowsky, W., L. Chambers. 1998. Near Global Sea Surface Temperature

Anomalies as Precictors of Australian Seasonal Rainfall. BMRC Research Report No 65,

Bureau of Meteorology Australia.

Epstein, E.S.,1991: On obtaining daily climatological values for monthly mean, J. Climate, 4, 365-368

Frederiksen, C.S. and R.C. Balgovind,1994: The Influence of the Indian Ocean/Indonesian SST Gradient on the Australian winter rainfall and circulation in an atmospheric GCM. Q. J. R. Meteorol. Soc., 120, 923-952.

Frederiksen, C.S. and J.S. Frederiksen, 1996 : A theoretical model of Australian Northwest Cloudband disturbances and Southern Hemisphere storm tracks : The role of SST anomalies. J. Atmos. Sci., 53, 1410-1432.

Frederiksen, C. S., D. P. Rowell , R.C. Balgovind and C. K. Folland, 1999: Multidecadal simulations of Australia rainfall variability: The role of SSTs. J. Climate. (In press).

Harzallah, A. and Sadourny, 1995: Internal versus SST-forced atmospheric variability as simulated by an atmospheric general circulation model. J. Climate, 8, 474-495.

Haykin, S. 1979. Nonlinear Methods of Spectral Analysis. Topics in Applied Physics

34. Springer-Verlag, Berlin.

Jones, D., 1999: The prediction of Australian land surface temperatures using near

global sea surface temperature patterns. BMRC Report. (In press)

Keppenne, C.L., M. Ghil. 1992. Adaptive Filtering and Prediction of the Southern

Oscillation Index. Journal of Geophysical Research 97(D18):20449-20454.

Kumar, A., M. P. Hoerling, M. Ji, A. Leetmaa and P. Sardeshmukh, 1996: Assessing a GCM’s suitability for making seasonal predictions. J. Climate, 9, 115-129.

Kuo, H.L., 1974: Further studies of the parameterization of the influence of cumulus convection on large-scale flow, J. Atmos. Sci., 31, 1232-1240.

Lorenz, E. N., 1973: The standard error of time-average estimates of climatic means. J. Appl. Meteor., 12, 1066-1069.

Madden, R. A., 1976: Estimates of the natural variability of time averaged sea level pressure. Mon. Wea. Rev., 104, 942-952.

Madden, R. A., 1981: A quantitative approach to long-range prediction. J. Geophys. Res., 86, 9817-9825.

Madden, R. A., 1983: Comments on "Natural variability and predictability", reply. Mon. Wea. Rev.,111, 586-589.

McAvaney, B.J. and R.A. Colman, 1993: The AMIP experiment: The BMRC AGCM

Configuration, BMRC Report, 38, 43pp.

McBride, J.L., and N. Nicholls, 1983 : Seasonal relationships between Australian

rainfall and the Southern Oscillation. Mon. Wea.Rev., 111, 1998-2004.

Nicholls, N., 1984 : Seasonal relationships between Australian rainfall and north Australian sea-surface temperature. Extended abstracts, Conference on Australian rainfall variability, Australian Academy of Science, Australian Bureau of Meteorology, 71-73.

Nicholls, N., 1988 : More on early ENSOs : Evidence from Australian documentary

sources. Bull. A M S, 69, 4-6.

Nicholls, N. 1989. Sea Surface Temperatures and Australian Winter Rainfall. Journal

of Climate 2:965-973.

Philander. S.G., 1990: El Nino, La Nina, and the Southern Oscillation, Academic Press,

New York, 293pp.

Pittock, A.B., 1975 : Climatic change and the patterns of variation in Australian rainfall. Search, 6, 498-504.

Pittock, A.B., 1984 : On the reality, stability, and usefulness of southern hemisphere teleconnections. Aust. Met. Mag., 32, 75-82.

Plaut, G., and R. Vautard, 1994: Spells of Low-Frequency Oscillations and Weather Regimes in the Northern Hemisphere. J. Atmos. Sci., 51, 210-236.

Potts, J.M., C.K. Folland, I.I. Jolliffe and D. Sexton, 1996: Revised "LEPS" scores for assessing climate model simulations and long-range forecasts, J. Climate, 9, 34-53.

Reynolds, R.W. and T.M. Smith, 1995: A high resolution global sea surface temperature climatology, J. Climate, 8, 1571-1583.

Ropelewski