Study area

Six basins were selected, located between 3° 30’ and 5° 8’ S and 79° 10’ and 80° 29’ W, which belong to the bi-national basins Catamayo-Chira and Puyango-Tumbes, located along the border between Ecuador and Peru.

The Catamayo-Chira Binational Basin has an extension of 17 199.19 km² and 817 968 inhabitants. The basin begins in the mountains of the Continental Divide (through Ecuador) and it ends in the Pacific Ocean (Peru). It goes through the Andes Mountain range downwards to the coast. It features tropical weather, diverse ecosystems, and various administrative/political systems. The geography of the basin is abrupt, with altitudinal variations between 0 m and 3700 masl and there are 11 life zones from tropical desert to mountain rainforest. The mean annual precipitation is about 800 mm, varying between 10 mm in the low-altitude areas to 1000 mm in the basin headwaters. Besides, 14 % of the basin surface is covered by humid forest, 41 % by dry forest, 30 % by pasture, 10 % by crops and 5 % by other uses.

The Puyango-Tumbes basin covers 5400 km² and it includes the provinces of El Oro and Loja in Ecuador and the department of Tumbes in Peru. The upper basin consists of mountainous terrain with steep slopes, which elevations between 500 m and 3700 m, and encloses wilderness areas, natural forest and crops. The upper basin is characterized by an important mining activity, especially in Calera and Amarillo. The lower basin has extensive plains and it is heavily cultivated, with rice fields; there are also large areas of dry forests, grasslands, and agricultural crops. The lower basin in the Peruvian portion suffers from frequent floods. Precipitation in the basin ranges from 200 mm to 1150 mm per year, with average temperatures between 13 °C and 25 °C. Three main tributaries compose the basin: the Calera river, the Yellow River, and the Pindo river. On average, the total volume of water produced by the basin (average annual mass) is approximately 3,400 million m³.

The availability of hydrometric information and geographic characteristics determined the choice of six sub-basins: the basins from the Amarillo and Pindo rivers pertaining to the PuyangoTumbes basin, and the basins from the Catamayo, Macará, Quiroz and Alamor rivers pertaining to the Catamayo-Chira basin (Table 1). The geographic location of the basins is shown in Figure 1.

Temez Model

The Temez Model (^{Temez, 1977}; ^{Estrela, 1999}; ^{Murillo and Navarro, 2011}) is a lumped hydrologic model that assumes that the soil profile is divided into an unsaturated upper zone and a saturated lower zone. The model behavior is similar to that of an underground reservoir that drains into the surface network.

Rainfall (P) is divided into evapotranspiration (ET) and excess (T). The excess is decomposed into a fraction that flows over the surface (E) (surface runoff), and into another fraction that infiltrates into the ground (I). The first fraction evacuates through the channel within the established period of time, whereas the second will remain in the underground reservoirs to be released at a later time. The excess is calculated through the following expression:

in wich:

P _t is the precipitation in the period between the instant i-1 to the instant i (mm); P ₀ is the runoff threshold that defines the height of precipitation, under which no runoff is produced (mm); T _t is the excess in the period between the instant i-1 to the instant i (mm); H _max is the maximum capacity of soil moisture (mm); H _i-1 is the soil moisture at the instant i-1 (mm); EP is the potential evapotranspiration from the instant i-1 (mm); and C is a parameter from the model.

The soil moisture, Hi, at the end of the period will result:

An actual evapotranspiration ER _t (mm) is produced and it is equal to:

Equation 6 shows that all the water available can be evapotranspired, with the potential evapotranspiration as an upper limit.

The model adopts an infiltration amount (I _t ) as a function of the excess, T _t and of the maximum infiltration parameter, I _max .

will be in mm.

The infiltration increases with the excess, being asymptotic for its higher values, having as a limit the I _max value.

This infiltration, I _i becomes recharge to the aquifer, R _i , whereas the rest of the excess (E _i = T _i I _i ) will become surface water. The model assumes that the time step in the unsaturated zone is lower than the simulation time step.

