Improvement of Hydroclimatic Projections over Southeast Spain by Applying a Novel RCM Ensemble Approach

Climate model outputs can be used as climate forcing for hydrological models to study the impact of climate change on the water cycle. This usually propagates cumulative uncertainties, transferring the errors from the climate models to the hydrological models. Then, methodologies are needed to evaluate the impact of climate change at basin scale by reducing the uncertainties involved in the modeling chain. The paper aims to assess the impact of climate change on the runoff, considering a novel approach to build a Regional Climate Model (RCM) ensemble as climate forcing for a parsimonious spatially distributed hydrological model. A semiarid basin of southeast of Spain was selected for the study. The RCM ensembles were built based on seasonal and annual variability of rainfall and temperature. If the runoff projections for 2021–2050 are compared to the 1961–1990 observed period, a significant decrease in runoff equal to −20% (p-value t-test 0.05) was projected. However, by changing the observed period to 1971–2000, a despicable change (2.5%) is identified. This fact demonstrates that trends based on short records are very sensitive to the beginning and end dates, due to the natural variability. Special attention should be paid to the selection of the period for impact studies.


Introduction
Numerous studies identified the Mediterranean region as one of the most vulnerable areas to climatic and anthropogenic changes and to water crisis [1][2][3][4].During the latter half of the 20th century, cumulative discharge from many mid-latitude rivers decreased by 60%, while high-latitude rivers experienced increased discharge [5].The impacts of climate change can be exacerbated when occurring in regions with imbalances between water demands and available water resources [1], such as the southeast of Spain.Moreover, climate change has a larger effect than human activities on decreases in runoff [6], so it must be studied carefully.
For understanding the plausible impacts of climate change on hydrological processes and water availability, the regional climate models (RCMs) and the global climate models (GCMs) are valuable tools, simulating precipitation, temperature and wind speed, among other variables.Nevertheless, the uncertainties in the climate projections, such as unknown future greenhouse gas emissions [7] and the simplified representation of processes, are still high [8,9].
Usually, two approaches are used to characterize the uncertainties of climate projections: perturbed physics ensembles [10] and multi-model ensembles generated with different initial and boundary conditions [11,12].On the one hand, some authors calculate multi-model ensembles as a simple mean of several models [13].On the other hand, the use of weighted averages where the model weight depends on some measure of performance, and several techniques to assess the performance of each model, are widely applied [14][15][16][17][18][19][20][21][22].Considering this last approach, the performance of two novel approaches to build RCM ensembles of monthly precipitation over Spain was evaluated [23].While [23] found the best ensemble approximation was based on the modified Reliability Ensemble average method [16,21] for weighting the ensemble RCM members.
Another problem with climate model output is the systematic error or bias between historical observations and climate model hindcast [24].Bias correction methods are increasingly applied in climate model studies and serve to overcome this error [24,25].The mean bias of the statistics of daily precipitation simulated by five RCMs was evaluated by [24].While [25] applied a bias correction of projected precipitation, provided by three RCMs and potential evapotranspiration projections, as input to a simple hydrological model.However, when many climate models are considered and the interdependency is important, other solutions should be taken into account.The building of RCM ensembles based on the performance to reproduce a historical period implicitly penalizes (less weight in the ponderation) those RCMs which exhibit higher bias in relation to the historical dataset.Therefore, it should be considered as an indirect method to take into account the bias or systematic errors, with the advantages of spatial assessment.
Due to both problems of uncertainty of climate projections, there are different strategies in reducing uncertainties (Figure 1).On the one hand, many researchers [25-29] investigate the impact of climate change on future water resources employing a hydrological model (HM) driven by one specific climate model (CM), as it is presented in Figure 1a.Then, a hydrological projection ensemble (PR ensemble ) could be built from the results of each simulation (Figure 1b).However, the effects of variability of precipitation and temperature from all available CMs have not been evaluated.
Water 2018, 10, 52 2 of 19 Usually, two approaches are used to characterize the uncertainties of climate projections: perturbed physics ensembles [10] and multi-model ensembles generated with different initial and boundary conditions [11,12].On the one hand, some authors calculate multi-model ensembles as a simple mean of several models [13].On the other hand, the use of weighted averages where the model weight depends on some measure of performance, and several techniques to assess the performance of each model, are widely applied [14][15][16][17][18][19][20][21][22].Considering this last approach, the performance of two novel approaches to build RCM ensembles of monthly precipitation over Spain was evaluated [23].While [23] found the best ensemble approximation was based on the modified Reliability Ensemble average method [16,21] for weighting the ensemble RCM members.
Another problem with climate model output is the systematic error or bias between historical observations and climate model hindcast [24].Bias correction methods are increasingly applied in climate model studies and serve to overcome this error [24,25].The mean bias of the statistics of daily precipitation simulated by five RCMs was evaluated by [24].While [25] applied a bias correction of projected precipitation, provided by three RCMs and potential evapotranspiration projections, as input to a simple hydrological model.However, when many climate models are considered and the interdependency is important, other solutions should be taken into account.The building of RCM ensembles based on the performance to reproduce a historical period implicitly penalizes (less weight in the ponderation) those RCMs which exhibit higher bias in relation to the historical dataset.Therefore, it should be considered as an indirect method to take into account the bias or systematic errors, with the advantages of spatial assessment.
Due to both problems of uncertainty of climate projections, there are different strategies in reducing uncertainties (Figure 1).On the one hand, many researchers [25-29] investigate the impact of climate change on future water resources employing a hydrological model (HM) driven by one specific climate model (CM), as it is presented in Figure 1a.Then, a hydrological projection ensemble (PRensemble) could be built from the results of each simulation (Figure 1b).However, the effects of variability of precipitation and temperature from all available CMs have not been evaluated.On the other hand, the assessment of climate change impacts on water resources using RCM and GCM ensembles and hydrological models has become an important challenge around the world [30][31][32][33][34] (Figure 1c).A multi-model average ensemble approach from 11 GCMs was used by [30] to drive a macroscale hydrology model in the Colorado River Basin.While [32] assessed actual evapotranspiration and streamflow with the Soil and Water Assessment Tool (SWAT) model in a subtropical Chinese basin for the 21st century.These authors used multi-model average ensemble approach from six GCMs (where each member is assigned by an equal probability of occurrence).An increase of the annual mean actual evapotranspiration ensemble and decrease of the annual mean streamflow compared to the baseline period  were projected [32].In these works, although attempts to reduce climatic uncertainties have been made using average climate On the other hand, the assessment of climate change impacts on water resources using RCM and GCM ensembles and hydrological models has become an important challenge around the world [30][31][32][33][34] (Figure 1c).A multi-model average ensemble approach from 11 GCMs was used by [30] to drive a macroscale hydrology model in the Colorado River Basin.While [32] assessed actual evapotranspiration and streamflow with the Soil and Water Assessment Tool (SWAT) model in a subtropical Chinese basin for the 21st century.These authors used multi-model average ensemble approach from six GCMs (where each member is assigned by an equal probability of occurrence).An increase of the annual mean actual evapotranspiration ensemble and decrease of the annual mean streamflow compared to the baseline period  were projected [32].In these works, although attempts to reduce climatic uncertainties have been made using average climate ensembles, the main weaknesses are related to the large hydrological uncertainties due to the parameterization of complex hydrological models.For example, both SWAT model used by [32,34], and Variable Infiltration Capacity (VIC) used by [30,33], present a large number of parameters and some of them are difficult to obtain.
However, there is no systematic approach available for reducing the combined uncertainty in streamflow projections from both climate projections and hydrological models [31].These authors proposed three different approaches, and concluded that reducing climate input uncertainty (combining climate models, Figure 1c) provides better results than multi-model streamflow forecasts (Figure 1b).Since the multi-model streamflow is obtained by combining the results from a single hydrological model driven by individual climate projections (Figure 1b).
In the present work, robust RCM ensembles of rainfall, minimum and maximum temperature were built based on the revisited REA method proposed by [24], giving different weights to RCMs according to their performance in simulating the yearly and seasonal present climate.In other words, the time series of the selected meteorological variables were built on seasonal and yearly timescales, then empirical Cumulative Density Functions (CDFs) were fitted.The reliability factors for each RCM were computed as the p-value of the two sample Smirnov-Kolmogorov tests (TSSK) between the RCM data and the observed dataset.Finally, the greater TSSK p-value defines the higher weight of the RCM in the final ensemble, because the final performance criterion is a combination of the reliability factors.
Then, with the aim to reduce the uncertainties in the hydroclimatic projections, a hydrological distributed model with only four parameters was forced by the climate ensemble (monthly precipitation and temperature) over the semiarid Segura River Basin (southeast of Spain).This basin presents the greater water scarcity of Spain, and it is usually affected by severe droughts.
In addition, the sensitivity of hydrometeorological projections to historical period selection has been evaluated.Two time periods have been considered for the analysis, the common 1961-1990 and 1971-2000.The latter one was affected by a strong El Niño event according to [35], and it was employed by several authors [32].The authors [36] found good predictability when El Niño-Southern Oscillation (ENSO) indexes were used to forecast dry events in the Mediterranean coast of Spain.Several authors summarize the influence of some climatic indexes over the Iberian Peninsula, and pointed out that the positive phase of North Atlantic Oscillation (NAO) index is related with below normal rainfall amounts over the Iberian Peninsula [37].Also, [37] indicated that the warming trends of climate would lead in more positive phases of NAO, which likely explain the precipitation negative trends over the Iberian Peninsula.
To the best of our knowledge, there are no published works considering distributed hydrological models driven by RCM ensembles based on seasonal and annual variability of meteorological variables, in semiarid water-stressed basins with strong problems relating to water scarcity and droughts.This type of approach constitutes a way to reduce the uncertainties involved in the climate and hydrological models, giving more support to the decision-making process.This is especially important in semiarid basins, such as the Segura River Basin and its headwaters basin, which have a high vulnerability to water stress and a fragile balance between water resources and demands.A particular case study was selected among the several headwaters of Segura River Basin: the Fuensanta Basin, which is one of the most important sources of water resources for the Murcia region in Spain.
The present work is organized into five sections.The description of the used dataset and a deeper description of the study area are presented in Section 2. Section 3 discusses the methodological aspects of the climate model ensemble building, hydrological model and calibration and validation procedures.The main results of the research are presented in Section 4, and finally the conclusions are presented in Section 5.

