CLIMATE CHANGE PROJECTIONS FOR WINTER STREAMFLOW IN DOURO RIVER

Climate change projections for the winter streamflow of the Douro River have been obtained for the period 2071-2099, using the Principal Component Regression (PCR) method. The winter streamflow time series (January to March) from eight stations distributed over the basin, covering the period 1950-2011, were used as predictand variables, while the principal components (PCs) of the winter (December to February) anomalies of sea level pressure (SLP) were used as predictors of the streamflow for the development of a statistical downscaling model. The period 19501995 was used for the calibration of the regression model, while 1996-2011 was used as validation period. The statistical downscaling model fitted from the observational SLP data was applied to the SLP outputs of three GCMs for the period 2071-2099, under the climate change scenarios RCP2.6, RCP4.5 and RCP8.5. The main result obtained is that all models and scenarios project a generalized decrease in the winter streamflow of the Douro River.


INTRODUCTION
Climate change has a direct and important impact on water resources.The effects on these resources will manifest not only in the variation of the quantity but also in the alteration of the quality and temporal distribution of the water (Moreno et al., 2005).
Although the flow and storage of water in the Earth's climate system are highly variable, it is clear that other changes apart from those due to natural variability are expected by the end of the current century (IPCC, 2013).The extent to which climate change affects water resources is evaluated through the changes observed in different variables (González-Zeas et al., 2010).The best example could be the proved increase in temperatures (Dasari et al., 2014), what will cause a net increase in rainfall, surface evaporation and sea level (Chang et al., 2015).
Spain is a country with a huge variety of water uses.Mainly its developed agriculture and its high hydroelectric potential have led to the development of a large amount of reservoirs, dams and underground water intakes (Moreno et al., 2005).And it is, by far, the agricultural activity the largest water consumer, representing on average 40% of the total water consumption in Europe, although in Mediterranean countries, like Spain, it represents 60-80% (Rasilla et al., 2013).Unfortunately, Spanish water resources have a high spatio-temporal irregularity in natural regime, especially when it is compared with the European average.For this reason, it has been necessary a marked human intervention in the hydrological cycle throughout the time, having profoundly altered its natural characteristics (Garrote de Marcos et al., 2008).
According to Falkenmark and Lindh (1976) a consumption of 20% of the total renewable water resources is considered as the limit of overexploitation of a system.In most of the basins of Spain water resources are overused.As shown in the Libro Blanco del Agua in Spain (2000), the Douro River is dangerously in the limit of overexploitation, with 14.175 m 3 of input and 2.929 hm 3 of intake, which makes a relationship of 21% of consumption.The precipitation regime in this basin, located over the north plateau of the Iberian Peninsula (IP), is characterized by a marked annual cycle, with rainy winters and very dry summer.
Natural water resources in the IP are used intensively, especially in agriculture, but also in other sectors, such as hydropower industry or tourism.Together with the increase in water demand, the observed decrease in water availability has intensified the situation of water stress in Spain.A country of this nature is very sensitive to declines that may present the inherent water resources to climate change, and, as mentioned before, adapting to the coming water stress conditions will demand economic, social and technical measures beyond the frontiers of each country.
In view of the above, the main objective of this work is to obtain climate change projections of the Douro River winter streamflow for the period 2071-2099.

DATA
The streamflow data base used in this work has been provided by the Center for Studies and Experimentation of Public Works (CEDEX).Time series from gauging stations and reservoirs with less than 5% of missing data for the period 1950-2011 have been considered.So, we have selected just eight stations whose gaps were filled by regression with well correlated neighboring stations.Table 1 shows a brief summary of the selected stations, while their localizations can be found in Figure 1.In the present study, we use the sea level pressure (SLP) averaged from December to February (DJF), as predictor variable for winter (JFM) Douro streamflow because the links between precipitation and atmospheric circulation tend to be the strongest in winter.Database used is the SLP monthly mean from the reanalysis data NCEP-NCAR.This SLP data set has a temporal coverage from 1950 to 2011, and a horizontal resolution of 2.5 o X 2.5 o .Also, the SLP outputs from three GCMs of the CMIP5 for the historic  and future periods (2071-2099) are used.The GCMs used are CESM1 (CAM5), IPSL-CM5A-MR, and MIROC5.All SLP data, reanalysis and GCMs, were chosen for a range covering the area from 20 o N to 90 o N and from 110 o O to 70 o E.
For computing the climate change projections of winter Douro streamflow, the RCPs scenarios 2.6, 4.5 and 8.5 have been chosen.

