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

Abstract

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 km²) 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.

1. 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].

Research on páramo hydrology is relatively new and has focused mainly on understanding hydrological processes such as evapotranspiration [13,14], interception [15], and temporal and spatial variability of precipitation [16,17]; hydrological functioning such as hydrological landscape controls [18,19] and water provenance and transit times [20,21,22]; impact of land use change [10,11,23]; and weather and climate [24,25,26]. Nevertheless, rainfall-runoff modeling (and the uncertainties related to model structure, input data and model parameters) has still been little studied, even though it is key for hydrological applications and water resources management.

Modeling of páramo catchments has concentrated on studying the impact of land use [27,28,29] and hypothesis testing [30,31]. Nevertheless, given the high rainfall variability found in this region [16,32], that can reach up to 25% volume differences between rain gauges within small catchments [33] and the scarcity of spatiotemporal hydrometeorological data, it is clear that the lack of rainfall monitoring can have a large impact on modeling [12].

Several studies have analyzed the impact of rainfall observations in model calibration and results (e.g., time to peak, peak flows, volume, performance, and parameters) [34,35,36,37,38,39]. Younger et al. [40] posits that synthetic rainfall perturbations, mostly located in the upper and lower catchment, induce changes in peaks and model parameters at the outlet, highlighting the importance of an adequate estimation of spatial precipitation. Faurés et al. [41] found differences of 2% to 65% in runoff volume if just one of five rain gauges was used in a small catchment. Similar conclusions were provided by Bárdossy and Das [42] and Xu et al. [43], showing that the performance of runoff simulation improves and the uncertainty is reduced by increasing the number of rain gauges until a given threshold. In order to overcome this measuring problem, Berne et al. [44] recommended spatial and temporal sampling for modeling purposes, making use of high quality and quantity of rainfall data, which is in most cases difficult to obtain. Nevertheless, it is not possible to generalize these results to other ecosystems, mainly due to the fact that the impact of precipitation on runoff modeling will be influenced by the characteristics of the catchments and storms [45,46].

Singh [47] concluded that the hydrological response of a catchment is related to the spatial and temporal rainfall variability, suggesting that the quality of hydrological modeling can be related to the capacity to measure this variability. Lobligeois et al. [48] tested this hypothesis using 181 catchments. It was found that those with higher rainfall variability needed a denser rain gauge network for improving runoff simulations. Finally, other studies have found that the quality of the model performance is indeed related to the quality of precipitation information [49,50].

These studies have clearly established the importance of analyzing the effect of rain gauge density for hydrological modeling [51]. However, none of the studies cited previously has been carried out in high mountainous regions. Indeed, only a few studies such as Arnaud et al. [52] and Girons, Lopez, and Seibert [53] have been carried out in high altitudes (with a maximum altitude of 2503 and above 5000 m a.s.l., respectively), leaving a gap in knowledge about this topic in mountain areas. Closing this gap will also have a positive impact on (i) operational hydrology, by providing solid recommendations to water managers about the design of rainfall observation networks in the region, a problem that has been previously identified by stakeholders [54], and (ii) the implementation and evaluation of watershed management programs as well as to confront the ecohydrological challenges of the region [55,56].

In this context, the objective of this study is to evaluate the impact of precipitation estimation in the hydrological simulation of a mountain catchment. For this purpose we implemented a dense rainfall network (12 rain gauges) in the Zhurucay Ecohydrological Observatory (7.53 km²), in southern Ecuador. The semidistributed conceptual rainfall-runoff HBV-light (Hydrologiska Byråns Vattenbalansavdelning) model was chosen because has been widely used with good results in several ecosystems including mountain ecosystems [53,57,58,59], and offers a good tradeoff between the amount of information required and the spatial representation. First, we identified the best model structure to simulate the total discharge and its subflow components (fast flow, interflow, and slow flow). Secondly, we evaluated model performance and parameter sensibility using six rainfall-monitoring scenarios using from 1 to 11 rain gauges. Finally, we explored one open question of Andean water managers. A common monitoring setting in mountain areas consists of having a single 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.

2. Study Area, Monitoring and Data Availability

2.1. 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² 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 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].

2.2. 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²). 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. Subcatchment aggregated area and vegetation cover.

Subcatchment	Aggregated Area (km2)	Vegetation Cover (%)
		Wetland	Tussock Grass
S1	0.20	15	85
S2	0.38	13	87
S3	0.38	18	82
S4	0.65	18	82
S5	1.40	17	83
S6	3.28	24	76
S7	1.22	65	35
S8	7.53	25	75

which was provided by Mosquera et al. [19].

3. 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.

3.1. 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. HBV-light model structures.

