Benefits from using combined dynamical-statistical downscaling approaches – lessons from a case study in the Mediterranean region

Various downscaling techniques have been developed to bridge the scale gap between global climate models (GCMs) and finer scales required to assess hydrological impacts of climate change. Such techniques may be grouped into two downscaling approaches: the deterministic dynamical downscaling (DD) and the statistical downscaling (SD). Although SD has been traditionally seen as an alternative to DD, recent works on statistical downscaling have aimed to combine the benefits of these two approaches. The overall objective of this study is to assess whether a DD processing performed before the SD permits to obtain more suitable climate scenarios for basin scale hydrological applications starting from GCM simulations. The case study presented here focuses on the Apulia region (South East of Italy, surface area about 20 000 km 2), characterised by a typical Mediterranean climate; the monthly cumulated precipitation and monthly mean of daily minimum and maximum temperature distribution were examined for the period 1953– 2000. The fifth-generation ECHAM model from the MaxPlanck-Institute for Meteorology was adopted as GCM. The DD was carried out with the Protheus system (ENEA), while the SD was performed through a monthly quantile-quantile correction. The SD resulted efficient in reducing the mean bias in the spatial distribution at both annual and seasonal scales, but it was not able to correct the miss-modelled nonstationary components of the GCM dynamics. The DD provided a partial correction by enhancing the spatial heterogeneity of trends and the long-term time evolution predicted by the GCM. The best results were obtained through the combination of both DD and SD approaches.


Introduction
Global climate models (GCMs) are the primary tool for understanding how global climate may change in the future.However, they currently do not provide reliable information on scales below about 200 km (Meehl et al., 2007).Hydrological processes typically occur at finer scales (Parry et al., 2007).Consequently, basin-scale assessments of climate change impacts usually produce large biases in the simulated hydrological processes whenever the raw output variables from a GCM are adopted (Mearns et al., 2003;Dibike and Coulibaly, 2007).Hence, to reliably assess hydrological impacts of climate change, higher resolution scenarios are required.
Various downscaling techniques have been developed to bridge this scale gap, and a number of papers have previously reviewed the downscaling concept (e.g., Hewitson and Crane, 1996;Wilby and Wigley, 1997;Xu, 1999;Fowler et al., 2007;Maraun et al., 2010).Two approaches to downscaling exist.Dynamical Downscaling (DD) nests a regional climate model (RCM) into a GCM to represent the atmospheric physics with a higher grid box resolution within a limited area of interest.Statistical Downscaling (SD) establishes statistical links between larger and local observed scale weather (Frìas et al., 2006).Traditionally, SD has been seen as an alternative to DD.With the increasing reliability and availability of RCM scenarios, recent works on statistical downscaling have aimed to combine the benefits of these two approaches (e.g., Wilby et al., 2004;Yoon et al., 2012;Vrac et al., 2012).
Published by Copernicus Publications on behalf of the European Geosciences Union.
Most of these studies have compared the performance of the two downscaling methods, but the use of different spatial domains, predictor variables and assessment criteria makes direct comparison of the relative performance difficult to achieve (Fowler et al., 2007).Moreover, studies that investigate more than one variable are rare (Dibike and Coulaby, 2005;Diaz-Nieto and Wilby, 2005;Khan et al., 2006) and few studies compare relative performances of DD and SD (Kidson and Thompson, 1998;Murphy, 1999;Hellström et al., 2001;Wilby et al., 2000;Haylock et al., 2006) or the performance of SD directly applied to GCM relative to the use of an intermediate DD (Hellström and Chen, 2003;Wood et al., 2004;Díez et al., 2005).According to these studies, the advantages of using a combined approach DD-SD are limited, although Wood et al. (2004) stress that an SD is mandatory in order to obtain downscaled scenarios useful for hydrologic simulation.More recently Vrac et al. (2012) highlight that each method brings independent benefits.
The main goal of the present study is to further investigate the advantages and limitations of using a DD as a preprocessing before SD to obtain monthly scenarios of temperature and precipitation.The adopted case-study concerns the Apulia region (SE of Italy) chosen for the availability of well-distributed long-term (collected from the middle of the past century) temperature and precipitation time series.The fifth-generation ECHAM model (Roeckner et al., 2003) and the Protheus system (Artale et al., 2010) have been selected as a state of the art GCM and DD, respectively.The widely used quantile mapping technique (Déqué, 2007) has been applied as SD (Hayhoe, 2010).
The performances of the complete GCM-DD-SD processing chain have been evaluated versus the performances of GCM-DD and GCM-SD alone using mainly correlation coefficients, distance measures such as root-mean-squarederror (RMSE), or explained variance (Fowler et al., 2007), although Busuioc et al. (2001) suggested that for climate change applications the more suitable downscaling model needs to be able to reproduce the low frequency variability.In this respect, the degree of non-stationarity between predictand and predictor has been considered by Hewitson and Crane (2006), while Benestad et al. (2007) and Fan et al. (2011) highlighted the difficulties in capturing longterm trends through downscaling when the GCM fails in this attempt.
In this context, we propose a methodology to evaluate the relative performance of the selected GCM, DD, SD and their combinations not only in terms of bias, but also in terms of time-variability, considering both the trend analysis and the non-stationarity.(3) SD applied directly to the GCM; (4) SD applied to the DD of the GCM; (ref) land observations.The spatial scale associated with each model is reported on the left.

