Natural Hazards and Earth System Sciences Evaluation of a statistical downscaling procedure for the estimation of climate change impacts on droughts

Despite uncertainties in future climates, there is considerable evidence that there will be substantial impacts on the environment and human interests. Climate change will affect the hydrology of a region through changes in the timing, amount, and form of precipitation, evaporation and transpiration rates, and soil moisture, which in turn affect also the drought characteristics in a region. Droughts are long-term phenomena affecting large regions causing significant damages both in human lives and economic losses. The most widely used approach in regional climate impact studies is to combine the output of the General Circulation Models (GCMs) with an impact model. The outputs of Global Circulation Model CGCMa2 were applied for two socioeconomic scenarios, namely, SRES A2 and SRES B2 for the assessment of climate change impact on droughts. In this study, a statistical downscaling method has been applied for monthly precipitation. The methodology is based on multiple regression of GCM predictant variables with observed precipitation developed in an earlier paper (Loukas et al., 2008) and the application of a stochastic timeseries model for precipitation residuals simulation (white noise). The methodology was developed for historical period (1960–1990) and validated against observed monthly precipitation for period 1990–2002 in Lake Karla watershed, Thessaly, Greece. The validation indicated the accuracy of the methodology and the uncertainties propagated by the downscaling procedure in the estimation of a meteorological drought index the Standardized Precipitation Index (SPI) at multiple timescales. Subsequently, monthly precipitation and SPI were estimated for two future periods 2020–2050 and 2070–2100. The results of the present study indicate the accuracy, reliability and uncertainty of the statistical downscaling method for the assessment of climate change on hydrological, agricultural and Correspondence to: A. Loukas (aloukas@civ.uth.gr) water resources droughts. Results show that climate change will have a major impact on droughts but the uncertainty introduced is quite large and is increasing as SPI timescale increases. Larger timescales of SPI, which, are used to monitor hydrological and water resources droughts, are more sensitive to climate change than smaller timescales, which, are used to monitor meteorological and agricultural droughts. Future drought predictions should be handled with caution and their uncertainty should always be evaluated as results demonstrate.


Introduction
Rainfall varies considerably over space and time.Agricultural and water resources systems have evolved in response to this variability, but in most regions of the world, rainfall variability continues to be a major source of risks that water resources managers face.Depending on spatial extent and persistence of drought, for example, entire communities and regions risk economic and food security problems.Research is being conducted to better understand climate variability, its impacts on agricultural and water resources systems, and how to reduce those risks through decisions and policies that consider climate variability.Nowadays anthropogenic climate change and its socioeconomic impacts are major concerns of mankind.Global surface temperature has been increased significantly during the last century and will continue to rise unless greenhouse gas emissions are drastically reduced (IPCC, 2007).Climate change effects are manifold and vary regionally, even locally, in their intensity, duration and areal extent.However, immediate damages to humans and their properties are not obviously caused by gradual changes in temperature or precipitation but mainly by socalled extreme events such as floods and droughts.The frequency and intensity of extreme events can be analysed with the use of long historical data series which are unavailable in Published by Copernicus Publications on behalf of the European Geosciences Union.
L. Vasiliades et al.: Downscaling climate change and droughts many parts of the world.Hence, coupled atmosphere-ocean general circulation models are suitable tools to simulate extreme events since there are able to generate long timeseries that can be used for model evaluation and also for analyses of possible future changes in extreme events.
However, there is a mismatch between the grid resolution of current climate models (generally hundreds of kilometers), and the resolution needed by environmental impacts models (typically ten kilometers or less).Techniques have been developed to downscale information from GCMs to regional scales.Downscaling is the process of transforming information from climate models at coarse resolutions to a fine spatial resolution.Downscaling is necessary, as the underlying processes described by the environmental impact models are very sensitive to local climate, and the drivers of local climate variations, such as topography, are not captured at coarse scales.There are two broad categories of downscaling: dynamic (which simulates physical processes at fine scales) and statistical (which transforms coarse-scale climate projections to a finer scale based on observed relationships between the climate at the two spatial resolutions) (IPCC, 2007).Dynamic downscaling, nesting a fine scale climate model in a coarse scale model, produces spatially complete fields of climate variables, thus preserving some spatial correlation as well as physically plausible relationships between variables.However, dynamic downscaling is very computationally intensive, making its use in impact studies limited, and essentially impossible for multi-decade simulations with different global climate models and/or multiple greenhouse gas emission scenarios.Thus, most impacts studies rely on some form of statistical downscaling, where variables of interest can be downscaled using historical observations.These relationships are empirical (i.e.calibrated from observations) and they are applied using the predictor fields from GCMs in order to construct scenarios.There are applications related criteria that contribute to an appropriate choice of downscaling method in a particular context (Mearns et al., 2004;Wilby et al., 2004).However, there are assumptions involved in both techniques which are difficult to verify a priori and contribute to the uncertainty of results (Giorgi et al., 2001).There has been extensive work developing and intercomparing statistical downscaling techniques for climate impact studies (Wilby and Wigley, 1997;Xu 1999;Giorgi et al., 2001;Varis et al., 2004;Xu et al., 2005;Fowler et al., 2007).
Most general circulation models predict a prominent change in precipitation (IPCC, 2007), supported by observations of precipitation trends (National Observatory of Athens, 2001) showing decreased winter precipitation and enhanced variability (IPCC, 2007).There is evidence that such changes are now reflected in low flows and hydrologic droughts (Hisdal et al., 2001).The frequency and severity of low flows has been extensively studied (Smakhtin, 2001).In contrast to various climate change drought studies of river discharge, limited studies of drought based on meteorological drought indices, which require considerably less input data when compared to weather, soil and land use information needed by meteorological, hydrologic, agrohydrologic and water management models, have been performed (i.e.Kothavala, 1999;Blenkinsop and Fowler, 2007;Loukas et al., 2007b;Mavromatis, 2007;Loukas et al., 2008;Dubrovski et al., 2009).This study, a continuation study of Loukas et al. (2008), examines, explicitly for the Lake Karla Watershed in Thessaly, Greece, whether the upward trend of droughts (IPCC, 2007;Weiss et al., 2007), as described above, could be depicted by a statistical downscaling method for precipitation using a global circulation model.It is used to reproduce present drought conditions; and, by reconstructing climatic records including climate and socioeconomic changes on future drought severities using two of the IPCC global emission scenarios, SRES A2 and SRES B2, to assess the uncertainty introduced to climate change impact studies on droughts.The methodology is based on multiple linear regression of GCM predictant variables with observed monthly precipitation, developed by Loukas et al. (2008) and extended in this study by a stochastic timeseries model component for regression residuals simulation (white noise).The methodology was developed for the base historical period  and validated against observed precipitation for the period 1990-2002.Subsequently, comparison of the Standardized Precipitation Index (SPI) timeseries calculated from observed and downscaled meteorological parameters will indicate the accuracy, reliability and uncertainty of the downscaling method for present and future climate conditions and the use of the downscaling method on climate impact studies in hydrology, agriculture and water resources.