The draining of the aquifer is modeled by an exponential function of the following type:

where Q _i is the discharge at the instant, i, α is the coefficient of the discharge branch of the aquifer that depends on the particular conditions of the basin being studied and t the time interval between instants i-1 and i.

The relation between the discharge Q _i and the volume V _i , stored in the aquifer is equal to:

The discharge for the infiltration is assumed to be concentrated in half of the time step; therefore, the law of the groundwater flow is:

in which R _i is the discharge to the aquifer in the period between i-1 and i, coinciding with infiltration, I _i .

The groundwater contribution throughout the period, A _SUBi .

The total contribution, A _Ti , is given by surface runoff (T _i I _i ) and the groundwater contribution:

There are four parameters of the model: Hmax, which is the maximum capacity of the soil moisture; Imax, maximum capacity of infiltration; C is the parameter of excess; a is the coefficient of the discharge of the aquifer. These four parameters are particular to each basin and are subject to calibration.

Soil and water assessment tool (SWAT)

SWAT is a continuous simulation model used to predict the impact of land management practices in the production of water, sediments, and nutrients in a watershed (^{DiLuzio et al., 2002}). Besides, it is a semi-distributed continuous model based on the equation of water balance in the land profile and it simulates the precipitation, infiltration, runoff, evapotranspiration, lateral flow, and percolation processes. The equation of water balance is:

where SW _t is the final soil water content (mm H2O), SW ₀ is the initial soil water content on day i (mm H2O), t is the time (days), Rday is the amount of precipitation on day i (mm H2O), Q _surf is the amount of surface runoff on day i (mm H2O), E _a is the amount of evapotranspiration on day i (mm H2O), w _seep is the amount of water entering the vadose zone from the soil profile on day i (mm H2O), and Q _gw is the amount of return flow on day i (mm H2O).

The SWAT model accomplishes a topographic division of the basin into sub-basins based on threshold of area; in turn, the sub-basins are subdivided into one or more homogeneous hydrological response units (HRU) that represent a unique combination of soil type and land use. The response to each HRU in water, sediments, nutrients, and pesticides are individually determined; then, they are added to a sub-basin level and moved toward the basin exit through its river network (^{Bouraoui et al., 2005}).

The surface runoff is estimated from daily precipitation data using the curve number methodology or Green-Ampt equation (^{Neitsch et al., 2002}). The evapotranspiration is determined applying the methodologies proposed by Hargreaves, PriestleyTaylor and Penman-Monteith (^{Neitsch et al., 2002}). A kinematic reservoir that considers variations in hydraulic conductivity, prevailing slope, and soil moisture is used to predict the lateral flow in each one of its layers. Groundwater is divided in two aquifer systems: 1) a non-confined aquifer, which is shallow and contributes to the backflow; 2) a confined aquifer, which is deep and disconnected from the system, unless groundwater pumping is being considered (^{Bouraoui, et al., 2005}). The sedimentation rate is estimated by applying the modified universal equation of soil loss (^{Neitsch et al., 2002}), using for that purpose the surface runoff, peak flow rate, erodibility of soil, slope length, its inclination, crop factor, and the management practice in the zone. ^{Neitsch et al. (2002)} explained the theoretical basis of the SWAT model.

Implementation, calibration, and validation of the model

For each basin, the lumped modeling was performed considering as input data the average spatial rainfall in the basin and the monthly potential evapotranspiration. For calibration purposes, the monthly average flows recorded in the downstream station of each basin included in Table 1 were considered. The precipitation data in 24 stations in Ecuador and in 19 in Peru stations as well as the temperature data in 10 stations in Ecuador and 4 stations in Peru were considered. The correspondence of the data series was verified through double-mass curve analysis. By regression analysis between stations that present geographical closeness, the same climatic regime, and the compatibility of records, data were homogenized to the common period 1970-2000 in the stations of the Catamayo-Chira basin, and 1965-1995 in the station of the Puyango-Tumbes basin. The Thornthwaite method (^{Thornthwaite and Mather, 1957}) was used to calculate potential evapotranspiration in each one of the 14 stations for which the temperature data was recorded.