Processing methods
In order to evaluate the relative performances of the DD and SD downscaling methods, the following four methods of data processing were compared with land observations: (1) direct output from the GCM control scenario (GCM); (2) DD applied to the GCM scenario (GCM-DD); (3) SD applied directly to the GCM scenario (GCM-SD); (4) SD applied to the DD of the GCM scenario (GCM-DD-SD).A spatial homogenisation through a Statistical Interpolation (SI) was performed before each comparison, as described below.Thus, data processing (1) to (4) refer to the SI performed on each processing output.Analogously, (ref) refers to the SI performed on the observations dataset.Data fluxes are schematised in Fig. 1.

Global circulation model
The global model simulations considered for this study are those produced by the ECHAM5/MPI-OM and included in the CMIP3 database (Roeckner et al., 2003;Marsland et al., 2002).In particular, the atmospheric component (ECHAM5) is run at spectral resolution T63, corresponding to approximately 200 km at mid-latitudes with 32 vertical levels.Many important topographic features are missing in the global model.For example, Dell 'Aquila et al. (2012) show that the Mediterranean area land-sea mask is an extremely approximated one and that the shape of the Italian peninsula cannot be conveniently captured at the adopted resolution.For ECHAM5 the maximum elevation over the Alps is 890 m at 13.370 • E, 47.560 • N (Eastern Alps), the maximum elevation over the Apennines is 155 m at 13.120 • E, 40.100 • N (Central Italy).
From the global model, we considered daily time series for rainfall, minimum and maximum temperature from the control simulation for the period 1950-2000.The data were then cumulated on a monthly basis for precipitation, whereas daily minimum and maximum temperature were averaged over the same time scale.

Dynamical downscaling
The dynamical downscaling was performed with the PROTHEUS system, an atmosphere-ocean regional climate model composed of the RegCM3 atmospheric regional model and the MITgcm ocean model.A detailed description of the coupled system is provided by Artale et al. (2010).The relevant aspects for the present study are that the atmospheric component RegCM3 is a 3-D, σ -coordinate, primitive equation, hydrostatic model.The dynamical downscaling was performed over an area ranging from 20 • N to 60 • N over the entire Mediterranean basin and produced adopting a uniform horizontal grid of 30 km horizontal resolution and 18σ -levels.For PROTHEUS the maximum elevation over the Alps is 2107 m at 10.516 • E, 46.773 • N (Eastern Alps), the maximum elevation over the Apennines is 930 m at 13.495 • E, 42.100 • N (Central Italy).
The dynamical downscaling of the ECHAM5/MPI-OM global scenarios considered in the present study was previously evaluated in Dell' Aquila et al. (2012), who emphasised the ability of the model in simulating critical aspects of the local climate.The fundamental improvements obtained with this modelling strategy are a partial reduction of the sea surface temperature bias produced in the driving global simulation and a better representation of the corresponding patterns.The dynamical downscaling tends to amplify the fluctuations of the sea-surface temperature seasonal cycle already present in the global driver, and to increase the frequency of large temperature anomalies (both warm and cold events).In particular, a more accurate description of complex orography surrounding the Mediterranean Sea, as well as of land surface processes, produces more organised patterns in the tendency of key impact indicators, such as the aridity index.Instead, the global driver produces extremely noisy results, difficult to interpret in the context of impact studies.
From the PROTHEUS simulations, we considered the cumulated rainfall and instantaneous temperature every six hours for the time period 1950-2000.The data were then cumulated on a monthly basis for precipitation, whereas daily minimum and maximum temperatures were averaged over the same time scale.