Study area and characteristics of droughts in the region
Lake Karla watershed is located in central Thessaly, Greece and is a plain region surrounded only by eastern high mountains (Fig. 1).It has an area of about 1171 km 2 .Elevation ranges from 50 m to more than 1900 m, and the mean elevation of the region is about 230 m.The plain is one of the most productive agricultural regions of Greece.The main crops cultivated in the plain area are cotton, wheat and maize whereas apple, apricot, cherry, olive trees and grapes are cultivated at the foothills of the eastern mountains.The climate is typical continental with cold and wet winters and hot and dry summers.Mean annual precipitation in Lake Karla watershed is about 560 mm and it is distributed unevenly in space and time.
In the Mediterranean Basin, and especially in Greece, the major methodological drawback for a long-term assessment of regional climate and its variability comes from the lack of suitable observations or simulated data.Global reanalyses databases have been created to overcome this obstacle (Kalnay et al., 1996;Gibson et al., 1997;Sotillo et al., 2005).However, the coarse spatial resolution of global reanalysis make these data sets inadequate tool to characterize regional, prevailing atmospheric conditions over areas where orography and land-sea contrasts are valuable (Morata et al., 2008).Processed monthly precipitation data from 12 precipitation stations for the period October 1960 to September 2002 were used (Fig. 1).The mean areal precipitation of Lake Karla watershed was estimated by the Thiessen polygon method modified by the precipitation gradient using the stations, which are within or in the vicinity of the watershed.Thessaly, and especially Lake Karla watershed, experienced severe, extreme and persistent droughts during the periods from mid to late 1970s, from late 1980s to early 1990s and the first years of 2000s (Loukas et al., 2007b).These three drought periods were quite remarkable and affected large areas.The first drought episode (1976)(1977) affected southern and western Europe, the second drought episode (1988)(1989)(1990)(1991) affected the whole Mediterranean Region with an estimated economic cost lager than 2.1 billion Euros, whereas the third drought episode (2000)(2001) affected Central Europe and the Balkans with total damage of 0.5 billion Euros (EEA, 2004).During these three periods the monthly and annual precipitation was significantly bellow normal in Thessaly.The prolonged and significant decrease of monthly and annual precipitation has a dramatic impact on natural vegetation, agricultural production and the water resources of the region (Loukas et al., 2007a).Large scale atmospheric circulation patterns affect the droughts over Greece and the Mediterranean basin, in general.In a recent study (Bordi et al., 2007) analysis of geopotential height anomaly of 500 mb indicated that a high positive anomaly over North-Eastern Europe is responsible for extended and severe droughts in Italy and Greece.These circulation patterns characterise mid-to high-latitude flow anomalies.These dipole-like geopotential anomalies characterize the large-scale circulation and produce long persistent droughts.Especially, the 1988-1991 drought episode has been observed during a high positive North Atlantic Oscillation (NAO) index (Xoplaki et al., 2000;Houssos and Bartzokas, 2006).During this period, the extension of the subtropical anticyclone of the Atlantic (Azores) up to central Mediterranean modified the tracks of the traveling depressions affecting precipitation in NW Greece.Furthermore, during this period, low pressure systems approached Greece mainly from the North, causing dry katabatic winds in NW Greece due to the NW-SE orientation of the Pindus mountain range, west of Thessaly.These atmospheric circulation patterns are considered typical for extreme dry periods and have been identified by many researchers (Xoplaki et al., 2000;Bartzokas et al., 2003).
Climate change with have a remarkable impact on future climate in Greece.Multimodel GCM experiments show a mean annual temperature increase of 0.5-1 • C for the Mediterranean region which is insensitive to the choice among Special Report on Emission Scenarios (SRES) for the period 2011-2030(IPCC, 2007)).In a recent study, where climate change impacts on temperature and precipitation were investigated on Greece using nine Regional Climate Models (RCMs), mean annual temperature will be increased by 3.7 • C and precipitation will be decreased by 15.8% for the period 2070-2100.The inter-annual variability of temperature will be increased in summer and reduced at winter, whereas summer precipitation variability for future climate is decreasing for the majority of the RCMs (Zanis et al., 2008).These pronounced changes in precipitation and temperature will have subsequent effects on droughts in the region.Loukas et al. (2007b), using the delta downscaling method of Global Circulation Model CGCMa2 (method of truncated means) on precipitation had assessed climate change impacts on drought impulses in the region of Thessaly.They found that future climate change would result in a significant increase in the number, severity and duration of drought events in Thessaly, which is evident even in the period 2020-2050.Drought events would be doubled and in some cases tripled by end of this century when using the socio-economic scenarios IS92a and SRES A2.Furthermore, in another study (Loukas et al., 2008), Annual Weighted Cumulative Drought Severity-Timescale-Frequency curves, which integrate the relationships between drought severity over the year, timescale and frequency and applied for the identification of various types of droughts, indicated that large increase in annual drought severity is expected towards the end of the century for SRES A2 and SRES B2 scenarios.

