Impact of Rain Gauges Distribution on the Runoff Simulation of a Small Mountain Catchment in Southern Ecuador

In places with high spatiotemporal rainfall variability, such as mountain regions, input data could be a large source of uncertainty in hydrological modeling. Here we evaluate the impact of rainfall estimation on runoff modeling in a small páramo catchment located in the Zhurucay Ecohydrological Observatory (7.53 km2) in the Ecuadorian Andes, using a network of 12 rain gauges. First, the HBV-light semidistributed model was analyzed in order to select the best model structure to represent the observed runoff and its subflow components. Then, we developed six rainfall monitoring scenarios to evaluate the impact of spatial rainfall estimation in model performance and parameters. Finally, we explored how a model calibrated with far-from-perfect rainfall estimation would perform using new improved rainfall data. Results show that while all model structures were able to represent the overall runoff, the standard model structure outperformed the others for simulating subflow components. Model performance (NSeff) was improved by increasing the quality of spatial rainfall estimation from 0.31 to 0.80 and from 0.14 to 0.73 for calibration and validation period, respectively. Finally, improved rainfall data enhanced the runoff simulation from a model calibrated with scarce rainfall data (NSeff 0.14) from 0.49 to 0.60. These results confirm that in mountain regions model uncertainty is highly related to spatial rainfall and, therefore, to the number and location of rain gauges.


Introduction
Páramo is a high-elevation tropical Andean ecosystem located in the upper belt of the northern Andes approximately from 3500 to 5000 m a.s.l. [1,2]. In southern Ecuador, it is characterized by the presence of tussock grasses, wetlands, and scarce patches of Polylepis sp. [3][4][5]. Like other mountain ecosystems worldwide recognized as water suppliers for downstream populations [6], the Andean páramo is the most important water source for Andean cities such as Quito, Bogota, Mérida, and Cuenca, mainly due to the high water retention capacity of its soils and the constant precipitation it receives throughout the year [7,8]. In Ecuador, it also provides water for some of the most important hydropower projects such as Paute Integral (2353 MW) and Coca-Codo Sinclair (1500 MW) and irrigation projects in the inter-Andean valley [9]. However, human activities such as grazing, cultivation, and pine plantations can alter their normal hydrological regulation [10,11]. Therefore, a good understanding of páramo hydrology is critical for present and future water resources development [12]. rain gauge station outside or in the outlet of the catchment area. Models are then calibrated using this configuration, evidently leading to low simulation efficiencies. However, as the catchment becomes more important (e.g., for water resources development), additional rain gauges are installed within the catchment. But, modelers will have to wait several years under this new monitoring system to recalibrate the model. This raises the question of how a model calibrated with a far-from-perfect rainfall estimation will perform using new improved rainfall data, i.e., without undergoing a recalibration process. Thus, we sought to answer this question.

Study Area
This study was carried out in the Zhurucay Ecohydrological Observatory located in the southern Ecuadorian Andes, near the continental divide, draining into the Pacific Ocean. The observatory consists of a nested catchment of 7.53 km 2 with elevation ranging from 3400 to 3900 m a.s.l. ( Figure  1). Mean annual precipitation is 1345 mm at 3780 m a.s.l. with weak seasonality. Rain mainly falls as drizzle and occurs almost daily [17]. Mean annual temperature is 6.0 °C [14]. The geomorphology of the catchment is glacial U-shaped. The average slope is 10°, although slopes up to 73° are found. The geology is compacted volcanic rock deposits formed during the last ice age [60]. In the northeastern of the catchment there is a ponded wetland at a flat hilltop. This structure drains towards the outlet of subcatchment S7 and occupies almost the entire subcatchment [21]. The main soil type in Zhurucay are Andosols (Ah horizon), covering 80% of the area and commonly located on the hillslopes. Histosols (H horizon) cover the 20% remaining area and are commonly located in flat areas where geomorphology allows water accumulation. These soils, formed in wetlands, are highly organic and are saturated most of the year [19]. Both soil types were The main soil type in Zhurucay are Andosols (Ah horizon), covering 80% of the area and commonly located on the hillslopes. Histosols (H horizon) cover the 20% remaining area and are commonly located in flat areas where geomorphology allows water accumulation. These soils, formed in wetlands, are highly organic and are saturated most of the year [19]. Both soil types were formed by volcanic ash accumulation; they are black, humic, acid, and rich in organic carbon and are located over a mineral horizon (C horizon) commonly rich in clay [21,61].
The vegetation within the catchment is highly correlated with soil type [19] and is typical of páramo grasslands. Tussock grasses grow in Andosols while cushion plants grow in Histosols wetlands [4,62], covering 70% and 25% of the area, respectively. The remaining 5% of the area is covered by polylepis trees and pine forest. Since the percentage of polylepis and pine forest is too low within each subcatchment, these vegetation covers were not included in the study. Further description of landscape characteristics are provided in Mosquera et al. [19].