Study Area
The variability of rainfall in the southeast of Spain can reach extreme irregularity.Rains are more frequent in autumn (between October and December).The study region corresponds to the Segura River Basin (SRB), with a surface area of 18,900 km 2 ; it is located in the semiarid southeast corner of the Iberian Peninsula (Figure 2).Average precipitation in the region ranges from 1000 mm/year in the headwater areas to 300 mm/year in the driest lowlands, while the potential average evapotranspiration is approximately 1400 mm/year.The main characteristic of the SRB is its water scarcity, considering the available water resources and the high agricultural and urban demands.The water exploitation index with a value of 237.2% [1], estimated from the mean annual runoff and mean annual total water demand, is the highest of all Spanish basins.The available water resource per inhabitant is only 442 m 3 /inhabitant/year for SRB, therefore the gap between water supply and demand is extremely high.In consequence, this basin is considered as one of the most water-stressed regions in the western Mediterranean basin.Consequently, the water transfers from other basins together with desalinization are considered the most attractive options to increase water availability in the basin.Moreover, recurrent and frequent drought periods (1980-1982, 1993-1995, and 2005-2008) have generated important crises and greatly complicated water management.
Once the RCM ensembles were obtained for temperature and rainfall projections over SRB, the impact analyses on hydrological projections were focused on the contribution of the basin to the Fuensanta River Reservoir (named Fuensanta Basin).Fuensanta Basin, located in the headwater of the Segura River Basin (Figure 2), covers an area of 1220.6 km 2 with an average altitude of 1263 m.This basin is one of the main sources of water for the whole Segura River Basin.Therefore, the strategic location of the Fuensanta Basin makes it one of the most suitable for studies of impacts on runoff generation, and thus for improving the understanding of the climate change impacts in all economic activities and livelihoods of the region.

