River Discharge Simulation in the High Andes of Southern Ecuador Using High-Resolution Radar Observations and Meteorological Station Data

Abstract

The prediction of river discharge using hydrological models (HMs) is of utmost importance, especially in basins that provide drinking water or serve as recreation areas, to mitigate damage to civil structures and to prevent the loss of human lives. Therefore, different HMs must be tested to determine their accuracy and usefulness as early warning tools, especially for extreme precipitation events. This study simulated the river discharge in an Andean watershed, for which the distributed HM Runoff Prediction Model (RPM) and the semi-distributed HM Hydrologic Modelling System (HEC-HMS) were applied. As precipitation input data for the RPM model, high-resolution radar observations were used, whereas the HEC-HMS model used the available meteorological station data. The obtained simulations were compared to measured discharges at the outlet of the watershed. The results highlighted the advantages of distributed HM (RPM) in combination with high-resolution radar images, which estimated accurately the discharges in magnitude and time. The statistical analysis showed good to very good accordance between observed and simulated discharge for the RPM model (R²: 0.85–0.92; NSE: 0.77–0.82), whereas for the HEC-HMS model accuracies were lower (R²: 0.68–0.86; NSE: 0.26–0.78). This was not only due to the application of means values for the watershed (HEC-HMS), but also to limited rain gauge information. Generally, station network density in tropical mountain regions is poor, for which reason the high spatiotemporal precipitation variability cannot be detected. For hydrological simulation and forecasting flash floods, as well as for environmental investigations and water resource management, meteorological radars are the better choice. The greater availability of cost-effective systems at the present time also reduces implementation and maintenance costs of dense meteorological station networks.

1. Introduction

Precipitation, along with topography, vegetation cover, and soil types, control the rain-runoff processes within a watershed [1,2,3]. On the one hand, precipitation duration, intensity, and spatial distribution are strongly influenced by local topography, wind speed, and direction, as well as air moisture [4,5], whereas natural vegetation generally regulates seasonal river discharge, generates good water quality, and reduces the splash effect of raindrops, mitigating negative effects upon soil erosion processes during extreme rainfall events [6,7]. On the other hand, soils can store large amounts of water, depending on organic matter content and specific soil parameters [8]. However, land use changes caused by anthropogenic activities can degrade soils, which reduces their water regulation capacity and increases the flood risk [9].

In this context, different hydrological models (HMs) have been developed to estimate river discharge due to precipitation events [10,11]. By means of these models it is possible to estimate the peak flows in magnitude and time, and, consequently, to reduce the vulnerability to natural disasters related to rain-runoff processes [12]. Furthermore, these models allow improvement of water resource management and the design of functional civil engineering structures, such as drinking water supply systems [13,14].

To execute HMs, different climatic and physical parameters of the basin as input data are necessary, in which precipitation is the most important because it determines the amount of water available for rain-runoff processes. The principal HMs are [15,16]: (i) stochastic models, which generate random predictions; (ii) conceptual models, which include the physical processes in the watershed, obtained by means of field measurements; and (iii) deterministic models. Deterministic models are divided into aggregate models, which consider the watershed as one unit (all values are averaged and assumed for the whole watershed); semi-distributed models, which divide the watershed into subunits and apply specific mean values for each of them; and distributed models, which consider the spatiotemporal distribution of climatic and physical parameters [17,18].

Currently, distributed HMs (DHMs) are recommended because they consider the spatiotemporal distribution of all input parameters [19,20]. Therefore, these models require grid data for a reliable discharge simulation, flood prediction, and, consequently, for human life loss prevention [18,21]. However, to obtain precise values respective to precipitation amounts and their spatiotemporal distribution is often complicated, because of the scarce information available in many parts of the world, especially in tropical high mountain watersheds where strong altitudinal gradients exist and instrumentation is generally poor due to difficulties in access [22,23,24]. Nonetheless, only with adequate precipitation data (spatially and temporally), can DHMs efficiently simulate the river discharge in high mountain watersheds [25,26,27].

Technology advances in recent decades, particularly the implementation of meteorological radar systems, can be an alternative for the collection of accurate precipitation data in high tropical mountains (e.g., [28,29,30]). Meteorological radars are especially useful in remote areas, where difficulties in access and the topography limit the operation of traditional monitoring networks. At these sites, meteorological radars can provide information in high spatiotemporal resolution, which improves the simulation of hydrological processes [18,31,32].

In general, meteorological radars emit microwaves in the form of electromagnetic pulses, which are back-scattered by precipitation particles inside the atmosphere (e.g., rain, snow, or hail [33]). The back-scattered electromagnetic waves are detected by the radar antenna and their reflectivity (Z) is registered, which can be converted to rain rate (R) by means of specific Z–R relationships for each rain type and region [34]. However, to estimate accurately precipitation amounts and their distribution (quantitative precipitation estimations; QPE), radar data calibration is necessary [35,36], especially in high mountains where obstacles (clutter), such as mountain peaks, obstruct the emitted pulses and do not leave valid information for these sites [22,31,37]. For correcting this shortcoming, different calibration methods based on robust algorithms exist, which generally also include limited information from surface rain gauges (e.g., [35,38,39]). Consequently, the collected data from meteorological radars are widely used in climatology, ecology, and hydrology [30,40]. However, their application in tropical zones is still scarce because of the high costs of this technology, which is the reason that developing countries generally do not possess these systems [29].

In high Andean watersheds, only isolated hydrological studies exist, due to lack of efficient monitoring networks [12], and the interaction of different mechanisms and processes in precipitation formation, which cause high spatiotemporal variability [22,34]. Therefore, precise precipitation amounts, and their distribution, are often not detectable by the poor station networks [41]. Nevertheless, hydrological studies in the Andes are of utmost importance because this information is vital for millions of people who benefit from the hydrological services of these watersheds [42,43].