Methodology
The aim of this study is to evaluate a statistical downscaling method for monthly precipitation and the subsequent estimation of climate change impacts on droughts.The downscaling method was developed using the outputs of the Canadian Centre for Climate Modeling Analysis General Circulation Model (CGCMa2) for the base historical period , validated against observed precipitation for the period 1990-2002, and used to estimate monthly precipitation timeseries for two future periods 2020-2050 and 2070-2100.The droughts have been assessed using the most commonly used drought index, the Standardized Precipitation Index (SPI).The SPI timeseries have been estimated at multiple timescales for the historical base period 1960-1990 for observed and downscaled monthly precipitation, validated for the period 1990-2002 for assessing drought severity classes, and used for evaluating future climate change impacts on droughts.The methodologies used in this study are presented in the next paragraphs.

Global circulation model
Global Circulation Models (GCMs) have been used to study the effects of the increasing concentration of carbon dioxide and the other greenhouse gases on the Earth's climate.These models link atmospheric processes with ocean and land surface processes and can be used to provide projections of the changes in temperature, precipitation and other climate variables in response to changes in greenhouse gas emissions.The second generation of GCMs (Manabe and Stouffer, 1996;Johns et al., 1997;Boer et al., 2000) is transient models assuming an increase of CO 2 equivalent concentration at a rate of 1% per annum from 1990 to 2100.In this study the gridpoint outputs from the second-generation Canadian Centre for Climate Modeling and Analysis GCM (CGCMa2) (Boer et al., 2000;Flato and Boer, 2001) and for two socio-economic development scenarios were used for the assessment of climate change impacts on monthly precipitation in Lake Karla watershed.The CGCMa2 is a spectral model with 10 atmospheric levels and has a resolution equivalent to 3.75 • of latitude by 3.75 • of longitude.The ocean component is based on the Geophysical Fluid Dynamics Laboratory MOM1.1 model and has a resolution of roughly 1.8 • of latitude by 1.8 • of longitude and 29 vertical levels.SRES A2 scenario assumes a strong, but regionally oriented economic growth and fragmented technological change with an emphasis on human wealth.It represents an high emissions scenario.The second scenario is the SRES B2 scenario which emphasizes the protection of the environment and social equity, but also relies on local solutions to economic, social, and environmental sustainability and represents a low emission scenario.These scenarios represent a world in which the differences between developed and developing countries remain strong.The two socio-economic scenarios used have been widely adopted as standard scenarios for use in climate change impact studies (IPCC, 2007).Scenario runs were taken over two time periods: a) 2020-2050 and b) 2070-2100.
The commonly used approach in climate change studies is to combine the output of the GCMs with an impact model.This approach is quite realistic although there are inherent uncertainties about the details of regional climate changes.These uncertainties stem from a number of sources, namely from uncertainties in GCM outputs, downscaling of GCM outputs and specification of the climate change scenarios.The major drawback of the current generation of GCMs is the limitation of their spatial resolution and the resolution of the output.Usually the output of GCMs is given for a much larger scale than the scale of even a large watershed.Interpolation techniques (McCabe and Wolock, 1999), statistical downscaling (Brandsma and Buishand, 1997;Wilby et al., 2002) and downscaling through coupling of GCM output and regional meteorological models (Giorgi et al., 2001) are methods that have been used to overcome the spatial resolution limitation of the GCMs.Uncertainty increases within and between every link of the approach.This uncertainty depends on: 1. quality of GCM simulations, regarding the predictor variables for downscaling (uncertainty of emission scenario included herein); 2. quality of downscaled scenarios, due to inhomogeneities in observed data and shortcomings of the technique applied; 3. quality and resolutions of the impact model(s), which are often strong simplifications of reality; and 4. errors in input data due to instrumentation and/or sample data error.
GCM uncertainty might be assessed by using different GCMs and by using Monte Carlo experiments with one GCM starting with different initial conditions.Uncertainty due to downscaling techniques might be assessed, e.g. by using different downscaling techniques or by varying parameterizations of the downscaling models.Likewise, uncertainties of impact models can be estimated by varying input parameters, taking into account, e.g.sampling errors.In this study, two types of uncertainty are addressed the downscaling technique and the impact model uncertainty.