Study Area
The variability of rainfall in the southeast of Spain can reach extreme irregularity.Rains are more frequent in autumn (between October and December).The study region corresponds to the Segura River Basin (SRB), with a surface area of 18,900 km 2 ; it is located in the semiarid southeast corner of the Iberian Peninsula (Figure 2).Average precipitation in the region ranges from 1000 mm/year in the headwater areas to 300 mm/year in the driest lowlands, while the potential average evapotranspiration is approximately 1400 mm/year.The main characteristic of the SRB is its water scarcity, considering the available water resources and the high agricultural and urban demands.The water exploitation index with a value of 237.2% [1], estimated from the mean annual runoff and mean annual total water demand, is the highest of all Spanish basins.The available water resource per inhabitant is only 442 m 3 /inhabitant/year for SRB, therefore the gap between water supply and demand is extremely high.In consequence, this basin is considered as one of the most water-stressed regions in the western Mediterranean basin.Consequently, the water transfers from other basins together with desalinization are considered the most attractive options to increase water availability in the basin.Moreover, recurrent and frequent drought periods (1980-1982, 1993-1995, and 2005-2008) have generated important crises and greatly complicated water management.
Once the RCM ensembles were obtained for temperature and rainfall projections over SRB, the impact analyses on hydrological projections were focused on the contribution of the basin to the Fuensanta River Reservoir (named Fuensanta Basin).Fuensanta Basin, located in the headwater of the Segura River Basin (Figure 2), covers an area of 1220.6 km 2 with an average altitude of 1263 m.This basin is one of the main sources of water for the whole Segura River Basin.Therefore, the strategic location of the Fuensanta Basin makes it one of the most suitable for studies of impacts on runoff generation, and thus for improving the understanding of the climate change impacts in all economic activities and livelihoods of the region.

Datasets
The high-resolution observational dataset of temperature and rainfall (1961-2007), named Spain02, was provided by [38].The regular grid presents 0.2° of spatial resolution.This dataset was extended in the Segura River Basin, for the period 2008-2013, considering the meteorological stations provided by several institutions such as SIAM (Sistema de Información Agraria de Murcia)

Datasets
The high-resolution observational dataset of temperature and rainfall , named Spain02, was provided by [38].The regular grid presents 0. Sixteen RCMs from ENSEMBLES European project [42], driven by different GCMs for scenario A1B [7], were used (Table 1).This dataset presents 0.25 • of spatial resolution for the 1961-2050 period.Considering the climate effects of SRES (Special Reports on Emission Scenarios) and RCP (Representative Concentration Pathways) scenarios, [43] indicate that the median of the projected warming to the end of the century under the SRES A1B scenario is found between RCP6 and RCP8.5 scenarios.Analyzing mean temperature and precipitation change, [44] found a very high spatial correlation between RCP8.5 and A1B results, with 0.82-0.97for temperature changes and 0.59-0.92for precipitation changes depending on the region for mid-century, and even stronger correlation for both parameters towards the end of the century.
The data of observed runoff for calibration in the Fuensanta Basin was provided by the Segura River Basin authority (CHS, Hydrographic Confederation of Segura River Basin).

