This study was conducted in the South American Andes, within the Catamayo–Chira transboundary river basin, which extends between 03°40′–05°07′ S and 79°05′–80°47′ W and ultimately discharges into the Pacific Ocean. The study focused on the upper portion of the basin, which drains into the reservoir formed by the Poechos Dam, which is located at 4°41′ S and 80°31’ W. Foucher et al. (2023) stated that the basin covers an area of approximately 13,565 km², with 54% of the total area situated in the province of Loja (Ecuador) and 46% situated in the department of Piura (Peru). As seen in Figure 1, the basin comprises the Chira system (1), the Alamor (2), and La Solana (3) sub-basins, which are located in the western and lower sectors, in proximity to the Poechos Reservoir, together with the Catamayo (4), Macará (5), and Quiroz (6) sub-basins, situated in the eastern and upper regions.

The basin’s elevation ranges from 81 meters above sea level (m a.s.l.) at the catchment outlet to 3,958 m a.s.l., rising from the Poechos Reservoir to the highest headwaters in the eastern part of the basin (Foucher et al., 2023). The rainy season occurs during the southern summer (December to May), when the prevailing tropical easterlies are frequently disrupted by low-level westerlies that transport moisture from the Pacific Ocean toward the study area. In contrast, the dry season extends from June to November and is dominated by persistent tropical easterlies that cause dry and sunny conditions. During this period, the elevated Eastern Cordillera acts as an orographic barrier, impeding the transport of moisture from the Amazon Basin (Ochoa et al., 2016).

According to the land-cover map generated by Foucher et al. (2023), the dominant land cover in the basin is rain-fed agriculture land, accounting for 35.2% (4,775 km²) of the total area. Dry forest covers 30.5% (4,142 km²) of the total area and is particularly well adapted to the rainfall regime and amounts that are characteristic of the study area. Herbaceous and/or shrubby vegetation occupies 22.4% (3,032 km²) of the area, primarily in areas where agricultural activity is limited. Montane forested areas, including both natural forests and plantations, cover 6.6% (896 km²) of the basin.

The Poechos Reservoir currently has an annual water demand of 2,343 hm³ (Ruiz Rosales, 2019). This volume is primarily managed by the Chira–Piura Special Project (PEChP), which operates both the reservoir and its associated hydraulic infrastructure (Seminario Gastelo, 2021). It supplies drinking water to approximately one million inhabitants in the cities of Piura, Sullana, Paita, and Talara (INEI, 2017). In addition, it irrigates 1,500 km² of agricultural land in the Chira and Piura valleys and provides water for two hydroelectric power plants located at the base of the dam and a plant situated at the end of the distribution channel to the neighboring Piura River Basin (Aguilar Villena, 2016).

Poechos began operating in 1976 with an original water storage capacity of 885 hm³ at its normal operating level of 112 m a.s.l. However, according to the bathymetric survey conducted at the end of 2018, 42 years after its inauguration, the reservoir’s capacity had decreased to just 365 hm³, a volume that is significantly lower than 437 hm³ projected for 2026 (Barboza Castillo, 2022; Foucher et al., 2023). In this context, Tote et al. (2011) and Aguilar Villena (2016) reported that nearly 49% of the reservoir’s sedimentation between 1976 and 2005 occurred during the El Niño events of 1983 and 1998, resulting in a loss of 150 hm³ in those 2 years alone. Reyes et al. (2019) pointed out that two other particularly wet years, 2008 and 2017, caused the loss of another 67 hm³.

The sediment accumulation (SA) in Poechos primarily originates from runoff and sheet water erosion occurring on the slopes of the Catamayo–Chira Basin (Arteaga et al., 2022). Additional contributions are linked to agricultural practices in the riverbeds of the tributary streams (Barboza Castillo, 2022). Furthermore, Foucher et al. (2023) observed that SA in the floodplains during dry periods is subsequently deposited in the reservoir during ENSO events.

The SSup from the Catamayo–Chira Basin was estimated using the SDR model from the InVEST 3.13.0 software (Natural Capital Project, 2022) to determine the annual mass and origin of the SA in the Poechos Reservoir. The annual soil-loss maps were generated using the revised universal soil-loss equation (RUSLE) proposed by Renard (1997), and each map was multiplied, on a per-cell basis, by its corresponding annual SDR map, following the workflow shown in Figure 2.