Monitoring and Data Availability
The Zhurucay Ecohydrological Observatory was established in 2009 with the purpose of studying the hydrological functioning of Andean mountain catchments. Over the years, it has been equipped with two automatic meteorological stations, a laser disdrometer, five permanent rain gauges, a nested hydrological network consisting of 10 weirs, a hillslope to study subsurface processes, and an environmental water quality monitoring system to study runoff flowpaths and runoff formation using isotopes and metals (in soils, streams, rainfall, and springs). Recent additions include an energy balance and eddy covariance system, and small-scale lysimeters.
At the beginning of 2013, the hydrological network was fully operational. For this experiment, the rain gauge network (five permanent rain gauges) was complemented with a total of 12 rain gauges: 11 were located inside the catchment and one was located 2 km downstream of the outlet. During the experiment, this rain gauge network was arguably the densest network in the Andes mountains at this spatial scale (1.46 rain gauges per km 2 ). The extended network operated from 2014 to 2016. Hence, daily data of precipitation, temperature, potential evapotranspiration, and runoff corresponding to the period October 2013 to October 2016 were selected for this study. These variables were used according to the requirements of the HBV-light model. Precipitation data for the first study year was obtained from the five permanent rain gauges, and for the remaining years from the 12 rain gauges. Temperature and meteorological variables to estimate potential evapotranspiration were acquired from a weather station at 3780 m a.s.l. Runoff data were obtained from V-notch weirs at the outlets of seven subcatchments (S1 to S7) and from one rectangular weir at the outlet of Zhurucay catchment (S8) which is the focus of this study. Figure 1 shows the distribution of sensors and subcatchments. Subcatchments aggregated area and vegetation cover are presented in Table 1 which was provided by Mosquera et al. [19].

Methodology
First, we calibrated and validated each model structure using all rain gauges located inside the catchment. Each model structure was evaluated and compared according to the ability to simulate total runoff and its subflow components in order to select the model structure with the highest performance. Then, the model structure selected previously was calibrated and validated using six rainfall monitoring scenarios to evaluate the impact of rainfall estimation on the simulated runoff and parameters. Finally, the calibrated parameters with the worst rainfall scenario were employed to run the model with the remaining scenarios to analyze the possibility to enhance runoff simulation by improving rainfall information. The methodological details are explained in the following subsections.

HBV-Light Model
The HBV-light model [63] is a semidistributed conceptual rainfall-runoff model based on the original HBV model developed by Bergström [64] and Bergström [57] at the Swedish Meteorological and Hydrological Institute (SMHI). This model can be distributed into different elevation and vegetation zones as well as subcatchments. HBV-light uses four routines to simulate runoff: (i) the snow routine represents snow accumulation and snow melt; (ii) the soil routine describes ground water recharge and actual evaporation as function of water storage for each elevation or vegetation zone; (iii) the response routine computes the runoff as function of water stored in reservoirs; and (iv) the routing routine simulates the routing of the runoff to the outlet of the catchment by a triangular weighting function. Before the routing routine, runoff from previous subcatchments are added to the generated runoff of the current subcatchment. A detailed description of the model can be found in Seibert and Vis [65].
The standard model consists of two serial reservoirs that receive water from the semidistributed soil routine. Storage in the upper soil reservoir (SUZ) simulates fast flow and interflow, representing the near surface and subsurface flow, respectively; while storage in the lower soil reservoir (SLZ) simulates slow flow, representing baseflow. Both reservoirs are connected by a percolation rate. In total, HBV-light offers 11 different structures of reservoirs varying from one to three reservoirs with varying degrees of spatial distribution according to elevation and vegetation zones. Further details of model structures are found in Uhlenbrook et al. [66].
For this study, the Zhurucay catchment was divided into eight subcatchments and two vegetation zones representing tussock grasses (grassland) and cushion plants (wetlands). Runoff for each subcatchment were simulated considering the contribution of previous subcatchments, and therefore, the runoff simulation of subcatchment S8 is the simulation of the whole Zhurucay catchment. Structures related to snow and very slow flow were not used since Zhurucay does not receive snow. Eight model structures were used in this study, as presented in Table 2.

