The Puyango-Tumbes basin has an area of 4,800 km², of which 60% belongs to the provinces of El Oro and Loja, located in the southeast of Ecuador, and 40% is in the Department of Tumbes, located in the north of Peru. The Puyango-Tumbes River originates at an altitude of 3,500 m above sea level in the Portovelo area, where it is called the Pindo River, then becomes the Puyango River, and finally the Tumbes River in Peruvian territory (Figure 1). Annually, the average annual precipitation is 1,200 mm with marked variations from 100 mm to 2,700 mm, and the temperature on the plains is 24.5°C and in the mountainous area is 22°C (UNDP, 2015).

Approximately 70% of the area has soil suitability for protection and/or restoration; almost two-thirds of the basin consists of fragile lands, whose natural conditions are not suitable for agricultural production, and their utilization implies severe risks of soil erosion and degradation (MAP and GIS, 2019).

The basin is of very high socioeconomic interest for both countries since it contains a large part of the population and their productive activities. However, the economy of these areas largely depends on water availability and, therefore, is vulnerable to poor management, overexploitation, and contamination of the resource, as well as the effects of variability and climate change (Jones et al., 2017).

This study used a daily time scale and sub-basin scheme derived from the digital elevation model (DEM). The sub-basin configuration through the TauDEM tool preserved natural channels and flow paths. Stream reaches were defined by choosing the minimum surface threshold automatically set by SWAT as a criterion to create them, thus obtaining a rigorous representation of the channels (Abbaspour et al., 2015a) (Figure 2A). Hydrological Response Units (HRUs) were obtained considering soil type, land use, and landscape slope (Zhang, 2014). Modeling with SWAT also considered daily climatic information, specifically precipitation, maximum and minimum temperatures for determining the main processes related to the hydrological cycle (Abbaspour et al., 2015a), and thus, with the calibrated parameters, served for application in future projections (Figure 4).