Statistical downscaling
According to Fowler et al. (2007), despite the multiplication of more sophisticated SD methods (as weather typing schemes or weather generators), simple statistical downscaling methods (regression models) seem to show similar performances in reproducing the mean climatological features, when compared with the more complex ones.The SD adopted in this study, is the quantile-quantile mapping correction technique method (Déqué, 2007); it has been selected as it is easy to implement, has a low computational cost and is, therefore, widely used for impact studies independently from the variables of interest.Over the monthly time series the quantile-quantile mapping performed at each station results in a global transfer function able to introduce downscaled patterns at the density scale of the land control network.Furthermore, the adopted SD by intrinsically assimilating the observations is able to reproduce the climatology with a small residual bias that can be interpreted as the intrinsic limit of the quantile mapping in projecting a modelled variable onto the distribution of the reference dataset.
The monthly dataset derived from GCM simulations and DD results were statistically downscaled versus the land stations.Each station was compared with the nearest model node.Quantiles were computed both for observations (predictor) and associated simulations (predictand) using a common uniform plotting position (Weibull, 1939;Makkonen, 2008).All values ranging between predictand quantiles corresponding to the plotting position p and p + 1 were then replaced by the predictor quantiles at p + 1. Predictand values lower and higher than the minimum and maximum observed quantiles were replaced by the minimum and maximum observed quantiles, respectively.In order to maximise the data samples' availability for the monthly quantile mapping, as well as the period of analysis, the same period was used for calibration and validation.An Anova test performed on the results from a 4-fold cross validation (Kohavi, 1995) and from the full calibration period confirms the null hypothesis that the quantile mapping performance is not significantly affected by the overlapping of calibration and validation periods.

Statistical interpolation
An ordinary kriging (Cressie, 1988), based on the covariance of the land observation data, was applied as SI (10 km grid) to the four data processing and to the land observations.The use of SI was motivated by the need to compare datasets having different spatial resolutions, including land observations.Monthly spatial covariances were estimated through the experimental semi-variogram of each month of the year.A cross-validation was performed through a leave-one-out cross-validation (LOOCV) (Hawkins et al., 2003) over the interpolated land observation data, in order to estimate the uncertainty introduced by the SI.The LOOCV was carried N. Guyennon et al.: Benefits from using combined dynamical-statistical downscaling approaches out using a single observation from the original sample as validation datum, and the remaining observations as training data.This procedure was repeated for each observation point.At each run, mean over time of the residues, defined as the difference between the SI estimation and the removed validation data, were computed.The obtained spatial mean values of the cross-validation residues were 0.187 mm, 0.019 and 0.023 • C for monthly precipitation, minimum and maximum temperature, respectively.

Indicators of performance
The following indicators of performance were used to evaluate the ability of each data processing in reproducing the land observed temperature and precipitation patterns.Hydrologists generally need the simulated data maintain features of the observed ones in terms of statistical moments.In view of the evaluation of climate change impacts, we consider importance key issues to evaluate the performances of the applied downscaling techniques in relation to the following objectives: (1) reducing the mean bias, (2) reproducing the observed non-stationarity and (3) reproducing the observed spatial heterogeneity of trends.
The measures adopted to evaluate the predictive performance are reported in the following sub-sections from Sects.2.2.1 to 2.2.3.In these sections, the monthly time series resulting from the SI of each data processing p and the land observation (ref), at each node n are referred to as SI p n and SI ref n , respectively.

Mean bias analysis
The mean bias is defined in each node n as: where the overbar stands for the mean over time, at a monthly scale.The spatial variability of mean bias values obtained from the four data processing methods was compared through the 5th, 25th, 50th, 75th and 95th percentile of the cumulative probability distribution over the 10 2 km 2 grid nodes resulting from the SI.The same elaboration was carried out after splitting residues into four seasonal sub dataset: winter (December, January, February), spring (March, April, May), summer (June, July, August) and autumn (September, October, November).

Non-stationarity analysis
The GCM adopted in this study cannot be considered a forecast product because it is not initialised with the observed state of the climate system at any given time.Therefore, an investigation of the time correlation between the model output and the observations reference is not viable.Instead, we require that the selected data processing is able to provide a sufficient description of the statistics of local climate, including the possible non-stationarity in the probability distribution of climate variables.We chose to analyse the nonstationarity in the distribution of climate variables by considering the evolution of quantiles of the corresponding probability distribution.The use of quantiles avoids assumptions on the shape of the probability distributions resulting from each processing method, thereby providing a more accurate detection of any possible change in the probability distribution of the variables of interest.Quantiles of each data processing were computed adopting the uniform plotting position suggested by Weibull (1939), recently confirmed by Makkonen (2008).The quantiles were computed for each data processing p, season s at each node n using a sliding time window centred on the year y and referred as Q The overall variability in the quantiles residues can be expressed through the mean-squared-error (MSE) which is intended here as a measure of the distribution variability for moving time windows: where L is the total number of plotting points.
The MSE can be disaggregated into the sum of the squared mean bias and the variance: (3) The seasonal time variation of the Mean of Quantiles was then computed for a given moving window centred at time y as: Here, the mean is referred to as the average performed over the L plotting points.The seasonal time variation of the quantiles variance was then computed as: The statistical meaning of the performance indicator defined by Eq. ( 2) is illustrated by the scheme in Fig. 2, where two quantile-quantile (q-q) plots between Q p n,s,y and Q p n,s at different times y and y + 1 are reported (black full line).The black dashed line in the q-q plot indicates the perfect stationarity, when the Q p n,s,y and Q p n,s have the same distribution.The mean bias in the quantile residues is indicated in Fig. 2 by a black dotted arrow.The grey full line schematises the q-q plot after removing the mean bias.The remaining error, expressed by the variance, is indicated on the plot by a greyshaded area.In the following, we will present the standard deviation Qstd p n,s,y in spite of the variance as it can be more intuitively expressed in the unit of the considered variable.
The seasonal time variation of the Standard Deviation of Quantiles were then computed as: The mean of quantiles Qmb Reporting the time variation of the quantiles mean and the standard deviations enables us to separately quantify the mean and the unbiased non-stationarity, i.e., the nonstationarity in the frequency of the events.