In Ecuador, hydrological studies are almost non-existent. Only two investigations have been conducted, of which one study simulated the runoff of the Zhurucay River (Province of Azuay), applying an aggregate HM (Nedbor–Afstromnings model; NAM). However, due to scarce and inaccurate precipitation data, model evaluation resulted in Nash–Sutcliffe efficiency coefficients (NSE) of 0.5, when comparing measured and simulated river discharges [10]. Another study was conducted by [20] in the San Francisco catchment (Province of Zamora-Chinchipe) using four semi-distributed models (HEC-HMS, CHIMP, SWAT, and LASCAM), one aggregate model (HBV-light), and one distributed model (HBV-N-D). Their analysis resulted in NSE coefficients of 0.60, 0.78, 0.79, 0.55, 0.73, and 0.73 respectively, although this basin is one of the best instrumented in southern Ecuador [44].

Currently, Ecuador relies on data from six meteorological radars, of which three are installed in Quito (the capital of Ecuador; urban radars) operated by the National Weather Service (Instituto Nacional de Meteorología e Hidrología; INAMHI). In southern Ecuador, three radar systems operate (LOXX, GUAXX, and CAXX), which monitor the precipitation in four provinces (Loja, El Oro, Zamora Chinchipe, and Azuay). These radars were installed by the RadarNet-Sur project between 2012 and 2015 [29], and are operated by the provincial government of Loja (radar GUAXX), the Climate Observatory of the Technical University of Loja (UTPL; radar LOXX), and the Public Municipal Company of Telecommunications, Drinking Water, Sewerage and Sanitation of Cuenca (ETAPA EP; radar CAXX). To date, RadarNet-Sur data have mainly been used for climate and ecological research, in addition to the development of methods for radar data calibration [22,29,39]; however, for hydrological applications this technology has not yet been applied.

The aim of the present study is to simulate the river discharge within an Andean watershed in southern Ecuador (El Carmen), using high-resolution radar observations and meteorological station data. Thereby, the utility of meteorological radars for hydrological, climatological, and ecological applications in tropical high mountains should be demonstrated, because meteorological station data is generally sparse or inexistent, especially in smaller watersheds. To reach this goal, two different HMs were tested, namely, the DHM Runoff Prediction Model (RPM), developed by the National Autonomous University of Mexico [45], and the semi-distributed HM Hydrologic Modelling System (HEC-HMS), developed by the US Army Corps of Engineers [46]. As precipitation input data for the RPM, the radar observations were used, whereas for the HEC-HMS, the available station data were employed. To evaluate the effectiveness of the models, three representative precipitation events for the rainy season were analyzed, and the simulated river discharges compared to the measured discharges at the outlet of the watershed. The El Carmen watershed was selected because it provides drinking water for the local and regional population [47], and is also used as a recreation area. Therefore, an accurate HM is of utmost importance to forecast river discharges and to provide an early warning tool for extreme rainfall events. This will also support the local water resource management to evaluate strategies for mitigating damage to civil engineering structures, such as those related to drinking water supply.

2. Materials and Methods

2.1. Study Area

The El Carmen watershed is located in the high Andes of southern Ecuador next to the city of Loja, between the coordinates UTM 17S 703883–706554 and 9550516–9553385. This watershed is part of the Eastern Cordillera of the Andes, which separates the Amazon Basin from the Inter-Andean valleys. The elevation ranges between 2315 m a.s.l. at the river outlet and 3417 m a.s.l. at the highest mountain peaks (Figure 1). The watershed has only an area of 5 km² but is of great hydrologic importance because it supplies drinking water for the city of Loja and its surroundings [7,47]. The middle and upper zone of the watershed is part of the Podocarpus National Park (PNP), where native vegetation predominates (tropical mountain forest and paramo). At the lower parts, intervened vegetation is more frequent (secondary forest, reforestation, and pastures), due to enhanced population pressure in recent decades (Figure 3b). Soils in the middle and upper zone are mainly Inceptisols with high organic matter contents, while in the lower parts Entisols prevail. Clay content is generally high in both soil types [7].

The climate in the study area is perhumid, characterized by strong altitudinal gradients of temperature and precipitation (e.g., [48]), caused by the local topography and predominant easterly winds [22,39]. However, during the rainy season (December to April) the easterly winds (trade winds) are frequently interrupted by westerly winds from the Pacific Ocean, generating strong convective storms over the Inter-Andean valleys [30,49]. A relatively dry season occurs between June and October, when the tropical easterlies are strongest and constant, but rainfall intensities are lower because the eastern Cordillera of the Andes forms a climatic barrier, which hampers the humidity transport into the watershed (Figure 1). May and November are transition months, where precipitation is variable [22]. The average annual precipitation in the El Carmen watershed varies between 1150 mm at the lower parts and 2150 mm at the mountain tops. Average annual temperature also depends on elevation, and ranges between 15.0 °C at the river outlet and 7.3 °C at the mountain peaks [7].

2.2. Data and Materials