Climatic information, specifically precipitation, was obtained from the gridded product RAIN4PE (Rain for Peru and Ecuador) (https://dataservices.gfz-potsdam.de), and maximum and minimum temperatures were obtained from the gridded product PISCO v2.0 for the period 1981–2015 (http://iridl.ldeo.columbia.edu).

Historical hydrometric information was obtained from the National Institute of Meteorology and Hydrology (INAMHI) of Ecuador and the National Service of Meteorology of Hydrology of Peru (SENAMHI) from three hydrometric stations “Pindo,” “Puyango,” and “El Tigre” for the period 1992–2015 (Figure 2C).

The SWAT model divided the slope into three categories, based on information obtained from the DEM (Figure 2A). The global soil map was obtained from the Harmonized World Soil Database v1.2 (HWSD) (Figure 2B).

In analyzing the LULC of the Puyango-Tumbes watershed over a 35-year period (1981–2015), two LULC scenarios from the MODIS product (https://ladsweb.modaps.eosdis.nasa.gov) were selected to assess the influence of future climate on streamflow and sediment yield: (a) an optimistic scenario characterized by reduced anthropogenic impacts (LULC_1985), and (b) a pessimistic scenario, representing the most realistic impacts (LULC_2015), considering the intensification of monoculture with pasture.

The configuration of the LULC scenarios in this research was based on an analysis of historical land-use trends, acknowledging that land-cover changes in the region are driven by local factors, including agricultural practices, land-use policies, and demographic pressures. This approach is justified by the fact that the SWAT model requires inputs that accurately reflect local conditions in order to perform precise simulations of hydrological behavior.

The model was calibrated and validated on a daily scale using daily flows along the Puyango-Tumbes River at three hydrometric stations: “Pindo,” “Puyango,” and “El Tigre” (Figure 2A). 28 years (1988–2015) were utilized for calibration and validation, including 4 years for model warm-up.

Parameter calibration was performed using the SWAT_CUP uncertainty procedures software (Abbaspour, 2015b). Based on previous studies (Cibin et al., 2010; Funk et al., 2015), 16 sensitive hydrological parameters were selected for analysis (Table 1).

A global sensitivity analysis was conducted to obtain both relative and absolute sensitivity. The Sequential Uncertainty Fitting algorithm, version No. 2 (Sufi-2) (Abbaspour et al., 2015a), was implemented to identify the most sensitive parameters according to the model’s response. In Sufi-2, parameter uncertainty, expressed as ranges (uniform distributions), accounts for all sources of uncertainty, such as uncertainty in input variables (e.g., rainfall), the conceptual model, parameters, and measurement data. The propagation of uncertainties in the parameters leads to uncertainties in the model’s output variable, expressed as 95% probability distributions, referred to as the 95% Prediction Uncertainty or 95 PPU.

These 95 PPU represent the outcomes of the models in a stochastic calibration approach. It is important to realize that we do not have a single signal representing the model results, but rather a set of good solutions expressed by the 95 PPU, generated by certain parameter ranges (Abbaspour et al., 2015a).

The calibrated parameters at the three hydrometric stations “Pindo,” “Puyango,” and “El Tigre” were used in the future flow projections. For the assessment of sediment load, the bottom solids discharge curve in liquid discharge fusion (Equation 1), developed by (Goyburo, 2017) for the El Tigre station, was employed. This involved characterizing the sediment evolution based on liquid flows during extreme events using the database of the Peruvian Geophysical Institute (IGP) and the total sediment transport during the sampling campaign (Jan-May/2016). It was determined that the methodology proposed by Rennie, Millar, and Church (2002) best approximates the observed data of bottom solid discharge, obtaining a correlation R = 0.43 between bottom solid discharge and observed liquid discharge. S y = 0.602 Q 0.938 (1) where S y : sediment bedload yield (tday⁻¹); Q : Flow (m³s⁻¹).

To generate long-term simulations of future climate, climate change scenarios were defined based on CMIP6 data, as they provide climate projections from multiple models based on alternative scenarios that are directly relevant to societal concerns regarding climate change mitigation, adaptation, or impacts (Riahi et al., 2016). The “Shared Socioeconomic Pathways” scenarios were used in this study:

SSP2-4.5 (this scenario represents the middle part of the range of future forcing pathways and updates the RCP4.5 pathway. SSP2 was chosen because its land use and aerosol pathways are not extreme compared to other SSPs) and SSP5-8.5 (this scenario represents the upper extreme of the range of future pathways, updating the RCP8.5 pathway). SSP5 was chosen for this forcing pathway because it is the only SSP scenario with emissions high enough to produce a radiative forcing of 8.5 W/m² in 2,100). The methodology applied in this study is based on the integration of 24 CMIP6 models, which allowed for the generation of more robust climate projections and minimized long-term uncertainties, especially in simulating critical variables such as precipitation and temperature (maximum and minimum) were downloaded for the near future (2035–2065) and far future (2070–2100) for the two scenarios. This was crucial for the hydrological analysis in the SWAT model, providing a solid foundation to assess the impacts of climate change on the Puyango–Tumbes River watershed under different radiative forcing scenarios (SSP2-4.5 and SSP5-8.5). By combining the mean of multiple models with delta change bias correction techniques applied to the 23 sub-watersheds, the reliability of the obtained results is ensured.

To correct CMIP6 data, bias correction was performed using the delta change method, which is an adjustment factor derived from simulated data and applied to the observed data series from 1981 to 2015. This procedure was applied to all records of the 23 obtained sub-catchments of the setup hydrological model, this method involves calculating an adjustment factor by dividing the averages of all months of the observed data by the averages of all months of the modeled data (Guo et al., 2022; Golian, El-Idrysy, and Stambuk, 2023). It is important to note that this method was also applied in the same manner for the 5th percentile (value below which 5% of the data lie) and 95th percentile (value below which 95% of the data lie) of the CMIP6 data.

Subsequently, a new series is obtained by multiplying the previously calculated adjustment factor with the modeled data series (Equation 2). ∝ j = x j o b s ¯ x j mod , (2) X i , j * = ∝ j X i , j mod

Where ∝ j = adjustment factor, x j o b s ¯ = observed data, x j mod = modeled data, and X i , j mod = scaled series.

The modeling of the Puyango-Tumbes basin was performed with the combination of land use and land cover scenarios (LULC_1981 and LULC_2015) and future projections (SSP2-4.5 and SSP5-8.5) (Figure 3).

According to the sensitivity analysis, six parameters significantly affect the modeled surface runoff in the Puyango-Tumbes basin. These parameters include the soil available water capacity (SOL_AWC), the average slope length (SLSUBBSN), and other parameters related to surface runoff (SURLAG) and groundwater (ALPHA_BF, GW_DELAYMM, GWQMN). These parameters were calibrated to fit the actual water balance based on literature information (Le et al., 2012). Calibration was executed using the Sequential Uncertainty Fitting version 2 (SUFI-2) algorithm (Arnold et al., 2012), a widely recognized method for automatic calibration and sensitivity analysis in SWAT models. Sensitive parameters affecting surface runoff and groundwater flow, such as soil available water capacity (SOL_AWC), average slope length (SLSUBBSN), and groundwater parameters (ALPHA_BF, GWQMN, GW_DELAYMM), were identified through global sensitivity analysis with SWAT-CUP. Defined parameter ranges enabled simulations that varied one parameter at a time, quantifying their sensitivity on model outputs. Subsequently, the SUFI-2 algorithm was applied to systematically adjust these parameters, aligning modeled outputs with observed water balance data, which ensured an accurate representation of hydrological processes in the Puyango-Tumbes basin and enhanced the reliability of water management decisions.

If compared to the default values, the ALPHA_BF parameter was increased in the model to facilitate the flow of water from the aquifer to the river, thus increasing the flow rate.

When compared to the default values, the ALPHA_BF parameter was increased in the model to facilitate the flow of water from the aquifer to the river, thereby increasing the base flow. Additionally, the GWQMN and GW_DELAY parameters were increased, leading to higher surface flows and consequently decreasing groundwater flow. The remaining parameters related to hydrological processes were calibrated to adjust both base and peak flows (Table 2).

This section presents a comparison of the SSP2-4.5 and SSP5-8.5 scenarios with observed precipitation data, regarding the spatial distribution of precipitation throughout the basin (Figure 4). Table 4 shows the comparison between observed and downscaled monthly average precipitation, with coefficient of determination (R²) values of 0.95, indicating a good correspondence between observed and downscaled data after bias correction by delta change. The analysis revealed that by the end of the 2065 decade, the annual mean precipitation is projected to increase by up to 6% for SSP2-4.5, with a similar trend for SSP5-8.5 showing an increase of 14% by 2065, indicating a rising trend. The projected temperature increases ranges from 2.8°C to 4.3°C and 3.7°C–7.4°C in the SSP2-4.5 and SSP5-8.5 scenarios, respectively.

The streamflow records in response to LULC change allowed the evaluation of sediment yield for the Puyango Tumbes river basin. The results indicate that sediment yield in the basin is directly related to precipitation, suggesting that this region has a rapid flow response. The highest sediment yield values occurred in 1985, as there was less anthropogenic influence. The lowest estimated sediment yield occurred in 2015, indicating a lower streamflow generation. Furthermore, the results show a variation in the standard deviation for streamflow and sediment yield data between LULC_1985 and LULC_2015 for Pindo equal to 4 m³ s⁻¹⁻ and 13%; Puyango equal to 33 m³ s⁻¹⁻ and 18%; and El Tigre equal to 21 m³ s⁻¹⁻ and 9%, respectively (Table 3).

As illustrated in Figures 5A, B, a comprehensive analysis of LULC transitions was conducted, revealing significant changes among various land cover classes during the study period. It is noteworthy that the forest, grassland, and savanna classes exhibited a high probability of persistence, exceeding 80%, indicating their stability within the basin. However, the grassland class emerged as the most dynamic, demonstrating a significant increase of 18.3%. This pronounced change suggests that grasslands are undergoing accelerated conversion, likely due to anthropogenic pressures and environmental changes.

These dynamics are crucial as they not only influence the current land cover but also have profound implications for future land management and ecological health in the region. The observed changes in grassland cover may contribute to alterations in hydrological cycles, habitat availability, and biodiversity, which in turn affect ecosystem services. Furthermore, the characteristics of these changes in land use and land cover are critical for future projections, as presented in Figures 7–9. Understanding these trends is essential for developing adaptive management strategies aimed at mitigating adverse impacts on the basin’s resources. This expanded analysis underscores the need for continuous monitoring and evaluation of land use dynamics to inform policies and enhance the sustainable management of resources in the Puyango-Tumbes River Basin.

Figure 6B shows the spatial distribution of sediment yield in the sub-basins. The results indicate that sediment production was higher in the sub-basins located in the upper part due to the predominance of areas with pasture cultivation and the influence of higher precipitation, resulting in increased streamflow (Figure 6A), contributing to accelerated erosion in these sub-basins. It was affirmed that the increase in river flow in the future depends mainly on increases in precipitation.

Table 4 shows the results of estimated streamflow and sediment yield for the Puyango-Tumbes River basin using future climate data, for the 5th percentile, mean, and 95th percentile. In this context, at the Pindo station, streamflow based on SSP2-4.5 and optimistic LULC increased relative to the baseline period (1981-2015) by 20% and 43%, respectively, in the periods 2035-2065 and 2070-2,100. According to SSP5-8.5, the expected streamflow exceeded the observed in the baseline period by 66% in the period 2035-2065, and by 86% in 2070–2,100 (Figure 7). The figure shows how the streamflow projections for future periods vary compared to the historical period, highlighting potential changes in flow patterns over time under different climate change scenarios.

At the Puyango station, a similar behavior to the Pindo station was observed. Streamflow based on SSP2-4.5 and optimistic LULC increased relative to the baseline period (1981-2015) by 17% and 29%, respectively, in the periods 2035-2065 and 2070-2,100. According to SSP5-8.5, the expected streamflow exceeded the observed in the baseline period by 32% in the period 2035-2065, and by 36% in 2070–2,100 (Table 4; Figure 8).

At the El Tigre station, streamflow based on SSP2-4.5 and optimistic LULC increased relative to the baseline period (1981-2015) by 11% and 40%, respectively, in the periods 2035-2065 and 2070-2100. According to SSP5-8.5, the expected streamflow exceeded the observed in the baseline period by 23% in the period 2035-2065, and by 52% in 2070–2100. In all three stations, this increase is due to the increase in estimated precipitation relative to SSP2-4.5 (Table 3; Figure 9). This figure provides a comparison of projected streamflow for the future and its difference from historical values, showing potential trends in flow as time progresses under different climatic scenarios.

In this study, the SWAT hydrological model, driven by climate simulations, estimated an overall increase in the flow of the Puyango-Tumbes River for both optimistic and pessimistic scenarios.

The results obtained showed that the simulated mean annual streamflow will mainly increase between 2035 and 2065 in SSP5-8.5, while in SSP2-4.5, it presented the lowest values. Between 2035-2065 and 2070-2100, there is less variation between the SSPs and scenarios, but also with an increase in values compared to the baseline. The results of the simulated annual streamflow showed an increase in the optimistic scenario in both climate scenarios (SSP2-4.5 and SSP5-8.5). This result is sensitive to the outcome of future climate projections (Santos et al., 2021; Silva et al., 2020) affirmed that the increase in river flow in the future depends mainly on increases in precipitation.

As evidenced in Table 4 and Figure 10, sediment production in the watershed showed the same trend as streamflow in both climate change scenarios. Projected sediment yield values for SSP2-4.5 exceeded those of the baseline by 18%–40% for the Pindo station, 16%–27% for the Puyango station, and 10%–37% for the El Tigre station. Regarding SSP5-8.5, the values were even higher, ranging from 61% to 79% in Pindo, 30%–33% in Puyango, and 21%–48% in El Tigre, for the different analysis periods. The pessimistic LULC scenario showed increases in sediment yield for all years from 72% to 113% for Pindo, 68%–86% for Puyango, and 45%–69% for El Tigre compared to SSP2-4.5. As for SSP5-8.5, the increase was greater, ranging from 122% to 145% for Pindo, 88%–101% for Puyango, and 53%–82% for El Tigre compared to the optimistic scenario, demonstrating that soil erosion in the region is more sensitive to climatic variability than streamflow, and that LULC changes can have serious impacts on the watershed. These results serve as an urgent call for the adoption of mitigating measures to prevent associated environmental problems and improve land and water resource management.

In the pessimistic scenario, at the Pindo station for SSP2-4.5, the flow rate exceeds the optimistic scenario by 11% in the period 2035-2065 and by 16% in the period 2070-2100. The flow rate results for the pessimistic scenario SSP5-8.5 show a slight increase of four%–5% compared to the optimistic scenario. At the Puyango station for SSP2-4.5, there is a slight elevation of the pessimistic scenario compared to the optimistic one by 1% in the period 2035-2065 and by 2% in the period 2,070–2,100. The flow rate in the pessimistic scenario SSP5-8.5 shows an increase of 1–5% compared to the optimistic one. Regarding the El Tigre station in SSP3-4.5, it surpasses the pessimistic scenario by 9% in the period 2,035–2,065 and by 2% in the period 2,070-2,100. The flow rate results for the pessimistic scenario SSP5-8.5 show an increase of 4% in the period 2,035–2,065 and 2% in the period 2,070–2,100.

The flow rate in the pessimistic scenario SSP5-8.5 shows an increase of between 1 and 5% compared to the optimistic one. Regarding the El Tigre station in SSP3-4.5, it exceeds the pessimistic scenario by 9% in the period 2,035-2,065 and by 2% in the period 2,070-2,100. The flow rate results for the pessimistic scenario SSP5-8.5 show an increase of 4% in the period 2,035-2,065 and 2% in the period 2,070-2,100.

The highest values of flow rate and sediment yield in the SSP2-4.5 and SSP5-8.5 scenarios occurred in 2,070–2,100, which is related to the precipitation pattern, as this period is the rainiest in these scenarios. The results suggest that in the near future (2,035-2,065), there may be negative impacts on the dynamics of the Puyango-Tumbes River basin, including floods, river sedimentation, affecting conservation and preservation efforts in the area.

Based on the SWAT model projections, the results show an increase in sediment yield for the analyzed LULC scenarios and climate change. These values suggest that responsible decision-makers should implement measures for flood management and avoid land use/land cover change transitions.