Establishment of the Baseline for the IWRM in the Ecuadorian Andean Basins: Land Use Change, Water Recharge, Meteorological Forecast and Hydrological Modeling

Abstract

This study was conducted in the Zamora Huayco (ZH) river basin, located in the inter-Andean region of southern Ecuador. The objective was to describe, through land use/land cover change (LUCC), the natural physical processes under current conditions and to project them to 2029. Moreover, temperature and precipitation forecasts were estimated to detail possible effects of climate change. Using remote sensing techniques, satellite images were processed to prepare a projection to 2029. Water recharge was estimated considering the effects of slope, groundcover, and soil texture. Flash floods were estimated using lumped models, concatenating the information to HEC RAS. Water availability was estimated with a semi-distributed hydrological model (SWAT). Precipitation and temperature data were forecasted using autoregressive and exponential smoothing models. Under the forecast, forest and shrub covers show a growth of 6.6%, water recharge projects an increase of 7.16%. Flood flows suffer a reduction of up to 16.54%, and the flow regime with a 90% of probability of exceedance is 1.85% (7.72 l/s) higher for 2029 than for the 2019 scenario, so an improvement in flow regulation is evident. Forecasts show an increase in average temperature of 0.11 °C and 15.63% in extreme rainfall by 2029. Therefore, intervention strategies in Andean basins should be supported by prospective studies that use these key variables of the system for an integrated management of water resources.

1. Introduction

Land use/land cover change (LUCC) impacts on various natural physical processes within watersheds have been extensively investigated worldwide [1,2]. The location of the Andean basins within intervention policies has also been the subject of analysis over the last several decades [3,4,5]. Changes and distribution of land use have an important influence on natural processes in a basin [1,2]. Water recharge, for example, shows sensitivity to different coverage due to the different infiltration rates that can occur under the same precipitation conditions [6,7], which is the reason why a change or alteration must be analyzed to evaluate its consequences. That is why water recharge was estimated in order to evaluate the effects that LUCC would have on it. In addition, the areas of the basin that favored it were spatially represented [8,9,10].

Another hydrological process that is influenced by LUCC are the flows associated with high intensity rains [11]; therefore, its analysis implies an element of judgment when making decisions at territorial level. The hydraulic simulation of flash floods was carried out to denote the impact of LUCC, generating information on flows and flood plains along the main channel [11] and giving technical support to the implementation of strategies associated with flood vulnerability reduction accompanied by land use planning policies and non-disruptive engineering practices [12].

Simulating the hydrological response of a basin under land use/land cover (LULC) scenarios is another key element for the integrated management of water resources [13], because it allows for the impact of LULC on the outlet flow of the basin [1,14,15]. This input is also the basis for the application of territorial policies and programs for the conservation of water sources.

Some studies have explored the relationship of LUCC with the hydrological response in basins. For example, for flow reproduction and calculation of sediment production in the Catamayo-Chira binational basin, the semi-distributed model called the Soil and Water Assessment Tool (SWAT) was implemented [16]. Similarly, the SWAT model was applied, obtaining satisfactory results at the Calumpang basin in Philippines [1].

In this study, the SWAT model is implemented for flow simulation, taking as scenarios: (a) a classified map of LULC in 2019, and (b) a projected map of LULC in 2029. With the generated hydrographs, flow duration curves (FDC) were prepared and the probability of exceedance of different flow values was analyzed for both scenarios.

Forecasted hydro-climatic information is scarce in the study region. In order to solve this need, the mean temperature and precipitation were forecasted through autoregressive and exponential smoothing models, widely used in annual series and meteorological variables [17,18,19,20]. Because these climatic data strongly influenced [21], their effects should be considered irreducible (even with the associated uncertainties), so intervention plans and adaptive strategies should be proposed to help mitigate these effects [22].

Our study focuses on the integrated water resources management (IWRM) in Andean basins. The ZH basin was used as a study case because of the interventions carried out there at the level of political institutions by decision-makers. Land use practices on the basin have been restricted in the last several years mainly due to municipal intervention policies for water sources preservation; a regional water fund operates in the area, ensuring that anthropogenic activities are not being carried out in areas considered as water recharge [23]. The basin has the also particularity of being located in a transition zone between a consolidated urban sector to a National Park [24].