METHODOLOGY
There are many different statistical techniques that have been applied in order to downscale climate data.In this work we follow the methodological scheme based on the Principal Component Regression (PCR), as in Palomino-Lemus et al. (2015).PCR is a method that can be used to overcome the problem of multicollinearity on predictor variables, because a multiple regression model needs predictors uncorrelated or not multicollinear (Wigena and Djuraidah, 2014).In summary, the PCR method has two parts: 1) obtaining streamflow predictors through a Principal Component Analysis (PCA) applied to the SLP field, and 2) the construction of a multiple regression model that use the selected predictors to simulate de streamflow (named the downscaling model).Finally, streamflow projections for the period 2071-2099 by applying the obtained downscaling model to the SLP data derived from the GCM simulations under the three chosen RCPs scenarios (RCP2.6,RCP4.5 and RCP8.5) are computed.In other to correct the models bias, the Delta method (Hijmans et al., 2005;Palomino-Lemus et al., 2015) was used, so the projected changes are calculated as the difference between the modeled streamflow using historical and RCPs output for each model.The difference between the mean values for present and future downscaled streamflow was evaluated by the Wilcoxon-Mann-Whitney rank sum test.

Spatio-temporal variability of SLP
For starting, the stability of the predictor field, winter SLP, is analyzed computing a Principal Component Analysis (PCA), used to identify the main variability patterns of winter SLP.This PCA obtains patterns (EOFs) that facilitate the comparison of the climate variability reproduced by the reanalysis data (NCEP) and the models.For this analysis, the first six modes of winter SLP variability identified by the PCA of the reanalysis data in the period 1950-2011, which explain 85.7% of the total variance, were retained.Figure 2 shows these spatial patterns (EOFs).EOF1 explains 46.48% of the variance for the winter SLP field, with a positive correlated center located over Greenland and another negative correlated center over the Mediterranean and the central North Atlantic.This spatial pattern clearly corresponds to the NAO.The second EOF, which explains 11.90% of winter SLP variance, exhibits a positive center located in the North Atlantic.EOF3 (10.86% of variance) presents a tripolar structure with a high negative correlation center in northern Europe.EOF4 (9.93% of variance) also shows a tripolar structure with a pattern with high negative loadings centered on the British Isles.EOF5 and EOF6 (3.54% and 3.02%, respectively) account for only 14.43% of the SLP variance and show different action centers over the study region with lower correlations.
The associated PC time series (not shown) present for all the cases an important interannual variability.For the PC1 decadal variability is also apparent.To explore the physical meaning of these spatial modes, correlations between the PC series and several teleconnection indices for the period 1951-2010 have been computed, and they are shown in Table 2.As expected, the time series of the first EOF is strongly linked with the NAO (r = -0.55),and it exhibits a significant increase in the year 2010.PC1 is also linked to both the EA (r = -0.33)and the EA-WR indices (r = -0.27).Also, the third PC series is linked with the SCAND index (r = 0.26).Figure 3 shows the spatial patterns of the SLP variability (EOFs) from the models MIROC5, CESM1 (CAM5) and IPSL-CM5A-MR, for the present time .For the first EOF, all models show similar structures and have almost the same percentage of explained variance than the EOF1 of NCEP, except for the MIROC5 model, whose structure is mainly the same but spreads to the west, being its variance a little bit lower.EOF2 is characterized by a positive center located in the North Atlantic for the NCEP data, being a negative center for the CESM1 (CAM5) and IPSL-CM5A-MR models, which explain a bit more of variance.The MIROC5 represents a structure that is further from the EOF of the reanalysis data, although its variance is higher.For EOF3, all models show similar structures to the EOF3 of reanalysis data, with the same percentage of explained variance.Regards the fourth EOF, it is the CESM1 (CAM5) model which is the nearest to represent the EOF4 from the reanalysis data, being the MIROC5 model which presents a more differentiated pattern to the NCEP data.Again, the variance explained is very similar.Finally, EOF5 and EOF6 represent different structure between them and from the NCEP data, and their variance is almost similar.
In order to select adequate predictors for winter Douro streamflow that can be used in a regression model, correlations between streamflow time series and the six leading PCs of the SLP from reanalysis data, for the period 1950-2010, have been calculated.PCs that show significant correlations at 95% confidence level with the streamflow time series in the Douro River were selected as predictor variables for the regression model.This result is shown in Table 3.Note that the first two variability modes and the fifth and sixth (also the fourth but just for the gauging station number 2082) of the winter SLP are significantly associated with winter Douro streamflow.Thus, only these five modes were considered for the development of the regression models.