Building PDF Ensemble
The ensemble probability density function (PDF) was built following the next steps.Firstly, the monthly time series of precipitation and temperature were obtained for each grid point from the RCMs and observed data.Secondly, the reliability factor (R) was calculated using a measure to compare the agreement among the cumulative density functions (CDFs) of each RCM and from observed data on 1961-1990.Thirdly, R was used to obtain the normalized reliability factor (Pm).Finally, the ensemble PDF was estimated using the normalized reliability factor Pm as the weighting factor.
The revisited REA method [23] was proposed to build the RCM ensembles of temperature and precipitation.The REA method has been applied in several studies [16,[45][46][47].This method is a flexible tool to assess the mean, uncertainty range, and the production of probabilistic outcomes of the ensemble climate change from both GCMs and RCMs.
The original REA method [15,16] calculates reliability weights (R i ) for different models as a product of two reliability factors: model performance criterion (R B ), which represents the capability of each model to represent observed climate; and convergence criterion (R D ), which evaluate the distance of a projected change to the ensemble average.Therefore, the reliability factor for model i is calculated as follows [15]: where the parameters m and n are the criterion weights.Due to the difficulty for assessing the convergence criterion, because the reference probability density functions (PDFs) for future climate are not known, several authors [21,23,48] eliminated R D from the definition of reliability factor.
In the present work, the estimation of the reliability factor R i is addressed by a novel approach based on seasonal and yearly cumulative density probabilities (CDFs), proposed by [23] as follows: where 1 ≤ i ≤ number of RCMs.
Equation ( 2) states that each model weight is given by the multiplication of five functional factors (R Winter , R Spring , R Summer , R Autumn and R Yearly ) to measure the model ability as a function of the model bias for reproducing the seasonal (winter, spring, summer and autumn) and annual mean variability of present climate.Furthermore, the bias analysis is based in the assessment of p-value from the well-known two sample Smirnov-Kolmogorov test [49], between RCM CDFs and observed data CDFs for the historical period .Therefore, the p-value is considered as a metric to compute the reliability factors R winter , R spring , R summer , R autumn and R yearly .Previous works successfully considered this metric to build RCM ensembles for the analysis of dry spell lengths [25,45].
The CDFs from the observational datasets represent the reference in this period.The inverse of the Weibull equation was used to compute the empirical quantiles of CDF.
The values of R i were estimated in the grid points defined over Segura River Basin with a 0.25 • of spatial resolution (Figure 2), and were interpolated applying the Kriging method for each model.The spatial distribution of the normalized reliability factor (Pm i ) was obtained for each RCM from the maps of R i , as follows [16]: where j = 1 to N, and N is the number of RCMs.The Pm i is the likelihood of each RCM.Then, based on Pm i as a weighting factor, the PDF ensembles (or CDF ensembles) was built on site 2021-2050 horizons, giving greater weight to RCMs with a high value of the Pm.The spatial interpolations methods named inverse distance weighting (IDW) and ordinary kriging presented a similar performance in the case of scarce spatial autocorrelation of the studied variable for complex terrain areas [50], because both are distance-based methods.The weights for the IDW are usually proportional to the inverse of the squared distance between the prediction point and the observation points.The precision of this geometric method is affected by the number of the closest neighboring samples.While kriging calculates the values of weights by estimating the spatial structure of the variable's distribution represented by a sample semivariogram.In the present case, the spatial representation of the studied variable obtained by the geostatistical method (kriging) was more realistic, therefore it was the selected method for the spatial interpolations.The assumption about the spatial correlation of the data implies that the closer the points are, the more similar are their values.Another relevant strength of kriging method is the ability to take into account the directional bias of data, which was expected because of the complex topography of studied region.Finally, the kriging method is a statistical method which computes the most likely value to interpolate, according to the spatial correlation which is considered trough the semivariogram.On the other hand, the IDW method does not take into account that spatial structure of data.

Hydrological Model and Its Calibration and Validation
A continuous, conceptual, and distributed hydrological model with few parameters, named the Témez model developed by [51], is applied.This model, which is widely applied in Spain [52], simulates monthly average flows in natural regime at any point of the drainage network.
The total runoff is produced from the rainfall not retained in the soil's upper layer, and not evaporated.The monthly model reproduces the essential water transport processes that occur in the different phases of the hydrological cycle, and establishes the principle of continuity and sharing laws and water transfers between different storages.
The meteorological input to the model corresponds to the datasets of precipitation (P) and potential evapotranspiration (PET).Information about sub-basins and maps of hydrogeological units are also required for model parameterization (H max : maximum soil holding capacity; I max : maximum infiltration rate; C: runoff coefficient; and α: aquifer discharge coefficient).As output, the model generates the spatial distributions of runoff, infiltration, soil water, and actual evapotranspiration.
In this study, before calibration, three parameters (H max , I max and α) were specified considering recommendations from [53].The parameter C was defined according to [54] guidelines.
Observed hydrometeorological datasets were considered to calibrate and validate the Témez model in the Fuensanta Basin.For calibrating the model, a representative period of six years (2000-2005) was used, where land cover was provided by CORINE Land Cover 2006 [55].During the time period 2000-2005, dry and wet years were observed with maximum and minimum flow rate fluctuations.In consequence, the parameters obtained in calibration phase are considered more robust and better able to represent both wet and dry years.
Calibration parameters were defined by minimizing the differences between observed runoff provided by the Segura River Basin authority (CHS) and simulated runoff by the model.The performance analysis was based on graphical and numerical goodness-of-fit, considering the Nash-Sutcliffe efficiency value (NSE) [56].The comparison between the observed and simulated hydrograph for the calibration phase in the period 2000-2005, at Fuensanta Reservoir Station (Figure 3a), demonstrates a good fit with a NSE value of 0.71.Then, the validation of the hydrological model was performed for a period of 30 years, from 1970 to 1999 (Figure 3b), with satisfactory results (NSE equal to 0.732).A limitation of hydrological models is the degree of predictability, that is, there is no assurance that a model obtains good performance when it is used under conditions different from those of the calibration period [57].In the present work, the selected calibration period was 5 years.Nevertheless, the validation process for a long period of 30 years demonstrates the good predictability of the model by the performance and reliability achieved [58].

Hydrological Model and Its Calibration and Validation
A continuous, conceptual, and distributed hydrological model with few parameters, named the Témez model developed by [51], is applied.This model, which is widely applied in Spain [52], simulates monthly average flows in natural regime at any point of the drainage network.
The total runoff is produced from the rainfall not retained in the soil's upper layer, and not evaporated.The monthly model reproduces the essential water transport processes that occur in the different phases of the hydrological cycle, and establishes the principle of continuity and sharing laws and water transfers between different storages.
The meteorological input to the model corresponds to the datasets of precipitation (P) and potential evapotranspiration (PET).Information about sub-basins and maps of hydrogeological units are also required for model parameterization (Hmax: maximum soil holding capacity; Imax: maximum infiltration rate; C: runoff coefficient; and α: aquifer discharge coefficient).As output, the model generates the spatial distributions of runoff, infiltration, soil water, and actual evapotranspiration.
In this study, before calibration, three parameters (Hmax, Imax and α) were specified considering recommendations from [53].The parameter C was defined according to [54] guidelines.
Observed hydrometeorological datasets were considered to calibrate and validate the Témez model in the Fuensanta Basin.For calibrating the model, a representative period of six years (2000-2005) was used, where land cover was provided by CORINE Land Cover 2006 [55].During the time period 2000-2005, dry and wet years were observed with maximum and minimum flow rate fluctuations.In consequence, the parameters obtained in calibration phase are considered more robust and better able to represent both wet and dry years.
Calibration parameters were defined by minimizing the differences between observed runoff provided by the Segura River Basin authority (CHS) and simulated runoff by the model.The performance analysis was based on graphical and numerical goodness-of-fit, considering the Nash-Sutcliffe efficiency value (NSE) [56].The comparison between the observed and simulated hydrograph for the calibration phase in the period 2000-2005, at Fuensanta Reservoir Station (Figure 3a), demonstrates a good fit with a NSE value of 0.71.Then, the validation of the hydrological model was performed for a period of 30 years, from 1970 to 1999 (Figure 3b), with satisfactory results (NSE equal to 0.732).A limitation of hydrological models is the degree of predictability, that is, there is no assurance that a model obtains good performance when it is used under conditions different from those of the calibration period [57].In the present work, the selected calibration period was 5 years.Nevertheless, the validation process for a long period of 30 years demonstrates the good predictability of the model by the performance and reliability achieved [58].

Assessment of Potential Evapotranspiration
Different PET input estimated by several methods (Penman-Monteith, Hargreaves-Samani, Jensen-Haise, and Hamon methods), can produce similar runoff simulations in both dry and humid regions [59].Therefore, the temperature-based Hargreaves-Samani model [60] can be considered as an alternative to the Penman-Monteith model [61] when only air temperature records are available to calculate PET.However, several studies [62,63] have shown that the Hargreaves-Samani model overestimates reference evapotranspiration (ET 0 ) values in humid regions and underestimates them in dry regions.For this reason, other studies have even proposed new parametric expressions to calculate the adjusted Hargreaves coefficient [64][65][66].
The spatial distribution of ET 0 is estimated from the observational dataset of temperature Spain02 and from the RCM ensembles of temperature, applying the method explained previously (Equation ( 2)), considering the Hargreaves model adjusted to the study area [64]: where ET 0 is the reference evapotranspiration (mm • day −1 ); R a is the water equivalent of the monthly-averaged daily extraterrestrial radiation (mm • day −1 ), obtained from tables based on the latitude and month according to [61] where 1 mm • day −1 = 0.408 × MJulio/m 2 /day; T max and T min are the monthly-averaged maximum and minimum values of daily air temperature ( • C); and T m is the monthly-averaged daily temperature, calculated as the average of T max and T min .
When ET is evaluated within the baseline balance of a basin, the concepts of reference evapotranspiration (ET 0 ) and potential evapotranspiration (PET) are equivalent.Formulas that were designed to calculate PET or ET 0 are used interchangeably, as shown [59,67].

RCM Ensembles
The kriging interpolation method was applied to obtain the maps of reliability factor (R i ). Figure 4 presents the results in Segura River Basin for R i for each of the RCMs in the case of the rainfall variable.The strengths and weaknesses of each RCM could be identified from the spatial distributions of values of R i , presented in Figure 4 and estimated with the Equation (2).Therefore, the RCMs with satisfactory skills (in blue) to simulate the observed data in the Segura River Basin are the following: C4IRCA3, CNRM/RM5.1,DMI/ECHAM5-r3, ETHZ/CLM and METNO/HadCM3Q0.The models with worse performance according to the obtained results (lower values of R i in Figure 4, in yellow), are DMI/ARPEGE, DMI/BCM, SMHI/BCM, SMHI/ECHAM5-r3, SMHI/HadCM3Q3 and UCLM/PROMES.
The CDFs of rainfall, minimum and maximum temperature for each RCM, and observational datasets, were obtained.Then, the CDF ensemble was built on each site of the study region, applying the methodology described in Section 3.1.Figure 5 presents a detail of results for seasonal and yearly rainfall on site 759 (Fuensanta Basin).The CDF ensemble in general captures well the seasonal and annual variability of observed rainfall, with more adjusted results in winter and summer (Figure 5).
In the case of maximum temperature CDFs on the same site (Figure 6), the goodness-of-fit between the CDF ensemble and Spain02 CDF is satisfactory.The best results correspond to the summer season and at yearly scale (0.794 and 0.522 p-values respectively).While for winter, spring and autumn the p-values are not satisfactory (0.134 winter, 0.366 spring, and 0.135 autumn).
For minimum temperature, the greatest differences between the Spain02 CDF and CDF ensemble occur in winter and at yearly scale (0.06 and 0.04 p-values), and the lowest differences in the summer season (0.634 p-value) (Figure 7).It is noteworthy (Figures 5-7) that the best performance varies for each analyzed season, site, and variable, and that there are instances where a single model performs better than the ensemble.However, this fact is not generalizable, and the ensemble approach in general reduces the bias compared to the single model approach.
A more exhaustive analysis of the sites and p-values for rainfall over Spain is presented in [23].The maps of RCM ensemble projections for mean annual rainfall (P), maximum temperature (Tmax), and minimum temperature (Tmin) for 2021-2050 period, and the differences (%) with respect to the 1961-1990 period, are presented in Figure 8.A decreasing gradient of P from west-east is expected for 2021-2050 (Figure 8a), with the highest P values corresponding to the Fuensanta Basin   The maps of RCM ensemble projections for mean annual rainfall (P), maximum temperature (Tmax), and minimum temperature (Tmin) for 2021-2050 period, and the differences (%) with respect to the 1961-1990 period, are presented in Figure 8.A decreasing gradient of P from west-east is expected for 2021-2050 (Figure 8a), with the highest P values corresponding to the Fuensanta Basin The maps of RCM ensemble projections for mean annual rainfall (P), maximum temperature (T max ), and minimum temperature (T min ) for 2021-2050 period, and the differences (%) with respect to the 1961-1990 period, are presented in Figure 8.A decreasing gradient of P from west-east is expected for 2021-2050 (Figure 8a), with the highest P values corresponding to the Fuensanta Basin (above 800 mm/year) and minimum values in the coast (less than 300 mm/year).The T max and T min (Figure 8c,e) have a spatial distribution with a strong increasing west-east gradient.Focusing on the difference maps (Figure 8b), decreases of P for 2021-2050 are expected in the middle east and some parts of the northwest, while increases (greater than 20%) in the southwest and north areas.On the other hand, the T max (Figure 8d) is predicted to increase by 10% to in the Fuensanta Basin, while an increase of 5% on average is expected for the rest of the basin.Finally, a general increase in the T min (Figure 8f), of between 10% and 20% is expected for most of the Segura River Basin, while in the Fuensanta Basin the highest increases are identified (greater than 60%).
(above 800 mm/year) and minimum values in the coast (less than 300 mm/year).The Tmax and Tmin (Figure 8c,e) have a spatial distribution with a strong increasing west-east gradient.Focusing on the difference maps (Figure 8b), decreases of P for 2021-2050 are expected in the middle east and some parts of the northwest, while increases (greater than 20%) in the southwest and north areas.On the other hand, the Tmax (Figure 8d) is predicted to increase by 10% to 20% in the Fuensanta Basin, while an increase of 5% on average is expected for the rest of the basin.Finally, a general increase in the Tmin (Figure 8f), of between 10% and 20% is expected for most of the Segura River Basin, while in the Fuensanta Basin the highest increases are identified (greater than 60%).

Trends for Meteorological and Hydrological Variables for 2021-2050 Horizon
Once calibrated and validated for the Fuensanta basin, the Témez hydrological model was used to obtain the hydrological projections driven by the RCM ensembles (rainfall and potential evapotranspiration) previously generated applying the methodology described in Section 3.1.The runoff projections were obtained for the period 1961-2050.
Figure 9 presents the projected anomalies of rainfall and runoff for the time period 2021-2050, related to the mean of the periods 1961-1990 and 1971-2000.The difference between the means from both time periods was estimated.The t-test was used to determine if the two means were significantly different from each other.In Table 2, the p-value of t-test is presented in parenthesis, identifying the significant differences (denoted by * in Table 2).If the mean annual runoff for the future (2021-2050) is compared with the historical dataset for 1961-1990 (Figure 9a), a plausible significant decrease (−20%) is expected (Table 2).However, from Figure 9a, the decrease of runoff is more clearly identified in the first part of the period.If the comparison is based on the difference between 2021-2050 and 1971-2000 period (Figure 9b), the expected trend of runoff would be positive (2.5% from Table 2).Analyzing the other variables in the hydrological cycle, the reduction of runoff for 2021-2050 (Table 2) in relation with the baseline period ) is justified in the reduction of rainfall and the increase of PET (≈4%).However, for 2021-2050 in relation to 1971-2000, an increase in rainfall over 11% and a similar increase of PET (≈3.7%) are expected.