The first results have been divided into two parts: this article, whose objective is to describe LUCC with some natural physical processes such as water recharge, flash floods, and water availability; temperature and precipitation forecasts have also been prepared to detail the possible effects of climate change to a horizon of 10 years. The second part is detailed in an article whose objective is to make a prospective analysis of water management in the ZH basin, including physical, socio-economic, and political-institutional variables, generating future scenarios and detailing strategies to reach a target horizon.

2. Materials and Methods

2.1. Study Area

The Zamora Huayco (ZH) river basin (3806 ha) is located in the inter-Andean region of the Loja province in southern Ecuador between the geographic coordinates 3°59′42″–4°04′03″ S and 79°11′54″–79°07′35″ W (Figure 1). The elevation of the basin ranges between 2120 and 3420 m asl, and its average slope is 0.65 m/m.

The basin’s climate is cold temperate mesothermal [25], characterized by an average annual temperature between 12 and 18 °C and average annual precipitation of 1047 mm. The wet season occurs from December to May and the dry season from June to November [26].

The basin has a predominantly forested vegetation. Since 1976, the forests of the basin have decreased by 19.3% [27]. The upper zone of the ZH basin is shared with the buffer zone of the Podocarpus National Park (PNP) [24]. Since the 1960s, natural vegetation near to PNP has been extensively removed to create pastures and farmland [27,28]. The main productive activities in the basin are agriculture and cattle raising [24,29]. In the ZH basin, there are two water catchments for potabilization that supply approximately 50% of the demand of the city of Loja with 450 l/s [30].

2.2. Land Use/Land Cover Change (LUCC)

Using remote sensing techniques, satellite images were processed to analyze LUCC. Three images were obtained: Landsat7 ETM (11 September 2009), Landsat8 OLI_TIRS (20 September 2017), and Sentinel 2B Level 1C (31 July 2019), which were ortho-rectified [31,32]. The satellite’s geometric coincidence with two temporarily invariable control points was verified, taking as reference an aerial photograph obtained from the SigTierras portal [33].

For the Landsat7 ETM and Landsat8 OLI_TIR images, the conversion to reflectance was performed to top of atmosphere (TOA) and to brightness temperature according to L. Congedo (2016). For the Sentinel 2B Level 1C image, the conversion was not performed as it is already scaled in TOA [32]. In the three images, the atmospheric correction Dark Object Subtraction (DOS) was performed [34]. In the three images, the topographic correction was also applied due to shading (CTS), using the Lambertian Cosine Model proposed by Teillet [35], for which it was necessary to determine the angle of incidence according to [36].

Radiometrically Terrain Corrected (RTC) data from DEM ALOS PALSAR, purchased from the Alaska Satellite Facility (ASF) Fairbanks, with a spatial resolution of 12.5 m, geocoded and radiometric corrected in its geographical extension [37], was used to carry out the CTS. LULC was obtained by supervised classification, using the method of the maximum likelihood classification algorithm based on Bayes’ theorem, of the three images previously processed. Three LULC categories of anthropogenic type and three of natural-physical type were classified.

For LUCC analysis, the procedure was analogous to that proposed by [38]. First, coverages obtained for the years 2009 and 2017 were analyzed; the quantitative change represented graphically in terms of gains and losses by categories was obtained.

As explanatory variables, maps were used showing the exchanges between the categories of grasslands and shrub vegetation, forest and bare soil, crops and grasslands, in addition to forest and crops. Additionally, a map containing the spatial trend of change was made, from forest to shrub vegetation, in order to generalize the pattern of change between these categories [39]. A ninth-degree polynomial order was considered.

Changes below 600 cells were ignored, in order to limit the transition model to 9 sub- transitions. The sub-transitions selected were: shrub vegetation to forest, bare soil to shrub vegetation, bare soil to grassland, bare soil to forest, grassland to shrub vegetation, grassland to forest, forest to shrub vegetation, forest to grassland, and grassland to crops.

The objective of the transition model is to create potential transition maps with a high degree of precision to execute the change model [38]. Transitions were grouped into a single sub-model; also, 6 explanatory variables were included.