Model Calibration and Validation
The HBV-light model was calibrated using the standard split sample model calibration procedure [67,68]. The first, second, and third year of data (October 2013-October 2016) were used for model warming up, calibration, and validation, respectively. A year of data for each calibration and validation period were considered sufficient because of they cover the hydrological cycle of the study catchment. For the warming up period, spatial rainfall was estimated from the five permanent rain gauges (R03, R06, R07, R08, and R12), by inverse distance weighting (IDW) interpolation. This interpolation method was used due to the reduced number of rain gauges available. For the calibration and validation periods (years two and three) rainfall was estimated from the 11 rain gauges located inside the catchment, interpolated by ordinary Kriging. Potential evapotranspiration was estimated by the Penman-Monteith method [69] which has been used previously in Zhurucay [14].
The Monte Carlo procedure established by the HBV-light model was used to select an optimal parameter set after performing 10,000 simulations. Table 3 lists all parameters used and their calibration range. Parameter sets were generated using random numbers from a uniform distribution. The recession coefficient for fast flow (K0) was not calibrated and was set to 0.99 day −1 (close to one day). Interflow (K1) and slow flow (K2) coefficients were calibrated between relative fast ranges. These considerations were taken due to the fast recession observed in these subcatchments. The simulation results for each set of parameters in the Monte Carlo procedure were optimized by selecting the Nash-Sutcliffe efficiency (NSeff) [70] as the objective function.

Evaluation of Model Structures
Model structures were evaluated in two steps. First, for each model structure, HBV-light was calibrated and validated based on NSeff. Nevertheless, other performance indexes such as Nash-Sutcliffe efficiency with logarithms (Log NSeff), coefficient of determination (R 2 ) and percentage annual difference (%ADiff) were additionally calculated to provide more information about the quality of the simulation [71,72]. NSeff has been widely used to evaluate the performance of hydrological models and it is sensitive to differences in the observed and simulated means and variances. Log NSeff is the NSeff using the logarithm of the runoff values reducing its sensitivity of extreme values. R 2 describes how much of the observed dispersion is explained by the prediction. %ADiff is the percentage difference of the total observed runoff against the total simulated runoff used to observe the under or overestimation in the water balance.
Second, to identify the most behavioral models, the simulated runoff generated from the different storages of the model structures (e.g., STZ, SUZ, and SLZ) representing fast flow (FF), interflow (IF), and slow flow (SF) were compared to the flow components derived from measured runoff using the Water Engineering Time Series PROcessing tool (WETSPRO) [73]. WETSPRO is a tool that extracts the subflow components of a hydrograph based on Chapman filter [74] and using its recession constants (K) and the average fraction of each subflow component over the total flow (w). This tool has been used successfully to separate flow components and to test models by multicriteria approach [75][76][77]. Since model structures M1 and M3 simulate only two flow components, e.g., slow flow and the combination of interflow and fast flow, WETSPRO-derived interflow and fast flow were added for an adequate comparison.

Rainfall Monitoring Scenarios
We selected six rainfall-monitoring scenarios. These scenarios ranged from using all the rainfall network to using a single rain gauge located outside the catchment; they are described in Table 4. The first scenario (referred to as E0) uses 11 rain gauges inside the catchment and therefore provides the best rainfall estimation. The remaining five monitoring scenarios (E1-E5) were selected based on common practices used by water managers in the region. They simulate an scarce network by reducing the number of gauges as follows: one configuration using three rain gauges located in the upper, middle, and lower catchment (R01, R06, and R11) and four configurations using only one rain gauge in each one, located in the upper (R01), middle (R06), lower (R11), and outside the catchment (R12). For scenarios based on one rain gauge, spatial rainfall was considered uniform throughout the catchment. For the three-rain-gauge scenario, spatial rainfall was interpolated by inverse distance weighting (IDW). This method was selected rather than ordinary kriging due to the low rain gauge density.