The matrix model of RUSLE used by Ochoa et al. (2015) for an Andean Basin is shown in Equation 1: A x , y = R x , y *   K x , y *   L S x , y * C x , y *   P x , y . (1)

Here, A _(x,y) is the estimated average annual soil loss by water erosion per unit area (t ha^-1 yr ^-1), R _(x,y) is the rainfall erosivity factor (MJ mm ha^-1 h^-1 yr^-1), K _(x,y) is the soil erodibility factor (t ha h MJ^-1 mm^-1 ha^-1), LS _(x,y) is the dimensionless factor reflecting the length and slope steepness, C _(x,y) is a dimensionless factor representing land use and land cover, and P _(x,y) is a dimensionless factor for soil conservation practices.

The R factor for each year was estimated using monthly precipitation data from the Peruvian interpolated data of the SENAMHI’s Climatological and Hydrological Observations (PISCOp V2.1) database (Aybar et al., 2020). Based on these data, the modified Fournier index (F) proposed by Arnoldus (1977) was calculated. This index evaluates the contribution of the wettest month relative to the total annual precipitation and is particularly suitable when only monthly or daily precipitation records are available. Its application in an Andean Basin in southern Ecuador by Ochoa et al. (2016) supports its relevance in similar climatic and topographic contexts. This method enables the spatial characterization of rainfall intensity, thereby facilitating the assessment of soil erosion risk and vulnerability (Equation 2). F = ∑ i = 1 12 p i 2 P . (2) Here, p _i is the precipitation of month i or of each month of the year (mm), and P is the average annual precipitation (mm). Once F was calculated, the rainfall erosivity factor R (MJ mm ha^-1 h^-1 yr^-1) was estimated with Equation 3 using the empirical equation proposed by Renard and Freimund (1994): R = 95.77 − 6.081   F + 0.4770   F 2 . (3)

Here, R is the R-factor in MJ mm ha^-1 h^-1 yr^-1 and F is the modified Fournier index proposed by Arnoldus (1977).

The K factor was determined using Equation 4 proposed by Paulet (1973) and applied for this basin by Arteaga et al. (2022). K = 0.010356 − 0.00378082 * C + 0.00232882 * S + 0.323545 * B d   /   105.33 . (4)

Here, K is the soil erodibility index in (t ha h); C = % clay, S = % silt, and Bd = bulk density (t · m^-3). The denominator accounts for the unit conversion from t ha · h t ^-1 mm ^-1 ha ^-1 to t ha h MJ ^-1 mm ^-1 ha ^-1. The soil texture and bulk density data used in this analysis correspond to the average values of the topsoil layer (0 cm–5 cm depth). These data have a spatial resolution of 500 m and were obtained from the SoilGrids database. They are available at https://www.isric.org/explore/soilgrids (Poggio et al., 2021).

The LS factor was calculated using the automated method for USLE developed by Desmet and Govers (1996), which is applicable to a two-dimensional surface. This calculation was based on a digital elevation model (DEM) with a 12.5-m resolution, which was downloaded from the ALOS PALSAR data portal https://search.asf.alaska.edu/#/.