ID	Structure a
M1	Two boxes (SUZ and SLZ, standard).
M2	Two boxes (SUZ and SLZ). UZL threshold.
M3	Two boxes (SUZ and SLZ). SUZ distributed.
M4	Two boxes (SUZ and SLZ). SUZ distributed and UZL threshold.
M5	Three boxes (STZ, SUZ and SLZ).
M6	Three boxes (STZ, SUZ and SLZ). STZ distributed.
M7	Three boxes (STZ, SUZ and SLZ). STZ and SUZ distributed.
M8	One box. UZL and PERC thresholds.

^a STZ = Storage in top zone; SUZ = Storage in upper zone; SLZ = Storage in lower zone; UZL = Threshold parameter above which overland flow is produced; PERC = Threshold parameter above which interflow is produced.

.

3.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. HBV-light model parameters and their range used in the Monte Carlo procedure.

Parameter	Description	Unit	Minimum	Maximum
Soil routinea
FC	Maximum soil moisture	mm	250	400
LP	Soil moisture (MS) above which actual evapotranspiration reaches potential evapotranspiration (MS/FC)	-	0.5	1
BETA	Relative contribution to runoff from rain or snowmelt	-	1	3
Response routine
K0 b	Recession coefficient for quick flow	day−1	0.999	0.999
K1	Recession coefficient for interflow	day−1	0.33	0.999
K2	Recession coefficient for baseflow	day−1	0.066	0.2
ALPHA c	Nonlinearity coefficient	-	0	1
UZL d	Threshold parameter for K0 outflow	mm	0	30
PERC e	Percolation from SUZ to SLZ	mm day−1	0	2
UZL f	Percolation from STZ to SUZ	mm day−1	0	30
PERC g	Threshold parameter for K1 outflow	mm	0	2
Routing routine
MAXBAS	Length of triangular weighting function	d	1	1.5

^a Parameters for each vegetation cover; ^b Parameter used only for structures M2, M4, M5, M6, M7, and M8; ^c Parameter used only for structures M1 and M3; ^d Parameter used only for structures M2, M4 and M8; ^e Parameter used only for structures M1, M2, M3, M4, M5, M6, and M7; ^f Parameter used only for structures M5, M6, and M7; ^g Parameter used only for structure M8; FC is the previous parameter and the expression defines the relationship to obtain LP.

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⁻¹ (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.

3.3. 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²) 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² 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.

3.4. 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. Rainfall monitoring scenarios.

ID	Rain Gauge(s)	Spatial Rainfall Estimation Method
E0	R01 to R11	Ordinary kriging
E1	R01, R06, R11	IDW
E2	R01	Punctual
E3	R06	Punctual
E4	R11	Punctual
E5	R12	Punctual

. 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.

3.5. 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² 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:

$$\mathrm{GORE}=1- \frac{{\sum }_{i=1}^n{\left(\sqrt{P_i^E}- \sqrt{P_i}\right)}^2}{{\sum }_{i=1}^n{\left(\sqrt{P_i}- \overline{\sqrt{P}}\right)}^2}$$

(1)

where, n is the number of time steps of the period, $P_i^E$ 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 $\overline{\sqrt{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.

3.6. 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:

$$\mathrm{BALANCE}\mathrm{index}= \frac{{\sum }_{i=1}^n{P}_i^E}{{\sum }_{i=1}^n{P}_i}$$

(2)

3.7. 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.

4. Results and Discussion

4.1. 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², 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].

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² 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² 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. Nash–Sutcliffe efficiency for simulated flow components from each model structure. For outlet subcatchment S8 and the average of S1–S8.

Model Structure	Flow Components S8			Average Flow Component (S1–S8)
	SF	IF	FF	SF	IF	FF
M1	0.42	0.70		0.32	0.72
M2	0.17	0.68	0.67	0.07	0.61	0.08
M3	0.33	0.67		0.37	0.71
M4	0.28	0.64	0.72	0.13	0.60	0.27
M5	0.39	0.48	0.66	0.27	0.52	0.16
M6	0.35	0.53	0.65	0.05	0.55	0.20
M7	0.36	0.62	0.66	0.19	0.55	0.09
M8	0.31	0.63	0.70	0.05	0.60	0.30

).

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. Optimal parameters from the best runoff simulation according to NSeff.

Routine	Parameter	Value
Soil routine 1 (Wetland)	FC_1	294.44
	LP_1	0.97
	BETA_1	2.79
Soil routine 2 (Tussock grass)	FC_2	254.11
	LP_2	1.00
	BETA_2	2.46
Response routine	PERC	1.90
	ALPHA	0.43
	K1	0.42
	K2	0.16
Routing routine	MAXBAS	1.00

. 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.

4.2. 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.

4.3. 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 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.

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.

4.4. 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.

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.

5. 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.