Evaluation of the Quality of Spatial Rainfall Estimation
The best rainfall estimation (E0) was compared to the estimations obtained from the different monitoring scenarios (E1 to E5) as a way to understand the impact of the degraded rain gauge network. For this purpose, we used an analysis of scatter plot, R 2 and the Goodness of Rainfall Estimation index (GORE) [50]. GORE was used to quantify the quality of the estimated rainfall time distribution. GORE is the transposition of NSeff in the precipitation domain, using the square root of the rainfall data to reduce the weight of extreme events and is expressed as: where, n is the number of time steps of the period, P E i is the estimated rainfall from rainfall monitoring scenario at time i, P i is the best estimated rainfall available from rainfall monitoring scenario at time i, and √ P is the mean of the best estimated rainfall available (considered as reference) over the study period. Like NSeff, the GORE index varies between −∞ and 1, where 1 represents that the estimated rainfall is equal to the reference. GORE will be smaller as the rainfall estimates become poorer.

Evaluation of the Impact of Spatial Rainfall Estimation on Model Calibration
To evaluate the effect of the reduction of precipitation information by operating a sparse precipitation network on hydrological simulation, the HBV-light structure selected in Section 3.3 was calibrated and validated using each of the rainfall scenarios in Section 3.4. The performance of the simulated runoff (using NSeff) and calibrated parameters were related to the quality of the estimated rainfall (GORE) used to run the model.
Besides GORE index, BALANCE index was used to quantify the overestimation or underestimation of the total estimated rainfall and its relation to the model performance (NSeff) [50]. BALANCE index is greater than 1 in the case of rainfall overestimation and smaller than 1 in the case of underestimation, and is expressed as:

Effect of Improved Rainfall Estimation on A Model Calibrated with Scarce Data
For this objective, the model calibrated with the worst rainfall scenario was run using the input rainfall of the estimates from the remaining rainfall scenarios in a step-wise fashion, i.e., by improving the spatial rainfall estimation step-by-step up to the reference scenario (E0). Then, the performances of runoff simulations (NSeff index) were analyzed regarding the quality of rainfall information (GORE index) used to run the model.

Evaluation of the Model Structures
The performances of the eight model structures, for each subcatchment and for calibration and validation periods are shown in Figure 2. Each line represents a model structure and the overlapping of these lines indicates that all model structures have very similar performances. Indexes values can be found in supplementary Tables S1 and S2. From the magnitudes of NSeff, LogNSeff, R 2 , and %ADiff, it is hard to identify a model structure that outperforms the others, as there is little difference among them. Similar results were found by Uhlenbrook et al. [66] using six structures of HBV model in a mountainous catchment in Germany where NSeff values varied between 0.825 and 0.876 in the calibration period. This suggests that model parameters compensate for the differences among structures. Additionally, similar results are observed for all subcatchments, which can be explained by two reasons. First, the reason explained previously and second, that physical differences between subcatchments are not large enough to produce an impact in specific model structures [78].
The average of NSeff for all subcatchments shows low values for slow flow (SF) and fast flow (FF), below 0.37 and 0.20, respectively. For most of the structures, NSeff values showed high variability between subcatchments; for example, for the structure M2, slow flow NSeff ranges from 0.52 to −0.01. This shows that while model parameters can compensate the overall performance of a given model structure, the simulated subflow components do not represent the physical functioning of the catchments, and therefore these models are not behavioral [79]. For the catchment S8, which is the outlet of Zhurucay and the focus of this study, the NSeff values from eight structures vary between 0.80 and 0.83 and 0.68 and 0.73 for calibration and validation, respectively. These results were similar than those obtained by Buytaert and Beven [30] in a páramo catchment of 2.53 km 2 with similar vegetation cover and altitude. The study evaluated nine model structures with NSeff values between 0.72 and 0.87 and 0.45 and 0.77 for calibration and validation, respectively; showing that HBV-light structures reached expected efficiencies for this ecosystem.
Considering these high NSeff values and the high R 2 values (over 0.82), one might consider that the eight structures can represent the overall runoff dynamics. However, lower values of LogNSeff compared to NSeff suggest that structures may have problems to simulate low flows [71]. Therefore, we compared the observed and simulated subflow components for each model structure. This was done for the S8 catchment and the average of S1-S8 (Table 5). The average of NSeff for all subcatchments shows low values for slow flow (SF) and fast flow (FF), below 0.37 and 0.20, respectively. For most of the structures, NSeff values showed high variability between subcatchments; for example, for the structure M2, slow flow NSeff ranges from 0.52 to −0.01. This shows that while model parameters can compensate the overall performance of a given model structure, the simulated subflow components do not represent the physical functioning of the catchments, and therefore these models are not behavioral [79].
For catchment S8, NSeff values for the eight model structures present similar results between subflow components, ranging from 0.17 to 0.42, 0.48 to 0.70, and 0.65 to 0.72 for SF, IF, and FF, respectively. These results indicate that SF was difficult to simulate; indeed, all model structures had NSeff values below 0.42 for this subflow. M1 was the model structure that obtained the highest performances (0.42 SF and 0.70 IF). Therefore, the M1 structure with two reservoirs that simulate slow flow and the combination of interflow and fast flow was chosen for the remaining of the study. The fact that the simplest structure showed the best performance in simulating all subflows suggests that the páramo ecosystem is hydrologically relatively simple [19], and that its water storage-release processes are mainly controlled by water moving laterally through the organic soils to the streams [22]. Thus, a two-reservoir model structure provided an accurate simulation of the observed runoff of the catchment.
The best simulation according to NSeff (0.8 and 0.72 for calibration and validation, respectively) using the M1 structure for the subcatchment S8 (Zhurucay catchment) is shown in Figure 3. In addition, the optimal parameters found for this simulation are listed in Table 6. Temporarily we can observe that in general the simulation fits well to the dry and humid periods. However, there is a slight subestimation of the simulation in the peaks.
The best simulation according to NSeff (0.8 and 0.72 for calibration and validation, respectively) using the M1 structure for the subcatchment S8 (Zhurucay catchment) is shown in Figure 3. In addition, the optimal parameters found for this simulation are listed in Table 6. Temporarily we can observe that in general the simulation fits well to the dry and humid periods. However, there is a slight subestimation of the simulation in the peaks.