Statistical downscaling method
Statistical downscaling is the process of building an empirical model: for a small-scale feature y, not adequately described in GCMs, and large-scale features x, well resolved.As predictands, y has been used as weather variables, such as monthly temperatures, and/or as monthly precipitation amounts.The predictor x has often been chosen as characteristics of the weather circulation.If the function F is linear, Eq. ( 1) becomes with ε drawn from a normal distribution with zero mean and standard deviation σ , and 0<a<1.The variations in ε are assumed to be independent from x.In this setting, the randomness in y stems from the randomness in ε.Hence, Eq. ( 1) must be understood as a stochastic equation (Von Storch, 1999).
Several statistical methods are applicable to the description of relationships between large-scale upper-air fields and local climate elements.Multiple linear regression (MLR), either based directly on grid point data or on principal components (PCs) of predictor fields, canonical correlation analysis (CCA), and non-linear methods such as multivariate splines and neural networks have been used most widely (Wilby and Wigley, 1997;Xu, 1999;Xu et al., 2005;Fowler et al., 2007).In this study, the GCM grid point outputs were downscaled using multiple regression equations between GCM predictor output variables and areal monthly precipitation.
Stepwise screening of gridpoint data was found to be the best statistical model among canonical correlation analysis, singular value decomposition, and multiple regression models on principal components (PCs) of predictor fields for downscaling daily temperature in Europe (Huth, 1999).The predictors used in such analyses should be: a) well simulated by the GCM, b) strongly correlated with the predictand variable (precipitation), and c) available.Using these criteria, six predictor grid variables were used, namely the mean sea level pressure (mslp), the mean 2 m wind speed (swa), the precipitation (pcp), the mean surface temperature (st), the 500 hPa geopotential height (gz500), and the geopotential thickness between 500 and 1000 hPa (gz500-1000).These are the most commonly used predictors in statistical downscaling of precipitation (IPCC, 2007).
A procedure based on forward selection stepwise regression technique and included testing with various linear and non-linear regression models was employed (Loukas et al., 2008).All of these models rely on homogeneous long timeseries of the target parameter on the local scale and one or several atmospheric predictors on the large-scale.A major limitation is the assumption that the relationships obtained under present conditions will also hold true under a changing climate.In this study, dummy variables (a set of twelve categorical variables assigned to the 12 months of the year) are used to account for the effect of the "month" on precipitation.The best regression downscaling model containing monthly dummy variables is expressed as (Loukas et al., 2008): , take the value of 0 and so on.However, the monthly downscaled precipitation (P MLR ) values will always have smaller variance than the local values (i.e.areal observed precipitation) (Von Storch, 1999).In many climate impact studies the variance of the downscaled timeseries should be the same with the variance of the observed values.To meet this requirement various methods have been proposed such as variance inflation (Karl et al., 1990;Huth, 1999), expanded downscaling (Burger, 1996), and randomization (Dehn and Duma, 1999).In this study, to preserve the variability of the observed series, the estimated precipitation was combined with the residual values of the regression.These can be viewed as a noise component, statistically independent of the large-scale climate.In the formula: with P = observed monthly precipitation, P MLR = monthly precipitation explained by multiple linear regression and P residual = residuals of MLR.If this operation is carried out on the estimated series of the regression fitting period (October 1960-September 1990), the result is the observed series.For the climate scenarios, P MLR is obtained by downscaling the GCM outputs while P residual remains unchanged.In this way, the problem of limited correlation between predictor and predictand variables may be tackled.However, in order to estimate the uncertainty of the downscaling method stochastic timeseries modelling was applied for the treatment of the residuals.Stochastic simulation of hydrologic timeseries such as precipitation is typically based on mathematical models.For this purpose a number of stochastic models have been suggested in literature (Salas, 1993;Hipel and McLeod, 1994).Using one type of model or another for a particular case at hand depends on several factors such as, physical and statistical characteristics of the process under consideration, data availability, the complexity of the system, and the overall purpose of the simulation study.Given the historical record, one would like the model to reproduce the historical statistics.This is why a standard step in hydrologic simulation studies is to determine the historical statistics.Once a model has been selected, the next step is to estimate the model parameters, then to test whether the model represents reasonably well the process under consideration, and finally to carry out the needed simulation study.  .Univariate stationary ARMA have been applied in the standardized monthly residuals.However, this procedure failed to reproduce the monthly residual correlations and periodic ARMA models have been fitted.The best fitted model, according to Akaike Information Criterion (AIC), was the Periodic Autoregressive Model of order four (4), PAR(4) defined as: where P res v,t respresents the monthly precipitation residual for year v and month(season) t; it is normally distributed with mean zero and variance σ 2 t (P res ); e v,t is the normally distributed and uncorrelated noise which has mean zero and variance σ 2 t (e); and ϕ 1,t . . . .ϕ 4,t are the monthly autoregressive parameters.Several statistics were estimated to evaluate Eq. ( 5) in simulating residual monthly precipitation for development and validation periods and then Eq. ( 5) was applied stochastically to generate 100 timeseries of the P residual .The calculated residual precipitation timeseries were added to the downscaled P MLR to reproduce the observed monthly precipitation pattern that used in drought estimation using a meteorological drought index the SPI.Finally, the developed MLR equation (Eq. 3) was used to downscale monthly GCM precipitation timeseries P MLR for the future periods 2020-2050, and 2070-2100, and then the precipitation residuals (P residual ) were added to P MLR using Eq. ( 4), assuming that the precipitation residual timeseries in the future have the same statistical characteristics of the historical period.Essentially, the residual precipitation timeseries for the future periods were the timeseries generated by Eq. ( 5) for the historical base period.