Principal component regression models for the streamflow
We used the principal component regression (PCR) method to build the forecasting model for the streamflow.For this analysis, we retained the first two variability modes of winter SLP field for all the stations, but also the fourth, fifth and sixth PCs appeared in some cases as significant.The training period 1950-1995 was used as calibration period, while the period 1996-2011 was used to verify the model.Table 4 presents a summary with the statistics obtained from the regression models, showing the significant correlations values, at 95% level, between the observed and predicted streamflow, and the root-mean-square error (RMSE), for both periods (calibration and validation).Table 4: Correlation (r) and root-mean-square error (RMSE) between the observed and predicted streamflow by the regression model for the calibration period, the validation period and the recalibration period.
In general, the regression models work properly, performing well during the training period, but a bit worse during the verification period, although it clearly fails when estimating peak values (not shown) as the high flows registered during the years 1979 and 2001, which correspond to very rainy years in the region.We may attribute to these years the principal cause of the high RMSE values obtained.
Once the downscaling model has been validated, we have recalibrated it using all data, that is, calibrating it again for the full period.The correlation coefficients between the observed and predicted values and their RMSE are also presented in Table 4.There is an improvement in the prediction withrespect to the verification period, although the model is not able to reproduce the extreme streamflow values of the last records.

Climate change projections
After having established the principal component regression model for the statistic downscaling of the streamflow, and in order to identify the potential impact of the climate change, we have applied it to the SLP data derived from the GCM simulations for both the present  and future (2071-2099) climate under the RCP2.6,RCP4.5 and RCP8.5 scenarios.
In general, for present climate (results not shown) we could say that CESM1 and IPSL-CM5A-MR models present a clear streamflow underestimation, whereas MI-ROC5 shows important overestimations.Furthermore, there are marked differences between the performance of the different models, being the IPSL-CM5A-MR which better represents the present conditions for most of the stations.Therefore, there are a considerable bias between the observational streamflow data and the present climate estimations from the models, which could be attributed to the tendency of the GCMs to show a more zonal pattern of the SLP than the NCEP data.This bias is taken into account when calculating the winter streamflow projections for the period 2071-2099.Thus, the streamflow differences have been obtained by subtracting the average model values in the present time to the projected future data, being the percentages relative to the present modeled average.The results obtained for the future are shown in Table 5.
The results of the projected changes for later this century show generalized decreases of the winter streamflow for all models and scenarios.As expected, the RCP2.6 scenario presents the lowest decreases, in contrast to RCP8.5, which shows the highest decreases in the winter streamflow for CESM1 and IPSL-CM5A-MR models, but being the RCP4.5 scenario which surprisingly projects the greatest changes for MIROC5.The stations 2002The stations , 2003The stations and 2029 suffer the biggest changes, due to the fact that are reservoirs and contain the biggest amount of water.
The p values obtained from the Wilcoxon-Mann-Whitney rank sum test is higher for the results obtained with the CESM1 and IPSL-CM5A-MR models, being the lowest values these which were obtained for the RCP8.5 scenario, so for this scenario the streamflow projected are significant for some locations.On the other hand, the MIROC5 model presents low p values in general, therefore, it shows the most significant changes.