Evaluation of the Estimated Rainfall from Rainfall Monitoring Scenarios
The comparison of the daily rainfall estimated from the five monitoring scenarios (E1-E5) against the best rainfall estimation for the catchment (E0) is shown in Figure 4. It is observed that the scenario E1 using three rain gauges distributed in the upper, middle, and lower catchment represents well the spatial rainfall. For scenarios that use only one rain gauge to represent the spatial rainfall, scenario E3 (placing the rain gauge in the middle of the catchment) gives the best estimation. This result is similar to results observed by Hrachowitz and Weiler [80], and suggests that the rainfall observation in the middle of the catchment is very similar to the average of the rainfall of the entire catchments. On the other hand, it is observed that a rain gauge located in the upper catchment (E3) has better results than a rain gauge located in the lower catchment (E4). The quality of the estimation is reduced significantly when using rainfall observations from the rain gauge at 2 km downstream of the outlet (E5), producing a large scatter and overestimation. These results indicate that in this mountain setting there is a large rainfall variability at short distances, in line with results found by Buytaert et al. [33] who indicated that páramo rainfall can be highly variable, even at short distances of 4 km, suggesting a strong orographic influence in rainfall.
Water 2018, 10, x FOR PEER REVIEW 12 of 19 Figure 4. Scatter plot, determination coefficient (R 2 ), and GORE of the estimated spatial rainfall of catchment S8 from monitoring scenarios E1-E5 (axis X) against the reference estimated rainfall from monitoring scenario E0 (axis Y).

Evaluation of the Impact of Rainfall Estimation on Hydrological Simulation
The impact of the rainfall estimation from the six monitoring scenarios on the performance of the simulated runoff is presented in Figure 5. This figure shows the model performance (NSeff) against the rainfall scenario (Figure 5a), the quality of the estimated rainfall (GORE, Figure 5b), and the over-or underestimation of rainfall (Balance index, Figure 5c).
As can be seen in Figure 5a, rainfall scenarios E0 and E1 produce very similar, good results, for both calibration and validation periods. This corroborates that the good rainfall quality of scenario E1 identified in Section 4.2, is also translated into a good runoff simulation. On the other hand, and as expected, the rainfall scenario with the worst rainfall quality (E5) produces the poorest runoff simulation (below 0.31 NSeff). Similar to results of Zhang and Han [81], we can observe in the rainfall scenario E5 that the model performance decreased drastically when the spatial rainfall variability is not considered. Model performance using rainfall scenarios with one rain gauge deteriorates when compared to the best spatial rainfall. And, although the validation period seems satisfactory for scenarios E2 and E3, the efficiency in the calibration period is significantly reduced. These results show how inadequate rainfall monitoring can produce high modeling uncertainty [36].
Additionally, we can observe in Figure 4 that E3 scenario shows better correlation than E4 with E0, however, in Figure 5, the E4 scenario gives better model performance for the calibration period compared with E3. The reason of the variation of the performance between calibration and validation periods can be due to annual variations that can be strongly affected by specific orographic factors.
Results also suggest that there is a direct relation between the rainfall quality (GORE) and the simulated runoff performance (NSeff), as can be seen in Figure 5b. Additionally, the slope of this relation is 0.88 indicating that a slight reduction in the quality of the estimated rainfall can produce a big reduction in the performance of the runoff simulation. This relation was also found in some Mediterranean catchments in France [49,50]. Nevertheless, the slope found in this study is more pronounced compared to the results of Andréassian et al. [50], which suggests a stronger influence of the quality of rainfall input on the runoff simulation in this mountain catchment.
Finally, the highest performance of the simulated runoff was obtained when the over-or underestimation of rainfall was minimal. According to the results of Andréassian et al. [50], it was expected that the performance of runoff simulation would increase when the Balance Index is closer Figure 4. Scatter plot, determination coefficient (R 2 ), and GORE of the estimated spatial rainfall of catchment S8 from monitoring scenarios E1-E5 (axis X) against the reference estimated rainfall from monitoring scenario E0 (axis Y).