Trend analysis
In the case of spatial heterogeneity in the observed trends, climate simulations either with or without downscaling should be able to resolve such a spatial variation.In order to quantify this ability, the annual Sen's slope (Sen, 1986) and the associated significance, through the Mann Kendall coefficients (Mann, 1945;Kendall, 1975) were computed over the whole study period at each node of the SI p n grid on the annual variables and referred as SS

Case study
The proposed methodology was applied to a meaningful case study located in Southern Italy, the Apulia Region, in which the climate and landscape features, including the water exploitation policy, represent a serious threat for water resources availability in the near future.The regional territory, with a total extension of 19 500 km 2 , is, in fact, mainly devoted to agriculture with more than 70 % of the total area occupied by cropped land which brought to a fast growing trend towards irrigation farming over the last four decades with a massive exploitation of groundwater resources.On the other hand, climate variables (rainfall in particular) exhibit a marked inter-annual variability, which makes water availability a worrying issue to the economic development and ecosystem conservation of the region (Portoghese et al., 2013).
Monthly observations from 77 temperature stations and 111 rainfall gauge stations covering the period 1950-2000 were used as land measurements.From the original dataset provided by the Apulia Hydrographic Service, only stations with less than 20 % of missing data were selected.Figure 3 shows the location of the temperature and precipitation stations, whose spatial density is about 1 per 2.76 × 10 2 km 2 and 1 per 1.91 × 10 2 km 2 , respectively.In the following, we will refer to the cumulated precipitation on a yearly and monthly basis as annual and monthly precipitation, respectively; the daily minimum temperature averaged over one year (one month) will be referred as annual (monthly) minimum temperature; similar definitions are adopted for annual and monthly maximum temperature.
Figure 4 shows the spatial distribution over the case study region through the associated spatial quantiles (5th, 25th, 50th, 75th and 95th quantiles) of the five time series resulting from the four data processing methods and the land observation (ref), for annual precipitation, minimum and maximum temperature.This gives an overview of the GCM mismatch in space and time, and the impacts of subsequent downscaling processes.

Results
In the next sections, the performances of the adopted downscaling methods and their combination are compared through the indicators presented in Sect.2.2.

Mean bias
The spatial variability of the mean bias M p n (Eq. 1) computed between (ref) and each of the data processing results is shown, in terms of percentiles (25th, 75th, 5th and 95th), in Fig. 5, while the numerical results are reported in Table 1. Figure 5 can be read as follows: the closer the mean bias to zero, the higher the ability of the data processing to reproduce the spatial mean condition for each variable; the narrower the distribution, the higher the ability of the data processing to reproduce the spatial heterogeneity of each variable.
The mean bias analysis highlights the poor ability of the adopted GCM to reproduce the spatial mean behaviour of precipitation in relation to the different seasons: a large overestimation is evident during winter and a large underestimation during summer (+15.4 and −20.5 mm, respectively) resulting in the low mean bias at annual scale (−2.3 mm).During spring and autumn, the GCM shows intermediate performances.Moreover the GCM's mean bias is associated with a large spatial heterogeneity, except in spring and summer.The application of the DD allows to reduce the mean bias in winter (−0.1 mm) and summer (+7.5 mm), but its performance degrades significantly in spring (+24.2 mm).The DD in general slightly reduces the spatial heterogeneity, while the SD is successful in reducing (by at least one order of magnitude) the mean bias and its variance, independently from the season.Finally, the combined DD-SD presents further improvements in reducing the mean bias and variance (0.14 mm).The Anova test carried out on the mean bias resulting from the data processing (3) and (4) confirms the significance of DD-SD improvement versus the SD alone in the spring, autumn and summer.
The GCM performance for minimum and maximum temperature is similar to those reported for precipitation.An overall mean bias is found for both variables (typically about ± 2 • C, respectively).The minimum temperature is systematically overestimated, while the maximum temperature is underestimated.The DD significantly reduces the mean bias for both variables (∼ 1 • C), except for the winter maximum temperature, but keeps almost unchanged the spatial heterogeneity.The SD reduces the annual and seasonal mean bias and its spatial heterogeneity by at least one order of magnitude (∼ 0.1 • C).Finally, also for temperature the combined DD-SD presents the best results (Table 1) with a further reduction in all the percentiles (0.07 and 0.08 • C for minimum and maximum temperature, respectively).Also in this case, the Anova test applied to the mean bias highlights the significance of DD-SD improvement versus the single SD for the annual, spring and summer minimum temperature, and for the annual and seasonal maxima.