The average precipitation and evapotranspiration for each one of the basins were calculated by Thiessen polygons (Figure 2).

In each one of the selected basins, manual calibration of four parameters (C, Hma, Imax, a) of the Temez model was carried out based on successive approximations. For this calibration, 70% of the measured records were used; the flows calculated by the model were compared to the recorded flows in the stations included in Table 1. In the calibration process, the parameters values were changed to no more than 12%. Once the model parameters were calibrated, they were validated, generating flows for a period covering the remaining 30% of the historical record. The goodness of fit between the calculated and observed values was measured by calculating the correlation coefficient (R), the Nash-Sutcliffe model efficiency coefficient (EF) (^{Nash and Sutcliffe, 1970}) and the root-mean-square error (RMSE). The correlation coefficient will measure the fit of the observed and calculated data to a straight line, the EF efficiency will measure the 1 to 1 relation of the observed data, and the RMSE will allow the evaluation of the average difference between observed and calculated values.

The application of the Temez model was performed by using the CHAC software which was developed by CEDEX of Spain. The Temez model may be found in http://hercules.cedex.es/ Chac/.

For comparison purposes, the SWAT model was implemented in each one of the sub-basins. A DEM from the Shuttle Radar Topography Mission (SRTM) of the study area was used for the delimitation of contributing basins and the calculation of the different morphometric parameters required by the model. Also, the information regarding the type of soil was gathered, and the soil parameters were determined by laboratory tests. Land use was extracted from a LANDSAT7 ETM+ image of October 2, 2011 through a supervised classification; in this extraction the minimum spectral angle (^{Chuvieco, 2002}; ^{Richards and Jia, 2006}) was applied. The land use parameters were taken from the SWAT model database, with minor modifications. An existing soil type map (^{Valarezo, 2007}) was validated and adapted for SWAT model implementation. The validation was done through perforations and pits in each of the edaphologic units identified. The infiltration test was performed in situ. Unaltered soil samples and general characteristics of soil horizons were taken. With all the collected data and the laboratory tests, the soil type parameters were determined for each of the edaphologic units identified.

The combination of soil type and land use allowed the determination of hydrological response units (HRU). The soil type information was used to characterize the non-confined aquifer. Additionally, precipitation and daily temperature data, as well as temperature parameters, solar radiation, wind speed, monthly relative humidity for each one of the stations were assembled to calculate evapotranspiration using the Penman Monteith method. The precipitation missing data were estimated by orthogonal correlation analysis between geographically neighboring stations with similar climatic conditions. In the case of temperature, an elevation-temperature equation was determined. A similar process was used for other climatic variables.

Before calibration, an analysis of model sensitivity was made by systematically changing the values of the model parameters and watching for the effect that such changes have on the obtained results, by selecting a reduced number of sensitive parameters of a higher relevance in the calculation and optimizing them in the calibration phase. The sensitive analysis was made using the module that the SWAT model has for this purpose, which combines the Latin hypercube sampling method with the one factor at a time (OAT) simulation (van Griensven et al. 2006). The Latin hypercube sampling method, in opposition to the conventional Monte Carlo method, makes a stratified sampling between the rank of possible values of each parameter; moreover, the OAT simulation assures that changes in the model can be attributed to an change made to the entrance variable in each simulation (van Griensven et al. 2006).

The SWAT model was calibrated manually, adjusting the values of the most sensitive parameters changing its initial values to no more than 12%, using 70% of the gathered records, and optimizing R, EF, and RMSE. The SWAT model parameters about soil loss were not considered in the calibration process, because the Temez model does not reproduce sediment losses, so the comparison between Temez model and SWAT model was made only with flow data. The validation was performed with the 30 % remaining from the record. Finally, a graphic comparison of calculated and observed flows was carried out. See ^{Oñate-Valdivieso and Bosque (2014)} for additional information on the implementation of the SWAT model for basins study.