Standardized precipitation index
Many indices have been used for the identification of more than one type of drought (Tate and Gustard, 2000;Keyantash and Dracup, 2002) and their categorization may not be appropriate, although it is widely used (Wilhite and Glantz, 1985;AMS, 2004).The Standardized Precipitation Index (SPI) has been developed by McKee and his associates (1993) for defining and monitoring droughts.It is used, among others, by the US Colorado Climate Center, the US Western Regional Climate Center, and the US National Drought Mitigation Center to monitor drought in the United States.The main advantage of the SPI is that can be calculated for multiple time-scales.This is very important because the timescale over which precipitation deficits accumulate functionally separates different types of drought (McKee et al., 1995) and, therefore, allows to quantify the natural lags between precipitation and other water usable sources such as river discharge, soil moisture and reservoir storage.Recent studies have used SPI as indicator of hydrological and water resources variables, like soil moisture, surface runoff and reservoir storage (Loukas and Vasiliades, 2005;Vicente-Serrano and Lopez-Moreno, 2005).The US National Drought Mitigation Center computes the SPI with five running time intervals, i.e. 1-, 3-, 6-, 9-, and 12-months, but the index is flexible with respect to the period chosen.This powerful feature can provide an overwhelming amount of information unless researchers have a clear idea of the desired intervals (Loukas and Vasiliades, 2005;Vicente-Serrano and Lopez-Moreno, 2005).
Computation of the SPI involves fitting a Gamma probability density function to a given frequency distribution of precipitation totals for a station, area or a watershed.The alpha and beta parameters of the Gamma probability density function are estimated for each station, for each timescale of interest (1, 3, 6, 9, 12 months, etc.), and for each month of the year.The Gamma distribution is defined by its probability density function: where α, β are the shape and scale parameters respectively, P is the precipitation amount and (α) is the gamma function.
Maximum likelihood solutions are used to optimally estimate α and β: α , and n is the number of observations.The resulting parameters are then used to find the cumulative probability of an observed precipitation event for the given month and timescale for the station in question.Since the gamma function is undefined for P =0 and a precipitation distribution may contain zeros, the cumulative probability becomes: where q is the probability of a zero and G(P ) the cumulative probability of the incomplete gamma function.If m is the number of zeros in a precipitation timeseries, then q can be estimated by m/n.The cumulative probability, H (P ), is then transformed to the standard normal random variable z with mean equal to zero and variance of one, which is the value of the SPI.Once standardized the strength of the anomaly is classified as set out in Table 1.This table also contains the corresponding probabilities of occurrence of each severity arising naturally from the Normal probability density function.Thus, at a given location for an individual month, moderate dry periods (SPI≤−1) have an occurrence probability of 15.9%, whereas extreme dry periods (SPI≤−2) have an event probability of 2.3%.Extreme values in the SPI will, by definition, occur with the same frequency at all locations.
Negative SPI values indicate droughts and positive SPI values denote wet weather conditions (Table 1).
In this study, areal monthly precipitation accumulations were used for the estimation of the monthly SPI for 1month, 3-month, 6-month, 9-month, 12-month, and 24month timescales, in the development period at Lake Karla watershed.The estimation error of the parameters calculation was tested using the 100 generated precipitation timeseries.Five evaluation statistics were used to assess the comparison of observed SPI timeseries and the use of median/or average parameter values of the generated 100 SPI timeseries in simulation the observed precipitation pattern.These statistics are the mean absolute error (MAE), the root mean square error (RMSE), the coefficient of efficiency (CE), the index of agreement (IA) and the persistence index (PI).MAE and RMSE are statistical parameters of residual error between observed and modelled timeseries datasets.MAE and RMSE account in real units the level of overall agreement between the observed and modelled datasets.They are non-negative metrics that have no upper bound and for a perfect model the result would be zero.MAE is unbiased where RMSE consists of a weighted measure of the error in which the largest deviations between the observed and modelled values contribute the most and subsequently is more sensitive than MAE.CE, IA, and PI are dimensionless coefficients that contrast model performance with accepted norms or standards (Dawson et al., 2007).The mathematical formulations of these statistics are: where, SPI i and SPI i−1 is the observed SPI value on month i and i−1, respectively, ∧ SPI i the simulated SPI value on month i, and SPI the average value of observed SPI for the simulation period (i=1 to n datapoints).CE and PI range from − ∝ to one, whereas, IA ranges from 0.0 (poor model) to 1.0 (perfect model).These three statistics record as a ratio the level of overall agreement between the observed and modelled datasets and represent improvements over the coefficient of determination (R 2 ) for model evaluation and forecasting purposes since they are sensitive to differences in the observed and modelled means and variances (Dawson et al., 2007).
The observed and generated precipitation timeseries were also used in the estimation of the SPI for present and future periods and assessed against the observed and simulated SPI timeseries in calculating the respective SPI classes.Because SPI is as a standardized drought index is designed to express drought conditions with respect to normal conditions at a given site.The result is that the range of the SPI is about the same for every meteorological station or watershed and/or for the study period represented by the input precipitation timeseries.Hence, this index cannot be used for comparison between-stations for the identification of drought magnitude, nor between different time periods, which, are essential in evaluating the potential effects of climate change.The parameters of the gamma distribution, α and β, are assumed unchanged in the future and their respective values for the historical period have been used.Several studies that assessing climate change impacts on drought indices have adopted the same technique (Loukas et al., 2007b;Loukas et al., 2008;Dubrovsky et al., 2009).

Results and discussion
The methods described above were used, firstly, to statically downscale the monthly areal precipitation using information from the second-generation Canadian Centre for Climate Modeling and Analysis GCM (CGCMa2) in Lake Karla watershed.The statistical downscaling method was developed for historical period , and validated against observed precipitation for the period 1990-2002.Furthermore, the SPI timeseries were estimated using the observed and the stochastically generated precipitation and compared for the development and validation periods.The historical period from October 1960 to September 1990 was considered as the base period for further analysis.Secondly, the future climate areal monthly precipitation timeseries were estimated from the outputs of CGCMa2 for two socio-economic scenarios using the statistical downscaling method.This procedure was repeated for the two future periods a) October 2020-September 2050, and b) October 2070 to September 2100.Finally, future climate timeseries of SPI for various timescales were calculated and the effect of climate change on droughts was assessed from the changes in the number of negative monthly SPI values by severity classes.