Trends for Meteorological and Hydrological Variables for 2021-2050 Horizon
Once calibrated and validated for the Fuensanta basin, the Témez hydrological model was used to obtain the hydrological projections driven by the RCM ensembles (rainfall and potential evapotranspiration) previously generated applying the methodology described in Section 3.1.The runoff projections were obtained for the period 1961-2050.
Figure 9 presents the projected anomalies of rainfall and runoff for the time period 2021-2050, related to the mean of the periods 1961-1990 and 1971-2000.The difference between the means from both time periods was estimated.The t-test was used to determine if the two means were significantly different from each other.In Table 2, the p-value of t-test is presented in parenthesis, identifying the significant differences (denoted by * in Table 2).If the mean annual runoff for the future (2021-2050) is compared with the historical dataset for 1961-1990 (Figure 9a), a plausible significant decrease (−20%) is expected (Table 2).However, from Figure 9a, the decrease of runoff is more clearly identified in the first part of the period.If the comparison is based on the difference between 2021-2050 and 1971-2000 period (Figure 9b), the expected trend of runoff would be positive (2.5% from Table 2).Analyzing the other variables in the hydrological cycle, the reduction of runoff for 2021-2050 (Table 2) in relation with the baseline period ) is justified in the reduction of rainfall and the increase of PET (≈4%).However, for 2021-2050 in relation to 1971-2000, an increase in rainfall over 11% and a similar increase of PET (≈3.7%) are expected.The variability of the simulated seasonal mean rainfall, PET, mean temperature and runoff for the Fuensanta Basin, in two observed (1961-1990 and 1971-2000) and one future (2021-2050) time periods, were comparatively analyzed (Figure 10).If the results for 2021-2050 are contrasted with the corresponding results for the observational periods 1961-1990 and 1971-2000, a lower range of variability of the mean of rainfall, PET, and temperature for the future, are expected.Nevertheless, in the case of future runoff, this trend is not detected.In Figure 10a, increases in the median rainfall (although not in the maximum value) are expected for all seasons in the period 2021-2050, in contrast with 1961-1990.However, these increases are greater in the case of contrast with the observed period 1971-2000.In the case of PET (Figure 10b), clear increases in the minimum, maximum, first, third quartile, and median are only expected in summer in 2021-2050 with respect  The variability of the simulated seasonal mean rainfall, PET, mean temperature and runoff for the Fuensanta Basin, in two observed (1961-1990 and 1971-2000) and one future (2021-2050) time periods, were comparatively analyzed (Figure 10).If the results for 2021-2050 are contrasted with the corresponding results for the observational periods 1961-1990 and 1971-2000, a lower range of variability of the mean of rainfall, PET, and temperature for the future, are expected.Nevertheless, in the case of future runoff, this trend is not detected.In Figure 10a, increases in the median rainfall (although not in the maximum value) are expected for all seasons in the period 2021-2050, in contrast with 1961-1990.However, these increases are greater in the case of contrast with the observed period 1971-2000.In the case of PET (Figure 10b), clear increases in the minimum, maximum, first, third quartile, and median are only expected in summer in 2021-2050 with respect to both observational periods.These results match the summer boxplot of the mean temperature (Figure 10c).In Figure 10d, contrasting the future period with the observational periods, a significant decrease in the runoff in autumn is identified.However, if the future scenario of runoff is contrasted with the 1971-2000 period, a clear increase in the median in winter, spring and summer is predicted.
to both observational periods.These results match the summer boxplot of the mean temperature (Figure 10c).In Figure 10d, contrasting the future period with the observational periods, a significant decrease in the runoff in autumn is identified.However, if the future scenario of runoff is contrasted with the 1971-2000 period, a clear increase in the median in winter, spring and summer is predicted.11a), important decreases of rainfall are identified because this historical period was wetter.Consequently, the decreases in runoff (Figure 11e) are higher and cover a wider area, in comparison with the contrast results from the time period 1971-2000 (Figure 11f).In this last case, the runoff decreases are only identified in the lower part of the basin.
The maps of PET projections relative to 1961-1990 (Figure 11c) and 1971-2000 (Figure 11d) historical periods are spatially consistent.A general increase of PET (4% on average) in the majority of the basin and a decrease in some small areas from 0 to 5% are expected.In the case of contrast with the period 1961-1990 (Figure 11a), important decreases of rainfall are identified because this historical period was wetter.Consequently, the decreases in runoff (Figure 11e) are higher and cover a wider area, in comparison with the contrast results from the time period 1971-2000 (Figure 11f).In this last case, the runoff decreases are only identified in the lower part of the basin.
The maps of PET projections relative to 1961-1990 (Figure 11c) and 1971-2000 (Figure 11d) historical periods are spatially consistent.A general increase of PET (4% on average) in the majority of the basin and a decrease in some small areas from 0 to 5% are expected.In summary, increases in mean temperature and PET are expected for the 2021-2050 period in contrast to the 1961-1990 historical dataset, especially in summer and autumn.Variations in precipitation are expected according to observed period and season.Finally, reductions in runoff, especially in autumn, are expected in the future with respect to both observed periods.However, with respect to the 1971-2000 period, increases in runoff for winter, summer and spring are projected.