Calibration of parameters for the Temez model

The calibrated parameters for the Temez model in each one of the basins are shown in Table 2.

The maximum capacity of soil moisture (H_max) allows the definition of the threshold runoff (P_o), i.e., the rainfall below which there is no runoff; thus, the excess (T_i) is defined and it is composed of the runoff and infiltration into the aquifers. Table 2 shows that the values of H_max surpass 150 mm in most of the cases, reaching a maximum value of 400 mm. The values of this parameter show a high storage capacity of the soils in the surveyed zone, which experience a considerable water deficit during most of the year. This situation can be found at the sub-basin of the Alamor river. The opposite can be observed at the subbasin of the Amarillo river, which is dry in its lower section but shows higher levels of rainfall throughout the year. For this reason, soils are more humid and have a lower maximum capacity of infiltration.

The parameter of excess C is used to determine the runoff threshold P _o and acts as a weighting factor of the difference between the maximum humidity of the soil and its humidity at a given moment. The values of parameter C (Table) 2 are lower than 1; thus, controlling a possible sub-estimation of the excess produced when considering relatively high humidity factors.

The actual infiltration I is directly proportional to the maximum infiltration I _max , the former being a fraction of the latter. Additionally, the values of this parameter, which are included in Table 2, show a level of maximum infiltration corresponding with H _max . Since this is an arid zone, high levels of infiltration may be expected.

The coefficient α reaches average values between 0.01 and 0.08. This parameter is related to the contributions of groundwater flow to the base flow of the basins analyzed. The higher values can be found in the basins with a larger area, since these basins extend across zones with different levels of contribution of groundwater. In higher zones, which are usually more humid and exhibit better conditions of plant growth, the contribution of groundwater is higher, unlike lower zones, where the opposite occurs. The calibrated parameter in the larger basins is representative of such conditions.

Evaluation of the Temez model

The parameters used to determine the goodness of fit between the observed and calculated flows for the period of validations are included in Table 3. The graphical relation between observed and calculated flows using the Temez model is shown in Figure 3.

In smaller basins (Amarillo, Pindo and Alamor), the model recreates the flows of the dry season better than those of the rainy season (Figure 3). As the size of the basin increases, an overestimation of the flows may be observed, and the correspondence between observed and calculated flows almost disappears. The overestimation of flows is due to the calculation of the flow as the product of the contribution (mm per month) multiplied by the area of the basin assuming that the levels of rainfall and evapotranspiration are uniform across the basin. This does not necessarily apply to medium-sized or large basins.

The Amarillo river basin shows a correlation coefficient of 0.64 and an efficiency EF of 0.63 (Table 3). These values confirm the existence of an acceptable adjustment between the observed and calculated values. The average quadratic error is 7.32 m³ s^-1, which, if compared to the average flow recorded (14.32 m³ s^-1), is excessive. In Figure 3, it is shown that the model does not accurately recreate peak flows, and the flows of the dry season are overestimated.

In the Pindo river basin, which approximately doubles the size of the Amarillo river basin, there was a correlation coefficient of 0.63 between the observed and calculated data, which is barely acceptable; besides, EF was 0.49, which is lower than that observed in the Amarillo river basin. The average quadratic error is 14.47 m³ s^-1, which, if compared to the average flow of the Pindo river (23.43 m³ s^-1), turns out to be excessive. The model overestimates flows. As shown in Figure 3, the best adjustment is produced in the flows of the dry season, and the peak flows shows the worst adjustment.

