Quantifying Chaos in the Atmosphere: the use of ensembles of General Circulation Model Integrations

Richard Washington
School of Geography, University of Oxford
Email: richard.washington@nta.geog.ox.ac.uk

Abstract

The atmosphere is known to be forced by a variety of energy sources, including radiation and heat fluxes emanating from the boundary layer associated with sea surface temperature anomalies. The atmosphere is also subject to internal variability which is essentially unforced and is thought to be a basic characteristics of fluids. Whereas much work has been done in quantifying the links between external forcing of the atmosphere and its long term response as well as the influence of boundary layer forcing in determining organised large scale modes of planetary scale circulation, less is known about the importance of internal variability or chaos in determining the evolution of weather and climate. In part this derives from the inability of scientists to run experiments on the weather. General Circulation Models (GCMs) now provide for this possibility. Multiple evolutions of the climate system may be computed in GCM simulations. Where these simulations are similar except for the conditions by which the model is initialised, the degree of departure in the evolution of climate from one model run to the next corresponds precisely to the degree of internal variability or chaos present in the model atmosphere. A methodology for quantifying this chaotic forcing is presented and is applied to century long intergrations of the U K Meteorological Office model HADAM2A..

Keywords: Chaos, internal variability, forced variability, GCMs, initial conditions, atmosphere

1. Introduction

As if in dual encouragement and warning to aspiring students, a colleague sports the following extract from Einstein's unpublished letters on her office door:

"As a young man, my fondest dream was to become a geographer. However, while working in the customs office I thought deeply about the matter and concluded it was far too difficult a subject. With some reluctance, I then turned to physics as a substitute."

Einstein's well known reluctance to embrace the notion of a cassino-like God in the guise of Quantum Theory is undoubtedly more widely appreciated than his caution of Physical Geography. But one seems forced to wonder whether he unconsciously coupled the difficulty he foresaw in Geography to his dislike of the non-deterministic route he painfully witnessed Physics taking, with the similar ancestry of determinism set against chance as the competing champions controlling process in the physical world. With or without Einstein, there can be no doubt that Physical Geography has done less to acquaint itself with the master-croupiers' personality than Physics has. This is not withstanding the playground which surrounds Physical Geography, be it along coastlines (Mandlebrot, 1967), in soil development (Phillips, 1993a), hillslopes (Phillips, 1993b) or in the atmosphere (Lorenz, 1963).

The atmosphere is known to be a chaotic system (Lorenz, 1963, 1993). We therefore face some stark choices when studying its ultimate character. Either we write down a system of equations, whose solutions behave chaotically and seem to emulate the climate system or we tackle the problem indirectly. This is essentially a choice of waiting for the mathematics to evolve, which could take a very long time, or we use what tools we currently have. This paper considers some simple approaches to representing the main attractors of the atmosphere from observed data. It moves on to suggest that the only indirect route we have at present to achieve the goal of quantifying chaos may be made through the use of ensembles of model integrations made over multidecadal time scales. The solution is not unique but probabilistic. An example of the forced and free response of the UK Meteorological Office model HADAM2A is given over the high seasons.

2. Phase Space Modes: Attractor Analogies?

At present we do not know the nature of the equations governing attractors in the atmosphere. We have no notion even of the number of equations there might be. What is possible, though, is the empirical analysis of atmospheric fields for which there is reasonable spatial and temporal data coverage. Some of the techniques used to understand the main modes present in such fields come close to defining the phase space patterns (or attractors) around which the atmosphere appears to vary. It should come as no surprise that Lorenz (1956) was the first to introduce one of the most frequently used statistical approaches to the study of climate, namely Empirical Orthogonal Functions (EOFs). This could be thought of as the separation between signal sub-space, controlled by dynamics, and the noise sub-space which is high-dimensional, though not unimportant to the individual realisations of any climate evolution. Such separation is founded on non-linear degeneration of variance from the leading modes downward.

A generalised partitioning of full phase space may be formalised as follows:

with t usually taken as time. The K modes and time coefficients describe the dynamics in the signal sub-space while the last term is the noise sub-space.

The EOF solution of an n by n correlation matrix, R, of an m by n data matrix, X, may be given as:

where l is the n by n identity matrix.

The eigenvectors of R may be obtained from solutions to the series of equations

where v_i are the eigenvectors.