Discussion and Conclusions
In this study, the average runoff decreases significantly for the 2021-2050 period in comparison with 1961-1990.The mean annual runoff reduction (20%) is mainly explained by the significant increase of mean annual evapotranspiration (4%), rather than the reduction of mean annual precipitation.However, if the 1971-2000 period is selected for comparison, a divergent trend (positive trend) is detected, although without significant differences.This fact demonstrates the high sensitivity of the results to the time period selected for comparison, where specific events could influence the results and conclusions.As an example, for 1971-2000, although two periods of severe drought (1980-1982 and 1993-1995) were observed, a strong El Niño occurred at the end of the period 1971-2000.
Therefore, in concordance with [37], trends based on short records are very sensitive to the beginning and end dates and do not in general reflect long-term climate trends.Moreover, due to the natural variability, special attention should be paid to the selection of the period for analysis of the impacts It should be highlighted that the negative impacts of climate change on water cycle will be exacerbated in semi-arid basins with water shortages, which will have important implications for water management.By building RCM ensembles of meteorological data that are used as climate forcing for a parsimonious hydrological model, the reliability of hydrological projections in a semi-arid basin is increased.
Summarizing, from the results the following conclusions can be drawn: • A new methodology for reducing the uncertainties involved in the modeling chain, from the RCMs to hydrological models, is presented.