Evaluation of the Impact of Rainfall Estimation on Hydrological Simulation
The impact of the rainfall estimation from the six monitoring scenarios on the performance of the simulated runoff is presented in Figure 5. This figure shows the model performance (NSeff) against the rainfall scenario (Figure 5a), the quality of the estimated rainfall (GORE, Figure 5b), and the over-or underestimation of rainfall (Balance index, Figure 5c). above this range we found uncorrelated values of NSeff and BI (i.e., while for a BI of 1.11 we found a very low NSeff of 0.31, for a higher BI of 1.19 we found a NSeff of 0.63). The sensibility of model parameters to rainfall quality is shown in Figure 6. ALPHA and BETA parameters are plotted against the GORE index. These two parameters were selected to illustrate their sensibility to rainfall quality.
The only parameter sensible to rainfall quality was ALPHA (Figure 6a). In the study of Andréassian et al. [50], the sensibility is expressed as the reduction of the variability of the parameter values as the rainfall quality increases until reaching an optimal parameter value. Our study indicated that the higher the GORE, the higher the ALPHA. On the other hand, Figure 6b shows BETA as example of a parameter showing no sensibility to rainfall quality. In contrast to the findings of Andréassian et al. [50], we found that the majority of parameters were not sensible to rainfall As can be seen in Figure 5a, rainfall scenarios E0 and E1 produce very similar, good results, for both calibration and validation periods. This corroborates that the good rainfall quality of scenario E1 identified in Section 4.2, is also translated into a good runoff simulation. On the other hand, and as expected, the rainfall scenario with the worst rainfall quality (E5) produces the poorest runoff simulation (below 0.31 NSeff). Similar to results of Zhang and Han [81], we can observe in the rainfall scenario E5 that the model performance decreased drastically when the spatial rainfall variability is not considered. Model performance using rainfall scenarios with one rain gauge deteriorates when compared to the best spatial rainfall. And, although the validation period seems satisfactory for scenarios E2 and E3, the efficiency in the calibration period is significantly reduced. These results show how inadequate rainfall monitoring can produce high modeling uncertainty [36].
Additionally, we can observe in Figure 4 that E3 scenario shows better correlation than E4 with E0, however, in Figure 5, the E4 scenario gives better model performance for the calibration period compared with E3. The reason of the variation of the performance between calibration and validation periods can be due to annual variations that can be strongly affected by specific orographic factors.
Results also suggest that there is a direct relation between the rainfall quality (GORE) and the simulated runoff performance (NSeff), as can be seen in Figure 5b. Additionally, the slope of this relation is 0.88 indicating that a slight reduction in the quality of the estimated rainfall can produce a big reduction in the performance of the runoff simulation. This relation was also found in some Mediterranean catchments in France [49,50]. Nevertheless, the slope found in this study is more pronounced compared to the results of Andréassian et al. [50], which suggests a stronger influence of the quality of rainfall input on the runoff simulation in this mountain catchment.
Finally, the highest performance of the simulated runoff was obtained when the over-or underestimation of rainfall was minimal. According to the results of Andréassian et al. [50], it was expected that the performance of runoff simulation would increase when the Balance Index is closer to 1. This was found true when the Balance Index ranged from 1.07 to 0.97 (Figure 5c). However, above this range we found uncorrelated values of NSeff and BI (i.e., while for a BI of 1.11 we found a very low NSeff of 0.31, for a higher BI of 1.19 we found a NSeff of 0.63).
The sensibility of model parameters to rainfall quality is shown in Figure 6. ALPHA and BETA parameters are plotted against the GORE index. These two parameters were selected to illustrate their sensibility to rainfall quality. very low NSeff of 0.31, for a higher BI of 1.19 we found a NSeff of 0.63). The sensibility of model parameters to rainfall quality is shown in Figure 6. ALPHA and BETA parameters are plotted against the GORE index. These two parameters were selected to illustrate their sensibility to rainfall quality.
The only parameter sensible to rainfall quality was ALPHA (Figure 6a). In the study of Andréassian et al. [50], the sensibility is expressed as the reduction of the variability of the parameter values as the rainfall quality increases until reaching an optimal parameter value. Our study indicated that the higher the GORE, the higher the ALPHA. On the other hand, Figure 6b shows BETA as example of a parameter showing no sensibility to rainfall quality. In contrast to the findings of Andréassian et al. [50], we found that the majority of parameters were not sensible to rainfall quality. This may be caused by the relative high number of parameters which allows an overadjustment to the input data compared to the three and six parameters used in Andréassian et al. [50]. Although, it is still an open question if the sensibility of parameters to the rainfall quality is a function of the number of parameters used by the model. Nevertheless, at least one parameter was found sensible, showing that despite this situation, rainfall quality may still have impact on the parameterization. This can be due to a correct estimation of precipitation which allows a better description of ALPHA parameter (which controls the amount of water that recharges the streams) to be adapted to the high water recharge of these soils.