Historical period
The use of statistical downscaling targeted in generating monthly timeseries of precipitation for evaluating climate change impacts on meteorological droughts.The downscaling method (Eq.4) is a combination of multiple linear regression (Eq. 3) and stochastic modeling of the residuals derived from MLR (Eq.5).The analysis results of the MLR has shown that the correlation coefficient, r, between the logarithmically transformed estimated downscaled monthly areal precipitation and the logarithmically transformed observed monthly basin-wide precipitation was equal to 0.66 for the development base period, 1960-1990, and 0.65 for the validation period, 1990-2002.The developed relationship has been found to be statistically significant at α=5% significance level using the t-test.These results are comparable with the results of previous studies on statistical monthly precipitation downscaling with more sophisticated methodologies (Dehn and Buma, 1999;Schoof and Pryor, 2001;Buishand et al., 2004;Tatli et al., 2004) and are better than the results obtained by Loukas et al. (2008) with the same MLR statistical downscaling method (Eq. 3) for Lake Karla watershed (r=0.55 for development period, 1960-1990, and 0.57 for validation period, 1990-1993).The latter is explained by the use, in this study, of a more dense precipitation network to estimate areal monthly precipitation for Lake Karla watershed (Fig. 1).However, this regression model (Eq.3) failed to reproduce the variance of precipitation, although simulated quite well the mean monthly precipitation for historical (Fig. 2) and validation (Fig. 3) periods.The same results were also valid in application of Eq. ( 3) to twelve hydrological homogeneous areas in Thessaly, Greece (Loukas et al., 2008).Stochastic timeseries theory (Eq.5) was applied for the treatment of the residuals to preserve the variance of observed monthly precipitation.The Least Squares (LS) method have been used to estimate the model parameters of the PAR(4) model by minimizing the sum of squares of the residuals with the moment estimates of model parameters taken as the initial values in the search algorithm.The skewness test on a month-by-month basis was applied for testing the normality and independence, respectively, of the e v,t (Salas et al., 1980).The skewness test of normality at 10% significance level was successful for all months expect for September where the hypothesis of normality was rejected.Application of the Portmanteau lack of fit test, based on the autocorrelation of the entire residual series, for testing the independence of the e v,t was successful at 5% significance level.After the successful fitting the PAR(4) model, 100 generated timeseries were produced to evaluate the results of the process, and the uncertainty introduced for historical period.Furthermore, 100 residual timeseries were generated to validate the methodology for the period 1990-2002.Application of the fitted PAR(4) model to monthly precipitation residuals was able to reproduce the historical statis- tical properties (mean, standard deviation, skewness coefficient, month-to-month correlations and autocorrelations) of the residual precipitation (results are not shown due to paper length limitations).The generated 100 P residual were added to downscaled P MLR to reproduce the monthly precipitation pattern for Lake Karla watershed (Eq.4). Figure 2 shows the statistical properties (mean and standard deviation) of the generated sample for the development period  whereas Fig. 3 for the validation period (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002).These figures show that the method is able to reproduce the statistical properties of historical monthly precipitation for historical period.
The SPI was then calculated for the development period for observed and generated monthly precipitation timeseries at multiple timescales (1, 3, 6, 9, 12, 24 months) and compared at validation period.The authors are aware that inconsistent conclusions could be obtained if smaller time lengths of precipitation record are involved in the SPI calculation (Wu et al., 2005).The longer the length of record used in the SPI calculation, the more reliable the SPI values will be, especially for long-time-scale SPI values (Guttman, 1994).to instability of the parameter estimates.If the parameter estimates have little confidence, then the resulting SPI values will also have little confidence.Therefore, operational use of SPI probably needs the longer length of record, because the shorter one is likely not to capture the 'signals' of climate variability (Wu et al., 2005).However, operational use is not the aim of this study but the testing of a statistical downscaling method for assessing climate change impacts on droughts.Model parameter convergence has been guaranteed comparing the historical SPI timeseries with the calculated SPI using the whole sample .This analysis had shown slightly differences in the timeseries when the whole sample is used in the estimation of SPI for all timescales.
Model parameters of the gamma distribution were estimated based on historical period.The results will be presented for SPI timescales of 3-, 9-, 24-months since there are representative of agricultural, hydrological and water resources drought, respectively (Loukas and Vasiliades, 2005).The α and β parameters of the gamma distribution (Eq.6) were assessed for each study timeseries individually.Model parameter convergence was tested using median and/or average parameter values of the generated 100 SPI timeseries in simulating the observed precipitation pattern.The use of average values for Gamma parameters of the 100 generated precipitation timeseries failed to reproduce the observed drought pattern for historical and validation periods (Fig. 4) due to the skewness of the accumulated precipitation which increases as timescale increases.However, the use of median values in α and β parameters of the gamma distribution is able to reproduce the historical temporal evolution of SPI timeseries for all timescales (Fig. 4.).Table 2 presents the evaluation statistics of the use of median parameters values in simulating the historical observed precipitation timeseries for development  and validation periods (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002).All evaluation statistics are considered satisfactorily www.nat-hazards-earth-syst-sci.net/9/879/2009/Nat.Hazards Earth Syst.Sci., 9, 879-894, 2009 for development and validation periods, respectively.However, the error between observed and median values simulated SPI slightly increases as timescale increases (Table 2).The temporal evolution of observed 3-month SPI timeseries (Fig. 4a) showed that Lake Karla watershed experienced frequent moderate and severe droughts (i.e.SPI≤−1) for all months of the year.As timescale increases, the monthly SPI timeseries were smoothed and drought duration of the identified drought events could easily be assessed (Fig. 4b and c).The temporal variation of observed and generated SPI timeseries was investigated, firstly, by comparing the number of months for which the SPI values for all timescales indicated moderate drought (−1.50<SPI≤−1.00),severe drought (−2.00<SPI≤−1.50), and extreme drought (SPI≤−2.00).Table 3 presents the number of dry months at various drought severity classes for development and validation periods.In general the statistical downscaling method is capable to simulate drought patterns since the generated timeseries produce similar number of total dry months with the observed timeseries for development period (average and median of Table 3).However, the allocation of the dry months in the respective drought classes is quite different.The downscaling method simulates satisfactorily the severe droughts but produces smaller moderate dry months and larger number of dry months for extreme drought.For validation period, the downscaling method gives larger number of dry months (Table 3).Moderate dry events are simulated with reasonable accuracy whereas overestimation of severe and extreme dry months is observed.This is attributed to monthly precipitation distribution of the validation period and especially the September precipitation variance.However, the observed drought pattern is in the range of the stochastic simulation.Furthermore, the superiority of the median statistical characteristic in depicting dry months is evident.Table 3 also indicates that about 15.4% of the time Lake Karla watershed experienced drought for the historical period for SPI 3-month.Stochastic simulation results for the same timescale show that Lake Karla experienced droughts, on average, for 15.8% of time with a range of 13.1% to 18.7%.Similar results are observed and for the other timescales.On the other hand, for validation period, stochastic simulation results show a small increase in the drought time and subsequently drought duration.For example, for 3-month timescale, the stochastic results show that 19.10%, with a range from 9.7% to 35.4%, of time is on drought conditions.Using the observed precipitation timeseries this percentage is 13.9%.This difference of 5% in the time percentage of dry months is maintained in all timescales and it is always within the range of stochastic simulation results.