•
The seasonal and annual variability of climate variables considered in the RCM ensembles are well captured by the proposed novel approach.

•
Taking into account the built RCM ensembles of rainfall and PET, different trends in runoff are detected depending on the comparison period.

•
Therefore, the hydroclimatic trends are very sensitive to small changes in the climate drivers and the selection of the baseline period.

Figure 2 .
Figure 2. Location of the Segura River Basin in the Iberian Peninsula, and the Fuensanta Basin with the selected sites (731 and 759) for the analyses.

Figure 2 .
Figure 2. Location of the Segura River Basin in the Iberian Peninsula, and the Fuensanta Basin with the selected sites (731 and 759) for the analyses.

Figure 5 .
Figure 5. Seasonal and yearly rainfall Cumulative Density Functions (CDFs) on site 759, located in the Fuensanta Basin, from observational dataset (in black), RCMs (in color) and RCMs ensembles (black dashed black), for the 1961-1990 period.

Figure 5 .
Figure 5. Seasonal and yearly rainfall Cumulative Density Functions (CDFs) on site 759, located in the Fuensanta Basin, from observational dataset (in black), RCMs (in color) and RCMs ensembles (black dashed black), for the 1961-1990 period.

Figure 5 .
Figure 5. Seasonal and yearly rainfall Cumulative Density Functions (CDFs) on site 759, located in the Fuensanta Basin, from observational dataset (in black), RCMs (in color) and RCMs ensembles (black dashed black), for the 1961-1990 period.