The first explanatory variable was the shrub vegetation cover contained in the 2009 map. This model was included as a dynamic component that is recalculated over time during the course of the prediction. The other explanatory variables were included as static components, corresponding to the spatial trend maps and exchanges between categories previously generated.

Physical restrictions of any kind have not been considered, because they are not necessary in a simple prediction [38]. For the transition model, the Multi-Layer Perceptron (MLP) algorithm has been applied, based on the Backpropagation (BP) algorithm, which is a supervised training algorithm [40]. For the change model, Markov chains were applied, generating a matrix with the probability that each category of land cover will change to any other category [38]. A validation was performed with the map generated for the year 2019.

To verify the image projected to 2019, a confusion matrix was performed and the kappa coefficient was determined. The confusion matrix shows the relationship between the reference data (2019 classified map) and the data to be evaluated (2019 projected map), constructing a matrix comparison of the classes verifying the overall accuracy of the projection by relating the number of points correctly assigned to the total [36]. The kappa coefficient evaluates the degree of agreement between categorical variables and takes into account the coincidences by randomness and by decision criteria; that is, it shows the degree of agreement that exists above random [41].

Similar to what was performed by [7,42], following the same conditions used in the 2009–2017 projection, the land use map for the year 2029 was generated, based on the 2009 and 2019 maps. The 2019 classified map containing LULC is called Scenario (1), and the 2029 projected map containing LULC is called Scenario (2).

2.3. Hydric Recharge Estimation

The multi-year average hydric recharge of the ZH basin was estimated, following the methodology set forth by [8] and also considering the infiltration criteria of [10], taking into account multi-year average precipitation, evapotranspiration (ET), basic soil infiltration rate, LULC, and terrain relief. The multi-year average rainfall was determined using data from La Argelia meteorological station over a 30-year period (1985 to 2015). Potential evapotranspiration ($E{T}_P$) was determined through Hargreaves and Samani [43], who considered precipitation over the same time period (30 years).

As there is a correlation between the crop coefficient ($K_c$) and the NDVI [44,45,46,47], $K_c$ was deducted by applying the expression [44]:

$$K_c=1.84 NDV{I}^2-1.03 NDVI+0.42$$

(1)

NDVI was obtained from the Landsat8 OLI_TIRS image (20 September 2017). The multi-year average value of the calculated ETP was multiplied by the $K_c$ value to estimate the actual ET [48].

A slope map was generated using the DEM ALOS PALSAR RTC. This map was classified into 6 different categories (limits of <15, 15, 30, 50, 70, and >70%), assigning a coefficient to each. $\left(k_p\right)$, known as the surface runoff coefficient, indicates the fraction that infiltrates due to the effect of slope [9]. Its value was taken from the table proposed by [8].

To each land use scenario (2019 and 2029) was assigned a value, named the land use coefficient $\left(k_v\right)$, which represents the fraction that infiltrates due to the effect of the vegetation cover [9]. Its values were extracted from the table proposed by [8].

The infiltration fraction related to soil texture $\left(k{f}_c\right)$ was estimated using the following expression [10]:

$$k{f}_c=0.267\ln \left(fc\right)-0.000154\left(fc\right)-0.723$$

(2)

where $fc$ is the soil basic infiltration rate in mm/d. This was estimated using the Green-Ampt infiltration model [49], from soil data sampled in the ZH basin by [27].

According to [27], the most common soil texture in the study area is loam, with a bulk density from 0.07–1.09 g/cm³, organic carbon percentage from 1.66 to 5.98%, field capacity around 25.34%, and hydraulic conductivity that varies between 4.6 and 8.9 mm/h in a saturated state.

The multi-year average recharge, for LULC Scenarios (1) and (2), was determined using the following expression [8,10]:

$$R=\left(Pma-E{T}_P·K_c\right)\left(k_p+k_v+k{f}_c\right)$$

(3)

where $Pma$ is the multi-year average precipitation corresponding to the La Argelia station. $E{T}_P$ is the potential precipitation determined through Hargreaves and Samani.

2.4. Flash Flood Risk Assessment