Climate and hydrology data for the study area were collected from various sources. The National University of Loja (UNL), together with the Regional Water Fund (FORAGUA), have operated three automatic hydrometeorological stations in El Carmen watershed since October 2015. The first is a meteorological station (Carmen_Met: 2550 m a.s.l, Campbell Scientific, 2009), which measures basic climatic variables within a 10-minute interval (temperature, precipitation, air humidity, solar radiation, wind, and pressure). The second is a rain gauge station (Carmen_Plu: 2354 m a.s.l), which comprises a TB6 rain gauge (Hydrological Services, 2015), also installed at the Carmen_Met station, and recording precipitation in intervals of 10 minutes with an accuracy of 0.1 mm. The third is a hydrological station (Carmen_Hid: 2350 m a.s.l, Campbell Scientific, 2013), which records the water level by means of a submersible pressure sensor (model CS451) in a 10-minute interval and accuracy of ±0.1%. The measured water level is transformed to discharge using a specific rating curve calibrated for the catchment (UNL). The hydrological station is installed at the river outlet near the drinking water collection plant (Figure 1).

Additional meteorological information was provided by UTPL, which has operated five automatic stations (Militar: 2033 m a.s.l; Jipiro: 2218 m a.s.l; Ventanas: 2816 m a.s.l; Técnico: 2377 m a.s.l; and Villonaco: 2950 m a.s.l.; Figure 1) in the valley of the city of Loja since 2011. These stations are Davis Wireless Vantage Pro2^TM Plus stations and contain sensors for all basic meteorological variables, as well as ultraviolet radiation. The information is stored in a datalogger in a 10-minute resolution [30]. In addition, there is the official automatic meteorological station for the city of Loja, operated by INAMHI (La Argelia: 2160 m a.s.l; Figure 1), which has provided historical data since 1964. The current automatic station is from Vaisala and stores data in a 1-minute interval. The rain gauge is from Texas Electronics (model TR-525M, Texas Electronics, 2015) and records precipitation with a 0.1 mm accuracy.

An X-band meteorological radar (local area weather radar; LAWR), manufactured by the Danish Hydraulic Institute (DHI; [50]), was installed in 2013 at the upper zone in the eastern mountain ridge next to the city of Loja (El Tiro; LOXX radar, 2850 m a.s.l.; Figure 1). The LAWR is based on the FURUNO 1525 Mk III (maritime radar) with a transmission frequency of 9410 ± 30 MHz and a bandwidth of 3 MHz. The pulse length is 1.2 μs with a repetition rate of 600 Hz. The radar emits a microwave beam within a weight of 0.92° horizontally and 10° vertically (up and down). The radar automatically generates 3 observations every 5 minutes with different ranges and resolutions (radius: 15 km, resolution: 100 × 100 m; radius: 30 km, resolution: 250 × 250 m; radius: 60 km, resolution 500 × 500 m [22,29,30]), which are transferred in real time to a server located at UTPL. For this study the highest resolution data (100 × 100 m) were used, because the El Carmen watershed lies within a 15 km radius and the DHMs prefer the highest resolution available for an accurate estimation of rainfall-runoff processes [46].

A digital elevation model (DEM) in a 30 × 30 m resolution was downloaded from the NASA database (Shuttle Radar Topography Mission, SRTM; [51]) and the errors (cells without information) were corrected using the PIT REMOVAL module of the Idrisi Selva software [52]. Subsequently, the resolution was re-sampled to obtain the same precision as the radar data (100 × 100 m; [53]). The resulting DEM was used for the calibration of the radar observations and to calculate the morphometric parameters of the watershed.

The vegetation map in a 3 × 3 m resolution was facilitated by the National Planning and Development Secretary of Ecuador (SENPLADES, [54]). It was generated by means of a supervised classification of 23 orthophotos from the year 2012. A soils map with a 30 × 30 m resolution was prepared by [7] based on 38 field measurements. The soil organic matter was determined by “humid oxidation” applying the Walkley–Black method [55], whereas permeability was estimated qualitatively in the field using the “Variable load and Structural class” method. The vegetation and soil maps were adjusted to obtain the same resolution as the radar data (100 × 100 m) by grid accumulation, in which the grids were continuously reclassified according to the dominant vegetation and soil type [56].

Two HMs were selected for runoff modelling in the El Carmen watershed: (i) the semi-distributed model HEC-HMS developed by the US Army Corps of Engineers for rain-runoff simulations, and (ii) the distributed Runoff Prediction Model (RPM) developed by the National Autonomous University of Mexico [45]. Both models are open access, which facilitates their use in developing countries, and are available online at [57] (HEC-HMS, Version 4.3) and [58] (RPM).

For the present study, three representative precipitation events for the rainy season were selected, specifically 5 March, 2016 (EVENT 1); 11–12 April, 2016 (EVENT 2); and 27–30 April, 2016 (EVENT 3). According to historical data of La Argelia meteorological station (INAMHI), 2016 was a typical year in terms of precipitation amounts.

2.3. Preparation of Input Data

By means of the DEM and the ArcGIS 10.3 Hydrology extension, the morphometric parameters of the watershed were calculated, including area, perimeter, length, and slope of the main stream, as well as the average watershed slope. In order to know the watershed susceptibility to flash floods [59], the form factor (Equation (1)) and the Gravelius coefficient (Equation (2)) were calculated using the following equations [60,61,62]:

$$IF=\frac{A_p}{L_a}$$

(1)

$$K_C=0.28\times \frac{P}{\sqrt{A}}$$

(2)

where IF is the form factor, Ap is the average width (km), La is the axial length (distance (km) in a straight line from the outlet to the highest point of the watershed), K_C is the Gravelius coefficient, P is the watershed perimeter (km), and A is the area of the watershed (km²).