Mean of quantiles
The analysis of the time evolution of quantiles provides clear information about the non-stationarity of local climate.Figure 6 shows the time variation of the mean of quantiles Qmb p s,y (Eq.4).The absolute value of Qmb p s,y indicates how much the quantile of each year (considering a 21 yr window) differs from the mean of quantiles computed over the whole period.In general, a flat signal (centred on 0 by construction) indicates that the considered variable is stationary along the analysed period.On the contrary, the non-stationarity could be detected by the presence of trends.The amplitude of the trend is directly expressed by the amplitude of the Qmb p s,y variation in the unit of the considered variable.To support these results, the p-value associated with a Mann Kendall test is computed over the whole period for each data processing method and for the (ref) dataset (Table 2).
The observed precipitation presents a negative trend in winter and spring, while summer and autumn are characterised by an initial increase in precipitation, followed by a negative trend and by a stationary period.In the case of model data, only in winter and at annual scale the negative trend results significant over the whole period (Table 2).The GCM correctly reproduces the annual behaviour, but underestimates the negative winter and spring trends, which result not significant for the Mann Kendall test (Table 2).Furthermore, the model does not reproduce the stationarity in the autumn observed during the 1990s.The SD has a negligible impact when combined with the climate model output.In the case of DD the global driver is expected to play a leading role in determining long-term tendencies (Dequé et al., 2007).However, the DD modulates the GCM output at the local scale, leading to a better representation of the observations in most of the cases (annual, winter and spring).
The observed minimum temperature presents positive trends in winter and spring, and a relative stationarity during the first half of the considered period, followed by a positive trend during the second half in summer and autumn.The observed annual time series of the minimum temperature is stationary until late 1970s, followed by a positive trend.All the observed trends, except in autumn, are significant over the whole period (Table 2).The GCM underestimates the positive trend and the associated significance.Benefits from the different downscaling methods are similar to those discussed for precipitation.
The observed maximum temperature is stationary in winter, whereas a negative trend is observed in spring, summer and autumn during the first half of the considered period; the second half is characterised by a positive trend.The GCM fails in reproducing both the annual and the seasonal non-stationarity, overestimating the positive trends and underestimating the negative ones.As for precipitation and minimum temperature, the DD properly modulates the GCM outcome, mainly enhancing the positive trends when it is already present in the GCM, and generating positive trends when the raw GCM output has stationary behaviour.Likewise, the SD has a negligible impact when combined with the DD.

Standard deviation of quantiles
The analysis of the time evolution of the variance of quantiles, hereinafter referred to as unbiased non-stationarity, describes the non-stationarity in the frequency of events of given intensity.Figure 7 shows the unbiased non-stationarity Qstd p s,y (Eq.6) used to describe the evolution of the standard deviation between the quantiles computed on a moving 21 yr window and those computed over the whole period.The absolute value of Qstd once the mean bias is removed.A flat signal indicates that the probability distribution of a variable is basically stationary throughout the analysed period and a different quantile distribution during each of the 21-yr window is a side effect of the subsampling.For example, this is the case for the GCM summer rainfall.The non-stationarity observed during the 1950s and during the 1990s may be affected by the variable size of the time windows shorter than 21 yr (grey rectangles), and will not be discussed.
The observed precipitation presents non-stationarity from the half to the late 1980s in autumn and at the annual scale.
Instead, the GCM reproduces correctly the observed pattern of the unbiased non-stationarity at the annual scale, in winter and autumn, whereas it slightly underestimates the results in spring.The DD enhances the non-stationarity simulated by the GCM in spring and summer and slightly reduces it in winter and autumn.This results in a better representation of the observed unbiased non-stationarity at annual scale, in spring and autumn.In particular, the results obtained for the summer unbiased non-stationarity suggest a key role for local processes at a spatial scale which is not well captured by the GCM.The SD has a low impact on the unbiased non-stationarity when it is directly applied to the GCM, except in summer when the non-stationarity is enhanced.Combined with the DD, the SD mostly reduces the non-stationarity when overestimated (spring and summer) and systematically lays between the underestimated GCM and the overestimated DD non-stationarity.In term of unbiased non-stationarity, the combined DD-SD presents a high covariance with the DD and a mean value similar to the SD results when applied directly to the GCM.
In the case of minimum temperature, non-stationarity is observed from the half of 1960s to the half 1970s in summer, and from mid 1970s to mid 1980s in winter.The GCM reproduces correctly the observed pattern of unbiased nonstationarity at annual scale in winter and spring, but results generally underestimated, except in autumn.In terms of relative impact of the downscaling, both the DD and the SD modulate the unbiased non-stationarity of the GCM mostly by increasing the standard deviation of quantiles.As for precipitation, the combined DD-SD presents a high covariance with the DD and an amplitude similar to the SD results when applied directly to the GCM.Compared to the precipitation, minimum temperature presents relatively low differences among data processing, except in spring after the 1970s, where the combined DD-SD better represents the reference.
For maximum temperature, non-stationarity is observed from early 1970s to early 1980s in summer, and from mid-1970s to mid-1980s in winter and spring.The GCM generally fails in reproducing the observed level of unbiased non-stationarity which results systematically underestimated.In terms of relative impact of the downscaling both the DD and the SD strongly modulate the GCM by systematically increasing the unbiased non-stationarity.In particular SD only and combined DD-SD lay always between the underestimated GCM and the overestimated DD signals.