An attractive aspect of EOFs is the geometrical orthogonality of the vectors and therefore their statistical independence. In essence they provide an indication of the leading patterns of independent variability. Figure 1 shows the correlation form of the leading eigenmode of satellite derived rainfall for January over the period 1979-1996. The EOF has been calculated on a grid 2.5 latitude by 2.5 longitude using the Xie and Arkin (1996) data. The eigenvector weights clearly point to an ENSO like structure with enhanced (suppressed) rainfall over Indonesia correlated with suppressed (enhanced) rainfall in the central Pacific. The pattern is coherent across most of the global tropics. From this perspective, ENSO can be seen as the main attractor of global tropical rainfall.

Figure 1: Leading mode of January satellite derived rainfall.

Association between an oceanic index of ENSO, here taken as the sea surface temperatures (SSTs) in the Nino 3 region in the central Pacific (5^oN-5^oS, 150^oW-90^oW), and global observed sea level pressure is shown as correlation coefficients for January, February and March for 1950-1993 in Figure 2. The sea level pressure data is the UKMO Global Mean Sea Level Pressure version 2.1F (GMSLP21F). This correlation map may be thought of as the linear association between the largest known mode of ocean surface temperature variability and the near surface atmosphere. While the modulus of the correlation coefficients is generally large, particularly in the Pacific, the correlation is nowhere perfect. Other modes in addition to internal (or unforced) variability account for a still larger portion of the variance in this field. This simple approach is able to indicate the variance of sea level pressure explained by Nino 3 SST (given as the square of the correlation coefficient), but it can say little else about the forced or free atmosphere. It would, however, be possible to partition the forced and free variance more exactly if we had several experiences of the twentieth century climate so that the extent of divergence in each of the trajectories could be quantified, given the same forcing from the Nino 3 region. Such requirements are obviously fanciful in the real world so that this approach to considering forced and free responses is severely limited.

Figure 2: Correlation between Nino3 SST and SLP JFM 1950-1993

2. Quantifying Forced and Free Manifolds

The atmosphere is known to be sensitive to initial conditions so that precisely the same forcing can produce hugely divergent evolutions of the climate system given infinitesimally small differences in initial conditions (Lorenz, 1963). The only opportunity we have to test this sensitivity is in a model whereby only one variable, here the initial conditions, may be arbitrarily changed while a forcing function remains constant. Multiple integrations of such an experiment, with each integration subjected to identical forcing yet each being uniquely defined by arbitrarily small differences in starting conditions of the model atmosphere, would yield a set a data for which a quantitative scheme could be used to partition the data variance into components essentially common to all integrations and variance unique to each integration. We might label these the Forced and Free Manifolds respectively.

Such ideas are a necessary extension of Lorenz's (1963) ideas, but have taken time to appear in the literature (Rowell, 1998; Ward and Navarra, 1997). The process of identifying the forced and free model response consists of separating the variance within and between the model integration for each grid box in the model for any particular field. The free response, corresponding to unforced variability (internal variability) may be measured by the variance of each observation from the ensemble mean for each year of the model integration. Rowell et al. (1995) specify this as follows:

where N is the number of years in the integration, n is the number of ensemble members and is the ensemble mean for the i^th year.

The ensemble mean variance is defined by:

where is the ensemble mean of all data,other terms being defined earlier.

Variance of the forced variance is then calculated as the difference between the ensemble mean variance and the internal variance. The total variance is the sum of the internal variance and the forced variance. Ratios between the internal and total variance and between the forced and total variance, when calculated for each of the model grid boxes for a particular field, yield the Free and the Forced Manifold respectively.

3. Model Formulation

This work uses the UKMO Hadley Centre atmosphere general circulation model (HADAM2A). HADAM2A has been successfully used in studies of anthropogenic climate impacts (Folland et al., 1988). Features include a resolution of 2.5^o latitude by 3.75^o longitude, with 19 hybrid layers in the vertical (5 representing boundary layer processes). Physical parameterisations include a gravity wave drag scheme and a radiation scheme that computes fluxes in several parts of the shortwave and longwave bands and which responds to cloud growth occuring dynamically in the model. Large scale cloud amount, convective cloud amount, liquid water concentrations and ice concentration are the prognostics used in the cloud scheme. A surface hydrology scheme includes parameterised surface moisture storage, surface and sub-surface runoff and a prescribed vegetation canopy that interacts with precipitation and runoff. A detailed description of the model physics is available in Phillips (1994).