The adjusted soil map based on [7] was reclassified according to its permeability group, using the information gathered during the field sampling campaign (Variable load and Structural class) and applying the hydrological group classification (Type A = Very High, B = Moderate, C = Low, and D = Very Low) of the US Soil Conservation Service (SCS) [63]. Subsequently, the vegetation map and soil map were combined to calculate the curve number (CN) at each grid cell within this basin. The CN was calculated using the HEC-GeoHMS module in ArcGis 10.3 software, which also establishes a uniform coordinate system to preserve the spatial properties of each parameter at each grid cell. The HEC-GeoHMS module offers two coordinate systems (SHG, standard hydrologic grid; and HRAP, Hydrologic Rainfall Analysis Project), of which the SHG was chosen as recommended by the Hydrologic Engineering Center for hydrological studies [64]. The conversion generates a text file (.txt), containing the coordinates and values for each grid cell as follows: (i) ShgX (X coordinate), (ii) ShgY (Y coordinate), (iii) Mod_Area, (iv) FlowLength (flow length), and (v) CN [65]. This information allows for watershed discretization using the ModClark method [3]. More detail on the different steps can be found in [45].

High spatiotemporal precipitation maps for the study area were generated by means of the LOXX radar data (Figure 1). For this, first, the 5-minute raw data were corrected and then calibrated by means of rainfall information from ground stations following the methodology proposed by [22,29,30]. This radar type (LAWR, X-band) has certain limitations which demand geometric and sensitivity corrections, due to the emitter’s degradation (magnetron; [49]). These corrections included: (i) attenuation correction and rectification of the decreasing sensitivity of the radar beam with increasing distance to the emitter by means of an empirical exponential function; (ii) detection of shaded or blocked cells by means of the DEM, and correction of these cells by proportionally scaling up the remaining beam fraction; and (iii) detection and extraction of clutter, using a map derived from radar data during periods without rain. The affected grids were filled by bilinear interpolation. Detailed information respective to the single correction steps can be found in [29].

After the geometric corrections of the 5-minute observations, they were added to obtain 10-minute data, because the rain gauges in the study area provide information for this time interval. Based on the corrected and summed radar observations, the path of the storm for each selected event could be determined, identifying the center of the storm and its displacement during each event (maximum reflectivity, [30]).

Then, the LAWR data, summed and corrected, were calibrated by means of the available rain gauge information, which were interpolated by ordinary kriging using the software Surfer to generate an isohyets map in the same spatial resolution as the radar images (100 × 100 m; [26]). Then, the radar reflectivity data and the isohyet map for the same time interval were related at all locations where a meteorological station is installed. The relation of these grids was interpolated by ordinary kriging to generate a relation surface with the same spatial resolution as the radar data [47]. Finally, the corrected and summed LAWR observations were calibrated for each 10-minute interval based on the relation surface, dividing the corrected radar reflectivity by the relation value of each grid cell. Additionally, a weight map was generated to make optimal use of the available station information. The weight map considers the horizontal and vertical location of each station and defines their areas of influence. This method assigns the measured precipitation amount to each grid cell where a station is located, and, with increasing distance to the station, the value slowly changes to the calibrated value of the radar data [22].

The high-resolution precipitation maps based on the radar observations were generated for the entire valley of the city of Loja (Figure 1), because the inclusion of more rain gauge information close to the watershed under study improves the accuracy of the radar calibration [39]. Subsequently, the study area was separated to obtain the spatiotemporal rainfall distribution within the El Carmen watershed for each selected event, which was used as an input for the DHM.

On the other hand, to generate precipitation maps for the semi-distributed HM, the Thiessen polygon method was applied [66], considering only the two rain gauge stations inside the catchment (Carmen_Met and Carmen_Plu; Figure 1) and the official weather station of INAMHI (Argelia), which is closest to the study catchment. The Thiessen polygon method draws straight lines between the stations and generates perpendicular mediatrices, which will be prolonged until they are cut with other neighboring lines. This also determines the area of influence of each meteorological station and assigns the measured station value to this area.

2.3.1. Distributed Hydrological Model (RPM)

The spatial CN and the high-resolution precipitation maps are the main input variables to execute the DHM selected for this study. The RPM, like other DHMs, requires raster information (grids), mostly with .txt extension, as data input. Therefore, the calibrated LAWR observations were converted into text format using the Watershed Modelling System software (WMS; student version [67]), whereas the CN text file (.txt) was generated by means of the Hec-GeoHMS tool as described previously. The precipitation and CN information for the RPM must be provided by two separated files, where the first file contains the hydrological parameters of the watershed (i.e., (i) ShgX (X coordinate), (ii) ShgY (Y coordinate), (iii) Mod_Area, (iv) FlowLength, and (v) CN), and the second file contains the high-resolution precipitation information for the established time interval (here: 10 min), as well as the respective coordinates of each grid cell. Furthermore, RPM uses an additional parameter, named “forgetting factor” (fx), which considers former rainfall events. The fx value ranges between 0 and 1, where 1 indicates higher runoff and lower evapotranspiration, because recent precipitations resulted in saturated soils [45].

Based on the precipitation amounts and the CN values at each grid cell, the effective precipitation (Pe) was calculated, applying the following equation (Equation (3); [46]):

$$P_e=\frac{{\left(P-\frac{5080}{CN}+50.8\right)}^2}{P+\frac{20320}{CN}-203.20}$$

(3)

where P_e is the effective precipitation (mm), P the measured precipitation (mm), and CN the curve number at each grid cell.

The transformation of Pe to runoff was carried out by the unit hydrograph method (UH) applying the modified Clark method ([3]; Figure 2). In general, the ModClark method consists of transiting the runoff generated at each grid cell to the watershed outlet by means of an isochronous map (curves with equal concentration time; tc) and a linear reservoir (K), which represents the storage effects at each grid cell [68]. To generate the UH, the two required parameters for each grid cell (concentration time (tc) and storage coefficient (K)) must be determined; tc was calculated using the Kirpich equation (Equation (4); [18]), because of its applicability for scarce information regions, whereas for K the suggested equation for practical purposes by [45] was applied (K = 0.8 * tc).