Figure 6 .
Figure 6.Seasonal and yearly maximum temperature CDFs on site 759, located in the Fuensanta Basin, from observational dataset (in black), RCMs (in color) and RCM ensembles (black dashed black) for the 1961-1990 period.

Figure 7 .
Figure 7. Seasonal and yearly minimum temperature CDFs on site 759, located in the Fuensanta Basin, from observational dataset (in black), RCMs (in color) and RCM ensembles (black dashed black), for the 1961-1990 period.

Figure 6 .
Figure 6.Seasonal and yearly maximum temperature CDFs on site 759, located in the Fuensanta Basin, from observational dataset (in black), RCMs (in color) and RCM ensembles (black dashed black) for the 1961-1990 period.

Figure 6 .
Figure 6.Seasonal and yearly maximum temperature CDFs on site 759, located in the Fuensanta Basin, from observational dataset (in black), RCMs (in color) and RCM ensembles (black dashed black) for the 1961-1990 period.

Figure 7 .
Figure 7. Seasonal and yearly minimum temperature CDFs on site 759, located in the Fuensanta Basin, from observational dataset (in black), RCMs (in color) and RCM ensembles (black dashed black), for the 1961-1990 period.

Figure 7 .
Figure 7. Seasonal and yearly minimum temperature CDFs on site 759, located in the Fuensanta Basin, from observational dataset (in black), RCMs (in color) and RCM ensembles (black dashed black), for the 1961-1990 period.

Figure 9 .
Figure 9. Projected anomalies 2021-2050 of rainfall and runoff related to the mean for: (a) the 1961-1990 and (b) the 1971-2000 period.

Figure 9 .
Figure 9. Projected anomalies 2021-2050 of rainfall and runoff related to the mean for: (a) the 1961-1990 and (b) the 1971-2000 period.

Figure 11
Figure 11 presents the spatial distributions of predicted changes in annual rainfall, potential evapotranspiration and runoff for 2021-2050 period, in relation to 1961-1990 and 1971-2000 historical periods.In the case of contrast with the period 1961-1990 (Figure11a), important decreases of rainfall are identified because this historical period was wetter.Consequently, the decreases in runoff (Figure11e) are higher and cover a wider area, in comparison with the contrast results from the time period 1971-2000 (Figure11f).In this last case, the runoff decreases are only identified in the lower part of the basin.The maps of PET projections relative to 1961-1990 (Figure11c) and 1971-2000 (Figure11d) historical periods are spatially consistent.A general increase of PET (4% on average) in the majority of the basin and a decrease in some small areas from 0 to 5% are expected.

Figure 11
Figure 11 presents the spatial distributions of predicted changes in annual rainfall, potential evapotranspiration and runoff for 2021-2050 period, in relation to 1961-1990 and 1971-2000 historical periods.In the case of contrast with the period 1961-1990 (Figure11a), important decreases of rainfall are identified because this historical period was wetter.Consequently, the decreases in runoff (Figure11e) are higher and cover a wider area, in comparison with the contrast results from the time period 1971-2000 (Figure11f).In this last case, the runoff decreases are only identified in the lower part of the basin.The maps of PET projections relative to 1961-1990 (Figure11c) and 1971-2000 (Figure11d) historical periods are spatially consistent.A general increase of PET (4% on average) in the majority of the basin and a decrease in some small areas from 0 to 5% are expected.

Table 1 .
Selected the regional climate models (RCMs) from ENSEMBLES project.

Table 2 .
Contrast of mean annual runoff, rainfall and PET for Fuensanta Basin.

Table 2 .
Contrast of mean annual runoff, rainfall and PET for Fuensanta Basin.