Trends analysis
The spatial distribution of the Sen's slope SS p n for the annual precipitation, minimum and maximum temperature are reported in Fig. 8.The variance of the spatial distribution of the Sen's slopes defined as Var p in Eq. ( 7) is reported in Table 3 for each data processing and (ref).
The trend slope in the observed annual precipitation presents a large spatial heterogeneity, with values ranging from −1.4 mm yr −1 in the central areas of the case study to −7.2 mm yr −1 in the North.Most of these trends are significant, except in the extreme South.Because of its low resolution, the GCM almost does not reveal any spatial heterogeneity (mean trend of −0.6 mm yr −1 ).In spite of a relatively small spatial heterogeneity in temperature, the DD modulates the GCM spatial rainfall trends, generating trends ranging from −1.0 mm yr −1 in the North to −3.3 mm yr −1 in the South.Thus, the DD is able to reproduce almost half of the observed variance (0.34 and 0.74 (mm yr −1 ) 2 , respectively).In general, the SD generates lower spatial variance than the DD (0.11 (mm yr −1 ) 2 ) as indicated by the trend slopes ranging from −0.2 to −2.1 mm yr −1 .The resulting covariance is of the same order of magnitude as the DD (6.5 % of the observed variance), but positive.Finally, the combination DD-SD presents the highest spatial variance (0.40 (mm yr −1 ) 2 ) with trend slopes ranging from −0.3 to −3.5 mm yr −1 .Neither the GCM nor further downscaling processes have shown significant trends in the annual values.
The observed trend slopes in annual minimum temperature present a large spatial heterogeneity, with values ranging from −0.024 • C yr −1 in the central areas of the case study to +0.045 • C yr −1 in the extreme North and in the extreme South.In general, significant positive trends are dominant in most of the study region.On the contrary, the GCM does not show any spatial heterogeneity, with a mean trend of +0.010 • C yr −1 .The DD does not modulate the spatial trends of the GCM leading to a mean trend of +0.015 • C yr −1 with almost no spatial variance (0.1 % of the observed variance), leading to a negative covariance of −0.3 % of the observed variance.The SD shows a spatial variance slightly higher than the DD (2.4 % of the observed variance) with trends ranging from +0.012 to +0.024 • C yr −1 , which are still far from a correct representation of the observed spatial heterogeneity.Finally, the combination DD-SD presents trend slopes slightly higher than the other downscaling (ranging from +0.016 to +0.024 • C yr −1 ) and the best results in terms of covariance with the reference (5 % of the observed variance), but still fails in representing correctly the (ref).In general, all the trend slopes resulting from the GCM and further downscaling are statistically significant.
The observed slope of annual maximum temperature shows spatial heterogeneity larger than the annual minimum temperature, with values ranging between −0.040 • C yr −1 in the South and +0.034 • C yr −1 in the centre and the North-West, both minimum (at South) and maximum (centre and North-West) slopes resulting significant.The GCM presents a mean trend of +0.012 • C yr −1 , not significant in the northern and central portions of the region.The DD slightly modulates the GCM spatial trends with slopes ranging from +0.017 to +0.025 • C yr −1 , though far from the spatial variance of the (ref) (0.8 % of the observed variance) and generating a positive spatial covariance of 1.1 % of the observed variance.The DD also enhances the GCM slopes significance in the entire region.The SD slightly increases the spatial variance compared to the DD (4.8 % of the observed variance) with trends ranging from +0.014 to +0.030 • C yr −1 , though with lower covariance with the reference than the DD (0.9 % of the observed variance).Finally, the combination DD-SD shows similar results to the DD, with slopes ranging from +0.016 to +0.025 • C yr −1 and significant in all grid boxes, but a negative covariance with the reference (−2.9 % of the observed variance) and, hence, a poor representation of observations.
The different behaviour of the DD temperature and rainfall in terms of their spatial heterogeneity is linked to two important aspects of the dynamical downscaling process.First, in DD the local topography can generate rainfall patterns even in the absence of significant temperature inhomogeneities at the surface.Second, local feedbacks affecting rainfall patterns and involving convective instability depend on soil moisture and, therefore, on the land use characteristics and land-sea contrast which can amplify irregular rainfall patterns even in the presence of weak temperature gradients (Schär et al., 1999).

