Robust stochastic seasonal precipitation scenarios

INTERNATIONAL JOURNAL OF CLIMATOLOGY Int. J. Climatol. 26: 2077–2095 (2006) Published online 5 May 2006 in Wiley InterScience (www.interscience.wiley.com) DOI: 10.1002/joc.1351

ROBUST STOCHASTIC SEASONAL PRECIPITATION SCENARIOS IOANNIS KIOUTSIOUKIS,* SPYRIDON RAPSOMANIKIS and REA LOUPA Democritus University of Xanthi, Environmental Engineering Department, Laboratory of Atmospheric Pollution and Pollution Control Engineering, Xanthi, Greece Received 10 May 2005 Revised 11 January 2006 Accepted 10 March 2006

ABSTRACT In this paper, a stochastic statistical forecasting methodology is employed for long-term predictions of winter precipitation over Greece. Lagged climatic indices and North Atlantic (NA) sea-level pressure (SLP) field are explored as potential predictors of the teleconnection. Rather than employing traditional stationary models, two dynamic regression-modelling schemes are analysed and validated and their parameter variation is interpreted. Dynamic regression models, in contrast to static (constant parameter) regression models, have time variable parameters (TVPs) evaluated through recursive optimisation and are suitable for analysis of non-stationary phenomena like most atmospheric processes. The analysis of the spectrum with non-stationary models points out that the most influential seasonal components of the winter precipitation anomalies have periods of 14 and 3.5 years, explain 40% of its variance, possess significant amplitude change and correlate significantly with the North Atlantic Oscillation Index Anomaly (NAOIA) and Southern Oscillation Index Anomaly, indicating their climatic origin. Furthermore, the forecasting skill of the dynamic models (R 2 = 0.71), in addition to reproducing the peaks, was found superior even to the hindcasting skill of the stationary model (R 2 = 0.55). Copyright  2006 Royal Meteorological Society. KEY WORDS:

non-stationary time series analysis; North Atlantic Oscillation; Southern Oscillation; winter precipitation anomaly; Greece

1. INTRODUCTION Seasonal climate prediction seeks to forecast the likely state of the climate several months in advance. Its scientific basis lies in the lower boundary conditions of the atmosphere (e.g. sea-surface temperatures (SST), land surface characteristics) that once perturbed, alter the likelihood of the occurrence of weather regimes on seasonal and long-term timescales (Goddard et al., 2001). Statistical and dynamical models have been utilized to extract these potentially predictable low-frequency signals. The dominant mode of seasonal climate variability on a global scale is the El Nino – Southern Oscillation (ENSO) phenomenon that arises from the strong ocean–atmosphere interactions internal to the tropical Pacific and overlying atmosphere. ENSO (Trenberth, 1997) seems to have a global influence beyond influencing tropical climate. The shifts in the location of the organised rainfall in the tropics and the associated latent heat release alters the heating patterns of the atmosphere that forces large-scale waves in the atmosphere that in turn establish teleconnections. In practice, most remote climate impacts (or teleconnections) arise through the propagation of SST disturbances that are excited in the tropical Pacific. In addition to teleconnections directly linked to SST changes, some arise from natural preferred modes of the atmosphere associated with the mean climate state and the land–sea distribution. The most prominent are the Pacific-North American (PNA) and the North Atlantic Oscillation (NAO). The NAO (Hurrell, 2001; Wanner et al., 2001; for a review) is one of the most dominant patterns of wintertime atmospheric circulation variability. Although NAO is evidently a mode of variability internal * Correspondence to: Ioannis Kioutsioukis, Democritus University of Xanthi, Environmental Engineering Department, Laboratory of Atmospheric Pollution and Pollution Control Engineering, Kimmeria Campus, GR-67100, Xanthi, Greece; e-mail: [email protected]

Copyright  2006 Royal Meteorological Society

2078

I. KIOUTSIOUKIS, S. RAPSOMANIKIS AND R. LOUPA