The C factor was calculated for each grid cell (x, y) using the annual land-use/land-cover maps corresponding to each year evaluated (2001–2017). For this purpose, Moderate Resolution Imaging Spectroradiometer (MODIS) imagery was processed using the Google Earth Engine platform. The MODIS Terra Vegetation Indices 16-Day Global collection was retrieved from the Earth Engine Data Catalog (https://developers.google.com/earth-engine/datasets/catalog/MODIS_061_MOD13Q1). Then, a seasonal mosaic was generated for each of the 17 years analyzed, using imagery from the dry season (May to December). The median composite was applied as it is robust to extreme outliers—typically caused by clouds and shadows in satellite data. The Soil-Adjusted Vegetation Index (SAVI) was estimated using Equation 5 to quantify the vegetation vigor: S A V I = N I R − R   1.5   /   N I R + R + 0.5 . (5)

Here, NIR and R refer to the near-infrared and red bands, respectively, of the MODIS satellite imagery. The SAVI was applied to the 17 annual mosaics to assess the vegetation cover vigor, thereby enabling the subsequent classification of the mosaics (Linares-Palomino et al., 2022).

Factor C maps were obtained based on the supervised classification on the 17 tiles, applying random forest (Svoboda et al., 2022). Two-thirds of the spectral seeds were used as training sites for the model, and the remaining third was used to validate the classification (Belgiu and Dragut, 2016). To determine the accuracy of the classification model, the Kappa statistical coefficient was applied, which is a metric to evaluate the agreement between two classifications in a random manner (Ali et al., 2024). The Kappa coefficient ranged from 0.75 to 0.80. Values above 0.67 are considered acceptable classifications as described by Carletta (1996).

Considering this supervised classification, the C factor implies the P factor of soil conservation according to Benavidez et al. (2018). The values of the C factor used are those proposed for Peru by Sabino Rojas et al. (2017): 0 for water bodies and moors, assuming that the erosion generated in them tends to be null; 0.003 for montane forest; 0.005 for TDF in hills and mountains with greater vigor; 0.05 for TDF in hills with less vigor; 0.09 for the urban area; 0.22 for scrubland; 0.25 for dry areas and pastures; and finally, 0.63 for the irrigated agricultural area, because this land-use type generates greater erosion in the studied basin (Arteaga et al., 2022).

As shown in the maps in Figure 3, the vegetation vigor and cover vary significantly from year to year, depending on precipitation. At the end of a dry year, areas dedicated to rain-fed agriculture typically exhibit bare soil. In regions dominated by dry forest—characterized by predominantly deciduous vegetation—foliage and understory cover are lost, leading to a decrease in vegetation vigor. In both cases, soil protection is reduced, increasing susceptibility to erosion.

The sediment delivery rate, SDR_(x,y), was determined by applying the SDR model described in the Natural Capital Project (2022). These results (SDR_(x,y)) were multiplied with the corresponding cells of the A _(x,y), and Equation 6 was used to estimate the sediment supplied (SSup_(x,y)) in t km^-2 from each cell of the basin map _(x,y). S S u p x , y = A x , y   *   S D R x , y . (6)

Thus, the flux of water supplied by the basin and deposited in the Poechos Reservoir was determined (Benavidez et al., 2018; Borselli et al., 2008). To estimate the total sediment supply SSup _(t) from the basin to the reservoir, the sediment contributions from all the cells in the basin, SSup _(x,y) , were added, as shown in Equation 7. S S u p t = ∑ i E x , y . (7)

Finally, a simple regression analysis was performed between the series of volume (hm³) values of SA and the series of mass (t) values of SSup from the basin, considering a soil bulk density of 1.22 t m^-3 in the basin, based on the cartographic data from https://www.isric.org/explore/soilgrids. To perform this analysis, a 95% confidence level was considered, and it was assumed that both variables coincide at the origin of the coordinate axis (Mejía et al., 2024). Zero SSup contributes with zero SA in the reservoir, and, therefore, β ₀ = 0.

Figure 4 shows that the highest values of the R factor (erosivity) in a specific year have a strong relationship with the amount of SSup estimated. The spatial distribution of erosivity by the sub-basin is provided in the figure in the supplementary material accompanying this article. Moreover, a clear relationship is observed in the annual course between SSup and SA, except for 2008, when the bathymetric measurement tripled the estimate. Specifically, this year saw the occurrence of a coastal El Niño phenomenon described by Farías de Reyes and Montero (2023). Because the occurrence of these coastal events was not yet understood, the event was underestimated, leading to poor reservoir management and high sediment retention. The coastal El Niño was first detected by the Geophysical Institute of Peru (IGP) in 2012 and became widely recognized in 2017 (Takahashi, 2017).

Figure 4 also confirms the trend line associated with the annual course between SA and SSup. The correlation coefficient r (p < 0.05) for the period 2001–2017 is 0.8424 and increases to 0.9403 when the 2008 ENSO data are excluded, indicating a strong positive correlation. Similarly, the adjusted coefficient of determination R ² is 0.6474, which may appear modest. However, when 2008 is excluded from the time series, R ² increases to 0.8174. This suggests that 81.74% of the variance in SA is explained by the independent variable SSup, supporting the use of the SSup variable to better understand soil erosion and sediment transport processes within the basin.

Using the SDR model, it was estimated that the basin supplied a median of 6.91 x 10⁶ t yr^-1 of sediment during the period 2001–2017, a value close to 6.53 x 10⁶ t yr^-1 of sediment accumulation measured through bathymetry over the same period. Additionally, the highest SSup to the Poechos Reservoir during the evaluated period occurred during the 2017 El Niño event, reaching 34.92 x 10⁶ t, which is approximately five times the median value (Figure 5).

Figure 6 presents the annual distribution of SSup by the sub-basins to the Poechos Reservoir for the period 2001–2017. The Catamayo, Macará, and Quiroz sub-basins located in the eastern part of the study area contributed 84.3% of the total SSup, with median values of 1.50, 2.34, and 1.98 x 10⁶ t yr^-1, respectively. These sub-basins are characterized by dryland agriculture on steep terrain, with slopes exceeding 25%, which results in high LS, C, and P factors. In contrast, the remaining sub-basins, namely, Alamor (0.65 x 10⁶ t yr^-1), La Solana (0.18 x 10⁶ t yr^-1), and Chira (0.05 x 10⁶ t yr^-1), accounted for 12.8% of SSup. These are located in the western portion of the study area, where the seasonally dry forest is confined to the mountainous zones.

An analysis of SSup from the sub-basins revealed that Catamayo, Macará, and Quiroz exhibited lower interannual variability in SSup than the La Solana and Alamor sub-basins throughout the study period. The eastern sub-basins contributed between 1 and 6 x 10⁶ t yr ^-1, whereas the western sub-basins supplied between 0 and 10.5 x 10⁶ t yr ^-1.

In 2017, during the occurrence of the so-called “coastal” El Niño event, the total SSup from the basin reached 34.92 x 10⁶ t. Each of the three eastern sub-basins tripled their median sediment contribution, reaching approximately 5.60 x 10⁶ t. Meanwhile, Alamor increased its sediment contribution tenfold to 6.68 x 10⁶ t and La Solana supplied 10.50 x 10⁶ t, which represents a 57-fold increase compared to its median value. During that year, approximately 52% of the SSup to the Poechos Reservoir originated from the western sub-basins, as illustrated in Figures 6, 7.

By dividing the annual median of the total sediment mass supplied from the basin, 6.91 x 10⁶ t yr^-1, by its area, the median sediment yield was estimated at 510 t km^-2 yr^-1, a value close to the tolerable soil loss rate of 600 t km^-2 yr^-1 (Gashaw et al., 2021). When analyzing the sediment yield per unit area, the Macará and Quiroz sub-basins exhibited median values of 828 and 632 t km^-2 yr^-1, respectively, and both exceeded the tolerable limit. They were followed by Alamor (547 t km^-2 yr^-1), Catamayo (358 t km^-2 yr^-1), La Solana (152 t km^-2 yr^-1), and the Chira system (46.3 t km^-2 yr^-1).

Figure 7 shows that the highest SSup occurred in 2017, when the basin supplied an estimated 2,578 t km^-2 in that year alone. This value is slightly higher than that reported by Rosas et al. (2023) for a neighboring basin. Notably, in 2017, the La Solana and Alamor sub-basins, where slopes (the LS factor) are the least pronounced in the study area, supplied 8,723 and 5,610 t km^-2, respectively. In contrast, the eastern sub-basins, that is, Macará, Quiroz, and Catamayo, which are characterized by steeper slopes, supplied 1,997, 1,822, and 1,306 t km^-2, respectively.

A detailed analysis of the sediment sources reveals extensive areas within the basin where SSup is recurrent. The highest contributions were estimated to originate from the basin region between the Quiroz and Macará sub-basins, particularly in the districts of Jililí, Montero, and Sícchez (Ayabaca Province, Peru), where seasonally dry forests and rain-fed agriculture dominate on slopes greater than 16% (Figure 7). Even during dry years, this area supplied more than 500 t km^-2·yr^-1, whereas in the wet years, it exceeded 5,000 t km^-2·yr^-1, reaching up to 20,000 t km^-2·yr^-1 in certain zones.

The Alamor sub-basin also contributed more than 2,000 t km^-2· yr^-1 during dry years. However, Figure 7 highlights specific areas within this sub-basin, particularly in the Cantons of Pindal and Célica (Loja Province, Ecuador), that exhibited much higher sediment yields. These zones are characterized by rain-fed maize cultivation (Ochoa et al., 2022). During wet years, sediment exports in this region reached 10,000 t km^-2·yr^-1 and peaked at 20,000 t km^-2·yr^-1 in 2017. Similarly, the lower part of the Catamayo sub-basin, which is also located in the Célica Canton, exported between 5,000 and 20,000 t km^-2·yr^-1 of sediment during wet years.