Precipitation records of 51 years (from 1964 to 2015) were prepared and the data were controlled, corrected, and revised; the maximum precipitation in 24 h for each year was ordered. Different distribution functions were evaluated. Their precision was estimated through accumulated error, the Nash–Sutcliffe coefficient (NSE) and root mean square error (RMSE). Using GIS, physical parameters of the basin were determined such as area, average slope, channel length, channel slope, etc. Concentration time was estimated by Kirpich, Clark, and Temez [50,51].

For subsequent calculations, the concentration time was considered equal to design storm time. The maximum intensity was determined with this time, taking into account the intensity, duration, and frequency equations based on the maximum rainfall in 24 h corresponding to La Argelia station [52].

The hydrological response to extreme events was evaluated through synthetic unit hydrographs of the concentrated models Ven Te Chow, Snyder, Triangular SCS, and Temez. The mean and standard deviation of these results were calculated to assess dispersion and choose a model whose response is close to the mean. Curve numbers (CN) for scenario (1) and scenario (2) was determined using its respective land cover maps, taking into account the soil data provided by [27] through GIS.

Information corresponding to paths and flows was input to the HEC RAS hydraulic model and the outputs obtained were returned to GIS [11,53]. For the hydraulic evaluation, the main channel was divided into 29 cross sections. The roughness coefficients of the main channel and the banks were determined with the procedure explained by the USGS [54,55].

2.5. Meteorological Forecast

Meteorological information was taken from La Argelia station records, analogously to the method undertaken by [56]. It was developed with the annual mean temperature values; in addition, we obtained the values of maximum precipitation in 24 h of each month, and they were averaged obtaining an annual mean. These parameters will be referred to as T_average and Pmax_average from now on in the document.

The maximum precipitation data in 24 h of each month were prepared, since, as indicated by [21], one of the effects to a near horizon due to the increase of temperatures in high Andean basins is the increase in high intensity rainfall. These two annual series were forecast for the year 2029 using the Integrated Autoregressive Moving Average (ARIMA) model, the Holt exponential smoothing model (double exponent), and the Holt–Winters model (triple exponent). These models seek to describe the future behavior of the variables in relation to their past values [18].

The ARIMA model was used to predict series of a single variable. These are optimal with data series without seasonal variation, so they are widely used in annual series and meteorological variables [18,19,20]. In the ARIMA model (p, d, q), p indicates the correlation between current values and their immediate past values (autoregressive component), d indicates the order of differentiation to be applied in the series so that it becomes seasonal, and q indicates the moving average component order [17,18,57].

The double exponential model (Holt) considers an exponential smoothing, an estimation of trend, and a forecast. It involves adjusting the trend of the series at the end of each period. Triple Exponential Smoothing (Holt–Winters) adds an additional seasonality component to Holt model [58].

Series were attempted to be temporarily consistent, from 1985 to 2015. For the ARIMA model and the exponential smoothing models, their parameters were estimated (autoregressive components, seasonality, moving average and functional transformations), selecting the appropriate ones according to the best fit [17,19,56,57]. In order to verify the precision of the model and the chosen parameters, statistical measures such as the mean absolute error (MAE) and the root mean square error (RMSE) were determined.

2.6. Water Availability Estimation

In order to estimate water availability, a semi-distributed hydrological modeling was carried out through SWAT, with LULC scenarios (1) and (2). With the flow values simulated, FDC were determined for each scenario. For the modelling, DEM ALOS PALSAR RTC was used and weather information was collected from La Argelia station. The model worked with 24 h precipitation and daily average values of maximum temperature, minimum temperature, relative humidity, and wind speed. Solar radiation was estimated applying the Hargreaves and Samani equations [43].

Monthly climatic parameters required by SWAT were calculated using the SWAT Weather Database tool [59], similar to what was done by [14]. Missing values were input as −99.0, so that SWAT can estimate data for that day [60]. A soil map was generated from the soil sampling data provided by [27], data was supplemented and adapted to SWAT requirements, and the $K_{USLE}$ variable was determined by applying Williams equations [60]. As these are soils with unique characteristics, they were added to the SWAT user-soil database table, analogously to the research carried out by [61].