to the atmosphere, there is an increasing evidence that the observed changes in the NAO may well be a response of the system to observed changes in SSTs and there are some indications that the warming of tropical oceans is a key part of this. Relations between tropical SSTs and extra-tropical climate that have been recently documented include links between tropical Pacific SST and the NAO (Hoerling et al., 2001) and the significant influence of El Nino on the climate of the North Atlantic (NA) region (Mathieu et al., 2004). Cassou and Terray (2001) presented evidence that the winter European atmospheric variability is connected to anomalous oceanic conditions in the North Atlantic and the Pacific. In addition, Sutton and Allen (1997) presented evolving NA SST anomaly patterns and suggested that their propagation can yield downstream predictability several years ahead. Consequently, since SST evolution is slow relative to characteristic atmospheric timescales (Hansen and Bezdek, 1996), NA presents considerable potential for predictability of European climate more than a year in advance (Rodwell et al., 1999; Kushnir, 1999). The Mediterranean is a small-scale coupled atmosphere ocean system with a short response time and with high estimated human-induced climate change impacts (Hulme et al., 1999; Palutikof et al., 1996). The understanding of the relevant processes and physical mechanisms and their link to large-scale climate will improve the management of climate-related risks in the region. Over the last years, several papers have addressed the relationship mechanisms between climatic variability in the Mediterranean and fluctuations in the global circulation. ENSO episodes and the NAO are among the key factors that influence the periodicity of droughts in the Mediterranean region (e.g. Haylock and Goodness, 2004; Eshel and Farrell, 2000; Hurrell, 1995; Fernandez et al., 2003). Both time series (ENSO and NAO) do not show a relatively simple cyclical behaviour; additionally, the well-established association between the NAO circulation mode and the surface climate of Europe has recently been shown to be non-stationary (Rodo et al., 1997) and asymmetric (Cassou et al., 2004). The period from October to March (hereafter winter) is the epoch with the highest precipitation amounts in the Mediterranean as a whole. Despite the large spatio-temporal variability of precipitation, a significant fraction of its variation can be explained by large-scale circulation changes at different heights as explored in recent studies. Eshel et al. (2000) studied the eastern Mediterranean (EM) winter rainfall variability in terms of subsidence anomalies associated with NA sea level pressure (SLP) anomalies. The linear statistical model trained from this teleconnection closely matched the observed precipitation anomalies (PRECA) in hindcast mode but its forecast skill was moderate. Dunkeloh and Jacobeit (2003) analysed the Mediterranean rainfall variability for winter as well as for the other seasons in terms of NA – European area geopotential heights (HGT) at 500 and 1000 hPa. The first set of canonical correlation patterns were explaining 30% of the precipitation variance and were depicting the Mediterranean Oscillation (MO) (Conte et al., 1989; Palutikof et al., 1996), which in addition correlates significantly with the Northern hemisphere modes of the NAO. Xoplaki et al. (2004) studied winter precipitation variability in the Mediterranean using surface (SLP, SST) as well as lower and upper troposphere variables (HGT: 850, 700, 500 and 300 hPa) in a canonical correlation analysis (CCA). They found that the SST does not improve the overall performance while a linear combination of the other five predictors explains 30% of the Mediterranean winter precipitation variability. The normalised time coefficient of the first canonical correlation pattern of precipitation correlates with the negative NAO and is responsible for the decadal and long-term variation in precipitation. Moreover, the inclusion of upper level predictors did not improve the winter predictability in the southeastern Mediterranean. Consequently, higher predictability is identified for the west and the northern part of the Mediterranean and lower for the south and the east. At a smaller scale, the amount and distribution of precipitation in Greece is highly irregular in both the spatial and temporal dimensions. Maheras et al. (2004) examined the rainfall variability over Greece in relation to the 500 hPa circulation types (14 in total) established on the high correlation between seasonal precipitation amounts (predictant) and the frequency of the cyclonic circulation types at 500 hPa (predictors). The multiple linear regression (MLR) models per station achieved an adjusted R 2 in the range 0.2–0.8 for winter (December–February) and much lower for the other seasons. The authors conclude that the unexplained variance is due to phenomena that could not be explained only by changes in the circulation-type frequency. Dynamic regression models, in contrast to static (constant parameter) regression models, have time variable parameters (TVPs) evaluated through recursive optimisation (Young et al. (2004); Young (1984); Young Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 2077–2095 (2006) DOI: 10.1002/joc

ROBUST STOCHASTIC SEASONAL PRECIPITATION SCENARIOS

2079

(1999) and the references therein). The approach has a Bayesian origin with the evolution of the stochastic (unobserved) state space (SS) variables assumed to be described by a generalised random walk (GRW) process. The hyperparameters of the GRW process are estimated from the data unless they are known a priori. The skill of the dynamic models has been validated in several different applications, ranging from environmental to economic systems (Young; 1998). In summary, there is enough physical evidence that ENSO influences the NA region and in turn NA affect the European and Mediterranean climate. Those teleconnections are based on well-documented physical schemes (Cullen and deMenocal, 2000; Lamb and Peppler, 1987). This work attempts to further understanding of seasonal precipitation variability over Greece from the perspective of the changes in large-scale dynamics. Dynamic regression models are employed to cope with the non-stationarity of the predictor–predictant relationship. In particular, the main objectives of the analysis are: (1) To analyse the winter PRECA over Greece for 1951–2002 and explore its principal modes of oscillation using non-stationary time series analysis. (2) To quantify the influence of NAO/ENSO in the winter PRECA for Greece. (3) To develop robust winter precipitation predictions, at least one year in advance, for Greece (an area with lower predictability than the northwest part and the whole Mediterranean on average) through the use of dynamic models and explore their effectiveness. 2. DATA ANALYSIS 2.1. Predictant Monthly precipitation data for Greece were derived from 26 stations of the National Meteorological Service network (Figure 1) covering the period from October 1950 to March 2002. The stations are quite homogenously distributed over the domain. The generation of the annual time series of winter (PRECA can be viewed as a three-step procedure. First, we construct for each station monthly PRECA (point estimate) by subtracting the climatological mean (for the particular station and month over the period 1951–2002) from the respective monthly value. Then, winter PRECA (point estimate) are calculated for each station by taking the average of the six-monthly PRECA (from October to March). Finally, the winter PRECA for the whole study area (area estimate) is constructed by averaging winter rainfall anomalies using the cell declustering technique (Isaaks and Srivastava, 1989). According to the approach, the entire area (domain 20–27 ° E, 35–41° N) is divided into rectangular regions (cells) and each sample (meteorological station) receives a weight inversely proportional to the number of samples that fall within the same cell. The spatial distribution of the mean winter precipitation field is shown in Figure 2(a). An uneven distribution is evident with the rainiest areas located in the western part of Greece and the driest in the central part of the domain. The detrended winter PRECA (i.e. the dependent variable in the analysis) time series (Figure 2(b)) correctly identified the winter 1962–1963 as the wettest of the whole period and the winter 1989–1990 as the driest (Maheras and Anagnostopoulou, 2003). The spatial distribution of the PRECA for the extreme years (Figure 2(c,d)) shows that the western part of Greece is affected mostly by extreme weather events (drought and excess rainfall). 2.2. Predictors 2.2.1. Fields. Two sources of monthly NA sea-level pressure anomaly (SLPA) data are used: (1) Atlas of Surface Marine Data (DaSilva et al., 1994) spanning from January 1948 to December 1994 with a spatial resolution of 1° over the area 90.5 ° W–0.5 ° W and 0.5 ° N–89.5 ° N, comprising 2627 (sea) grid points. (2) Climate Data Assimilation System (NCEP/NCAR Reanalysis) spanning from January 1948 to the present with a spatial resolution of 2.5° (Kalnay et al., 1996) over the area 87.5 ° W–2.5 ° W and 2.5 ° N–87.5 ° N, comprising 382 (sea) grid points. Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 2077–2095 (2006) DOI: 10.1002/joc

2080

I. KIOUTSIOUKIS, S. RAPSOMANIKIS AND R. LOUPA 42

41

KAVALA

FLORINA KOZANI

40

LIMNOS

IOANNINA CORFU

LARISSA NEAGHIALOS

ARTA

39

ALEXANDROUPOLIS

THESSALONIKI

MYTILINI

Latitude

AGRINIO

38

PATRA

ELEYSINA ATHENS

ARGOSTOLI

SAMOS

TRIPOLI KALAMATA

37

KOS KITHIRA

36

RHODES RETHYMNO HERAKLIO

35

34 19

IERAPETRA

20

21

22

23

24

25

26

27

28

29

Longitude

Figure 1. Distribution of the meteorological stations over Greece

The steps to generate the input time series for the large-scale 2D fields (Eshel et al., 2000) are briefly reported here: (1) Noise filtering. A 3-month average of the 2D monthly SLPA dataset is created. Thus, for example, the December map of each year is replaced with the mean value of November–December–January. (2) Correlation pattern. The prognostic skill of the NA-SLPA over the detrended PRECA is investigated for lead times ranging from 1 to 18 (lagged teleconnection). Correlation maps are constructed between the PRECA time series and the detrended NA-SLPA time series at each grid point and for all considered lead times. Only locally significant correlations (p < 0.05) are reserved with the rest set to zero. (3) Field significance. The significance of each correlation map (18 in total) as a whole is examined through Monte Carlo analysis (Livezy and Chen, 1983). We fit an AutoRegressive Moving Average (ARMA) model to the PRECA time series and generate 1000 different realisations of it using its statistical properties (spectral properties and variance). The areal fraction of the domain with locally significant correlations is calculated for each realisation (and lead time) and the cumulative probability density function (cdf) is constructed for each lead. The 95th percentile of the cdf is taken as the field significance threshold value for that lead. (4) Time series projection. The detrended monthly (3-month average) NA-SLPA map of each year is projected on the correlation map Cl of the leads (l) that passed the significance tests yielding the annual lagged SLPA time series (the bar indicates the 3-month average SLPA and the superscript T the transpose of the matrix): T

SLPAl (t) = SLPAl Cl

(1)

We assume here that the 3D matrices have been reshaped into their 2D equivalents, i.e. the (X,Y,T) matrix with dimension NX*NY*NT has been reshaped to (S,T) with dimension (NX*NY)*NT. Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 2077–2095 (2006) DOI: 10.1002/joc

2081

ROBUST STOCHASTIC SEASONAL PRECIPITATION SCENARIOS 42 41 0 300 40 0 60000 7 800

0 50

60 0

400 300

38

0 50

700

37

500

400

600

Latitude

39

300

50

40

60

36

0

500

35 34 19

20

21

(a)

22

23

24 25 Longitude

26

27

28

29

40

Winter Rainfall Anomaly (mm/month)

30 20 10 0 −10 −20 −30 −40 1950 1954 1958 1962 1966 1970 1974 1978 1982 1986 1990 1994 1998 2002 (b)

Calendar Year

Figure 2. (a) Geographical distribution of winter (October–March) precipitation amounts (mm). (b) Winter rainfall anomaly time series (mm/month). (c) Geographical distribution of winter (October–March) PRECA (mm/month) for the wettest year (1962–1963). (d) Geographical distribution of winter (October–March) PRECA (mm/month) for the driest year (1989–19990)

2.2.2. Indices. In addition to the 2D fields, the following monthly climatic indices were exploited: • North Atlantic Oscillation Index Anomaly (NAOIA) (Climate Prediction Centre: http://www.cpc.ncep.noaa. gov/data/teledoc/telecontents.html). Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 2077–2095 (2006) DOI: 10.1002/joc

2082

I. KIOUTSIOUKIS, S. RAPSOMANIKIS AND R. LOUPA 42

120

39

40

60

38

80

40

00

1

20

60

37

80

Latitude

20

80

60

100

40

20

40

0 20

40

41

60 40

40

20

36 40

0

35 34 19

20

21

22

23

(c)

24 25 Longitude

26

27

28

29

28

29

42

−60

−8 0

39 38

−40

37

−40

0

−6

Latitude

−2 0

−60

−40

40

−4 0

41

0 −4

0

−8 0

−6

−2

0

−20

36 35 34 19

20

21

(d)

22

23

24 25 Longitude

26

27

Figure 2. (Continued)

• Southern Oscillation Index Anomaly (SOIA) (Climate Prediction Centre). • Tropical North Atlantic Index (TNAI) (anomaly of the average of the monthly SST from 5.5 ° N to 23.5 ° N and 15 ° W to 57.5 ° W) (Climate Diagnostic Center, http://www.cdc.noaa.gov/index.html). The selection of the input time series from the climatic indices is explored through correlation analysis for lead times up to 10 months. Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 2077–2095 (2006) DOI: 10.1002/joc

ROBUST STOCHASTIC SEASONAL PRECIPITATION SCENARIOS

2083

3. MODELLING SCHEMES The PRECA forecasts are established on inputs such as projected NA-SLPA and climatic indices at various lead times. The input–output mapping is calculated using dynamic regression models with time variable parameter such as dynamic linear regression (DLR), dynamic harmonic regression (DHR) and dynamic auto-regressive exogenous (DARX) variables. Traditional MLR models are also calculated for comparative purposes. The approach to TVP estimation in this paper is established on the basis of a multiple input–single output unobserved component s (UC) model of the general form: yt = Tt + St + f (ut ) + et

(2)

The subscript t symbolises the discrete-time of the observed values t = 1, 2, . . . , N (sample size). Tt is the trend element, St is the periodic (seasonal/cyclical) component, u denotes the vector of the auxiliary variables used for forecasting and finally et is the discrete-time white noise (∼N (0, σ 2 )) element. Equation (2) describes the observation equation of the SS model. The stochastic evolution of the state vector is assumed to be expressed by a GRW process: + Gi ηit I = 1, 2, . . . , k (analogous to the total number of TVPs) xti = Fi xt−1 i

(3)

with ηit the zero mean white noise vector (independent from the observation noise et ) with covariance Qi and Fi =

!

α 0

" ! " β δ 0 , Gi = γ 0 ε

(4)

The optimisation of the GRW model hyperparameters (elements of Fi , Gi and the noise variance ratio (NVR) defined as Qi /σ 2 ) is achieved in a recursive context accomplished by several approaches (maximum likelihood (ML) optimisation, minimisation of the sum of squares of the recursive multiple-step-ahead prediction error, frequency domain optimisation) (Young, 1999). Fixed interval smoothing (FIS) estimates (derived from Kalman’s work on optimal SS filter theory in the time domain; Kalman, 1960; Young and Pedregal, 1999) of the TVPs are then calculated recursively under the assumption that the stochastic SS state vector parameters vary as a GRW process (the term ‘fixed interval’ refers to the interval covered by the total sample size N ). Applications of the full UC model are generally limited due to problems in the simultaneous estimate of all the components. Indeed, special cases of the general UC model have been proven valuable practically (Young, 1999): (1) The dynamic linear regression (DLR) model is yt = Tt +

m # i=1

bit uti + et

(5)

where bit are the regression coefficients and uti the ‘regressors’ (auxiliary input variables). Note that when the regression coefficients are constant we obtain the normal regression model. (2) The dynamic harmonic regression (DHR) model is yt = Tt +

Rs # i=1

ait cos(ωi t) + bit sin(ωi t) + et

(6)

where the summation represents the variable amplitude periodic (seasonal or cyclical) term. Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 2077–2095 (2006) DOI: 10.1002/joc

2084

I. KIOUTSIOUKIS, S. RAPSOMANIKIS AND R. LOUPA

(3) The dynamic auto-regressive exogenous variables (DARX) model is yt = Tt +

n # i=1

−ait yt−i +

#

bit xt−δ−1 + et

(7)

where δ is a pure time delay. In DARX, in addition to auxiliary variables, past values of the output variable are also used. For DLR and DARX, the hyperparameters (NVR, elements of Fi and Gi ) are optimised via ML based on prediction error decomposition. The approach, although established on a strong theoretical basis, has some practical disadvantages such as the intense dependence on the length of the series, the shape of the likelihood surface (e.g. multiple minima), etc. (Young, 1999). To counter these limitations, DHR model utilizes a special form of optimisation in the frequency domain; specifically, optimisation is based on fitting the DHR model pseudo-spectrum to the logarithm of the AutoRegressive (AR) spectrum. The approach, although it is not as general as ML, leads to a considerably better-defined optimum in the objective function-hyperparameter space with consequent advantages to DHR modelling (convergence time and the number of parameters that can be optimised simultaneously). The final step in the analysis using linear dynamic modelling techniques would be the physical explanation of the nature of the parameter variation in terms of other variables. The analysis can then be further extended with the use of more complex dynamic models arriving at fully non-linear stochastic SS models. 4. RESULTS 4.1. Hindcast with DHR – analysis of the PRECA time series The identification of the frequency values in the DHR model is accomplished with reference to the AR spectrum. The Akaike information criterion (AIC) (Akaike, 1974) pointed out the best model order as AR(20) followed by AR(14). The AR(14) spectrum is selected in the analysis because it represents the spectral properties while being the most parsimonious. Figure 3(a) demonstrates the AR(14) spectrum; the peaks appear at periods close to 14, 7, 3.5 and 2 years indicating seasonal/cyclical oscillations with a fundamental period of 14 years with most significant periods at 14 and 3.5 years. Periodogram analysis of the North Atlantic Oscillation Index (NAOI) demonstrate a 7.5-year cyclical component (Werner and Schonwiese, 2002) for winter. Sutton and Allen (1997) identified the existence of decadal fluctuations (regular period of 12–14 years) in NA SST. Furthermore, typical ENSO cycles are between 3 and 4 years (Schneider and Schonwiese, 1989). Hence, PRECA cycles at first seem to be related to typical NAO and ENSO cycles. The analysis of whether this coincidence is random or not is investigated in the next paragraphs. We therefore built a DHR model with the five harmonics estimated from the data and a trend component. The 11 TVPs (Tt , ait , bit , i = 1, 2, 3, 4, 5) are modelled as random walk (RW) processes (the NVR values are as follows NVR(Tt ) = 0.2692, NVR(a1t ) = NVR(b1t ) = 0.0677, NVR(a2t ) = NVR(b2t ) = 0.00004, NVR(a3t ) = NVR(b3t ) = 0.0498, NVR(a4t ) = NVR(b4t ) = 0.00003 and NVR(a5t ) = NVR(b5t ) = 0.0023). The estimates of the TVPs for the particular interval (1951–2002) are then obtained in a recursive manner (Equations (6) – in its SS form – and 3) using a two-step (prediction-correction) version of Kalman filtering followed by an optimal smoothing procedure (FIS). The hindcasted PRECA is shown in Figure 3(b) and the seasonal components in Figure 3(c). The hindcast skill of the DHR model in terms of the coefficient of determination is R 2 = 0.94. The negative precipitation trend is a fact for the second half of the twentieth century in Greece (e.g. Amanatidis et al., 1993; Giorgi, 2002). We estimate this to 0.57 mm/winter/year as identified from the nonlinear trend component. Three harmonics have essential contribution to the total seasonal component. The fundamental harmonic (14 years) possess an amplitude increase and explains 22% of the PRECA variance. In addition, it correlates negatively with NAOIA (r = −0.41, p < 0.02). The 3.5 years harmonic exerts an Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 2077–2095 (2006) DOI: 10.1002/joc

2085

ROBUST STOCHASTIC SEASONAL PRECIPITATION SCENARIOS 4 3.5 3

log10(P)

2.5 2 1.5 1 0.5 0 −0.5

20

10

6.67

(a)

5 4 3.33 2.86 Period (samples/cycle)

2.5

2.22

50

Winter Rainfall Anomaly (mm/month)

40 30 20 10 0 −10 −20 Simulated PRECA (DHR) Observed PRECA Simulated Trend Standard Error

−30 −40 −50 1950

1955

1960

(b)

1965

1970

1975

1980

1985

1990

1995

2000

Calendar Year

Figure 3. (a) The AR(14) spectrum of the PRECA time series, (b) Comparison of DHR model output with measured data, (c) the TVPs

amplitude oscillation, explains 17% of the PRECA variance and correlates with SOIA (r = −0.52, p < 0.04). Finally, we observe that the local minima of the fundamental harmonic coincide with the plateau (at zero amplitude) of the 2 years harmonic, which explains 8% of the PRECA variance, demonstrating an interaction between the two oscillations. The origin of this mode is unknown (correlates weakly with all indices, including the Mediterranean Oscillation Index) and probably expresses the superposition-interaction of the existing forcings. Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 2077–2095 (2006) DOI: 10.1002/joc

2086

2y Period

2.5y Period

3.7y Period

6.7y Period

14y Period

I. KIOUTSIOUKIS, S. RAPSOMANIKIS AND R. LOUPA 20 0 −20 1950 20

1955

1960

1965

1970

1975

1980

1985

1990

1995

2000

1955

1960

1965

1970

1975

1980

1985

1990

1995

2000

1955

1960

1965

1970

1975

1980

1985

1990

1995

2000

1955

1960

1965

1970

1975

1980

1985

1990

1995

2000

1955

1960

1965

1970

1975

1980

1985

1990

1995

2000

0 −20 1950 20 0 −20 1950 20 0 −20 1950 20 0 −20 1950

(c)

Calendar Year Figure 3. (Continued)

Consequently, the principal cycles of the PRECA time series appear at periods of 14 and 3.5 years, explain together 40% of its variance and in addition they correlate significantly with NAOIA and SOIA, which show similar principal cycles. 4.2. Hindcast with MLR and DLR 4.2.1. Stationary models. Before analysing the detrended input–output time series for the period 1951–2002, we compare the results obtained using the two different SLPA (input) datasets for the overlapping period (1951–1992). The whole analysis is performed, according to Section 2.2.1, and MLR models are built based on the highest skill SLPA predictors. The two different input datasets revealed almost identical results as shown in Table I. These results also agree with the findings of Eshel et al. (2000), whose PRECA time series was averaged over the whole EM. The analysis is now performed to the full dataset (1951–2002). The leads that passed the field significance test (we focus on long-term forecast and hence we do not consider leads prior to 8 months) as well as the t-test of the lead’s hindcast skill in univariate mode are then combined in a backward stepwise MLR analysis. The best models in terms of R 2 , mean square error (MSE) and Malow’s Cp (Weisberg, 1985) are presented in Table II. The best bivariate model captures only 33% of the PRECA variability (the result was 49% for the 42-year dataset), whereas the overall best model has three predictors and achieves a Pearson correlation coefficient of 0.63. The reason for this notable change in the model skill, when data from the last decade are added, is shown in Figure 4. Weaker correlation fields accompanied with smaller fraction of the domain with significant correlations produce regressor time series with less significance. From the physical point of view, Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 2077–2095 (2006) DOI: 10.1002/joc

ROBUST STOCHASTIC SEASONAL PRECIPITATION SCENARIOS

2087

Table I. Hindcast skill (Pearson correlation between observed and modelled PRECA) of MLR models trained with lagged (by 9, 15, 16 and 18 months) SLPA predictors from different resolution SLPA datasets (I:DaSilva and II: CDAS) over 1951–1992. All correlations are t-test significant at p < 10−5 Predictors (lagged SLPA)

Pearson correlation (SLPA from DaSilva)

Pearson correlation (SLPA from CDAS)

0.70 0.71

0.70 0.68

(9, 16) (15, 18)

Table II. Hindcast skill of MLR models trained with different combinations of lagged (by 10, 12, 15 and 18 months) SLPA predictors from CDAS SLPA dataset (1951–2002). All correlations are t-test significant at p < 10−4 . MSE is the mean square error Predictors (12, (10, (15, (10, (12, (10, (10,

18) 18) 18) 12, 18) 15, 18) 15, 18) 12, 15, 18)

R2 0.33 0.30 0.29 0.35 0.40 0.35 0.41

MSE 119.3 124.8 126.8 115.7 106.8 115.3 105.8

Mallow’s Cp 7.0 9.5 10.4 7.4 3.5 7.2 5

those differences reflect the probabilistic nature of the teleconnection relationship arising from fluctuations in atmospheric circulation that are inherently unpredictable in some cases. From the mathematical point of view, there exist no hindcast univariate correlation significant at the 0.001 level (there were two at the 1951–1992 dataset: leads 9 and 16); the ratio of the percentage of the domain with significant points, calculated as 1951–2002 over 1951–1992 is 0.54 for the 9 month lead and 0.65 for the 16 month lead. We now explore the relationship of the PRECA with the climatic indices (Table III). Significant correlations are identified for lead times of 8–9 months. In other words, the winter PRECA over Greece is influenced by anomalies in the climatic indices (calculated based on SLP and SST anomalies) of the previous spring. This result matches the earlier finding based on NA-SLPA. Significance level is p < 0.001 for the highlighted NAOIA and SOIA and p < 0.01 for the highlighted TNAI. Thus, we will train an MLR model utilising the lagged by 8–9 months climatic indices (NAOIA MAY, SOIA APR and TNAI MAY). The MLR model achieved R 2 = 0.42 (MSE = 104.1) slightly better than the one based on the SLPA (12, 15, 18). Using all six inputs for training, the upper bound of R 2 in MLR mode reaches 0.55 (MSE = 80.3). 4.2.2. Non-stationary models. The spatial and temporal non-stationarity of the correlation field influenced considerably the hindcast skill of the MLR models trained with SLPA predictors. The same result was obtained for the MLR trained with climatic indices as predictors. The DHR model (Figure 3b), established on a probabilistic theoretical basis, tackled with the probabilistic nature of the teleconnections and achieved a hindcast score of R 2 = 0.94. We now explore the behaviour of another dynamic model trained with auxiliary variables, that is, the DLR model. The DLR model is built using as regressors the set – NAOIA MAY, SOIA APR, TNAI MAY, SLPA12, SLPA15 and SLPA18, i.e. the winter PRECA is analysed in terms of SLPA variations of the previous winter (12 m), the previous early autumn (15 m) and the previous early summer (18 m) as well as with variations in the climatic indices of the previous spring. The six TVPs (bit , i = 1, 2, 3, 4, 5 Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 2077–2095 (2006) DOI: 10.1002/joc

2088

I. KIOUTSIOUKIS, S. RAPSOMANIKIS AND R. LOUPA Lead Time 9 months 0.3

0.4

0.3

.3 −0

20 −80

−60 −40 −20 Lead Time 10 months

40

−0.3

0. −0.4 0.4 .5 3 −0

20

−0.3

−80

−60

Latitude

80 0.4

60 40

−40

−80

−20

Latitude

80

60

−0.3

40 0.3

−80

−0.3 −0.4

20

Lead Time 12 months

0.3

20

40

−60 −40 −20 Lead Time 10 months

60

80

3

0.3

3

Latitude

80

−0.3

−0.5 −0 .4−0.3

0.

60

−0.4

40

−80

Latitude

0.40.3

0.3

20

Latitude

60

Lead Time 9 months 80

−0.

Latitude

80

−60

0.3

60 40

−0.3

−60 −40 Longitude

−80

−20

−60 −40 Longitude

80

−0.5 −0.3 −0.4

60

−0.3 4 0.

40

Latitude

Latitude

−20

Lead Time 15 months

Lead Time 15 months

20

60

−0.3

40

0.3

20 −80

−60

−40

−20

−80

−60

−40

−20

Lead Time 16 months

Lead Time 16 months 80

60

Latitude

80

−0.3

40

0.

4

Latitude

−20

20

80

−80

−60

−40

60 40 20

0.3

−0.3 −0.4

20

−20

0.3 −80

−60

−40

−20

Lead Time 18 months

Lead Time 18 months 80

40 0.4

−80

0.4

−60 −40 Longitude

40 20

0.3

20

60

−20

0.3

0.3

0.3

60

Latitude

80 Latitude

−40

Lead Time 12 months

−80

0.3

−60

−40

−20

Longitude

Figure 4. Correlation patterns between PRECA and SLPA from CDAS-1 dataset over 1951–1992 (left) and 1951–2002 (right) for lead times 9, 10, 12, 15, 16 and 18 months. This figure is available in colour online at www.interscience.wiley.com/ijoc

Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 2077–2095 (2006) DOI: 10.1002/joc

ROBUST STOCHASTIC SEASONAL PRECIPITATION SCENARIOS

2089

Table III. Pearson correlation between PRECA and climatic indices. Correlations are significant at p < 0.01

J F M A M J J A S O N D

NAOIA

SOIA

TNAI

– – – – −0.47 – – –

– – −0.38 −0.48 – – – – – – – –

– – 0.33 0.38 0.41 0.36 – – – – 0.41 0.37

– –

Values in bold are significant at p < 0.001

and 6) are modelled as both random walk and integrated random walk (IRW) processes (The NVR’s obtained after ML unconstrained optimisation for RW (IRW) are given as follows NVR(NAOIA) = 7.29 × 10−3 (1.46 × 10−4 ), NVR(SOIA) = 7.88 × 10−9 (3.36 × 10−13 ), NVR(TNAI) = 4.00 × 10−7 (7.87 × 10−8 ), NVR(SLPA12) = 1.43 × 10−7 (3.12 × 10−9 ), NVR(SLPA15) = 3.76 × 10−6 (4.43 × 10−27 ), NVR(SLPA18) = 1.48 × 10−6 (1.57 × 10−26 )). FIS estimates of the TVPs are shown in Figure 5. The coefficient of model determination is R 2 = 0.71 for DLR and R 2 = 0.55 for the constant parameter MLR shown dotted; DLR model explain better with the large deviations from the mean. The statistical significance of the estimated TVPs is examined through the properties of the normalised recursive residuals. If the assumptions regarding the dynamic regression problem are satisfied, those quantities should be a zero mean, serially uncorrelated sequence of random variables with changing variance (Young, 1984). The residuals were found statistically independent of the input functions as validated from the lag-correlation of the autocorrelation function and the cross-correlation function. The best models built on two, three, four, five and six regressors are also given in Table IV. The SLPA auxiliary variables have higher variability than the climatic variables. This variation is considered in the dynamic models resulting in lower MSE, even when the correlation coefficient gives indistinguishable results. If we constrain the NVRs to higher values (NVR = 0.01) than the ones estimated after optimisation, we obtain a model that ‘fits’ the PRECA time series (R 2 = 0.9888). The negative effect is less smooth variations (Figure 5(b)) of the regression coefficients (especially for the coefficients of SLPA). Those significant parameter variations then suggest the use of higher order non-linear models. However, such an analysis is prohibitive due to the limited length of the time series. Alternatively, if possible, one should explain the variation of the TVPs in terms of known variables (to improve the forecast skill of the model). We have demonstrated that the analysis with dynamic models is superior in comparison with the constant parameter models. However, their forecast skill must be evaluated before any inferences are made. 4.3. Forecast The ability of the models (MLR, DLR and DHR) to forecast one year in advance the PRECA is now examined. The 52-year time series due to its autocorrelation has approximately 24 degrees of freedom; hence a forecast for one degree of freedom should be applied to 2.2 years. We apply the omission approach and sequentially withdraw 3 years of data (50 cases) and repeat the whole procedure presented in Section 2 for the censored dataset. Then we test the skill of both the constant parameter models and the dynamic models to correctly forecast over the gap in the input data. The observed time series is compared with the mean forecast Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 2077–2095 (2006) DOI: 10.1002/joc

2090

I. KIOUTSIOUKIS, S. RAPSOMANIKIS AND R. LOUPA 40

Winter Rainfall Anomaly (mm/month)

30 20 10 0 −10 −20 −30 −40 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 2000 Calendar Year

(a)

10

1960

1970

1980

1990

−5 −10 1950 10

1960

1970

1980

1990

1960

1970

1980

1990

bSLPA12

1980

1990

−0.2 1950 0.2

1960

1970

1980

1990

2000

0.1 1960

1970 1980 1990 Calendar Year

2000

1980

1990

2000

1960

1970

1980

1990

2000

1960

1970

1980

1990

2000

1960

1970

1980

1990

2000

1960

1970

1980

1990

2000

1970 1980 1990 Calendar Year

2000

0 −0.2 1950 0.5

2000 bSLPA15

1970

0

0 1950 (b)

1960

bSLPA18

bSLPA18

bSLPA15

0 1950 0.2

1970

0 −20 1950 0.2

2000 bSLPA12

−10 1950 0.05

1960

−5 −10 1950 20

2000

0

0 −10 1950 0

2000 bSOIA

bSOIA

−10 1950 0

bTNA

bNAOA

0

bTNA

bNAOA

10

0 −0.5 1950 1 0 −1 1950

1960

Figure 5. (a) Comparison of DLR model output (RW: full line, IRW: full line with dots) with measured data (circles) for unconstrained optimisation. The output of the constant parameter model is shown as dotted. (b) The corresponding TVPs, for unconstrained optimisation (left) and for TVPs constrained to equal values (right) (dashed-dotted line: RW, full line: IRW, dotted line: MLR coefficient) Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 2077–2095 (2006) DOI: 10.1002/joc

2091

ROBUST STOCHASTIC SEASONAL PRECIPITATION SCENARIOS

Table IV. Hindcast skill of MLR and DLR (unconstrained optimisation) models trained with lagged SLPA (by 12, 15 and 18 months) predictors and climate indices (1951–2002). Notation: Climate Group (C1: NAOIA MAY, C2: SOIA APR and C3: TNAI MAY), SLPA Group (S1: SLPA12, S2: SLPA15 and S3: SLPA18). Correlations are t-test significant at p < 10−4 . CC is the Pearson Correlation Predictors

C(1,2,3) S(1,2,3) C(1,2) C(1,2)S(2) C(1,2)S(1,2) C(1,2)S(1,2,3) C(1,2,3)S(1,2,3)

MLR

DLR

CC

MSE

CC

MSE

0.65 0.63 0.61 0.67 0.71 0.73 0.74

104.1 106.8 112.8 97.1 87.9 82.6 80.3

0.66 0.76 0.63 0.74 0.77 0.84 0.84

100.1 77.4 108.6 80.1 71.9 54.1 53.1

Values in bold are significant at p < 10−5 ). CC is the Pearson correlation and MSE the mean square error.

Table V. Forecast skill of MLR and DLR (unconstrained optimisation). Notation: see Table IV. Correlations are t-test significant at p < 10−3 Predictors

MLR

C(1,2,3) S(1,2,3) C(1,2)S(2) C(1,2)S(1,2) C(1,2,3)S(1,2,3) Value with

∗

DLR

CC

MSE

CC

MSE

0.57 0.32(∗ ) 0.53 0.56 0.54

122.1 171.2 130.7 123.4 130.1

0.60 0.47(∗∗ ) 0.61 0.64 0.64

115.5 142.0 111.7 106.5 107.3

is significant at p < 0.03; value with

∗∗

is significant at p < 0.01).

time series, calculated as the mean of the forecasts for each year (there exist one forecast for 1951 and 2002, two forecasts for 1952 and 2001 and three forecasts for all other years). Table V demonstrates the results for DLR/MLR models (unconstrained optimisation). Forecasting using only climate indices is more robust than forecasting based on SLPAs. In forecast mode, the projected SLPA time series is recomputed each time after the sequential withholding, repeating all the steps given in Section 2 (noise filtering, correlation pattern, field significance and time series projection). The recalculation of the correlation field for each censored dataset assigns high variability to the SLPA time series (i.e. the SLPA time series is very sensitive to the correlation field) resulting in higher forecast skill to the model utilising only climate indices and lower to the model built only on SLPAs (also recall that higher significance correlations with PRECA were observed for climatic predictors than for SLPA predictors in the full dataset), although in hindcast mode their skill was equal (however, the forecast results are also alike when random withholding is applied). For comparative reasons, the MLR forecast built only on SLPAs for the period 1951–1992 achieved r = 0.51. Mixed inputs improve the skill only in DLR mode, while in MLR mode the addition of SLPA regressors does not always give better forecasts. Figure 6 demonstrates results for DLR/MLR using all six regressors. DLR models explain more of the PRECA variability in comparison with the stationary models, although they require additional information (forecasting of the regression coefficients variation for the withheld sample). The scatterplots of the forecast Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 2077–2095 (2006) DOI: 10.1002/joc

2092 40

40

30

30

20

20

10

10

0 −10

−30

−30 −20

0 DLR Forecast

20

40

−40 −40

−20

0 MLR Forecast 0.5

0.4

0.4

0.4

0.3

0.3

0.3

0.2

0

SRCTNA

0.5

SRCSOIA

0.5

0.2 0.1

0.1

0

0

200 400 Forecast MSE

0

0

200 400 Forecast MSE

0.4

0.4

0.4

0.3 0.2 0.1

0.1 0

SRCSLPA18

0.5

0.2

0

200 400 Forecast MSE

0

40

0.1

0.5

0.3

20

0.2

0.5

SRCSLPA15

SRCSLPA12

−10 −20

(a)

(b)

0

−20

−40 −40

SRCNAOA

Observed

Observed

I. KIOUTSIOUKIS, S. RAPSOMANIKIS AND R. LOUPA

0

200 400 Forecast MSE

0

200 400 Forecast MSE

0.3 0.2 0.1

0

200 400 Forecast MSE

0

Figure 6. (a) Forecast scores versus observations for DLR model (left) and MLR (right) model. (b) Forecast error (MSE) versus standardised regression coefficients

error (MSE) in terms of the standardised regression coefficients (Figure 6(b)) shows clearly the principal importance of the NAOI input. The nominal values (hindcasting mode) are shown as dotted. Significant improvement of the results by means of constrained optimisation, as seen before, is highly improbable due to the higher variability of the regression coefficients. The analysis in now repeated for the DHR model. Figure 7 illustrates the DHR forecasts after sequential 3-year withdrawal and their associated standard errors, with the NVRs modelled as RW (top) and IRW Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 2077–2095 (2006) DOI: 10.1002/joc

ROBUST STOCHASTIC SEASONAL PRECIPITATION SCENARIOS

2093

80

Winter Rainfall Anomaly (mm/month)

60 40 20 0 −20 −40 −60 −80 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 2010 (a) Calendar Year 80

Winter Rainfall Anomaly (mm/month)

60 40 20 0 −20 −40 −60 −80 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 2010 (b)

Calendar Year

Figure 7. Forecast scores for DHR model, for TVPs modelled as RW (top) and IRW (bottom)

(bottom). The Pearson correlation between forecasts and observations is 0.77 and 0.84, for RW and IRW, respectively, while the respective MSE is 114.2 and 79.7. IRW statistics are better, however, those forecasts are associated with higher standard errors. Forecasting is less complicated when using DHR than DLR because the regressors in the former are wellknown periodic functions (sine and cosine) that can be easily identified inside or outside the sample. The DHR results are also quite unaltered even for larger segments of missing data. The high robustness of the DHR model allows us to force the model to forecast beyond the end of the series. An 8 years forecast (years Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 2077–2095 (2006) DOI: 10.1002/joc

2094

I. KIOUTSIOUKIS, S. RAPSOMANIKIS AND R. LOUPA

2003–2010) is also presented in Figure 7. Here we assume than an 8-year period is reasonably short for retaining the same dynamics in the system (e.g. the human-induced pressure is non-radically different). We could expect, given the uncertainty expressed by the standard errors, the winter during 2004–2005 will be drier in comparison with the climatological mean.

5. CONCLUSIONS The relationship between winter PRECA over Greece and disturbances linked to SLP and SST is investigated with the aim to develop robust seasonal precipitation predictions. Teleconnections between NA and EM as well as between tropical and extra-tropical climate form the scientific basis for the predictability. Initially, the data are validated against published results for the EM over the period 1951–1992. Subsequently, fluctuations in the EM – NA teleconnection were identified and quantified, which limited radically the applicability of traditional MLR models. The effectiveness of two (DLR and DHR) sophisticated dynamic regression models was then examined as the most natural choice to manage non-stationary processes. Significant parameter variation was found in the case of DLR that was growing with increasing input lag. However, the development of non-linear state dependent regression models was prohibited from the time series length. Even with such an inherent limitation, DLR achieved more than 10% higher R 2 than traditional regression models in forecast mode, with the latter failing systematically in predicting high variations from the mean. On the other hand, the analysis of the PRECA time series through the DHR model pointed out some fundamental properties of the spectrum like its harmonic terms and their associated seasonal components. It was found that the amplitude of the most influential seasonal components with periods of 14 and 3.5 years is significantly varying over time and specifically increasing in the former and oscillating in the latter; additionally, the respective terms correlate significantly with NAOIA and SOIA, indicating their climatic origin. Finally, the forecasting capabilities of DHR were examined. Model predictions appeared very robust when missing data were included in the analysis reaching a forecast skill (r = 0.84) higher even than the hindcast skill of the traditional regression models. The dynamic model with TVPs reflecting the inherent variability of the time series produced the best results. This is explained partly by the periodicity of the regression coefficients in the latter as well as to the higher total number of TVPs. Given the uncertainty arising from the internal variability of the atmosphere and the variations in the strength and spatial distribution of SST anomalies, the dynamic scheme provided robust forecasts accompanied with uncertainty estimates. The forecast, although in aggregated form, could be always used in the direction of reduced risk/cost or increased production by, e.g. farmers. Its influence is currently investigated in the skill of a stochastic daily downscaling scheme for the same region. ACKNOWLEDGEMENTS

The authors wish to thank Prof. P. Young and C. Taylor for providing an extended license of the CAPTAIN toolbox (Young et al., 2004) and Associate Prof. D. Melas for the fruitful discussions. The authors also wish to thank the two anonymous reviewers for their significant comments on the earlier version of this manuscript. REFERENCES Akaike H. 1974. A new look at the statistical model identification. IEEE Transactions on Automatic Control 19: 716–723. Amanatidis G, Paliatsos A, Repapis C, Barttzis J. 1993. Decreasing precipitation trend in the Marathon area, Greece. International Journal of Climatology 13: 191–201. Cassou C, Terray L. 2001. Dual Influence of Atlantic and Pacific SST anomalies on the North Atlantic/Europe winter climate. Geophysical Research Letters 28: 3195–3198. Cassou C, Terray L, Hurrel J, Deser C. 2004. North Atlantic winter climate regimes: Spatial asymmetry, stationarity with time, and oceanic forcing. Journal of Climate 17: 1055–1068. Conte M, Giuffrida A, Tedesco S. 1989. The Mediterranean Oscillation. Impact on Precipitation and Hydrology in Italy. Conference on Climate, Water. Publications of the Academy of Finland: Helsinki; 121–137. Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 2077–2095 (2006) DOI: 10.1002/joc

ROBUST STOCHASTIC SEASONAL PRECIPITATION SCENARIOS

2095

Cullen HM, deMenocal PB. 2000. North Atlantic influence on tigris-euphrates streamflow. International Journal of Climatology 20: 853–863. DaSilva A, Young A, Leritus S. 1994. Atlas of Surface Marine Data. Algorithms and Procedures. NOAA Atlas NESDIS 6, U.S , Vol. 1. Department of Commerce: Washington, DC. Dunkeloh A, Jacobeit J. 2003. Circulation dynamics of Mediterranean precipitation variability 1948–98. International Journal of Climatology 23: 1843–1866. Eshel G, Farrell B. 2000. Mechanisms of eastern Mediterranean rainfall variability. Journal of the Atmospheric Sciences 57: 3219–3232. Eshel G, Cane M, Farrell B. 2000. Forecasting eastern Mediterranean Droughts. Monthly Weather Review 128: 3618–3630. Fernandez J, Saenz J, Zorita E. 2003. Analysis of wintertime atmospheric moisture transport and its variability over Southern Europe in the NCEP-Reanalyses. Climate Research 23: 195–215. Giorgi F. 2002. Variability and trends of sub-continental scale surface climate in the twentieth century. Part I: observations. Climate Dynamics 18: 675–691. Goddard L, Mason SJ, Zebiak SE, Ropelewski CF, Basher R, Cane MA. 2001. Current approaches to seasonal-to-interannual climate predictions. International Journal of Climatology 21: 1111–1152. Hansen DV, Bezdek HE. 1996. On the nature of decadal anomalies in North Atlantic sea surface temperature. Journal of Geophysical Research 101: 8749–8758. Haylock M, Goodness C. 2004. Interannual variability of European extreme winter rainfall and links with mean large-scale circulation. International Journal of Climatology 24: 759–776. Hoerling MP, Hurrell JW, Xu T. 2001. Tropical origins for recent North Atlantic climate change. Science 292: 90–92. Hulme M, Barrow E, Arnell N, Harrison P, Johns T, Downing T. 1999. Relative impacts of human-induced climate change and natural variability. Nature 397: 688–691. Hurrell JW. 1995. Decadal trends in the North Atlantic oscillation: Regional temperatures and precipitation. Science 269: 676–679. Hurrell JW. 2001. The North Atlantic Oscillation. Science 291: 603–605. Isaaks EH, Srivastava RH. 1989. Applied Geostatistics. Oxford University Press: Oxford. Kalman R. 1960. A new approach to linear filtering and prediction problems. ASME Transactions Journal Basic Engineering 83D: 95–108. Kalnay E, Kanamitsu M, Kistler R, Collins W, Deaven D, Gandin L, Iredell M, Saha S, White G, Woollen J, Zhu Y, Chelliah M, Ebisuzaki W, Higgins W, Janowiak J, Mo KC, Ropelewski C, Wang J, Leetmaa A, Reynolds R, Roy J, Dennis J. 1996. The NMC/NCAR 40-year reanalysis project. Bulletin of the American Meteorological Society 77: 437–471. Kushnir Y. 1999. Europe’s winter prospects. Nature 398: 289–291. Lamb PJ, Peppler RA. 1987. North Atlantic oscillation: Concept and application. Bulletin of the American Meteorological Society 68: 1218–1225. Livezy R, Chen W. 1983. Statistical field significance and its determination by Monte-Carlo techniques. Monthly Weather Review 111: 46–59. Maheras P, Anagnostopoulou C. 2003. Circulation types and their influence on the interannual variability and precipitation changes in Greece. In Mediterranean Climate, Variability and Trends. Springer-Verlag: Berlin, 215–239. Maheras P, Tolika K, Anagnostopoulou C, Vafiadis M, Patrikas I, Flocas H. 2004. On the relationship between circulation types and changes in rainfall variability in Greece. International Journal of Climatology 24: 1695–1712. Mathieu PP, Sutton RT, Dong B, Collins M. 2004. Predictability of winter climate over the North Atlantic European region during ENSO events. Journal of Climate 17: 1953–1974. Palutikof JP, Conte M, Casimiro Mendes J, Goodess CM, Espirito Santo F. 1996. Climate and climate change. In Mediterranean Desertification and Land Use, Brandt CJ, Thornes JB (eds). John Wiley and Sons: London. Rodo X, Baert E, Comin F. 1997. Variations in seasonal rainfall in southern Europe during the present century: relationships with the North Atlantic oscillation and the El Nino-southern oscillation. Climate Dynamics 13: 275–284. Rodwell MJ, Rowell DP, Follan CK. 1999. Oceanic forcing of the wintertime North Atlantic oscillation and European climate. Nature 398: 320–323. Schneider U, Schonwiese C. 1989. Some statistical characteristics of El Nino/southern oscillation and North Atlantic oscillation indices. Atmosfera 2: 167–180. Sutton RT, Allen MR. 1997. Decadal predictability of North Atlantic sea surface temperature and climate. Nature 388: 563–567. Trenberth K. 1997. The definition of El Ni˜no. Bulletin of the American Meteorological Society 78: 2771–2777. Wanner H, Bronnimann S, Casty C, Gyalistras D, Luterbacher J, Schmutz C, Stephenson D, Xoplaki E. 2001. North Atlantic oscillation–concepts and studies. Surveys in Geophysics 22: 321–382. Weisberg S. 1985. Applied Linear Regression. Wiley & Sons: New York. Werner A, Schonwiese C. 2002. A statistical analysis of the North Atlantic oscillation and its impact on European temperature. The Global Atmosphere and Ocean System 8: 293–306. Xoplaki E, Gonzalez-Rouco J, Luterbacher J, Wanner H. 2004. Wet season Mediterranean precipitation variability: influence of largescale dynamics and trends. Climate Dynamics 23: 63–78. Young P. 1984. Recursive Estimation and Time Series Analysis. Springer-Verlag: Berlin. Young P. 1998. Data-based mechanistic modelling of environmental, ecological, economic and engineering systems. Environmental Modelling & Software. 13: 105–122. Young P. 1999. Nonstationary time series analysis and forecasting. Progress in Environmental Science 1: 3–48. Young P, Pedregal D. 1999. Recursive and en-bloc approaches to signal extraction. Journal of Applied Statistics 26: 103–128. Young PC, Taylor CJ, Tych W, Pedregal DJ, McKenna PG. The Captain Toolbox. Centre for Research on Environmental Systems and Statistics. Lancaster University: 2004. (www.es.lancs.ac.uk/cres/captain).

Copyright  2006 Royal Meteorological Society

Int. J. Climatol. 26: 2077–2095 (2006) DOI: 10.1002/joc