Effect of Improved Rainfall Estimation in A Calibrated Model
In this section, we analyzed the performance of the runoff simulation obtained from a calibrated model with far-from-perfect rainfall using new improved rainfall estimations. In this way, the model The only parameter sensible to rainfall quality was ALPHA (Figure 6a). In the study of Andréassian et al. [50], the sensibility is expressed as the reduction of the variability of the parameter values as the rainfall quality increases until reaching an optimal parameter value. Our study indicated that the higher the GORE, the higher the ALPHA. On the other hand, Figure 6b shows BETA as example of a parameter showing no sensibility to rainfall quality. In contrast to the findings of Andréassian et al. [50], we found that the majority of parameters were not sensible to rainfall quality. This may be caused by the relative high number of parameters which allows an overadjustment to the input data compared to the three and six parameters used in Andréassian et al. [50]. Although, it is still an open question if the sensibility of parameters to the rainfall quality is a function of the number of parameters used by the model. Nevertheless, at least one parameter was found sensible, showing that despite this situation, rainfall quality may still have impact on the parameterization. This can be due to a correct estimation of precipitation which allows a better description of ALPHA parameter (which controls the amount of water that recharges the streams) to be adapted to the high water recharge of these soils.

Effect of Improved Rainfall Estimation in A Calibrated Model
In this section, we analyzed the performance of the runoff simulation obtained from a calibrated model with far-from-perfect rainfall using new improved rainfall estimations. In this way, the model calibrated with input rainfall from scenario E5 was run for the validation period with the remaining scenarios. Figure 7 shows the NSeff index against the rainfall monitoring scenario and the GORE index. calibrated with input rainfall from scenario E5 was run for the validation period with the remaining scenarios. Figure 7 shows the NSeff index against the rainfall monitoring scenario and the GORE index. Rainfall scenarios E4 to E0 produce a considerably better runoff simulation compared to scenario E5 (Figure 7a). Furthermore, NSeff values obtained from these scenarios were very similar and increased to 0.49-0.60 from an original value of 0.14. Additionally, in Figure 7b it is observed that high values of GORE are related to high values of NSeff, although, it was not found a direct relation between these two indexes.
In this way, better rainfall estimations produced an improvement in model performance. Similar results were found by Bárdossy and Das [42], who used 10 and 20 rain gauges for calibration and validation, respectively. Nevertheless, the improvement in the model efficiency was higher in our case. In this way, any improvement in rainfall monitoring will produce a big gain in modeling efficiency, even when using a model calibrated with far-from-perfect rainfall. Additionally, the fact that bad and good efficiencies can be obtained for the same parameters suggests that for this mountain ecosystem the rainfall input could be the biggest source of model uncertainty.