Discussion
Some considerations on the intrinsic limitations due to the use of a single case study, a single GCM, a single DD and a single SD method may help to better contextualise the results obtained through the proposed indicators of performance.
The uncertainty introduced by the choice of the driving GCM, as well of the dynamical and/or statistical downscaling technique, was recently evaluated using ensemble run (Chen et al., 2006;Déqué et al., 2007;Fowler et al., 2007).It has been concluded that the choice of driving-GCM generally provides the largest source of uncertainty both for the RCM (Déqué et al., 2007) and the SD scenarios (Cavazos and Hewiston, 2005;Wilby and Wigley, 1997;Fowler et al., 2007).
Within this framework, the presented results are indicative of the relative role of each downscaling processing rather than of their absolute performances, which depends to a large degree on the quality of the driving GCM.In the presence of complex orography and land-sea contrast the DD approach considered in this study produces physically coherent patterns in the tendency of key impact indicators, which is a desired characteristic for the generation of usable climate scenarios (Dell'Aquila et al., 2012).Therefore, although the DD is certainly not sufficient for improving the quality of climate scenarios (e.g., the regional climate model may have its own deficiencies), it does appear to be a necessary pre-condition for the generation of climate scenarios that are usable for impact modelling in those cases when local dynamics (e.g., the mesoscale, between 100 and 1000 km) and feedbacks (e.g., interactions in the soil-vegetation-atmosphere system) are poorly represented in a GCM.This aspect is of particular relevance for the Mediterranean area, where recent studies suggest that the DD of global simulations does improve specific aspects of the modelling of regional climate (Dubois et al., 2011;Dell'Aquila et al., 2012;Gualdi et al., 2012).
The presented results reveal that by improving aspects of the local dynamics, the DD modulates the quantile distributions produced by the GCM, thus, improving the GCM predictive performance (Fig. 6).Although the GCM obviously plays a major role in reproducing long-term climate scenarios, the DD, when properly tuned to represent the local environmental conditions, is able to deviate significantly from the GCM behaviour at the interannual scale (Fig. 7).For example, simplified conceptual models demonstrate that the feedbacks between surface hydrology and the local energy budget support a large variability of the soil-vegetationatmosphere system, especially in water-limited areas where transitions between wet/cool and dry/hot conditions are possible (Baudena et al., 2008).In particular, past studies performed with the same atmospheric model adopted here shown that the variability of maximum surface temperature is sensitive to changes in the land-cover characteristics (e.g., Anav et al., 2010).Therefore, a more accurate representation of land-sea contrast and land cover characteristics adopted for the DD is expected to produce significant deviations from the GCM and the amplification of its unbiased nonstationarity shown in Fig. 7. Downstream SD may then be employed as a further correction of residual bias: although the quantile distribution of the SD follows essentially the background variability of the large-scale driving climate, regardless of its provenance from the GCM or from the RCM output, SD can be considered as a practical tool for removing the model's systematic bias.
The need of using the whole chain processing is confirmed by the fact that in this study GCM-DD-SD shows always the lowest mean bias compared to the reference observation network.A direct consequence of the combined benefits of a DD-SD processing of the GCM scenarios is that spatial heterogeneities in long-term climate fluctuations are captured only when both DD and SD are included in the processing chain.It is worthy to note that, as the driving control climate is not initialised with observations, the patterns shown in Fig. 8 is not expected to closely follow the observed ones.Nevertheless, especially for the case of temperature, none of the standalone approaches, either DD or SD, produce significant spatial heterogeneity, which starts to be detectable only in the case of the combined processing (corresponding to about 50 % of the observed heterogeneity for precipitation).