Climate change impacts on droughts
The statistical downscaling method (Eq.4) was applied to generate future monthly precipitation timeseries for the future periods 2020-2050, and 2070-2100.Using Eq. ( 3) of multiple linear regression, future monthly precipitation timeseries were produced and added to the same residual precipitation timeseries (Eq.5) of the development period, assuming that the timeseries of the residuals will remain un- changed.For the climate scenarios SRES A2 and SRES B2, P MLR is obtained by downscaling the GCM outputs while P residual , generated by PAR(4) timeseries model (Eq.5), assumed to remain unchanged in the future.In this way, the problem of limited correlation between predictor and predictand variables has been tackled.This procedure applied for the stochastically generated 100 residual timeseries as well as for the observed rainfall where the residuals remained unchanged, mentioned as deterministic approach in Tables 4  and 5.In this way a direct comparison is feasible in order to evaluate the uncertainty introduced in the downscaling method.Application of the downscaling precipitation method for the two socio-economic scenarios, namely SRES A2 and SRES B2, indicated that the scenario SRES A2 was the most severe.The larger reduction for average annual precipitation was observed for the period 2070-2100 as expected.Monthly precipitation was, in general, reduced for all months and time periods except for September when average monthly precipitation was increased.The average annual precipitation reduces about 3.9% and 3.4% for period 2020-2050, and about 13.5% and 8.5% for period 2070-2100 for SRES A2 and SRES B2, respectively.Subsequent decrease is observed for monthly precipitation in all months except for September.The average September precipitation increased by 12.7% and 13.2% for SRES A2 and SRES B2, respectively, in the period 2020-2050.Smaller increase was observed in the period 2070-2100, 7.5% for SRES A2 and 10.9% for SRES B2.
The future truncated precipitation timeseries and the 100 future generated monthly precipitation timeseries for the two socio-economic scenarios and the two future periods were used for the estimation of SPI timeseries at multiple timescales.The parameters of the gamma distribution, α and β, were assumed unchanged in the future and their respective values for the historical period have been used.Figure 5 shows the comparison of the 24-month SPI, where the larger differences are identified, between present and future climatic scenarios and the two climate study periods for the observed historical precipitation.Similar pattern but with smaller deviations are identified for shorter timescales.The same principle is obvious and in the generated SPI timeseries.Figure 5 indicates that the drought severity and duration would be increased in the future periods.Overall, the number of dry months (SPI≤−1.00)will increase for the two scenarios.Table 4 shows the number of dry months for the two scenarios in the period 2020-2050.The dry months are also categorized with the severity of drought classes.These results indicated that the number of dry months will be increased for the two study scenarios, with SRES A2 being the most severe and SRES B2 the most conservative scenario.According to this analysis, the total number of dry months will be increased by 45% for the deterministic case (observed future) for SPI 3-month for both socio-economic scenarios (Table 4).For the stochastically generated timeseries the increase is from 0% to 71% with a median increase of 38% and 36% for SRES A2 and SRES B2, respectively.The largest increase (about 800%), compared to the base historical period, is observed for the extreme dry months (SPI≤−2.00)for the two scenarios.Similar results have been found and for the other timescales of SPI (Table 4).The increase in the number of dry months is rising at larger timescales.For example, in 24-month timescale, the total number of dry months in the stochastic approach will be increased by about 88% with a range from 33% to 153% for scenario SRES A2.In the scenario SRES B2 this increase is from 23% to 146% with an average of 84%.Overall, when comparing the deterministic and the stochastic approach, the stochastic approach gives more conservative results than the deterministic (Table 4).
Table 4 also indicates that about 22.4% of the time Lake Karla watershed will be experienced droughts using the deterministic approach for the 3-month timescale and SRES A2 and SRES B2 scenarios.Stochastic simulation results for the same timescale show that Lake Karla experienced droughts 21.2% of time, on average, with a range of 15.6% to 26.3% for SRES A2 scenario, and 20.9% of time, on average, with range of 15.4% to 25.1% for SRES B2.Similar results are observed and for the other timescales.Again the stochastic approach gives more conservative results than the deterministic approach.It should be mentioned that the increase in time percentage of dry months is higher for larger timescales than the smaller timescales.For example for 24-month timescale the increase in time percentage of dry months is on average 14.8% (range: 5.6%-25.8%)and 14.2% (range 3.9%-24.6%)for SRES A2 and SRES B2, respectively, when compared with the historical period 1960-1990.This result is very important since larger timescales are used to monitor hydrological and water resources droughts and smaller timescales to monitor meteorological and agricultural droughts.Hence, water resources seem to be more vulnerable to climate change.
The increase in the number of dry months is rising for the period 2070-2100 again at larger SPI timescales.For SPI-24 , the total number of dry months in the stochastic approach will be increased by 284% with a range from 137% to 377% for scenario SRES A2.For the scenario SRES B2 this increase is from 95% to 270% with a median increase of 182% (Table 5).Overall, when comparing the deterministic and the stochastic approach, the stochastic approach gives more conservative results than the deterministic (Table 5).
Table 5 also indicates that 37.2% and 31.1% of the time Lake Karla watershed will be experienced droughts using the deterministic approach for the SPI-3 and SRES A2 and SRES B2 scenarios, respectively.Stochastic simulation results for the same timescale show that Lake Karla watershed will experience droughts for SRES A2 and SRES B2 scenarios for 32.1% (25.7% to 37.2%) and 26.2% (19.8% to 31.6%) of time, respectively.The distribution of the time percentage in the dry months is uniform distributed in the drought classes.That means that moderate, severe and extreme droughts will be expected with the same frequency.Similar results are observed for the other SPI timescales.The stochastic approach gives more conservative results than the deterministic approach as in the period 2020-2050.It should be mentioned that the increase in time percentage of dry months is higher for larger timescales than the smaller timescales of SPI.For example, for SPI-24, the increase in time percentage of dry months, when compared with the historical period 1960-1990, is on median average 48% (range: 23%-64%) and 31% (range 16%-46%) for SRES A2 and SRES B2 scenarios, respectively.That means that 65% and 48% of the 360 months in period 2070-2100 will be on drought conditions for SRES A2 and SRES B2, respectively.
The above results are similar with the results of a recent study (Loukas et al., 2007b) in which climate change impacts on drought impulses in the region of Thessaly had been assessed.That study (Loukas et al., 2007b) had used the delta method of downscaling (method of truncated means) for evaluating climate change impacts on droughts.In the simplest "delta" or "perturbation" downscaling method, where differences between the control and future GCM simulations are applied to baseline observations by simply adding or scaling the mean climatic change factors, assumes a constant bias through time and no consideration is made for changes in the variability of descriptors with climate change (Evans and Schreider, 2002;Diaz-Nieto and Wilby, 2005).The present study extends the statistical downscaling method of Loukas et al. (2008) with a stochastic component to account the variability of descriptors with climate change, and to evaluate the uncertainty introduced on climate change impact on droughts.Results show that climate change will have a major impact on various types of droughts but with high uncertainty.