The LULC maps of 2019 and 2029 contained in scenarios (1) and (2) respectively were used for analysis. Coverages of forest, shrub vegetation, grassland, bare soil, agriculture, and urban life were concatenated to SWAT database coverages considering similarity between characteristics and parameters, similar to what was done by [16].

The slopes range considered for the definition of the hydrological response units (HRU) were selected based on those established by [8], which respond to an adaptation of the considerations made by [9] considering water infiltration easiness into the ground due to slope, similar to what was performed by [62].

Simulation was performed considering a daily periodicity. The number of years of omission (NYSKIP) was entered, which is necessary to train the model; in this case 5 years was input as NYSKIP. As a general rule, a minimum of 3 NYSKIP is an acceptable heating period [1]. The simulated flow begins in 1991 and ends in 2015.

In order to analyze the occurrence and frequency of the flow at the exit of the basin, and to be able to predict its availability, FDC were generated for each scenario [63,64]. This theoretical curve has been determined from modeled daily average data.

3. Results

3.1. LULC and LUCC

First, when analyzing the gains and losses by categories, it was evident that the most important changes between 2009 and 2017 occurred in relation to forest and grassland cover. When analyzing these categories individually, in the forest cover the most significant exchanges occurred with the categories of shrub vegetation, bare soil, and crops; for grassland cover, the categories with the most important exchanges were crops and shrub vegetation.

It is for this reason that the explanatory variables, and the projection’s transition model to year 2019 based on the years 2009 and 2017, included maps with exchanges between categories of grasslands and shrub vegetation, forest and bare soil, crops and grasslands, in addition to forest and crops. Moreover, because of the evident pattern of change between forest and shrub vegetation, this was included in the transition model as a map of spatial trend of change.

The spatial trend of change map was generated with a ninth-degree polynomial order, because lower polynomial orders do not significantly influence the projection. In addition, with a ninth-degree polynomial order, a higher value was reached in the confusion matrix overall accuracy. This map showed greater sensitivity within the transition model.

The transition model used for the projection to 2019 reached an accuracy rate of 78.47% (2009–2017) after being executed with the MLP algorithm and under a threshold of 10,000 iterations. An accuracy rate of about 80% is acceptable [65]. Nevertheless, other investigations support using a transition model with an accuracy rate greater than 70% for land cover prediction [38,66]. The confusion matrix, having as a reference parameter the supervised classification for 2019 and its projected data; delivered an overall accuracy of 77.03%, close to the value obtained [1].

The kappa index, analogously to the confusion matrix, provided an overall value of 0.6813. According to [67], it has a moderate concordance; for [68], it holds a substantial precision evaluation; for [69], it is considered a fair-to-good concordance; and for [70], it has a degree of agreement that is close to very good; this interpretive scale is detailed by [41]. Based on [41], the projection is considered acceptable.

Under the same conditions used in the 2019 projections with the maps of the years 2009 and 2017, the land use map for the year 2029 was generated, based on the classified maps of the years 2009 and 2019. Applying the MLP algorithm and under a 10,000 iterations threshold, an accuracy rate of 73.76% was achieved in the transition model of the map projected to 2029. Land cover classified maps for (a) 2009, (b) 2017, and (c) 2019 are presented in Figure 2. Figure 3 shows the projected maps (MLP) for (a) 2019 and (b) 2029.

Map (c) in Figure 2 represents LULC of scenario (1); map (b) in Figure 3 contains LULC of scenario (2). When analyzing the coverage extension, it is evident that the forest area is the most abundant, occupying in the basin 54.81% of its total extension for the year 2019, while on projection for 2029 it occupies 58.72%, so an increase of 3.91% is evident. Another coverage that shows an increase is shrub vegetation, occupying in 2019 an extension of 11.02%, and in the projection for 2029 of 13.71%, denoting an increase of 2.69%. Forest and shrub covers have a projected growth of 6.6%.

Bare soil, grasslands, and crop coverages have presented a reduction of 2.35%, 3.27%, and 0.99% respectively in their projection to 2029, in relation to 2019. Urban extension coverage has not shown changes. These results are presented in