Conclusions
The present study was designed to assess the impact of rainfall estimation on hydrological modeling using six rainfall monitoring scenarios in a small headwater páramo catchment. This was achieved by installing a dense network of 12 rain gauges in the Zhurucay Ecohydrological Observatory in southern Ecuador, and creating rainfall scenarios by withdrawing a given number of rain gauges.
This study showed that all model structures of HBV-light model can represent total runoff. However, the capacity to simulate subflow components strongly varied between structures. The simplest model structure M1 (standard model) had the highest performance to represent subflow components, although all model structures have problems properly representing slow flow.
The research has also shown that having good spatial rainfall measurements is essential to achieve good modeling results in mountainous areas. We can conclude that a limited number of rain gauges can produce acceptable modeling performance. However, it strongly depends on the location of the rain gauges. In this case, three rain gauges in the upper, middle, and lower catchment worked well, but this has to be confirmed in other páramo catchments in order to generalize this knowledge. Furthermore, a better description of the spatial rainfall of the catchment not only enhances the runoff Rainfall scenarios E4 to E0 produce a considerably better runoff simulation compared to scenario E5 (Figure 7a). Furthermore, NSeff values obtained from these scenarios were very similar and increased to 0.49-0.60 from an original value of 0.14. Additionally, in Figure 7b it is observed that high values of GORE are related to high values of NSeff, although, it was not found a direct relation between these two indexes.
In this way, better rainfall estimations produced an improvement in model performance. Similar results were found by Bárdossy and Das [42], who used 10 and 20 rain gauges for calibration and validation, respectively. Nevertheless, the improvement in the model efficiency was higher in our case. In this way, any improvement in rainfall monitoring will produce a big gain in modeling efficiency, even when using a model calibrated with far-from-perfect rainfall. Additionally, the fact that bad and good efficiencies can be obtained for the same parameters suggests that for this mountain ecosystem the rainfall input could be the biggest source of model uncertainty.

Conclusions
The present study was designed to assess the impact of rainfall estimation on hydrological modeling using six rainfall monitoring scenarios in a small headwater páramo catchment. This was achieved by installing a dense network of 12 rain gauges in the Zhurucay Ecohydrological Observatory in southern Ecuador, and creating rainfall scenarios by withdrawing a given number of rain gauges.
This study showed that all model structures of HBV-light model can represent total runoff. However, the capacity to simulate subflow components strongly varied between structures. The simplest model structure M1 (standard model) had the highest performance to represent subflow components, although all model structures have problems properly representing slow flow.
The research has also shown that having good spatial rainfall measurements is essential to achieve good modeling results in mountainous areas. We can conclude that a limited number of rain gauges can produce acceptable modeling performance. However, it strongly depends on the location of the rain gauges. In this case, three rain gauges in the upper, middle, and lower catchment worked well, but this has to be confirmed in other páramo catchments in order to generalize this knowledge. Furthermore, a better description of the spatial rainfall of the catchment not only enhances the runoff simulation but also the possibility to select an optimal parameter value.
Another significant finding that emerges from this study is that a calibrated model based on far-from-perfect rainfall estimates can produce acceptable runoff simulations when was run with new improved rainfall data, something that can be highly valuable by water managers.
Despite that rain gauges are considered the most trusted source of rainfall information, they are subject to errors. However, differences between rain gauges found here were larger (GORE from 0.41 to 0.98) than those caused by measurement errors. The effect of factors such as wind, elevation, topographic location, and mechanical errors on rainfall measuring were out of the scope of this study. Nevertheless, authors encourage the study of such factors mainly due to the lack of information about rainfall measuring on páramo ecosystem.
Overall, this study has showed that rainfall input could be the largest source of model uncertainty for this mountain ecosystem. Therefore, our findings have provided a first insight of the importance of rainfall monitoring for hydrological modeling in páramo catchments, which are the main water supply for millions of people.
Author Contributions: A.S. analyzed the data and wrote the manuscript. R.C. designed and supervised the study and provided a critical revision of the manuscript.
Funding: This research was part of the project "Identificación de los procesos hidrometereológicos que desencadenan crecidas extremas en la ciudad de Cuenca" funded by the Dirección de Investigación of Universidad de Cuenca (DIUC) and the Empresa Pública Municipal de Telecomunicaciones, Agua Potable, Alcantarillado y Saneamiento de Cuenca (ETAPA-EP).