In the Alamor river basin (Table 3) there was a relatively good correlation between calculated and observed values (R=0.72); however, there is a low value of EF (0.35), with an average quadratic error of 9.11 m³ s^-1 that surpasses the average observed flow of 7.3 m³ s^-1. It could be expected that with a high correlation coefficient RMSE should be lower, but as shown in Figure 3, the model recreates the trend; therefore, the value of the correlation coefficient is satisfactory. In addition, the recorded peak flows are overestimated to a great extent by the model and this has an influence on RMSE. The EF exhibits a more realistic value of 0.35, which, although it is not optimal (EF=1), it gives an idea of the model’s performance, which is, to a certain extent, acceptable.

The Catamayo river basin station (Table 3) presents the worst results, which means that there is no correlation between observed and calculated flows. A negative value of the Nash-Sutcliffe efficiency is obtained; it can be considered that the mean of the observed values is a better predictor than the model itself.

The Quiroz river basin shows acceptable values of the correlation coefficient (R=0.69) and the NashSutcliffe efficiency (EF=0.62), although the excessive average quadratic error, which duplicates the average observed flow, shows the overestimation of the calculated flows.

The Macara river basin has the best correlation coefficient (R=0.82); however, the low value of the Nash-Sutcliffe efficiency (EF=0.27) implies that the model cannot recreate the hydrological behavior of the basin. In addition, the average quadratic error almost doubles the average observed flow by overestimating flows in general.

The study area is mainly mountainous, with irregular rainfall and since the average number of rainfall station is only one every 270 km², and it is assumed that the existing network of meteorological stations in the area does not have sufficient density to accurately reflect the spatial variability of precipitation. Additionally, historical records of precipitation, temperature and especially flow, showed significant information gaps that, to some extent, point to problems in the management of meteorological and hydrological stations in the area.

Like precipitation, slope, soil type, and land use present a significant spatial variation. Therefore, considering a single lumped value, representative of each condition may be strictly valid only for very uniform and small basins.

Comparison with SWAT model

Meteorological data was the most difficult to achieve in the study area, because the number of weather stations is low, only 14. Information about wind speed, solar radiation and relative humidity was the scarcest.

The sensitivity analysis reported the 35 most sensitive parameters, which were classified in a scale of 35 levels and 1 corresponds to the most sensitive parameter. Table 4 shows the first 10 parameters reported by sensitivity analysis and their influence levels in the model.

As shown in Table 4, the model turned out to be sensitive to the curve number CN2 because this parameter is used to separate effective precipitation from total precipitation. Therefore, there is a direct effect on the runoff calculation; furthermore, the available water capacity SOL_AWC determines surface runoff since it causes higher or lower infiltration. The Baseflow recession constant ALPHA_ BF will regulate the contributions of groundwater to the flow. The moist soil albedo SOL_ALB and the soil evaporation compensation coefficient ESCO regulated losses by evapotranspiration. An important influence was observed in the average slope SLOPE and the surface runoff lag coefficient SURLAG, thus directly affecting the concentration time of each subbasin and varying the temporal occurrence of peak flows. The soil loss parameters were not considered

A graphical analysis of the relationship between the observed flow and calculated flow with the SWAT model is presented in Figure 4. As seen in that figure, the SWAT model provides a better representation of trends and better fit with the observed values, even when considering the peak flows of the wet season and the low flows in the dry season.

Comparative analysis of SWAT/ Temez model performance by validation parameters is summarized in Table 5. There it is observed that in all cases the SWAT model presents better correlation coefficient, reaching in most cases values higher that 80 %. The Nash-Sutcliffe efficiency is significantly higher when calculated with the values obtained with the SWAT model. Similarly, the RMSE was significantly smaller on values calculated with the SWAT model.

As a semi-distributed SWAT model, a better performance in all of basins may be expected. However, even though in the smaller basins, the Temez model performance is lower to that of the SWAT model, which could be acceptable considering the lower amount of information required for implementation of the Temez model.

Considering that in many places the lack of information is a reality, the Temez model may be a viable option when applied to a relatively small watershed. Additionally, the Temez model could be used as a basis for a semi-distributed model applicable to larger watersheds.