$$t_c=0.000325\times \left(\frac{L^{0.77}}{{\overline{S}}^{0.385}}\right)$$

(4)

where tc is the concentration time (hours), L the main stream length (m), and $\overline{S}$ the average watershed slope (%).

The base flow in the El Carmen basin was determined by means of the recession constant method, using the WETSPRO software (Water Engineering Time Series PROcessing Tool [70]). The software applies a linear reservoir model, based on the measured field data, and determines the watershed flows (base, subsurface, and surface flow), applying Chapman’s numerical filtering technique [71]. The function is based on the general slow-pass filter equation and assumes an exponential recession to sub-flows.

Finally, RPM was manually calibrated (method: trial and error [46]), adjusting the parameters tc and fx. In general, this type of calibration consists in manually modifying the most sensitive parameters for hydrological response, using an iterative and repetitive process until the best adjustments between measured and simulated flows are reached. The adjusted parameters were applied to all selected events to provide an early warning tool for the El Carmen watershed.

2.3.2. Semi-Distributed Hydrological Model (HEC-HMS)

HEC-HMS 4.3 software has four principal components: (i) an input parameter analysis; (ii) a graphic interface for illustrations such as the hietogram; (iii) a database, where analyzed data are stored and managed; and (iv) an interface to present the model results. The software provides different methods for runoff simulation [45]. For the present study, the SCS curve number (SCSCN) method was selected because it provides accurate results, despite its simplicity, comparable to other more complex methods [72,73]. Furthermore, this method can be applied for different environments and multiple applications, such as forecasting or early warning systems [74].

The SCSCN method needs two parameters as inputs, specifically the weighted CN of the watershed (CN_P) and the initial abstraction (Ia). CN_P was calculated using Equation (5) [75], and Ia was obtained by applying Equations (6) and (7) [76,77]. However, to calculate Ia, the factor of the original equation (Ia = 0.2 × S) was changed to 0.05, because hydrological studies executed in small catchments and at individual hillslopes recommended this value [78,79].

$$C{N}_P=\frac{1}{A_T}\left(A_{VC1}\times C{N}_1+A_{VC2}\times C{N}_2+\cdots +A_{VCN}\times C{N}_N\right)$$

(5)

where CN_P is the weighted CN for all basins; A_T is the total area of the watershed (km²); A_VC is the area covered by a specific vegetation type [km²]; and CN_1-N is the curve number value of a specific vegetation type.

$$S=254\times \left(\frac{100}{C{N}_P}-1\right)$$

(6)

$$I_a=0.05\times S$$

(7)

where CN_P is the weighted CN and S is the maximum initial infiltration (mm).

To calculate P_e and transform it into runoff, the SCS unit hydrograph method (transform method) was applied. This method requires one additional parameter, namely, lag time (T_LAG) [80], which can be calculated using Equation (8) [81]:

$$T_{LAG}=2.587\times L^{0.8}\times \frac{{\left(\frac{1000}{C{N}_P}-9\right)}^{0.7}}{\left(1900\times {\overline{S}}^{0.5}\right)}$$

(8)

where T_LAG is the lag time (hours); L is the length of the main river (m); CN_P is the weighted CN; and $\overline{S}$ is the average watershed slope (%).

Then, the base flow was determined as described in the previous section (RPM model). To adjust the HEC-HMS model, the automatic optimization module of the HEC-HMS 4.3 software was used, in addition to the manual calibration method described before (trial and error method [46]). However, the most sensitive parameters for hydrological response in the HEC-HMS model are CN, Ia, and T_LAG [45], which were adapted, and finally, the same adjusted values were applied to all three precipitation events, to facilitate an early warning tool for El Carmen watershed.

2.4. Validation

Two types of validation were used to evaluate results and efficiency of both model simulations: (i) graphical, which implies a visual analysis of the results by comparing graphs between measured and simulated flows, and (ii) numerical, which measures the model adjustment degree using different indices [82]. The indices applied for the present study were coefficient of determination (R²) and Nash–Sutcliffe efficiency coefficients (NSE) [83]. R² measures the degree of correlation between two variables (measured flows vs. simulated flows). The range of R² varies between 0 and 1, where 1 indicates a perfect correlation and 0 represents no correlation between variables. NSE is the most widely used tool to evaluate the efficiency of hydrological simulations due to its robustness and effectiveness (e.g., [20,84]). NSE coefficients range between −∞ and 1, where 1 indicates optimal simulation and negative values indicates low adjustment between measured and simulated flows [10].

3. Results

The determination of morphometric variables of the El Carmen watershed, based on DEM (100 × 100 m), resulted in an area of 4.70 km² and a perimeter of 12.17 km. The calculated mean slope of the main stream was 26%, while the average slope of the watershed was 70%. The calculated watershed form factor (IF) and Gravelius coefficient (KC) were 0.55 and 1.61, respectively, which indicated that the watershed is not susceptible to great-magnitude floods because of its elongation [59,61]. However, the calculated concentration time (tc) of the watershed was relatively short (approximately 20 minutes), due to very steep slopes, its extension, and the moderate to low soil permeability caused by high clay contents (hydrological groups B and C) [63,85]. Therefore, the hydrological response of the El Carmen watershed to precipitation events is generally fast, but extremely high runoff magnitudes cannot be expected [86,87].

3.1. RPM Model

