Hydroclimatic E ﬀ ects of a Hydropower Reservoir in a Tropical Hydrological Basin

: The consequent change in land cover from vegetation to water surface after inundation is the most obvious impact attributed to the impoundment of reservoirs and dam construction. However, river regulation also alters the magnitude and variability of water and energy ﬂuxes and local climatic parameters. Studies in Mediterranean, temperate and boreal hydrological basins, and even a global-scale study, have found a simultaneous decrease in the variation of runo ﬀ and increase in the mean evaporative ratio after impoundment. The aim here is to study the existence of these e ﬀ ects on a regulated tropical basin in Colombia with long-term data, as such studies in tropical regions are scarce. As expected, we observed a decrease in the long-term coe ﬃ cient of variation of runo ﬀ of 33% that can be attributed to the impoundment of the reservoir. However, we did not ﬁnd important changes in precipitation or the expected increasing evaporative ratio-e ﬀ ect from the impoundment of the reservoir, founding for the latter rather a decrease. This may be due to the humid conditions of the region where actual evapotranspiration is already close to its potential or to other land cover changes that decrease evapotranspiration during the studied period. Our study shows that the e ﬀ ects from impounded reservoirs in tropical regulated basins may di ﬀ er from those found in other climatic regions.


Introduction
The water cycle is not only a generator of ecosystem services and a critical agent of change in the Earth system; it is also a system that is subject to intense human interference. The resilience of water systems regulates the state of the ecosystems and biomes through interaction and feedback of the ecological functions related to water [1]. Out of the many important factors, climate variability and changes in land use/land cover (LULC) appear to be the main drivers of changes in these functions and in the availability of fresh water on land [2]. The storage and redistribution of water for agriculture and energy production, by irrigation and impoundment of water in artificial reservoirs is also known to induce important changes in the availability and variability of water [3][4][5][6] by increasing evaporation from the artificial reservoirs and regulating the flow of water [7][8][9]. These effects induce significant physical and ecological impacts on dependent terrestrial and aquatic ecosystems [10,11]. The effects of reservoirs have also been identified in variables of the hydrological cycle other than runoff (R). For example, Degu et al. [12] found that the presence of large reservoirs increases the extreme daily precipitation (P) (on 50th, 90th, 95th, and 99th percentiles) in the vicinity of the reservoirs in Mediterranean and arid climates. Hossain [13] found an increase in extreme precipitation after the hydroelectric generating capacity of 12.0 GW [23,24]. According to the Colombian Mining and Energy Planning Unit, the energetic demand of the country is projected to increase between 105 to 147% for the year 2050 in relation to that of 2010 [25]. Furthermore, 125 hydropower projects are in the pre-feasibility stage, that would add about 5600 MW to the existing installed capacity [26]. Since the hydroclimatic impacts of flow regulation from hydropower in Colombia and the corresponding effects on the climate, aquatic and terrestrial ecosystems have not been studied up to date, judging from the minimum amount of scientific publications in peer-reviewed journals, the construction of all these hydroelectric projects raises concern among environmental authorities, the communities in the areas of influence of such projects and its owners.
The Prado River is located in the central mountainous of Colombia. The Prado River originates in the Eastern Cordillera at an altitude of 2100 m above sea level (masl), and it flows from east to west. The longitude of the main river channel is approximately 96 km, presenting an average slope of 2% and discharging into the Magdalena River at an approximate elevation of 300 masl (Figure 1a). The basin has a mean monthly temperature (T) of 23.4 • C, oscillating between 28.1 and 17.5 • C, and P of 1895 mm year −1 . Most of the monthly runoff is generated in the upstream subbasins of the Negro and Cunday Rivers. The Prado Hydroelectric Project has the fourth largest reservoir in Colombia, 42 km 2 ( Figure 1a). The construction of the hydroelectric project started in 1961 and was followed by the filling of the reservoir and beginning of operation in 1973. The reservoir is multiuse, providing services for electricity, irrigation, aquaculture, fishing, recreation, and tourism, and is fed by the Cunday and Negro Rivers, as well as other smaller tributaries. The hydroelectric project comprises a large dam of 56 m in height and an associated water impounded reservoir that covers 2.6% of the area of the upstream basin. The Prado Hydroelectric Project has a volume of 1100 × 10 6 m 3 used to generate electric power with a nominal generation capacity of 60 MW that feeds into the Colombian National Electricity Grid. The reservoir also supplies water to the 2634-ha Prado River Irrigation and Drainage District located downstream of the dam. The area irrigated by the project consists mainly of rice (87%) and corn, grass for fodder and fruit trees. In addition, the reservoir supplies drinking water to 40,000 people in the Central section of the Colombian Andes [27]. According to CORINE Land Cover Colombia [28], the basin is covered in 50% by grasslands, 23% forests and 12% shrubs and bushes ( Figure 1b). The Cunday River sub-basin is mainly covered by grasslands (68%), shrubs and bushes (18%), and crops (4%) and the San Pablo River sub-basin by forests (60%), grasslands (29%) and shrubs and bushes (7%).
Sustainability 2020, 12, x FOR PEER REVIEW 3 of 19 the year 2050 in relation to that of 2010 [25]. Furthermore, 125 hydropower projects are in the prefeasibility stage, that would add about 5600 MW to the existing installed capacity [26]. Since the hydroclimatic impacts of flow regulation from hydropower in Colombia and the corresponding effects on the climate, aquatic and terrestrial ecosystems have not been studied up to date, judging from the minimum amount of scientific publications in peer-reviewed journals, the construction of all these hydroelectric projects raises concern among environmental authorities, the communities in the areas of influence of such projects and its owners. The Prado River is located in the central mountainous of Colombia. The Prado River originates in the Eastern Cordillera at an altitude of 2100 m above sea level (masl), and it flows from east to west. The longitude of the main river channel is approximately 96 km, presenting an average slope of 2% and discharging into the Magdalena River at an approximate elevation of 300 masl (Figure 1a). The basin has a mean monthly temperature (T) of 23.4 °C, oscillating between 28.1 and 17.5 °C, and P of 1895 mm year −1 . Most of the monthly runoff is generated in the upstream subbasins of the Negro and Cunday Rivers. The Prado Hydroelectric Project has the fourth largest reservoir in Colombia, 42 km 2 (Figure 1a). The construction of the hydroelectric project started in 1961 and was followed by the filling of the reservoir and beginning of operation in 1973. The reservoir is multiuse, providing services for electricity, irrigation, aquaculture, fishing, recreation, and tourism, and is fed by the Cunday and Negro Rivers, as well as other smaller tributaries. The hydroelectric project comprises a large dam of 56 m in height and an associated water impounded reservoir that covers 2.6% of the area of the upstream basin. The Prado Hydroelectric Project has a volume of 1100 × 10 6 m 3 used to generate electric power with a nominal generation capacity of 60 MW that feeds into the Colombian National Electricity Grid. The reservoir also supplies water to the 2634-ha Prado River Irrigation and Drainage District located downstream of the dam. The area irrigated by the project consists mainly of rice (87%) and corn, grass for fodder and fruit trees. In addition, the reservoir supplies drinking water to 40,000 people in the Central section of the Colombian Andes [27]. According to CORINE Land Cover Colombia [28], the basin is covered in 50% by grasslands, 23% forests and 12% shrubs and bushes (Figure 1b). The Cunday River sub-basin is mainly covered by grasslands (68%), shrubs and bushes (18%), and crops (4%) and the San Pablo River sub-basin by forests (60%), grasslands (29%) and shrubs and bushes (7%).
(a) (b) Figure 1. (a) Location of the Prado River Basin, its sub-basins, the reservoir, and precipitation, streamflow and meteorological stations (T: Temperature) (background map sourced from OpenStreetMap contributors [29]. (b) Percentage breakdown by land cover in the Prado River Basin and its sub-basins, based on Colombia's CORINE land cover [30].

Figure 1.
(a) Location of the Prado River Basin, its sub-basins, the reservoir, and precipitation, streamflow and meteorological stations (T: Temperature) (background map sourced from OpenStreetMap contributors [29]. (b) Percentage breakdown by land cover in the Prado River Basin and its sub-basins, based on Colombia's CORINE land cover [30].

Hydroclimatic Datasets and Processing
We collected data on discharge, precipitation and temperature in the region. We used the three discharge monitoring stations that create the Prado River Basin and its two unregulated upstream subbasins: Boquerón (main Prado Basin) and subbasins Cunday, and San Pablo, respectively, as shown in Figure 1a. The first one is located immediately downstream of the location of the dam with upstream drainage area of 1630 km 2 , and the second and third stations have upstream drainage areas of 131 and 174 km 2 , respectively. We converted the series from streamflow in m 3 /s to runoff (R) in mm/month or/year by dividing it by the drainage basin area. We used the calendar year for our analysis as we consider it convenient due to the tropical location of the hydrological basin. The period of runoff data availability differed among the three stations: Boquerón , Cunday (1954Cunday ( -2014 and San Pablo . We delineated the basin and its two sub-basins with the digital elevation model (3-arc second resolution) obtained from the EarthExplorer platform, which is hosted by the U.S. Geological Survey [31]. The drainage network shape file and the location of the three discharge monitoring stations were obtained from the web platform DHIME (in Spanish, Datos HIdrológicos y MEteorológicos) of the Colombian Institute of Hydrology, Meteorology and Environmental Studies (in Spanish, Instituto de Hidrología, Metereología y Estudios Ambientales -IDEAM) [32].
Out of the 41 meteorological stations available within and in the surroundings of the Prado River Basin, we selected the eight stations with the most complete monthly precipitation (P) data, that is, those stations that had less than 10% of monthly missing data during the period of research (i.e., 1954-2014). This information was interpolated spatially (cell size 300 × 300 m) using the inverse distance weighting interpolation method (IDW) with second power [33]. The spatially-interpolated and weighted monthly data of P was then used to calculate a monthly P average of all cells within each sub-basin of the Prado drainage basin.
The monthly T data from seven stations were obtained from IDEAM [32]. The stations had on average 33 years of monthly data with less than 5% of the data missing. One station located outside of the Prado River Basin but close to it (i.e., Guamo-21185030 with an altitude of 360 masl) had data throughout the period 1957 to 2013, with not more than 3% of missing data (Figure 1a). The missing values from 1954 to 1956 and 2014 in the time series of T in Guamo station were filled with T data from NASA Global Land Data Assimilation System (GLDAS) [34]. The GLDAS-2.0 dataset has two components: one forced entirely with the Princeton meteorological forcing data (hereafter, GLDAS-2.0) [35], and the other forced with a combination of model and observation forcing data sets (hereafter, GLDAS-2.1). The daily products of GLDAS have a 0.25 degree resolution and currently cover the period from January 1948 to December 2014 [36]. The results of this process were used to calculate T at the mean elevation of the Prado River Basin and its sub-basins by means of the thermal gradient (6.1 • C/km; [37]).

Post Processing Methods: Separation of Effects on Hydroclimatic Variables
There are several driver-effect separation techniques in the literature (see review in Jaramillo [38]). In this study, we combined two simple frameworks for separating the effects of landscape and climatic drivers on evapotranspiration [39,40] and the variability of runoff [41]. We basically compared these hydrological parameters before and after the construction of the dam, and between the regulated Prado river downstream and the unregulated upstream sub-basins Cunday and San Pablo, and complemented these analyses with statistical significance tests ( Table 1). The overarching framework of this study is shown in Figure 2. Furthermore, we compared the mean of the hydrological variables for the pre-dam (1954)(1955)(1956)(1957)(1958)(1959)(1960)(1961)(1962)(1963)(1964)(1965)(1966)(1967)(1968)(1969)(1970)(1971)(1972) and post-dam (1973-2014) periods by means of the non-parametric Wilcoxon's test [42]. If the resulting p-value of the Wilcoxon's test was less than the 0.05 significance level, the null hypothesis H 0 (i.e., nonidentical populations) was rejected, concluding that there was a statistically significant difference between the monthly medians in the populations before and after the construction of the dam (H 1 ).

Statistical Hypothesis Test Variables Purpose
Wilcoxon P and R Test the changes in the inter-annual monthly variability (i.e., for January, February… December), before and after the construction of the dam (P1 and P2) CVR and CVP Test any significant changes in the medians of the periods before and after the construction of the dam (P1 and P2) AETwb/P and AETclim/P

Before and After Changes in Precipitation and Runoff
In order to identify any potential effect of the reservoir on regional and local P, we studied the evolution of monthly P in the three rain gauge stations that were closest to the stream gauge stations in each sub-basin and located upstream of the dam (i.e., Aco (Prado Basin), Cunday (Cunday subbasin) and Villarrica (San Pablo sub-basin)). We assessed changes between the periods with available hydroclimatic data before construction (Period 1-P1) 1954-1972 and after construction (Period 2-P2) 1973-2014, where the value for each period was calculated as the arithmetic mean of all annual (or monthly) values with complete data within each period (See Equation (1) for the case of precipitation). We performed the same procedure for multiannual monthly and annual runoff averages.
Separating the effects of an impounded reservoir on precipitation can only be done with complex modelling such as that of moisture tracking, either by the Lagrangian (e.g., [43]) or Eulerian approaches (e.g., [44,45]) or by a spatial sample with several reservoirs and stations. Since the scarcity of data in the region does not allow for this analysis, we limited ourselves to the comparison of monthly total precipitation before and after the impoundment of the reservoir and the multicomparison of these differences across the regulated basin and its subbasins. Methodological framework. We used this methodology to separate or identify any potential effect of the impoundment of the reservoir on the hydrological fluxes or their temporal variability. The time series of the different hydroclimatic variables were split in two periods, P 1 (before dam) and P 2 (after dam). Table 1. Statistical tests and hydroclimatic variables to which they were applied, to detect any effect of the impoundment of the Prado reservoir.

Statistical Hypothesis Test Variables Purpose
Wilcoxon P and R Test the changes in the inter-annual monthly variability (i.e., for January, February . . . December), before and after the construction of the dam (P 1 and P 2 ) CV R and CV P Test any significant changes in the medians of the periods before and after the construction of the dam (P 1 and P 2 ) AET wb /P and AET clim /P

Before and After Changes in Precipitation and Runoff
In order to identify any potential effect of the reservoir on regional and local P, we studied the evolution of monthly P in the three rain gauge stations that were closest to the stream gauge stations in each sub-basin and located upstream of the dam (i.e., Aco (Prado Basin), Cunday (Cunday sub-basin) and Villarrica (San Pablo sub-basin)). We assessed changes between the periods with available hydroclimatic data before construction (Period 1-P 1 ) 1954-1972 and after construction (Period 2-P 2 ) 1973-2014, where the value for each period was calculated as the arithmetic mean of all annual (or monthly) values with complete data within each period (See Equation (1) for the case of precipitation). We performed the same procedure for multiannual monthly and annual runoff averages.
Separating the effects of an impounded reservoir on precipitation can only be done with complex modelling such as that of moisture tracking, either by the Lagrangian (e.g., [43]) or Eulerian approaches (e.g., [44,45]) or by a spatial sample with several reservoirs and stations. Since the scarcity of data in the region does not allow for this analysis, we limited ourselves to the comparison of monthly total precipitation before and after the impoundment of the reservoir and the multi-comparison of these differences across the regulated basin and its subbasins.

Driver-Effect Separation on Runoff Variability
We studied the evolution of the inter-annual monthly variability of P and R throughout the available period by calculating the coefficient of variation of precipitation (CV P ) and runoff (CV R ), Sustainability 2020, 12, 6795 6 of 18 respectively, in order to detect any changes in intra-annual variability of these variables before and after the construction of the hydroelectric project. We assumed that the coefficient of variation of precipitation can be used as a proxy of the climatic or natural variability of runoff at the interannual scale, following Jaramillo and Destouni [41]. The effect of regulation on the variability of runoff can then be obtained by comparing: (1) CV R before and after the impoundment of the reservoir, and (2) the relationship between CV R and CV P before and after impoundment, and (3) 1 and 2 between the regulated hydrological basin and the two unregulated subbasins. The coefficient of variation is the standard deviation of the 12 monthly values divided by their annual mean. As such, a time series is more homogeneous for low values of the coefficient of variation. In order to reduce the high inter-annual variability of R and P for the time series analysis and identify more clearly any long-term changes, we applied a ten-year moving window filter to each of the time series. We further tested other calculations of CV R and CV P by using a twelve-month CV R and also a sixty-month CV R , and then moving the CV calculation each month.

Separation on Evapotranspiration
Regional evapotranspiration data can be taken from satellite observations, although the available periods are generally too short (i.e., since 2000) [46]. Otherwise, the water mass balance can be used to estimate the AET of a basin by using mean basin-scale P and R over long periods of time and assuming a negligible change in water storage within the basin over the period of analysis, providing an overview of the net effect of all drivers of evapotranspiration (Equation (2)) [47][48][49]. However, the water balance method by itself does not provide additional information about the drivers of change in evapotranspiration due to the large number of interacting drivers in a given basin [50]. We refer to this actual evapotranspiration calculation as AET wb (Equation (2)).
The Budyko empirical equation expresses the evaporative ratio (i.e., AET/P) in terms of potential evapotranspiration (PET) over precipitation, also named the aridity index (i.e., PET/P). It can be used to obtain a climatic estimate of the annual evapotranspiration, which in this case is solely dependent on PET and P calculations and data, respectively ( [51,52]) (Equation (3)). We refer to this second estimate ( [51]) of AET as the climate-driven AET, AET clim .
For the use of Equation (3), we initially used a PET raster product from the Climatic Research Unit (CRU-TS 4.01), at a spatial resolution of 0.5 • (UEA CRU et al., 2017). This last set of PET data is a variant of the Penman-Monteith FAO method (FAO: Food and Agriculture Organization) (Ekström et al. [53], based on Allen et al. [54]), which uses average temperature T min , T max , vapor pressure, cloud cover, and wind speed, based on fixed weather conditions [55]. We also estimated PET based on other temperature and radiation-dependent equations, for which information was available (Table 2; Hargreaves-Samani [56], Hamon [57], Blaney-Criddle [58] and Langbein [59]). For the case of the Hargreaves-Samani, radiation data were computed based on latitude and month. Table 2. Empirical methods used to estimate the PET in the Prado River Basin and its sub-basins. T mean : mean temperature ( • C); T min : minimum temperature ( • C); T max : maximum temperature ( • C); Ra: incident solar radiation above the atmosphere (mm d −1 ); N: theoretical maximum daily insolation (hours); p: is the ratio of total daytime hours out of total daytime hours of the year; k: monthly consumptive use coefficient (k = 1); Ta: mean annual temperature; Rn: net radiation (MJ m −2 d −1 ); G: ground heat flux (MJ m −2 d −1 ); T: average temperature at two meters high ( • C); U 2 : average wind velocity (or estimate based on U 10 ) at 2 m high (m s −1 ); U 10 : wind velocity at 10 m high (m s −1 ); (e a -e d ): vapor pressure deficit in relation to saturation (kPa); ∆: slope of the saturation vapor pressure curve (kPa • C −1 ); γ: psychometric constant (kPa • C −1 ); 900: coefficient for the reference crop (kJ −1 kg K d −1 ); 0.34: wind coefficient for the reference crop (s m −1 ).

PET Model Equations
Hargreaves-Samani Lastly, we quantified the change in the residual evaporative ratio between the two periods, ∆(AET/P) r , as shown in Equation (4), since other studies have found that a positive ∆(AET/P) r can be an indicator of the effect of water use within a given basin on the evaporative ratio, specifically by increasing evapotranspiration related to irrigation and/or impoundment of reservoirs [18,19,60]. Hence, ∆(AET/P) r represents the component of the change in the evaporative ratio that is not explained by variability in potential evapotranspiration and/or precipitation, and as such represents the change component of any other driver of change related to changes in the surface such as changes in vegetation, water use, or soil characteristics, that may change the rates of evapotranspiration in the basin. In this case, ∆(AET clim /P) and ∆(AET wb /P) were calculated in the same manner as in Equation (1).
The nature of the change on ∆(AET wb /P) and ∆(AET clim /P) is based on the parameters used for their calculation. A positive ∆(AET wb /P) can really relate to anything that increases AET and decreases R (Equation (2)), either climatic or stemming from changes in land cover, vegetation, land use, etc., or based in this framework, anything that affects ∆(AET wb /P) or ∆(AET/P) r . That is why it is important to also calculate ∆(AET clim /P), the latter relating only to the effect of climatic parameters on the evaporative ratio, more specifically P and PET. It is worth noting that the calculation of PET depends in turn on several other climatic parameters such as relative humidity, wind speed, temperature and radiation, depending on the model used (See Table 2). Hence, a positive ∆(AET clim /P) can relate to an increase in PET-by either an increase in wind speed, in temperature, in radiation, or a decrease in P (see the structure of Equation (3)).
The change in residual can explain the driving effect of human activities on evapotranspiration change, or in our case, change in the evaporative ratio of a basin. For instance, if the ∆(AET/P) r is positive, it implies that ∆(AET wb /P) is larger than ∆(AET clim /P), and as such, the driver of the residual must have increased AET, considering that P is the same for all estimates. Drivers that can then increase AET include a change in land cover related to vegetation (i.e., grassland to forest or non-irrigated to irrigated crops systems) or to water (i.e., change of grassland to water surface when constructing Sustainability 2020, 12, 6795 8 of 18 an impounded reservoir). Conversely, a negative ∆(AET/P) r can then point to decreasing AET from deforestation due to less biomass and vegetation coverage.

Detecting Changes in Precipitation and Runoff
The time series analysis shows an inter-annual variation pattern that is similar for P and R, both on basin (Prado measured at Boquerón station) and sub-basin scales (measured at Cunday and San Pablo stations) (Figure 3). Furthermore, a simultaneous decrease in P and R occurs after the time of commissioning of the project in 1973. If the continuous decrease in P and R would have occurred only in the Boquerón station and not in the two sub-basins located upstream, it could be indicative of the effects of the commissioning of the project. However, since the decrease is more regional than local, also occurring in the two sub-basins located upstream of the Prado reservoir, it is more likely that the decrease in P and a large part of the decrease in R are driven by climatic variability. During the studied period, the sub-basin of San Pablo reported the lowest annual mean of P (1676 mm year −1 ), in comparison to Cunday sub-basin and Prado Basin (1885 mm year −1 and 2021 mm year −1 , respectively). The mean annual runoff for the Cunday sub-basin (635 mm year −1 ) was also less than the value calculated for San Pablo (860 mm year −1 ), meaning that although it might rain more in Cunday than in San Pablo, the amount of precipitation that partitions into runoff is less. constructing an impounded reservoir). Conversely, a negative Δ(AET/P)r can then point to decreasing AET from deforestation due to less biomass and vegetation coverage.

Detecting Changes in Precipitation and Runoff
The time series analysis shows an inter-annual variation pattern that is similar for P and R, both on basin (Prado measured at Boquerón station) and sub-basin scales (measured at Cunday and San Pablo stations) (Figure 3). Furthermore, a simultaneous decrease in P and R occurs after the time of commissioning of the project in 1973. If the continuous decrease in P and R would have occurred only in the Boquerón station and not in the two sub-basins located upstream, it could be indicative of the effects of the commissioning of the project. However, since the decrease is more regional than local, also occurring in the two sub-basins located upstream of the Prado reservoir, it is more likely that the decrease in P and a large part of the decrease in R are driven by climatic variability. During the studied period, the sub-basin of San Pablo reported the lowest annual mean of P (1676 mm year −1 ), in comparison to Cunday sub-basin and Prado Basin (1885 mm year −1 and 2021 mm year −1 , respectively). The mean annual runoff for the Cunday sub-basin (635 mm year −1 ) was also less than the value calculated for San Pablo (860 mm year −1 ), meaning that although it might rain more in Cunday than in San Pablo, the amount of precipitation that partitions into runoff is less. Regarding the annual variability of precipitation, the coefficient of variation of precipitation (CVP) over Prado, Cunday and San Pablo basins followed the same temporal pattern; an increase around 1968 and a stabilization, with minor oscillations, after this period (Figure 4) without any clear difference among the three basins or indication of the effect of the reservoir in precipitation variability. Rather, the precipitation variability in the region is generally associated with the interactions between the intertropical confluence zone (i.e., between 4° and 7° latitude), orography and local circulation, which affect the volumes of precipitated water [61]. As opposed to CVP, the overall trend of CVR in Boquerón (Prado Basin) is markedly different from that of Cunday and San Pablo and quite different before and after the hydroelectric power plant began operations. In the period before 1973, all stations show a similar CVR, but after 1973, CVR in Boquerón (Prado Basin) decreases while that of San Pablo and Cunday sub-basins increase. In addition, CVR after 1973 oscillates in a lower level in Boquerón than in the two sub-basins, meaning that the differences in R between the twelve months have decreased by the flow regulation of the dam. This can be also seen Regarding the annual variability of precipitation, the coefficient of variation of precipitation (CV P ) over Prado, Cunday and San Pablo basins followed the same temporal pattern; an increase around 1968 and a stabilization, with minor oscillations, after this period (Figure 4) without any clear difference among the three basins or indication of the effect of the reservoir in precipitation variability. Rather, the precipitation variability in the region is generally associated with the interactions between the intertropical confluence zone (i.e., between 4 • and 7 • latitude), orography and local circulation, which affect the volumes of precipitated water [61]. As opposed to CV P , the overall trend of CV R in Boquerón (Prado Basin) is markedly different from that of Cunday and San Pablo and quite different before and after the hydroelectric power plant began operations. In the period before 1973, all stations show a similar CV R , but after 1973, CV R in Boquerón (Prado Basin) decreases while that of San Pablo and Cunday sub-basins increase. In addition, CV R after 1973 oscillates in a lower level in Boquerón than in the two sub-basins, meaning that the differences in R between the twelve months have decreased by the flow regulation of the dam. This can be also seen for the variability of total monthly runoff ( Figure 5). These results are similar regardless of the methods of computation of CV (See Supplementary Materials, Figures A1 and A2).      Comparison between the multi-annual monthly averages of P (first row) and R (second row) in P1 (gray) and P2 (blue and red) in the Prado River Basin (Boquerón, upper section) and its Cunday and San Pablo sub-basins (left to right, respectively). The confidence intervals were calculated using the non-parametric bootstrap bias-corrected and accelerated method [62] with 2000 bootstrap replicas for each of the months within the study periods. Diamond points represented any significant changes in the medians of the periods before and after the construction of the dam (P1 and P2) by means of the Wilcoxon test (p < 0.05).

Figure 5.
Comparison between the multi-annual monthly averages of P (first row) and R (second row) in P 1 (gray) and P 2 (blue and red) in the Prado River Basin (Boquerón, upper section) and its Cunday and San Pablo sub-basins (left to right, respectively). The confidence intervals were calculated using the non-parametric bootstrap bias-corrected and accelerated method [62] with 2000 bootstrap replicas for each of the months within the study periods. Diamond points represented any significant changes in the medians of the periods before and after the construction of the dam (P 1 and P 2 ) by means of the Wilcoxon test (p < 0.05).
To test whether the CV P and CV R were significantly similar in the basin and its sub-basins, the Wilcoxon contrast method was used (Table 3) between the periods before and after the construction of the dam. We found that the only significant difference (p < 0.05; Wilcoxon) among basins existed for the median values of CV R for the period after the construction of the dam between Cunday and Boquerón (Prado Basin) and San Pablo and Boquerón (Prado Basin) stations, pointing to an effect of flow regulation. For the case of CV P , significant differences (p < 0.05; Wilcoxon) were found between Prado and Cunday in the period after the construction of the dam. The significant differences of CV P between the two upstream sub-basins are possibly due to the difference in precipitation patterns among them according to their specific orographic characteristics. However, the magnitude and direction of the CV P 's from Prado and San Pablo were significantly different (p < 0.05). This result can be related to monthly variability of P in San Pablo has been different from values observed in Cunday before and after the construction of the dam. To validate the differences between the CV P 's of San Pablo with Prado and Cunday, the contrast tests of the MK, Sen's slope and Levene, were applied on a monthly scale to rule out possible effects of the relative variability of the CV. In general, we estimate that the impoundment of the reservoir has decreased CV R by 33%, after removing the change in natural variability as observed in the development of CV R from the period before the dam to the period after. Table 3. Results of the Wilcoxon significance test to determine the significance of the differences in the means of CV P and CV R between the Prado River Basin and its sub-basins, before and after the hydroelectric power plant began operations. * indicates a p-value less than the significance level (α) of 0.05. An analysis of the average monthly R and P in the Prado Basin and its sub-basins (Cunday and San Pablo) gives more insight into the reasons behind the changes of CV R and monthly R occurring in Boquerón station after the regulation of the Prado River in 1973 ( Figure 5). Firstly, after the construction of the dam and during the dry period stretching from June to September, the regulation of the dam produced more runoff in Boquerón (Prado basin) than the actual amount of precipitation falling within its upstream basin. This only occurs in the Boquerón (Prado basin) station and is not observed in the other two discharge stations corresponding to the sub-basins upstream. Secondly, while the seasonal pattern between monthly P and R is consistent before and after the construction of the dam in the upstream Cunday sub-basin (i.e., they are in phase), the seasonal pattern of R is consistent with that of P only before the construction of the dam in the Prado Basin. This explains how the Prado reservoir retains water during the wet season of March to May and October to December to maximize the reservoir capacity for hydropower generation, to later release it in the dry period from June to August. Furthermore, the monthly confidence intervals (CI) of P ( Figure 5) show a reduction in the magnitude and variance of this variable during P 2 (blue) with respect to P 1 (gray) during the wettest months (April and November) and in the dry month of June. Although the reduction in P in these months could be attributed to the construction of the dam, the fact that it occurs also over the upstream subbasins suggests that these decreases are more related to natural climatic variability.

Main Hydroclimatic Changes in the Basins
Now, when studying simultaneously the changes in the means of R, P, and AET between P 1 and P 2 for the Prado basin and the two sub-basins, we found that P had decreased in the three basins. In the sub-basins Cunday and San Pablo, the decrease in P partitioned mostly into a decrease in R than a decrease in AET wb (Figure 6). For the case of the large Prado basin including the reservoir, the decrease in P of 120 mm year −1 partitioned mostly into a decrease in AET wb of 150 mm year −1 and an increase in R of 30 mm year −1 . Contrary to basin-wise analysis in regulated basins that have found an increase in AET wb after regulation, the Prado regulated basin experienced a decrease in AET wb , highlighting the potential role of other drivers of AET change that tend to reduce evapotranspiration or the negligible effect of the reservoir on evapotranspiration at the basin scale (See Figure A4).

Main Hydroclimatic Changes in the Basins
Now, when studying simultaneously the changes in the means of R, P, and AET between P1 and P2 for the Prado basin and the two sub-basins, we found that P had decreased in the three basins. In the sub-basins Cunday and San Pablo, the decrease in P partitioned mostly into a decrease in R than a decrease in AETwb ( Figure 6). For the case of the large Prado basin including the reservoir, the decrease in P of 120 mm year −1 partitioned mostly into a decrease in AETwb of 150 mm year −1 and an increase in R of 30 mm year −1 . Contrary to basin-wise analysis in regulated basins that have found an increase in AETwb after regulation, the Prado regulated basin experienced a decrease in AETwb, highlighting the potential role of other drivers of AET change that tend to reduce evapotranspiration or the negligible effect of the reservoir on evapotranspiration at the basin scale (See Figure A4). To verify if these basin-scale mean annual changes may be related to the impoundment of the reservoir and not to precipitation and/or PET changes, we analyzed the time series of AETwb/P and AETclim/P and calculated the residual ∆(AET/P)r (Figure 7). The AETwb/P in Prado has been generally below the mean and CI values of AETclim/P in both periods. On the other hand, in the Cunday subbasin, the AETwb/P has been over the AETclim/P throughout the analysis period, but both the magnitude and the direction of the tendencies are significantly similar before and after 1973 (p < 0.05). Although the rates and variability of P in Cunday are similar to those observed in San Pablo and Prado, there is a greater evapotranspiration in Cunday (see Figure A3), possibly due to the higher mean temperature (i.e., 25 °C in Cunday and 21 °C in San Pablo).
The residual change ∆(AET/P)r gives information on the other possible governing drivers of evapotranspiration change in the basin, with the sign and the magnitude of the change providing information on the nature of the main drivers different from PET or P such as other climatic factors, land or water use [63]. As ∆(AET/P)r is negative in the Prado River basin, there appears to be no effect of the impoundment of the reservoir on evapotranspiration, or at least is shadowed by other effects. Furthermore, ∆(AET/P)r is positive in its upstream unregulated subbasin, highlighting even more the lack of the expected signal from the reservoir. To verify if these basin-scale mean annual changes may be related to the impoundment of the reservoir and not to precipitation and/or PET changes, we analyzed the time series of AET wb /P and AET clim /P and calculated the residual ∆(AET/P) r (Figure 7). The AET wb /P in Prado has been generally below the mean and CI values of AET clim /P in both periods. On the other hand, in the Cunday sub-basin, the AET wb /P has been over the AET clim /P throughout the analysis period, but both the magnitude and the direction of the tendencies are significantly similar before and after 1973 (p < 0.05). Although the rates and variability of P in Cunday are similar to those observed in San Pablo and Prado, there is a greater evapotranspiration in Cunday (see Figure A3), possibly due to the higher mean temperature (i.e., 25 • C in Cunday and 21 • C in San Pablo).
The residual change ∆(AET/P) r gives information on the other possible governing drivers of evapotranspiration change in the basin, with the sign and the magnitude of the change providing information on the nature of the main drivers different from PET or P such as other climatic factors, land or water use [63]. As ∆(AET/P) r is negative in the Prado River basin, there appears to be no effect of the impoundment of the reservoir on evapotranspiration, or at least is shadowed by other effects. Furthermore, ∆(AET/P) r is positive in its upstream unregulated subbasin, highlighting even more the lack of the expected signal from the reservoir. Sustainability 2020, 12, x FOR PEER REVIEW 12 of 19 Figure 7. Distribution of the changes in the evaporative ratio, Δ(AETwb/P) before and after reservoir impoundment and the climatic Δ(AETclim/P) and residual Δ(AET/P)r effects over the Prado River Basin and its sub-basins. The vertical lines (black color) are upper and lower limits of confidence intervals (CIs) from the bootstrapping method after accounting for various methods of calculation.

Discussion
Since the long-term changes in the variables here studied may be driven by a combination of inter-annual and intra-annual climate variability and other factors such as deforestation [64], reforestation [63], rainfed and irrigated agricultural development [65] and impoundment of water in reservoirs [66], it is difficult to attribute these changes to a sole driver. However, changes in the hydrologic variables before and after the construction of the dam could suggest or at least point to a possible influence of the project on the basin's long-term hydroclimate. Such methodology has been used to study the effects of forestry, rainfed agriculture and hydropower development on the hydroclimate of a given basin (e.g., [18,19,60,63]).
The current results show a decrease in CVR by 33% that is explained by the regulation of the reservoir. This effect was not observed in the sub-basins upstream from the project. Regarding the change in precipitation and runoff observed in the basin and its subbasins, if the decreasing trend of precipitation and runoff continue in the Prado River Basin in the future, the operation of the hydroelectric project might be affected when compared to its current conditions, although more information is needed regarding the operation of the dam that in this case is confidential and not distributed by the operating company.
Furthermore, the need to release runoff resulting from the wet season (i.e., March, April, October, and November) to keep the storage for flood control means that the vulnerability to water security as ecosystem service downstream during the dry season will increase. In order to compensate this imbalance between the availability of water during the wet and dry months, the dam administrators must retain more water than what is naturally available in the rivers during the wet season, and release a greater amount of water stored during the dry season, as evidenced by Ehsani et al. [67] in regulated dams in the northeastern United States.
Hydropower production developments have been reported to co-occur with a AETwb/P increase and a CVR decrease [18,19,60]. However, the novelty of our study in a tropical regulated basin relies in the fact that only changes in CVR were here found, and not the increase in evaporative ratio expected from regulated reservoirs as seen in other Mediterranean, boreal and temperate regulated Figure 7. Distribution of the changes in the evaporative ratio, ∆(AET wb /P) before and after reservoir impoundment and the climatic ∆(AET clim /P) and residual ∆(AET/P) r effects over the Prado River Basin and its sub-basins. The vertical lines (black color) are upper and lower limits of confidence intervals (CIs) from the bootstrapping method after accounting for various methods of calculation.

Discussion
Since the long-term changes in the variables here studied may be driven by a combination of inter-annual and intra-annual climate variability and other factors such as deforestation [64], reforestation [63], rainfed and irrigated agricultural development [65] and impoundment of water in reservoirs [66], it is difficult to attribute these changes to a sole driver. However, changes in the hydrologic variables before and after the construction of the dam could suggest or at least point to a possible influence of the project on the basin's long-term hydroclimate. Such methodology has been used to study the effects of forestry, rainfed agriculture and hydropower development on the hydroclimate of a given basin (e.g., [18,19,60,63]).
The current results show a decrease in CV R by 33% that is explained by the regulation of the reservoir. This effect was not observed in the sub-basins upstream from the project. Regarding the change in precipitation and runoff observed in the basin and its subbasins, if the decreasing trend of precipitation and runoff continue in the Prado River Basin in the future, the operation of the hydroelectric project might be affected when compared to its current conditions, although more information is needed regarding the operation of the dam that in this case is confidential and not distributed by the operating company.
Furthermore, the need to release runoff resulting from the wet season (i.e., March, April, October, and November) to keep the storage for flood control means that the vulnerability to water security as ecosystem service downstream during the dry season will increase. In order to compensate this imbalance between the availability of water during the wet and dry months, the dam administrators must retain more water than what is naturally available in the rivers during the wet season, and release a greater amount of water stored during the dry season, as evidenced by Ehsani et al. [67] in regulated dams in the northeastern United States.
Hydropower production developments have been reported to co-occur with a AET wb /P increase and a CV R decrease [18,19,60]. However, the novelty of our study in a tropical regulated basin relies in the fact that only changes in CV R were here found, and not the increase in evaporative ratio expected from regulated reservoirs as seen in other Mediterranean, boreal and temperate regulated basins. This may be explained by the fact that the construction of only one reservoir within the basin may not be sufficient to evidence an effect on net evapotranspiration from flow regulation from the entire basin. On the other hand, a physical interpretation of the finding is related to the sensible and latent heat fluxes for open water bodies and vegetated covers. In tropical and humid regions, forested areas exhibit comparable moisture fluxes due to transpiration than evaporation from open water bodies [12]. Therefore, flooding or clearing of forest by an artificial lake are unlikely to create a distinctly different local climate. However, in semi-arid, Tundra, humid continental and Mediterranean regions, several studies [4,18,19,40] have demonstrated that the reservoirs can to add sufficiently more moisture than the sparsely vegetated surroundings, resulting in spatial gradients of water vapor flux in these types of climate. Although the reservoir in the Prado River basin has the characteristics of a large dam, the monthly precipitation magnitude and variability in the Prado basin and its subbasins Cunday and San Pablo appear not to be affected in this case.
Another possible uncertainty might be related to the limited availability of climatic data on wind speed and radiation, variables needed to calculate a more precise value of PET. The lack of specific humidity, radiation, moisture and wind, data for the region and during the time of construction of the dam does not allow us to specify which may be the most relevant climatic parameter. Not accounting for these factors based on the lack of measured data may also be linked to the decrease in relative actual evapotranspiration. Other driver-effects occurring in the basin may also be reasons why an increase in (AET/P) r is not observed. Furthermore, since deforestation in the basin should decrease (AET/P) r as here found, [68], deforestation may also be a factor changing evapotranspiration and the evaporative ratio before and after impoundment.

Conclusions
For the regulated tropical river basin here studied, we found a decrease in the long-term coefficient of variation of runoff as seen in other regulated basins around the world. However, we did not find changes in precipitation (both monthly and annual scales) or the expected increasing evaporative ratio-effect (annual scale) from the impoundment of the reservoir, founding for the latter rather a decrease. This may be due to the humid conditions of the region where actual evapotranspiration is already close to its potential or to other land cover changes that decrease evapotranspiration during the studied period. Our study shows that the effects from impounded reservoirs in tropical regulated basins may differ from those found in other climatic regions.

Appendix A
Sustainability 2020, 12, x FOR PEER REVIEW 14 of 19 Appendix A Figure A1. Coefficient of variation of monthly runoff (CVR; red) and precipitation (CVP; blue) and moving window (MW) using twelve months to calculate CV for runoff (MV-CVR; orange) and precipitation (MV-CVP; darkblue) in the Prado River Basin (Boquerón) and its two sub-basins (Cunday and San Pablo). Figure A2. Coefficient of variation of monthly runoff (CVR; red) and precipitation (CVP; blue) and moving window (MW) using 60 months to calculate CV for runoff (MV-CVR; orange) and precipitation (MV-CVP; darkblue) in the Prado River Basin (Boquerón) and its two sub-basins (Cunday and San Pablo). Figure A1. Coefficient of variation of monthly runoff (CV R ; red) and precipitation (CV P ; blue) and moving window (MW) using twelve months to calculate CV for runoff (MV-CV R ; orange) and precipitation (MV-CV P ; darkblue) in the Prado River Basin (Boquerón) and its two sub-basins (Cunday and San Pablo). Appendix A Figure A1. Coefficient of variation of monthly runoff (CVR; red) and precipitation (CVP; blue) and moving window (MW) using twelve months to calculate CV for runoff (MV-CVR; orange) and precipitation (MV-CVP; darkblue) in the Prado River Basin (Boquerón) and its two sub-basins (Cunday and San Pablo). Figure A2. Coefficient of variation of monthly runoff (CVR; red) and precipitation (CVP; blue) and moving window (MW) using 60 months to calculate CV for runoff (MV-CVR; orange) and precipitation (MV-CVP; darkblue) in the Prado River Basin (Boquerón) and its two sub-basins (Cunday and San Pablo). Figure A2. Coefficient of variation of monthly runoff (CV R ; red) and precipitation (CV P ; blue) and moving window (MW) using 60 months to calculate CV for runoff (MV-CV R ; orange) and precipitation (MV-CV P ; darkblue) in the Prado River Basin (Boquerón) and its two sub-basins (Cunday and San Pablo). Figure A3. Annual potential evapotranspiration (PET) estimated by five different models and its 95% confidence intervals (CI) from 1954 to 2014. The CIs were calculated using the non-parametric bootstrap bias-corrected and accelerated [62] with 2000 bootstrap replicas for each of the months within the study periods. Figure A4. Changes in the water balance-based actual evapotranspiration (AETwb) and estimated evapotranspiration based on a yearly climate approach (AETclim) from 1954 to 2014, presented through a 10-year moving average in the Prado River Basin and its sub-basins. The gray dotted line shows the year in which the hydroelectric plant began operations (1973). The red lines are the CIs of AETclim and were calculated using the non-parametric bootstrap bias-corrected and accelerated [62] with 2000 bootstrap replicas for each of the months within the study periods.

Appendix B
We compared the PET from Langbein model with pan evaporation (PE) measures from the Guamo meteorological station (see Figure 1a). The PE records are available from 1965 to 2012 with 22% of missing data on a monthly scale. Therefore, we considered the 21 years without monthly missing data of PE and aggregate them into an annual scale and finally compared them with PET. The mean interannual values of PE and PET were 1714 and 1604 mm year −1 , respectively, and the standard deviation among these variables was around 155 mm year −1 . Figure A3. Annual potential evapotranspiration (PET) estimated by five different models and its 95% confidence intervals (CI) from 1954 to 2014. The CIs were calculated using the non-parametric bootstrap bias-corrected and accelerated [62] with 2000 bootstrap replicas for each of the months within the study periods.  Figure A3. Annual potential evapotranspiration (PET) estimated by five different models and its 95% confidence intervals (CI) from 1954 to 2014. The CIs were calculated using the non-parametric bootstrap bias-corrected and accelerated [62] with 2000 bootstrap replicas for each of the months within the study periods. Figure A4. Changes in the water balance-based actual evapotranspiration (AETwb) and estimated evapotranspiration based on a yearly climate approach (AETclim) from 1954 to 2014, presented through a 10-year moving average in the Prado River Basin and its sub-basins. The gray dotted line shows the year in which the hydroelectric plant began operations (1973). The red lines are the CIs of AETclim and were calculated using the non-parametric bootstrap bias-corrected and accelerated [62] with 2000 bootstrap replicas for each of the months within the study periods.

Appendix B
We compared the PET from Langbein model with pan evaporation (PE) measures from the Guamo meteorological station (see Figure 1a). The PE records are available from 1965 to 2012 with 22% of missing data on a monthly scale. Therefore, we considered the 21 years without monthly missing data of PE and aggregate them into an annual scale and finally compared them with PET. The mean interannual values of PE and PET were 1714 and 1604 mm year −1 , respectively, and the standard deviation among these variables was around 155 mm year −1 . Figure A4. Changes in the water balance-based actual evapotranspiration (AET wb ) and estimated evapotranspiration based on a yearly climate approach (AET clim ) from 1954 to 2014, presented through a 10-year moving average in the Prado River Basin and its sub-basins. The gray dotted line shows the year in which the hydroelectric plant began operations (1973). The red lines are the CIs of AET clim and were calculated using the non-parametric bootstrap bias-corrected and accelerated [62] with 2000 bootstrap replicas for each of the months within the study periods.

Appendix B
We compared the PET from Langbein model with pan evaporation (PE) measures from the Guamo meteorological station (see Figure 1a). The PE records are available from 1965 to 2012 with 22% of missing data on a monthly scale. Therefore, we considered the 21 years without monthly missing data of PE and aggregate them into an annual scale and finally compared them with PET. The mean interannual values of PE and PET were 1714 and 1604 mm year −1 , respectively, and the standard deviation among these variables was around 155 mm year −1 .