Table 1. Coverage categories extension (%) in relation to 2019 and 2029.

Coverage	2019	2029
Forest	54.81	58.72
Shrub vegetation	11.02	13.71
Grassland	25.72	22.45
Crops	1.63	0.64
Bare soil	2.8	0.45
Urban	4.03	4.03

.

3.2. Hydric Recharge Analysis

Multi-year average precipitation obtained for the period under analysis is 1047.8 mm. ETP was equal to 1173.88 mm. The average $K_c$ coefficient of the basin is 0.293. Based on these, the multi-year average recharge estimated for 2019 in the basin is 614.05 mm; for 2029 it is 618.45 mm. This increase is associated with forest and shrub vegetation cover gain. It has also been determined that the highest recharge values are above the 2350 m asl level. This area has mostly forest and shrub vegetation cover. The soil is loam texture mostly, and in some areas, there are also high values of infiltration rates (>16 mm/h). Spatial distribution of hydric recharge results in the basin is presented in Figure 4.

3.3. Floodplains

A series of maximum precipitation in 24 h of each year was evaluated with different distribution functions using spreadsheets provided by [71]. Precision was estimated through accumulated error, the Nash–Sutcliffe coefficient (NSE), and root mean square error (RMSE). These results are detailed in

Table 2. Analyzed distribution functions and their adjustment statistics for the series of multiannual maximum precipitation in 24 h.

Distribution Functions	Accumulated Error	NSE	RSME
Gamma 2 parameters	8.754	0.99	1.21
Pearson Type III	168.38	−3.33	25.77
Exponential	195.86	−4.86	27.16
Nash	10.70	0.98	1.63
Normal	11.33	0.98	2.20
LogNormal	9.19	0.99	1.27
Gumbel	12.42	0.97	1.72

.

Gamma of 2 parameters was established as the distribution function that best fits the data series. Using this function, the maximum precipitation in 24 h for 3 return periods (Tr) (50, 100, and 500 years) were determined.

Using LULC maps and soil sampling data provided by [27], CN was determined, being 61.66 and 60.96 for the years 2019 and 2029 respectively. There is a projected decrease to the year 2029 due to the increase in shrub and forest vegetation.

Concentration times estimated by the equations provided by Kirpich, Clark, and Temez were 0.53, 1.46, and 2.71 h respectively. Temez was chosen because it adjusts better with the morphological and high Andean conditions of the ZH basin [72].

After calculating mean and standard deviation of peak flow rates obtained with synthetic unit hydrographs, it was determined that the model with results closest to the mean is the Ven Te Chow, which was therefore selected to determine the flood flow rates.

Table 3. Return periods, estimated maximum rainfall in 24 h, intensity for a storm time equal to concentration time, and flow rates for 2019 and 2029 coverages.

Tr (Years)	50	100	500
P 24 h (mm)	69.04	73.38	82.72
I Tc (mm/h)	16.24	17.26	19.46
2019 Q (m3/s)	3.59	5.29	9.92
2029 Q (m3/s)	3.00	4.55	8.87

summarizes obtained values.

The channel was divided into 29 cross sections. The roughness coefficients of the main channel and banks were 0.04 and 0.07 (s/m^1/3) respectively.

Table 4. Average distance from the center of the channel to floodplain margins along the channel, flooding-susceptible areas, and hydraulic results for each scenario.

Scenario	Tr (Years)	Flooding-Susceptible Areas (ha)	Floodplain Margin (m)	d (m)	v (m/s)
2019	50	28.80	33.46	1.53	0.41
	100	29.47	33.47	1.87	0.45
	500	30.86	35.65	2.05	0.58
2029	50	28.67	33.42	1.42	0.39
	100	29.29	33.45	1.61	0.44
	500	30.65	35.64	1.97	0.56

details a summary of the hydraulic results, the floodplains, and areas susceptible to flooding.

3.4. Meteorological Forecast

Statistical measures such as MAE, RMSE, and R² were determined to verify the model’s precision and the chosen parameters. For T average, adjustment statistics are presented in

Table 5. Annual average temperature, models, and adjustment statistics.