Figure 3a shows that soils of hydrological group B were found mainly at the lower part, covering about 10% of the watershed. These soils have loamy to loamy-clay texture, and therefore have moderate permeability [63]. The rest of the watershed was classified as hydrological group C (90%), because of the higher clay content (loamy-clay), indicating lower permeability [7]. Predominant land cover in the watershed (Figure 3b) is natural vegetation (tropical mountain forest and paramo [88]), especially at mid and higher elevations. These woody vegetation types create innumerable channels and small compartments in the subsoil which improve infiltration rates [24]. In addition, the soils at higher elevations have greater content of organic matter, due to the colder climate conditions, which delays decomposition [89,90], but improves soil water retention capacity and, subsequently, hydrological regulation [17]. Therefore, it can be expected that these parts mainly contribute to the discharge by means of subsurface and lateral flows (in the slope direction), except during saturated soil conditions (abundant former rainfalls) when principally superficial flows are produced [24,91]. In the lower zone of the watershed, hydrological group B soils are common, due to erosion processes at steeper slopes and the accumulation of the material at less inclined areas near the valley bottom [90]. However, the notable anthropogenic impact, particularly the clearing of the natural vegetation to create pasture land for livestock breeding, leads to lower soil depths (Entisols) and soil compaction (cattle) [92], which provokes changes in the hydro–physical soil properties. Therefore, superficial flows should be the main contributors to runoff at these areas [7].

The combination of soil type, hydrological condition (permeability), and vegetation cover allowed for the determination of CN at each grid cell (spatial CN; Figure 3c). The calculated CN values varied between 63 and 79, which indicates a moderate to low soil permeability and a rapid response to precipitation events (runoff formation), especially for the lower parts. The highest CN values were calculated for anthropogenic intervention areas (pasture), where lower soil depths and soil compaction reduced the infiltration capacities (CN of approximately 79). However, 80% of the watershed surface showed CN values of approximately 70 because of the natural forest cover, which improves the infiltration capacities. At the ridges, where paramo predominates, the conditions were even better (CN of approximately 63), mainly due to a thicker organic layer, which increases the soil porosity [92,93].

The precipitation distribution and amounts of the three selected events are shown in Figure 4. Storm direction was different during the events (black arrows), in which easterly atmospheric flows were determined for EVENT 1 and 3 (Figure 4a,c), whereas westerly flows prevailed for EVENT 2 (Figure 4b). During EVENTS 1 and 3, humidity was transported by advection from the Amazon Basin, but only moderate to low rainfall amounts were observed (maximum 10–25 mm). This is due to the local topography, because the eastern cordillera of the Andes forms a climatic barrier under these atmospheric conditions (easterly winds), which impedes the humidity transport into this basin [30]. Consequently, maximum precipitation amounts are shown at the eastern ridges (east–north-east), decreasing to the west (barrier effect; Figure 4a,c). However, the amount of humidity that reaches the watershed also depends on wind speed (average: EVENT 1 = 3.7 m/s; EVENT 3 = 7.0 m/s), because stronger winds make the humidity transport from the Amazon Basin more effective, for which reason precipitation amounts are higher during EVENT 3 [4,22].

EVENT 2 presented strong convective rainfall, especially in the lower part of this watershed (Figure 4b). As mentioned before, the formation of convective storms is common during the rainy season (December–April) over the Inter-Andean valleys, such as the valley of the city of Loja, especially during episodes of westerly winds [39]. During EVENT 2, the convective air was additionally forced orographically (mechanical lifting) by the eastern mountain chain, for which reason maximum precipitations (44 mm) were observed on the lower part of the watershed decreasing with altitude (orographic precipitation [5]). Under this atmospheric condition, the eastern mountain ridge of El Carmen watershed also forms a climatic barrier, reducing the humidity content of air towards higher elevations and preventing the humidity transport further to the east (Amazon Basin).

By means of the high-resolution precipitation maps (100 × 100 m, every 10 minutes) generated by means of the radar data, the spatial CN (Figure 4c), and the morphometric parameters, the RPM model was executed. After the first model run, without any parameter calibration (optimization), the results were compared to the measured discharge at the hydrological station (Carmen_Hid), which indicated good correlations between the simulated and the observed discharges (R²: 0.65–0.81; NSE: 0.59–0.71). The lowest correlation was obtained for EVENT 1, whereas the highest correlation was found for EVENT 2, which was due to the higher precipitation amounts, increasing surface runoff, which facilitates model simulation [82,94]. However, to improve the simulation efficiency, the parameters tc and fx were calibrated manually. The parameter tc was finally adjusted to 35 minutes, because of the predominant natural vegetation in the catchment, especially at the middle and upper parts, which improves infiltration and hydrological regulation [95]. The adjustment of tc also changed the parameter K, which increased from 0.26 to 0.46. Furthermore, the fx-parameter was optimized, increasing it slightly from 0.78 to 0.90, which generates greater surface runoff. This adjustment was realized because soil water contents are generally high (near saturation) in tropical mountain regions, especially under natural vegetation cover (primary forest and paramo), caused by the perhumid climate and the high storage capacity of the soils [43].