Conclusions
The present study aimed to further investigate the advantages and limitations of using a DD as a pre-processing before SD to obtain monthly scenarios of temperature and precipitation, to be used for hydrological simulation at local and/or basin scale.
Even if the study is limited by the singularity of the case study and the adopted models (GCM, DD and SD), the results are of general usefulness for climate impact modelling at the local and basin scale.
The sizeable effect of the DD on the description of nonstationarity of local climate, especially in the case of rainfall and of maximum temperature, highlights the key role of local processes (including the triggering of convection and the surface energy balance) in characterising the local climate.Our analysis suggests that SD is a necessary step in the processing of climate simulation for obtaining reliable statistics at the local scale.For example, the quantile-quantile transform is confirmed as one of the best SD tools for the removal model bias from meteorological variables.However, an explicit modelling of the physical system at a sufficiently high resolution (hence, the DD) appears a necessary pre-condition to a skillful SD, especially during the seasons in which local processes have a larger control on local fluctuations of climate.In particular, DD plays a key role in characterising the spatial distribution of trends.Moreover, only the DD is able to modulate the inter-annual variability simulated by the GCM by enhancing the role of local feedbacks, for example in the soil-vegetation-atmosphere system.However, it is worthy to note that for the GCM scenarios considered in this study, the correction introduced by the DD is not sufficient to reproduce the observed trends.
The resulting complementarity of the two downscaling techniques suggests that the combined DD-SD is a suitable choice for the generation of weather scenarios for impact modelling.In fact, the combined DD-SD presents the best results, both in terms of mean bias and spatial distribution of trends by retaining the improvements obtained by the DD in terms of climate non-stationarity as well as the intrinsic assimilation of observation datasets.

Fig. 1 .
Fig. 1.Methodological framework representing the adopted methods of data processing.The arrows indicate the data fluxes, while models (GCM and relevant downscaling) and land observations are shown with ellipses.The Statistical Interpolation (SI) is represented by a dashed rectangle.Data processing resulting from the data flux are referred as: (1) GCM; (2) DD applied to GCM;(3) SD applied directly to the GCM; (4) SD applied to the DD of the GCM; (ref) land observations.The spatial scale associated with each model is reported on the left.
,y were then compared with quantiles of the same node computed over the whole period Q p n,s .The Q p n,s were computed using the same plotting position as for the associated Q p n,s,y .The non-stationarity in the climate is then revealed by the time variation of the residues between the quantile Q p n,s,y and the quantile Q p n,s .Similarly we defined the quantiles of the reference as Q ref n,s,y and Q ref n,s .The ability of each data processing to reproduce the observed nonstationarity is revealed by the comparison with the analogue time variation of the residues computed for the reference dataset

Fig. 2 .
Fig. 2. Decomposition of non-stationarity in the quantile-quantile plots.Full black line indicates quantile-quantile plot.Dot black arrow indicates associated mean bias and the grey full line the unbiased quantile-quantile plot.The remaining variance is indicated by the grey surface.
p n,s,y and the standard deviation of quantiles Qstd p n,s,y were further averaged over the n grid nodes of the SI and indicated as Qmb p s,y and Qstd p s,y , respectively.
p n .The SS p n spatial distribution variance of each data processing was computed as an indicator of the spatial heterogeneity of the trend amplitude:

Fig. 3 .
Fig. 3. Location of Apulia region.The hydrological domain area is delimited by a grey full line.Locations of the temperature and precipitation sampling stations are shown with grey and black full circle, respectively.GCM nodes are shown with black stars and DD nodes are shown with black crosses.The grid boxes associated with GCM and DD nodes are delimited by black full line.

Fig. 6 .
Fig. 6.Annual and seasonal mean of quantiles variation between the quantiles computed on a 21 yr window and on the whole period.Precipitation, minimum and maximum temperature analysis for land observations (ref), GCM (1), GCM-DD (2), GCM-SD (3) and GCM-DD-SD (4).Grey rectangles represent values computed with less than 21 yr.

Fig. 7 .
Fig. 7. Annual and seasonal standard deviation of quantiles between the quantiles computed on a 21 yr window and on the whole period.Precipitation, minimum and maximum temperature analysis for land observations (ref), GCM (1), GCM-DD (2), GCM-SD (3) and GCM-DD-SD (4).Grey rectangles represent values computed with less than 21 yr.

Fig. 8 .
Fig. 8. Spatial distribution of Sen's slopes resulting from each data processing for annual cumulated precipitation, minimum temperature and maximum temperature.Grid boxes marked with stars are those in which the estimated trend is statistically significant.

Table 1 .
Annual and seasonal spatial distribution percentiles of the mean bias.
The bold values indicate the minimum bias among data processing.

Table 2 .
Annual and seasonal Mann Kendall p-value over the period 1953-2000.
The bold values indicate p-value within the 95 % of confidence.

Table 3 .
Variance of Annual Sen's slope spatial distribution.