Variable	Model	Transformation	R2	RMSE	MAE
T average	ARIMA (3,0,1)	None	0.44	0.09	0.06
T average	WINTERS	None	0.32	0.08	0.06
T average	HOLT	None	0.30	0.10	0.07

, and for Pmax average in

Table 6. Maximum rainfall recorded in 24 h of each month averaged by year, models, and adjustment statistics.

Variable	Model	Transformation	R2	RMSE	MAE
Pmax average	ARIMA (1,0,0)	Natural logarithm	0.12	3.31	2.63
Pmax average	WINTERS	None	0.06	2.93	2.67
Pmax average	HOLT	None	0.02	3.46	2.91

. According to the obtained results, forecast Pmax average precision is lower compared to that obtained for T average.

In Figure 5, the ARIMA (3,0,1) model and exponential smoothing models are closely similar, both indicating an 0.11 °C increase in the average temperature for the year 2029.

In Figure 6, models keep similarities in forecasts, but given the stochastic nature of rainfall [73], they do not reach a high adjustment value. Models show an increasing trend, according to the ARIMA model (1,0,0), toward an increase of 15.63% in the year 2029, taking the year 2015 as a reference. In contrast, other researchers point out that in high Andean zones of Ecuador there is an incremental tendency in annual rainfall [74].

The temperature increase in Figure 5 is related to the increase in precipitation in Figure 6. Temperature peaks precede the precipitation peaks and when the precipitation presents higher values the temperature shows a decrease with a convex curve. In the future, as temperature increases, the peaks of precipitation would also increase, resulting in storms of greater intensity and therefore an increased risk of extreme hydrological events such as floods.

In Figure 5 and Figure 6, the forecast final values are close to each other. After the training set, those have an upward growth, following the series’ central tendency. This is because the ARIMA model and the exponential smoothing models are adjusting to series that have no periodicity.

3.5. Basin Hydrological Response

Simulated hydrographs show a seasonality associated with 24 h precipitation. The 2019 scenario reaches higher values in peak flows and lower values during dry season. The scenario for 2029 shows better flow regulation by obtaining lower values at peak flows and higher in the dry season in relation to the 2019 scenario. The hydrological response of the basin is shown in Figure 7.

Surface runoff (annual average) decreased by 1.37%, due to the increase in forest and shrub vegetation areas and the decrease in agricultural cover. The flow duration curves (FDC) (Figure 8) show that, for the 2019 coverage scenario, there is a flow of 0.3228 m³/s with a 90% probability of exceedance, while for the 2029 scenario it is 0.328 m³/s. The flow regime with a 90% probability of exceedance is 1.85% (7.72 l/s) higher for the 2029 than it is for the 2019. This trend has found up to a 1% probability of exceedance, reaching an increase of 1.05%.

On the other hand, with a probability of exceedance of 0.01%, the 2019 coverage scenario has a flow of 2323 m³/s. For the 2029 scenario, it is 2298 m³/s; that is, there is a decrease of 1.07%, a trend that is maintained for exceedance probabilities of less than 1%. These results are related to the behavior of extreme flood events flows analyzed in Section 4.3.

4. Discussion

4.1. LULC and LUCC Analysis

Land use practices such as agriculture and livestock have been limited in recent years mainly due to municipal intervention policies for the preservation of water sources. Actions taken have allowed a forest recovery, which is evidenced in the projected increase in its coverage. Additionally, [75] indicates that the population of the area will decrease at a rate greater than 1.56%, which makes sense with the non-existent growth of urban cover.

4.2. Hydric Recharge Analysis

Forest and shrub vegetation cover gain increases hydric recharge rate since it supposes a greater fraction of rain that is intercepted by vegetation. Areas where this increase is greater are located above 2350 m asl, which is also the level where a catchwater is located.

All areas related to water supply are protected by a mercantile escrow subscribed to by local government and the Regional Water Fund FORAGUA. One of FORAGUA’s main activities is the purchase of lands that are located within zones identified as water recharge zones (WRA). Most of the purchased land was declared as a municipal reserve. In those zones, any exploitation is forbidden [23]. That is why forest species have the potential to extend their borders and improve the water regulation in the basin.

4.3. Floodplains Analysis