The final RPM simulations for the three selected events are presented in Figure 5, where the simulated hydrographs (red dotted line) are compared to the measured discharges at the Carmen_Hid station (black line). In general, discharge increases were simulated precisely for all three events. However, peak flows were better predicted for EVENT 2 and EVENT 3, where magnitudes were only slightly underestimated by about 8% (measured discharge: EVENT 2 = 2.54 m³ s⁻¹, EVENT 3 = 1.80 m³ s⁻¹; RPM: EVENT 2 = 2.37 m³ s⁻¹, EVENT 3 = 1.65 m³ s⁻¹), with timely synchronization. In contrast, peak discharge for EVENT 1 was simulated less accurately for both, magnitude (measured discharge: EVENT 1 = 1.64 m³ s⁻¹; RPM: EVENT 1 = 1.37 m³ s⁻¹) and time, which indicated that soil saturation (adjusted fx-factor) was still underestimated for this event. Usually it should be expected that low precipitation amounts during a rainfall event can, at least, partially infiltrate into the soils if natural vegetation cover exists, because of the improved infiltration and storage capacities [43,95]. However, the damp and cold climate in the study area, especially at higher elevations [96], causes soil saturation during most of the year, reducing infiltration rates, which are additionally reinforced by soil texture (hydrological groups B and C). Therefore, the simulated discharge peak for EVENT 1 was delayed and lower (error: 16%). In general, if precipitation amounts are higher, direct surface runoff must be assumed, because of soil saturation or with regard to the specific infiltration rates of a soil type, which is why peak flow simulation for EVENT 2 and 3 were more accurate [94].

The discharge recessions were simulated to be faster than observed for all events, although an exponential recession was presumed [71], which generally provides better results because soil water content after a precipitation event is considered [97]. As [98] and [99] clarified, an HM which adequately simulates discharge recessions does not exist, because precise information respective to soil water content after a specific precipitation event is necessary.

Nevertheless, the manual parameter calibration (t_c and fx), whose final values were used for all three selected events, improved the simulations notably. After the optimization, good to very good correlations between the simulated and the observed discharges were obtained (R²: 0.85–0.92; NSE: 0.77–0.82). This indicated that the applied DHM (RPM), in combination with high-resolution meteorological radar data, is suitable for hydrological forecasts in high tropical mountain watersheds.

3.2. HEC-HMS Model

Semi-distributed HMs, such as HEC-HMS, divide the watershed into subunits and apply specific mean values for each [18,100]. The calculation of the mean values for El Carmen watershed resulted in an CN_P of 70.20, an Ia of 4.90, and a T_LAG of 24.18. As mentioned before, for this HM the precipitation distribution and amounts were calculated by means of the available station data, applying the Thiessen polygon method [66]. This method divides the watershed into different sectors and assigns the measured precipitation amounts of the respective station to its whole area of influence (Carmen_Plu: EVENT 1: 10 mm, EVENT 2: 29 mm, EVENT 3: 25 mm; Carmen_Met: EVENT 1:7 mm, EVENT 2: 43 mm, EVENT 3: 20 mm; Figure 6).

By means of these input parameters, the HEC-HMS model was executed, first without calibration. The results showed only low to satisfactory simulations with an R² between 0.38 and 0.67, and an NSE between 0.05 and 0.56. Similar to the DHM applied (RPM), EVENT 1 was less accurate, while EVENT 2 showed the best results, which could be traced back to the higher precipitation amounts during this event, because higher rainfall amounts generally lead to higher runoff and faster hydrological response [94]. However, the HEC-HMS software includes an automatic optimization tool, which tries to optimize the input parameters for each event selected. Nonetheless, to provide an early warning tool for the El Carmen watershed, the input parameters were afterwards calibrated manually, and the same parameter values applied for all selected events. Finally, CN_P was slightly increased from 70.20 to 70.40, and Ia from 4.70 to 4.90, whereas T_LAG was reduced from 24.18 to 23.00, which prevented pronounced over- or under-estimation of the discharge peaks for all selected events. The optimized simulations are shown in Figure 7.

The discharge increase was overestimated for EVENT 1 and EVENT 2, which might be caused by the reduced T_LAG, which creates faster hydrological response [80]. However, discharge peaks were timely, but clearly underestimated for EVENT 1 (measured discharge: EVENT 1 = 1.64 m³ s⁻¹; HEC-HMS: EVENT 1 = 1.10 m³ s⁻¹) and overestimated for EVENT 2 (measured discharge: EVENT 2 = 2.54 m³ s⁻¹; HEC-HMS: EVENT 2 = 3.05 m³ s⁻¹). Although the discharge increase was better simulated for EVENT 3, the discharge peak was also clearly underestimated (measured discharge: EVENT 3 = 1.80 m³ s⁻¹; HEC-HMS: EVENT 3 = 1.21 m³ s⁻¹), which resulted in an overall error for the modeled discharge peaks between 20% (EVENT 2) and 30% (EVENT 1 and 3). As for the RPM model, the discharge recessions were simulated to be faster than observed for all selected events, because an exponential recession was presumed [97]. However, the measured discharge regression in the study catchment is slower, due to the predominant natural vegetation cover, which improves the hydrological regulation of the watershed [95].

Optimization of this model improved the results notably, and finally, satisfactory to good correlations were obtained when comparing simulated and measured discharges (R²: 0.68–0.86; NSE: 0.26–0.78). Better correlations were obtained for EVENT 2 and 3, when precipitation amounts were higher, which implies that the HEC-HMS model is especially appropriate for forecasting flash floods in poor instrumented watersheds.

4. Discussion

Comparing the two HMs applied in this study, it was found that the distributed RPM model in combination with high-resolution radar data is more precise than the semi-distributed HEC-HMS model using the available station data, because the spatiotemporal distribution of all input parameters is considered [19]. The assignment of mean values for the whole watershed (HEC-HMS) leads to higher inaccuracies in the simulation results, specifically if precipitation data is sparse and inexact, or if precipitation amounts are low [32]. This is illustrated in

Table 1. Statistical analysis and used precipitation input data for the selected events after HM optimization, comparing simulated and measured discharges, including the coefficient of determination (R²) and the Nash–Sutcliffe efficiency coefficient (NSE).