Several critical assumptions relating to the model performance are necessary for the method adopted here to be successful in quantifying forced and free modes. These include:

1. That a 'mean' or quasi-equilibrium state exists both in nature and in the model.
2. That the model converges to the correct solution. Whether this solution is regarded as correct cannot be determined analytically.

HADAM2A is an improved version of the model submitted by the UKMO as part of the AMIP (Atmospheric Model Intercomparison Project) (Gates, 1992) which identified the UKMO model as amongst the best four models currently available. More recent reviews and model intercomparisons have again shown the performance of HADAM model to be comparable with the best available models (Srinivasan et al., 1995; Airey et al., 1996; Gates et al., 1996; Joubert, 1997).The climatologies of HADAM1 (the version submitted to AMIP) and HADAM2 are compared in Gallani (1994).

Rowell (1996) has evaluated the HADAM1 model climatology in some detail. Comments here are restricted to precipitation and sea level pressure as these are observed variables which are available for the length of the model integration. Spatial correlations between observed and model 1961-1990 climatologies (mean and interannual standard deviations) area weighted over 90^oN to 80^oS at 5^o by 5^o resolution for sea level pressure using the Basnett-Parker SLP and over 70^oN to 60^oS at 5^o by 7.5^o resolution for precipitation using Legates-Willmott (1990) and Hulme (1994) are shown in Table 1. This corresponds with the best observed and best modelled periods (given the better quality of the SST forcing towards the end of the twentieth century). On the whole correlation coefficients are encouragingly high and very significant at greater than 99% level for one tailed tests. Improvements to HADAM1 in the form of corrections to the detrainment, downdraft and evaporation schemes together with improvements made to the horizontal diffusion parameters have seen a substantial improvement in HADAM2 (Gallani, 1994). The correlations shown below are therefore an underestimate of the HADAM2's performance.

The structure of southern hemisphere model fields from simulations by the Atmospheric Model Intercomparison Project, including HADAM1, have been scrutinised by Joubert (1997). This study shows that HADAM1 captures the zonally asymmetric flow at 500 hPa better than any of the major AMIP climate models in both January and July.

Table 1: Spatial correlations between observed and model (HADAM1) 1961-1990 climatologies (mean and interannual standard deviations) area weighted over 90^oN to 80^oS at 5^o by 5^o resolution for sea level pressure using the Basnett-Parker (B-P) SLP and over 70^oN to 60^oS at 5^o by 7.5^o resolution for precipitation using Legates-Willmott (L-W) and Hulme (From Rowell, 1996).

Field	Statistic	Verification Data	DJF	MAM	JJA	SON
SLP	Mean	B-P	0.96	0.97	0.97	0.97
SLP	Std Dev	B-P	0.90	0.82	0.77	0.74
Precip	Mean	L-W	0.65	0.65	0.63	0.66
Precip	Mean	Hulme	0.78	0.67	0.70	0.77
Precip	Std Dev	Hulme	0.81	0.75	0.78	0.77

Improvements made in HADAM2A have not removed all important differences between the observed and modelled climate. Remaining problems in modelled sea level pressure (Gallani, 1994) include the following:

1. The Aleutian Low is too deep from September to May and extends too far east over Alaska.
2. The low pressure trough to the east of the Philippines is too persistent in September to October.
3. Subtropical highs in the southern hemisphere are too intense in all seasons and the southern hemisphere high latitude low pressure belt is too intense in all seasons.

HADAM1 deficiencies were improved greatly by altering model diffusion, so that excessive rainfall in the south Pacific and southern Indian Ocean and deficient rainfall over South America during DJF were corrected in HADAM2. Likewise HADAM1 in JJA was too wet in west Pacific and monsoon regions of SE Asia. This has been corrected in HADAM2. There is still too much rain in East India between March and May with the monsoons starting too early in the year. Importantly the model Hadley Circulation is stronger by 10% than most recent NWP analyses.

3.1 Model Experiments