DISCUSSION
Climate change projections for winter Douro River streamflow under the scenarios RCP2.6,RCP4.5 and RCP8.5, for the period 2071-2099, using the outputs of three GCMs (MIROC5, CESM1 and IPSL-CM5A-MR) have been computed by using the PCR method.The PCA applied to the winter SLP reanalysis data from the NCEP showed the existence of six leading variability modes, which explain 85.7% of the total variance.We analyzed the association between these resulting principal components and the streamflow series in order to use the SLP PC series as predictor variables in a multiple regression model.
The statistical models built for the eight stations showed in general a good representation both for the calibration and validation periods, although they clearly failed at estimating peak values which correspond to very rainy years.After this step, we recalibrated it again for the full period.Finally, we applied the model data to our regression model, obtaining the streamflow projections.
The general conclusions are the following: Table 5: Absolute (hm 3 ) and percentual (%) differences between observed and downscaled GCM values for the future period (2070-299), for each GCM and RCP.Significant changes at 95% condidence level are in bold.
• The NAO is the teleconnection pattern that mainly controls the winter precipitation variability in the west part of the IP, and therefore it is directly reflected in the winter seasonal flow of the Douro River.This pattern is mainly represented by the first SLP EOF.
• For the downscaling model, there is an improvement in the regression model during the verification period after having recalibrated it for the full period.However, the model fails at estimating peak values related to very rainy years.
• CESM1 and IPSL-CM5A-MR models present a clear underestimation of the present streamflow, whereas MIROC5 shows important overestimations.Thus, that bias has effects on the results of the future projections.
• The Global Circulation Models MIROC5, CESM1 and IPSL-CM5A-MR clearly show generalized decreases of the winter Douro River streamflow for the period 2071-2099, under the scenarios RCP2.6,RCP4.5 and RCP8.5, being for this scenario more significant, particularly for MIROC5 model.
It should be noted that the downscaling method used in this work can only represent the changes in the streamflow, which are principally derived from the changes in the rain, linked to changes in the atmospheric circulation as we only used one predictor field (SLP) in the downscaling model.In this sense, this work represents a first step to the analysis of the climate change impact on IP streamflow.Subsequent research is necessary because changes in the streamflow could be also linked to other predictor fields such as temperature.

Fig 1 :
Fig 1: Localization of the streamflow series analyzed in this work within the Spanish part of the Douro Basin.Gauging stations are shown in red and reservoirs in green.

Fig 2 :
Fig 2: Factor loading of the first six leading Empirical Orthogonal Factors of winter SLP over the study area for the reanalysis data.Variance explained by each mode is indicated in parentheses.

Fig 3 :
Fig 3: EOF loading factors of winter SLP for the three GCMs considered in the period 1951-2005.Variance explained by each mode is indicated in parentheses.

Table 1 :
Basic characteristics of the gauging stations (dark gray) and reservoirs (light gray)used in this work.The average annual streamflow was obtained from CEDEX.

Table 2 :
Correlations between the six leading PCs of the winter SLP reanalysis data and the teleconnection indices NAO, EA, SCAND and EA-WR, for the period 2051-2010.Significant correlations at 95% confidence level are shown in bold.

Table 3 :
Correlations between the PCs of winter SLP reanalysis data and the streamflow.Significant correlations at the 95% confidence level are shown in bold.