HM	RPM			HEC-HMS
EVENT	R2	NSE	Precipitation Input	R2	NSE	Precipitation Input
1	0.85	0.80	Radar observations	0.68	0.26	Station data
2	0.92	0.82	Radar observations	0.86	0.78	Station data
3	0.85	0.77	Radar observations	0.81	0.67	Station data

, where the efficiencies after model calibration and optimization of both HMs (RPM vs. HEC-HMS) are shown, comparing the simulated discharges with the measured discharges.

To depict rainfall amounts and their distribution accurately inside a watershed by means of meteorological station data, dense networks are necessary. However, installation and operation are often complicated in developing countries, especially in tropical high mountains, due to the high costs and the difficulties in access [12,23,88]. Even if sufficient meteorological station data is available, the rainfall information must be interpolated or extrapolated, which leads to inaccuracies in the generated precipitation maps. This is confirmed by [20], who applied the HEC-HMS model and the DHM HBV-N-D in the San Francisco catchment, one of the best instrumented watersheds in southern Ecuador. Based on interpolated meteorological station data (inverse distance weighting; IDW), they obtained an NSE of 0.60 (HEC-HMS) and 0.73 (HBV-N-D), which is comparable to the HEC-HMS simulation efficiencies of the present study (Thiessen polygon method). Another study, executed in the Zhurucay micro-catchment in Ecuador, applied the HBV-light semi-distributed HM to analyze the effects of rainfall estimation based on rain gauge measurement respective to model efficiency. In addition to lower model accuracies compared to this study (NSE between 0.49 and 0.60), they concluded that a higher number of rain gauge stations significantly improves HM simulations [101]. The higher model accuracies of the applied semi-distributed HM in this study (HEC-HMS) can be explained by the smaller size of the El Carmen watershed, where at least two rain gauges are installed. For smaller watersheds, fewer point measurements are necessary to accurately estimate precipitation distribution and amounts. However, even with this rain gauge information, the precipitation distribution in the study area during the three selected events could not be depicted precisely, because notable under- or over-estimations of the discharge peaks were obtained (20%–30%), although they were accurate with respect to time.

Therefore, for precise discharge simulations by means of HMs, dense meteorological station networks are necessary, which is often impossible in tropical mountain catchments, which is why meteorological radar has an advantage. Despite the high installation and operation cost of sophisticated radar systems (C-Band or S-Band; e.g., [102,103]), recently, cost-effective X-band technologies have become available (e.g., [29,104]), which are also affordable for developing countries. The advantages of meteorological radars over point measurements (rain gauges) are that these systems detect the precipitation over larger areas in high spatiotemporal resolution, and for radar data calibration only isolated ground information at strategic points is required [22,30]. This permits the generation of precise precipitation maps for the whole area of radar coverage, even for unequipped or isolated watersheds [39]. This implies that river discharge in a tropical mountain catchment can only be modeled satisfactorily by means of meteorological radar data.

Nonetheless, for a precise discharge simulation, an adequate spatial and temporal resolution of the radar observations, and, consequently, of the deviated precipitation maps, is also necessary [26,27]. This is confirmed by the study of [84], who stated that the HM efficiency depends on the spatiotemporal resolution of the radar images. They simulated the discharge of the Amacuzac River basin (approximately 9000 km²) using daily information of a meteorological radar installed at the “Cerro Catedral” (México) with a spatial resolution of 27 × 27 km and applying the DHM CEQUEAU [105]. However, their model efficiency was only satisfactory with an NSE of 0.6, because of the low resolution of the radar images. Better model efficiencies were obtained by [46], who used hourly precipitation maps derived from the Next Generation Radar (NEXRAD; e.g., [106]) of the United States with a spatial resolution of 1 × 1km to estimate the discharge peaks in the Escondido River basin (3240 km²) in northern Mexico. They also applied a DHM, particularly the semi-distributed HEC-HMS, which was converted to a DHM by means of the Modified Clark transformation method [3]. Due to the relatively high image resolution with respect to the catchment area, they obtained an NSE of 0.96 for both simulated events.

In summary, meteorological radars provide accurate area-wide precipitation information in high spatiotemporal resolution for HM simulations, as well as for isolated or unequipped tropical mountain catchments. However, the spatiotemporal resolution of the radar images must be fine enough to accurately predict the river discharges. Depending on the watershed size, resolutions of 1 × 1 km or higher are required, especially for small watersheds in high tropical mountains where strong altitudinal gradients prevail [19].

5. Conclusions

For hydrological, climatological, and ecological investigations, as well as for water resource management (e.g., drinking water supply for the population), high-resolution radar observations are the best alternative to accurately determine precipitation amounts and their distribution in isolated or poorly equipped watersheds, because only ground information at strategic points is necessary to generate precise precipitation maps for the whole area of radar coverage. At the present time, cost-effective meteorological radar systems (mainly X-Band) are affordable for developing countries, which do not only provide accurate areawide precipitation information but also reduce implementation and maintenance costs of dense meteorological station networks. Furthermore, isolated storms are often not detected by traditional meteorological station networks, especially in tropical mountain regions, due to the high spatiotemporal precipitation variability, which is why radar systems are the better option.

This was shown by the present study, in which the application of a DHM in combination with high-resolution radar observations significantly improved the hydrological simulations. The DHM applied (RPM) resulted in a very good model efficiency for the El Carmen watershed, whereas the semi-distributed HM (HEC-HMS) only provided acceptable to satisfactory simulations, in which the latter holds particularly true for precipitation events with intense rainfall (here: EVENT 2). Therefore, only if sufficient ground measurements inside a watershed are available, can semi-distributed HMs be used as a tool for the implementation of an early warning system, especially for forecasting flash floods, because discharge peaks are better predicted in terms of magnitude and time for extreme precipitation events.