The suite of experiments analysed in this paper consist of an ensemble of four integrations forced throughout with observed SST, sea ice, changing equivalent carbon dioxide and the changing tropospheric sulphate aerosol data set parameterised as in Mitchell et al. (1995). Each integration runs from October 1903 to December 1994 making these integrations the longest of their kind yet undertaken. Ensemble members differ in terms of arbitrarily chosen initial conditions from recent UKMO operational analyses. The first two months of the experiment were discarded to allow the model atmosphere to reach near equilibrium with the boundary layer forcing after imposition of the arbitrary initial atmospheric conditions. SST data were interpolated to the model grid and linearly interpolated to 5-day means to ensure a smooth evolution of the boundary layer.

3.2 Sea Surface Temperature (SST) Data

The SST data used to force HADAM2A derive from the GISST (Global Sea Ice and Sea Surface Temperature) data set. GISST is the result of developments covering some two decades of research at the United Kingdom Meteorological Office (UKMO). Since these data set generation schemes are thoroughly documented elsewhere (Bottomley et al., 1990; Colman, 1992; Parker et al., 1995; Folland and Parker, 1995; Rayner et al., 1996) only the important details are covered here.

The primary data source for the GISST data set is ship observations with the addition of buoy data in recent years (Bottomley et al., 1990; Colman, 1992; Parker and Folland, 1995). Gaps in the data set were filled with data from the Comprehensive Ocean Atmosphere Data Set (COADS) (Woodruff et al., 1987). Observations averaged to form 1 degree by 1 degree pentad (5 day) means were converted to anomalies by subtracting a background average on a 1 degree by 1 degree pentad resolution. The monthly pentad anomalies were winsorised, though only if at least four 1 degree by 1 degree anomalies were being averaged. These data were then averaged to create 5 degree by 5 degree anomalies. Corrections to the pre-1942 SST (Folland and Parker, 1995) anomalies were necessary to compensate for biases resulting from the extensive use of uninsulated canvas and semi-insulated wooden buckets dependant on month and location. Regional corrections of errors resulting from truncation, which induced falsely low temperatures, have also been made to Japanese North Pacific SST data in the 1930s.

SST data used to force AGCMs must be complete. Extensive work based on eigenvector reconstruction of early twentieth century SST data, necessary for complete coverage, has been undertaken by Rayner et al. (1996). An early version of GISST2 has been used. These data give a considerably better representation of ENSO related SST variations before 1960 than in the first version of GISST (Parker et al., 1995).

Estimates of the Forced and Free Manifold modes were undertaken using NAG routine on the School of Geography SUN Enterprise 450. Graphics were produced using GrADS (Grid Analysis Display System) software. Model integrations were performed on the UKMO Cray supercomputers.

4. Results

Figure 3 shows the Forced Manifold of December to February (DJF) model rainfall. The statistics have been calculated over the complete period of model integration (1903-1994). The most striking result to emerge is the extremely high component of the Forced Manifold in the tropics but the extremely low values found in the midlatitudes. In essence this indicates that tropical atmosphere, as represented by rainfall which usefully integrates many atmospheric processes over several vertical layers in the model atmosphere, responds to SST forcing in a similar manner regardless of the initial conditions. The component of internal variability or chaos is therefore very small. The opposite is true in the midlatitudes. Here the atmosphere with which the model was initialised is still important many decades later, despite the forcing of the boundary layer by SSTs. While the tropical atmosphere is predictable on the basis of SST forcing, the midlatitude atmosphere, on the whole, is not. This result is not new (Charney and Shukla, 1981; Palmer and Anderson, 1994). On the whole, predictability is lower over the tropical and subtropical continents than over the oceans at the same latitude.

Figure 3: Model SST forced DJF rainfall 1904-1993 all ensemble members.

Figure 4 shows a detail of the DJF Forced Manifold for the tropical Atlantic region. The shift from forced to free response across the Amazon is abrupt. In the region of the Atlantic tropical-temperate cloud band in particular, forcing by SST is not important. Regions thought to be influenced by the dynamics of this feature, such as southern Africa (Harrison, 1986; Tyson, 1986) may therefore experience a larger component of chaos than typical for that latitude.

Figure 4: Model SST forced DJF rainfall - Atlantic detail. 1904-1993 all ensemble members.

Figure 5 shows the equivalent result for June to August (JJA). The general pattern of large forcing in the tropics against internally variability in the midlatitudes is apparent. An important observation is the limited size of the Forced Manifold in the Asian Monsoon region. it is clear that the Asian Monsoon constitutes an important agent of internal variability in the model atmosphere.

Figure 5: Model SST forced JJA rainfall 1904-1993 all ensemble members.