Conclusions
This study illustrated that Lake Karla watershed experienced frequent moderate and severe droughts during the period 1960-1990 and future climate change would result in a significant increase in drought severity.The outputs of CGCMa2 model have been employed to statistically downscale monthly precipitation, to account the uncertainty of the downscaling method and to estimate future precipitation timeseries for the periods of 2020-2050 and 2070-2100.The choice of a GCM is not a key component since all available CCM underestimate annual precipitation by about 20% for the Mediterranean basin (IPCC, 2007).The developed statistical downscaling procedure could be applied in any other available GCM or in multimodel GCMs, to simulate GCM uncertainty.In this study the uncertainty of a GCM is included in the stochastic simulation of the residuals.A lot of studies have proved that with a statistical or dynamical downscaling method, it could be explained a large proportion of the monthly observed precipitation for present climatic conditions (i.e.Dehn and Buma, 1999;Schoof and Pryor, 2001;Buishand et al., 2004;Tatli et al., 2004, Loukas et al., 2008).Future work is to test the developed methodology with the results of Regional Climate Models (PRUDENCE, ENSEM-BLES EU), which, for Central Eastern Greece underestimate annual precipitation by about 33±19% (Zanis et al., 2008).The monthly and annual precipitation future projections of the present study lie within the range of respective calculations from various GCMs (IPCC, 2007) and are comparable with the results for the period 2070-2100 and SRES A2 scenario of the nine RCMs used for simulating future annual precipitation in Central Eastern Greece, (−13.5%, this study; −15.8%,Zanis et al., 2008).
The historical, generated and downscaled future period precipitation timeseries were used for the estimation of SPI for various timescales and for two socio-economic scenarios SRES A2 and SRES B2.SPI calculation parameters were based according to the present climate and used to calculate SPI timeseries for the future climate.The authors believe that the scientific community must start a debate on whether the parameters of meteorological drought indices should be only limited to present climate.If the full length (historical and future climate precipitation timeseries) are used in the calculation procedure of SPI then the effects of climate change on droughts would have smaller effects than the results presented in the manuscript.However, because of www.nat-hazards-earth-syst-sci.net/9/879/2009/Nat.Hazards Earth Syst.Sci., 9, 879-894, 2009 the drier future climate, present drought episodes could not be easily identified, and thus SPI is loosing its operational character in predicting and quantifying drought episodes.
The results of the present study indicate the accuracy, reliability and uncertainty of the statistical downscaling method for present and future climate conditions and the suitability of the downscaling method for the assessment of climate change on hydrological, agricultural and water resources droughts.Results show that climate change will have a major impact on droughts but the uncertainty introduced is quite large.Larger timescales of SPI which are used to monitor hydrological and water resources droughts are more sensitive to climate change than smaller timescales which are used to monitor meteorological and agricultural droughts.These results indicated that climate change would largely affect drought severity and subsequently the design of future water resources projects (e.g.reservoirs).However, the uncertainty in the future drought episodes should always be accounted for since the range of the drought episodes is quite large especially for larger timescales.Future drought episodes due to climate change should be handled with caution and always with their respective ranges as this analysis demonstrates.

Figure 1 .
Figure 1.Study area, database and digital elevation model of Lake Karla watershed.Fig. 1. Study area, database and digital elevation model of Lake Karla watershed.
3) where P MLR is the logarithmically transformed monthly precipitation, b 1 , b 2 , b 3 , . . ., b 12 are the monthly weighing dummy variables, a 1 , a 2 , a 3 , . . ., a 12 are regression coefficients, and c is the regression constant.Dummy variables, b 1 −b 12 , are assigned binary values, 0 or 1, depending on the month in which precipitation is referred.For example, if month is October, then, b 1 takes the value of 1 and all the other dummy variables, b 2 −b 12 , take the value of 0. Similarly, if month is November, then, b 1 takes the value of 0, b 2 takes the value of 1 and all the other dummy variables, b 3 −b 12

Figure 2 .Fig. 2 .
Figure 2. Statistical properties of the statistical downscaling procedure for development period 1960-1990 for a) average monthly precipitation and b) standard deviation of monthly precipitation (Note: Box-Whisker plots and average refers to stochastic simulation results, whereas MLR are the results obtained from multiple linear regression)

Figure 3 .Fig. 3 .
Figure 3. Statistical properties of the statistical downscaling procedure for validation period 1990-2002 for a) average monthly precipitation and b) standard deviation of monthly precipitation (Note: Box-Whisker plots and average refers to stochastic simulation results, whereas MLR are the results obtained from multiple linear regression)

Figure 4 .Fig. 4 .
Figure 4. Comparison of observed SPI timeseries with simulated SPI timeseries using median and average values of the Gamma distribution parameters, α and β, of stochastically generated precipitation timeseries in reproducing the observed historical precipitation for timescales: a) 3-month, b) 9-month and c) 24-month.

Table 1 .
Drought classification by SPI values and corresponding event probabilities.

Table 2 .
Evaluation statistics between observed SPI timeseries and simulated SPI timeseries using median values of the Gamma distribution parameters, α and β, of stochastically generated precipitation timeseries in reproducing the observed historical precipitation for development(1960-1990) and validation periods (1990-2002).

Table 3 .
Total numbers of dry months (SPI≤−1) for various drought severity classes at development and validation periods

Table 4 .
Total numbers of dry months (SPI≤−1) for various drought severity classes at period 2020-2050 for the two socioeconomic scenarios.

Table 5 .
Total numbers of dry months (SPI≤−1) for various drought severity classes at period 2070-2100 for the two socioeconomic scenarios.