According to Table 4 and Figure 9, floodplain margins and flooding-susceptible areas have a limited increase for return periods greater than 50 years. In addition, a reduction is expected by 2029. A remarkable number of rural buildings, agricultural land, and orchards are located on the riverside within the margin of the floodplains. This problem occurs because of fertility and ease of access to the areas.

4.4. Meteorological Forecast Analysis

Average temperature increments indicate an increase in terrestrial radiation (reflected solar) due to the presence of clouds (mainly before the first hours of radiation) and greenhouse gases (GHG) [76,77].

Adaptations in the basin, due to climate change, should be channeled around the increase in average temperature and the increase in high intensity rainfall. Being temperature linked to atmospheric pressure and amount of water vapor in the air, as it increases, the air rises by convection and gets colder, producing cloudiness, and eventually precipitation [76,77]. As there is an increasing trend in maximum precipitation, there is a probable increase in peak flow rates of flash floods.

Forecasts show a tendency to increase in maximum rainfall in 24 h, so a scenario considering its possible effects is appropriate to achieve assertive management of hydrological risks [22].

Precipitation forecasts do not provide reliable estimations of future patterns at the local scale. However, they point to an increasing trend. [22] also reports uncertainties in the precipitation forecast, noting that these are only worthy when the impacts of climate change are considered in their scenarios.

4.5. Basin Hydrological Response Analysis

FDC comparison shows an improvement on the water regulation capacity. Forest and shrub vegetation cover increments influenced the surface runoff and water recharge, and therefore the flow regulation. The greatest contribution to base flow was caused by lateral and underground flows. According to the model, the annual average of both contributions to the base flow in the 2019 scenario was 455.35 mm, while for the 2029 scenario it was 459.67 mm. These results support the premise that a mainly wooded area limits the source material of surface runoff, and therefore favors water recharge [78]. Forest and shrub vegetation trap sediment and delay runoff between adjacent zones, facilitating percolation and drainage of water through the soil, increasing base flow [1,78].

4.6. Actual Context of the ZH Basin and the Need for IWRM

These particular models comply with the initial stage of IWRM diagnosis. In a complementary investigation, these results are expected to be used in a prospective study of future scenarios and to develop a strategic plan to achieve a desired future, particularly by developing territorial prospective and implementation strategy stages [3]. A detailed description of each model used is given in Table S1 in Supplementary Materials.

In addition, the management of this particular basin is of vital importance for Loja city due to the fact that it provides about 50% of the water for human consumption [30] not considering its agricultural use, so, there is demographic pressure around the hydric resource [79]. The basin has been intervened with in terms of conservation, with highly restrictive regulations related to human activities. The local government has also bought lands closed to catchments to protect water supply, however, wealthy owners still retain some properties focused on areas in the upper part of the basin [80]; therefore, the conditions of change of use and land cover are atypical and deserve to be studied, in addition to their effect on the different natural physical processes associated with water. The meteorological forecast was developed in order to contextualize the increase in temperature and extreme rainfall and to raise the need for adaptation measures to climate change.

5. Conclusions

Forest and shrub covers show an increase for the 2029 scenario, which has a strong influence on natural physical processes in the basin. Hydric recharge suggests an increase in 2029, with the areas that favor it mostly located above 2350 m asl. The drinking water catchment points above this altitude, so it must be ensured that agricultural practices be carried out below this level in the basin.

Peak flow rates of flash floods show a reduction for the future LULC scenario. Flood margins do not show a significant decrease; however, on average they are located 36 m from the course of the river, so agricultural practices and human settlements should be limited to outside that margin to decrease their vulnerability.

Forecasts made up with the three models considered (Holt, Winters, and ARIMA) show an increase in average temperature and extreme rainfall (even with the associated uncertainties). On this issue, management strategies must deal with unpredictability and a correct capacity to respond to possible conditions of climate change.

Considering that the flow regime with a 90% probability of exceedance is greater for scenario (2) than for scenario (1), it is evident that basin recovery is due to an improvement on flow regulation related to increase of forest and shrub vegetation cover.

Intervention strategies in Andean basins can be supported in prospective studies using these key system variables focused on the integrated management of water resources.