4.1 Internal Variability on multi-decadal time scales

The operation of the Forced and Free Manifold may be estimated on multi-decadal time scales by removing periodicities on time scales of less than ten years within each of the 7008 model grid boxes for each of the four ensemble members. This was undertaken by means of an Integrated Random Walk Kalman filter developed by Young et al. (1991).

The filtering technique fits a smooth line (called the trend component) through a time series. At time point k, the smooth line value is the level of the trend, lk, and the slope of the smooth line is sk. It is assumed that the slope varies through time as a random walk process. Since it is the slope of the smooth line sk that is modelled, the method is referred to as an Integrated Random Walk (IRW). Ward (1994) has noted that modelling the slope means that the IRW model is less vulnerable to large swings resulting from outliers in the observations than many simpler methods.

In the equation below it can be shown that the values of lk and sk define the state of the model, or the 'state space' of the system. The general equation for the evolution of the state of the model parameters is:

where Xk is a vector containing the values of the model parameters at time point k. F is the transition matrix and G the input matrix, which determine how the model parameters (X) and white noise input (h) yield the new estimate of the state vector. The constituent matrices for the IRW model is given by:

The two parameters that define the state of the model are given by:

which show that the actual level of the trend component lk is estimated as the level at the time point k-1, plus the implied change in the level based on the slope estimated for the time point k-1 (Ward, 1994).

The Integrated Random Walk observation equation is given by:

yk = (1 0 ) +ek

An arbitrary observation, yk consists of a residual ek and the level of the trend lk. The ratio of the variance hs to e controls the proximity of the smooth line fit to the last observation and is called the noise variance ratio (NVR) (Young et al., 1991).

The 50% amplitude of the filter's response function (F50) is given by:

Figure 6 shows the Forced Manifold on multidecadal time scales for the JJA season. Comparison of Figure 5 and 6 reveals that considerably less internal variability is present in certain parts of the midlatitudes on multidecadal time scales.

Figure 6: Model SST forced LF JJA rainfall 1904-1993 all ensemble members.

While much remains to be done in terms of assessing this method's sensitivity to the length of model integration and the number of ensemble members needed, the method is a step forward in our ability to quantify the chaotic component of the atmosphere.

5. Conclusions

This paper has reviewed some simple ways of quantifying chaos in the climate system. The direct route to achieving this is by writing down the equations which govern the atmospheric attractors. Dealing with the problem indirectly through the observational record is inherently problematic. Only 'analogous' solutions are available. Resorting to modelling allows a plausible scheme to be introduced whereby the probabilistic solution rather than the unique solution is achievable. Multiple runs of a climate model allow the nature of the Forced and the Free Manifold to be explored.

Acknowledgments

The author is very grateful for useful discussions with Chris Folland and David Rowell at the Hadley Centre of the UKMO. The idea of long integrations of an AGCM forced with SST belongs to them, as does the method of separating forced and free variance.

References

Airey, M. J., M. Hulme and T. C. Johns, 1996: Evaluation of simulations of terrestrial precipitation in UK Met. Office/Hadley Centre climate change experiments. Geophysical Research Letters, 23(13), 1657-1660.

Bottomley, M., C. K. Folland, J. Hsiung, R. E. Newell and D. E. Parker, 1990: Global Ocean Surface Temperature Atlas. HMSO, London.

Charney, J.G. and Shukla, J. 1981: Predictability of monsoons, in Monsoon Dynamics, eds J. Lighthill and R.P. Pearce, CUP, Cambridge, 99-109.

Colman, A.W. 1992: Development of worldwide marine data eigenvectors since 1985, CRTN, 28, Hadley Centre.

Folland, C. K. and D.E. Parker, 1995: Correction of instrumental biases in historical sea surface temperature data. Quarterly Journal of the Royal Meteorological Society, 121, 319-367.

Folland, C.K., D.M.H. Sexton, D.J. Karoly, C.E. Johnson, D.P. Rowell and D.E. Parker, 1998: Influences of anthropogenic and oceanic forcing on recent climate change. Geophysical Research Letters, 25, 353-356.

Gallani, M.L. 1994: Comparison between 1st and 2nd versions of the Climate Model, Hadley Centre. Internal Note no 58.

Gates, W.L. 1992: AMIP: The Atmospheric Model Intercomparison Project. Bulletin of the American Meteorological Society, 73, 1962-1970.

Gates, W.L., A.Henderson-Sellers, G.J. Boer, C.K. Folland, A. Kitoh, B.J. McAvaney, F. Semazzi, N. Smith, A.J. Weaver and Q-C. Zeng, 1996: Climate models - evaluation. In: Climate change 1995. The science of climate change, (eds) Houghton, J.T., L.G. Meira Filho, B.A. Callander, N. Harris, A. Kattenberg and K.Maskell, Cambridge University Press, Cambridge.

Harrison, M. S. J. 1986: A synoptic climatology of South African rainfall variations. PhD Thesis, University of Witwatersrand, South Africa.

Joubert, A. 1997: Simulation by the Atmospheric Model Intercomparison Project of the atmospheric circulation over southern Africa. International Journal of Climatology, 17, 1129-1154.

Lorenz, E.N. 1956: Empirical Orthogonal Functions and Statistical Weather Prediction. Scientific Report No 1, Contract AF19 (604) 1566, Met Dept, MIT.

Lorenz, E.N. 1963: Deterministic non-periodic flow, Journal of Atmospheric Science, 20, 130-141.

Mandlebrot, B.B. 1967: How long is the coastline of Britain? Statistical self-similarity and fractal dimensions. Science 56, 636-38.

Mitchell, J. F. B., T.C. Johns, J.M. Gregory and S.F.B. Tett, 1995: Climate response to increasing levels of greenhouse gases and sulphate aerosols. Nature, 376(6540), 501-504.

Palmer, T.N and Anderson, D.L.T. 1994: The prospects for seasonal forecasting - a review paper, Quarterly Journal of the Royal Meteorological Society 120, 755-793.

Parker, D. E., C. K. Folland and M. Jackson, 1995: Marine surface temperature: observed variations and data requirements. Climatic Change, 31(2-4), 559-600.

Phillips J. D. 1993a: Chaotic evolution of some Coastal Plain soils. Physical Geography 14(6), pp 566-580.

Phillips J. D. 1993b: Instability and chaos in hillslope evolution. American Journal of Science 293(1), pp 25-48.

Phillips, T.J. 1994: A summary documentation of the AMIP models. PCMDI report 18, Univ California, Lawrence Livermore Lab.

Rayner, N. A., E. B. Horton, D. E. Parker, C. K. Folland and R. B. Hackett, 1996: Version 2.2 of the global sea-ice and sea surface temperature data set, 1903-1994. CRTN 74, Hadley Centre.

Rowell, D.P. Folland, C.K., Maskell, K. and Ward, M.N. 1995: Variability of summer rainfall over tropical north Africa (1906-92): Observations and modelling, Quarterly Journal of the Royal Meteorological Society, 121, 669-704.

Rowell, D. P. 1996: Using an ensemble of multi-decadal GCM simulations to assess potential seasonal predictability. CRTN, 69, Hadley Centre, UKMO.

Rowell, D.P. 1998: Assessing potential seasonal predictability with an ensemble of multidecadal GCM simulations, Journal of Climate, 11(2), 109-120.

Srinavasan, G., M. Hulme and C.G. Jones, 1995: An evaluation of the spatial and interannual variability of tropical precipitation as simulated by GCMs. Geophysical Research Letters, 22(16), 2139-2142.

Tyson, P. D. 1986: Climatic change and variability in southern Africa. Oxford University Press, Cape Town.

Xie, P. and P. A. Arkin, 1996: Analyses of global monthly precipitation using gauge observations, satellite estimates, and numerical model predictions. Journal of Climate, 9(4), 840-858.

Ward, M. N. 1994: Tropical North African rainfall and worldwide monthly to multi-decadal climate variations. PhD thesis, Reading University, 313pp.

Ward, M.N. and Navarra, A. 1997: Pattern analysis of SST-forced variability in ensemble GCM simulations: Examples over Europe and the tropical Pacific, Journal of Climate, 10(9), 2210-2220.

Woodruff, S. D., R. J. Sutz, R. L. Jenne and P. M. Steurer, 1987: A comprehensive ocean-atmosphere data set. Bulletin American Meteorological Society, 68, 1239-1250.

Young, P. C., C. N. Ng, K. Lane, D. E. Parker, 1991: Recursive forecasting, smoothing and seasonal adjustment of non-stationary environmental data. Journal of Forecasting, 10, 57-89.