2.1. Study area

The study area is located in the western Cordillera and has an altitude between 1438m and 4427m (Figure 1). The 283km² Pangor catchment is delimited to the north by Mt. Chimborazo (6310 m), to the west by the Chimbo River and by the cities of Pallatanga to the south and Riobamba to the east. The mean slope value in the catchment is 25°, with 50% of the area being above this value. These result in steep topography with narrow and deeply incised river valleys. The Pangor river valley crosses the Palatanga right-lateral strike-slip fault, partly visible in the field at elevations above 3600m (Figure 1). Fault slip rates are estimated at 2.9 to 4.6mm.yr¹ (Winter et al., 1993).

Figure 1: Geomorphic setting of the Pangor catchment in the western Cordillera, with location of the two pluviometric (yellow circle) and one hydrological (blue circle) station. The two sub-catchments, Virgen Yacu and Panza, where the landslide mobilization rates were quantified are delineated in black. The sampling locations for 10Be are shown by red hexagons. Figura 1: Localización de la cuenca del Rio Pangór en la Cordillera Occidental de los Andes, con la ubicación de las estaciones pluviométricas (amarillo) y la estación hidrométrica. Las dos sub-cuencas, Virgen Yacu y Panza, por las cuales la tasa de producción de sedimentos por deslizamientos flue cuantificada, están delimitadas en el mapa. La localización de los sitios de muestreo para 10Be está demarcada con un hexágono rojo.

The climate is equatorial mesothermic with a semi-humid to humid regime (Molina et al, 2015), with precipitation increasing with altitude. Intra-annual variations are characterised by a dry season from June to September and a bimodal rainy season from October to May. Temperatures stay relatively constant throughout the year with a mean annual temperature of 8°C to 20°C depending on altitude. Two pluviometric and one hydrological station from the Instituto Nacional de Meteorología e hidrología (2009) have exceptionally long time series of precipitation and discharge data for this region. The mean annual precipitation at the Juan de Velasco station (3109 m a.s.l., Figure 1) is about 1400 mm (1970-2009), with high inter-annual variability, from 475mm (2002) to 3700mm (1994). At the Chimbo DJ Pangor station (1452 m a.s.l., Figure 1), the mean annual precipitation was only about 1000 mm.

Soils are developed in volcanic and meta-sedimentary rocks, and are classified as Andisols, Histosols and Mollisols following the USDA soil taxonomy. They are characterised by high organic matter contents and high fertility when undisturbed (Henry et al, 2013; Molina et al., 2019). The land use pattern is the result of several decades of anthropogenic land use changes: at mid and low altitudes, one can observe a patchwork of small agricultural plots, remnants of subalpine cloud forest, and patches of abandoned land with regeneration of shrubs. Patches of montane cloud forest are preserved on steep slopes with low accessibility. At high altitude, páramo grasslands are the dominant vegetation system, although plantation forests with exotic species are rapidly expanding since the 2000s (Balthazar et al, 2015).

2.2. Landslide-related sediment flux

A multi-temporal landslide inventory (1963, 1977, 1989, 2010) was established for two representative subcatchments, Virgen Yacu (11 km²) and Panza (30 km²), located in the middle and lower part of the Pangor catchment (Figure 1). For 1963, 1977 and 1989, landslides were identified and digitised in ArcGIS on orthorectified aerial photographs at 1:60000 scale following the method described in Guns & Vanacker (2014). For 2010, the landslide inventory was established using World View II satellite images with a resolution of 0.5m. To avoid overestimation of landslide numbers, the digital inventory only includes landslides with fresh morphological appearance and poor to no vegetation cover. Relying on field evidence of rapid recolonisation of the landslide scars by vegetation, we assumed that the landslides occurred shortly before the date of image acquisition. The area affected by landslides is here defined as the total disturbed area. Given that it was not always possible to distinguish the erosion and deposition areas on the aerial photographs, the total landslide area might be slightly overestimated.

Landslide volumes were derived from an empirical relationship between the landslide area and volume following (Brunetti et al., 2009). For 32 landslides, the landslide length, width and depth were measured in the field in 2010 and 2011. To represent the landslide form, and also to partly compensate for area overestimation, we decided not to approximate landslides by rectangle parallelepipeds, as many authors e.g. Brunetti et al., (2009) and Guzzetti et al. (2009) do, but instead by hemi-ellipsoids. For the 32 landslides for which depth (d), length (l) and width (w) was known, we calculated landslide volume (V) by applying the hemi-ellipsoid formulae (Eq.1):

An empirical relationship can then be obtained between the landslide volume (V) and area (A):

As conversion from log(V) to V creates a systematic bias, half of the residual variance has to be incorporated in the right side of Eq. 2 to correct for volume underprediction. The corrected empirical power-law scaling linking landslide volume (V) to its area (A) is the following:

where α and ϒ are parameter estimates and σ is the residual standard error. This empirical power-law equation (Eq.3) allowed us to estimate the volume of the individual landslides. The total volume of material displaced by the landslides (V_T) is obtained by summing up the values of all individual landslides.

Following Guzzetti et al (2009), the landslide mobilisation rate, φL (mm.yr^–1), can be calculated as

Where V_T is the total volume (m³) of material displaced by landslides over a period of time t (year) and for a given catchment area, A_C (m²). In this study, we use the landslide mobilisation rate as an approximation of the landslide erosion rate, as we do not have information on the proportion of landslide debris that is temporally stored on the slope and the tributary channels. We assume that the landslide material mobilised on the hillslopes is efficiently delivered so that the landslide mass is entirely evacuated within 3 years after the landslide occurrence (Guns & Vanacker, 2013; Clark et al., 2016). This assumption is reasonable: for the 283 landslides in the inventory, the median distance from the landslide runout zone to the river network is less than 75m (Figure 2) pointing to a short and direct connection with the river network.

Figure 2: Landslide inventory containing 283 landslides. On the left: Landslide frequency-area distribution of the landslides observed over the period (1963-2010), in the Virgen Yacu and Panza sub-catchments. On the right: Distance between run-out zone and the river network for all landslides located in (semi-)natural and anthropogenic environments. Figura 2: Inventario de 283 deslizamientos. A la izquierda: la relación entre la frecuencia y el área de afectación de los eventos registrados durante la época 1963-2010, en las sub-cuencas de Virgen Yacu y Panza. A la derecha: distancia entre la base del deslizamiento cubierta por el material acumulado y el cauce del río para los deslizamientos localizados en áreas (semi-)naturales y antropogénicas.

2.3. Sediment-gauging-based average sediment yield

The average sediment yield at the outlet of the Pangor catchment was derived from the sediment rating curve and flow frequency distribution using techniques described in (Crowder & Knapp, 2005). For the hydrological station Pangor AJ Chimbo (INAMHI, Figure 1), data on daily water discharge (m³.s^–1) is available for the period between 1/01/1974 and 31/12/2009 (except for the year 1997). Notwithstanding the exceptionally rich hydrological dataset, sediment data from INAMHI is scarce with only 16 data points on suspended sediment concentration (SSC) during this time period. To overcome this limitation, an intensive data collection was organised in 2011-2012, covering the period between 20/07/2011 and 2/08/2012. During one hydrological year, a local farmer recorded the river stage at the hydrological station every day, and took a 500 ml water sample once a week during the dry season and twice a week during the rainy season. In total, 93 water samples were collected and analysed for suspended sediment concentrations at the University of Cuenca, Ecuador.

By using regression techniques, we derived the sediment rating curve from the plot of the log(suspended sediment concentration) against log(daily water discharge). As conversion from log(SSC) to SSC creates a systematic bias that leads to sediment load underestimation (Ferguson, 1986), half of the residual variance was added to the estimator to correct for underprediction. The corrected sediment rating curve is calculated following Eq. 5

where SSC is the suspended sediment concentration (kg.m^–3), Q is the water discharge (m³.s^–1) and a and b are coefficients determined by least squares regression.

The daily discharge data, recorded between 1974 and 2009, are the input data for the flow frequency distribution. The discharge values were classified into equal-size intervals of log₁₀(Q). The frequency of observed discharges in each class interval was counted and converted in days per year. The results were then plotted as flow frequency (in d.yr^–1) on the Y-axis against log₁₀(Q) class intervals on the X-axis. Taking the discharge values as the class midpoints, a mathematical adjustment was performed to find the best flow frequency distribution. The return periods for particular discharges were computed from the flow frequency distribution.

For each discharge class defined above, we estimated the suspended sediment load corresponding to the mid-class discharge value using Eq. 5 and multiplied it by the percentage of time that discharge values of that particular class are predicted to occur per year according to the flow frequency distribution. The total suspended load is then obtained by summing the data over the flow duration curve intervals. By dividing the suspended load by the catchment’s area, the average annual suspended sediment yield is derived. In this paper, the suspended sediment yield is expressed in t.km^–2.yr^–1. The advantage of this approach is that it allows us to include discharge values with long recurrence intervals that span beyond the period of observations.

2.4. Cosmogenic 10Be-based erosion rates

Steady state catchment-average denudation rates were derived from the blank-corrected ¹⁰Be concentrations in river sediments. In-situ produced cosmogenic ¹⁰Be denudation rates integrate over time scales of 10³ to 10⁴ years, and dampen short-term anthropogenic or climatic fluctuations in denudation rates (Vanacker et al., 2007b; Vanacker et al., 2015). The in-situ ¹⁰Be concentration at the terrestrial surface is inversely proportional to the denudation rate, and can be assessed from:

Where P₀ is the production rate of ¹⁰Be in quartz (atoms g^–1 yr^–1)⁰, λ the nuclide decay constant (yr^–1), ε the denudation rate (cm yr^–1), ρ the density of the material (g cm^–3) and Λ the cosmic ray attenuation length (g cm^–2). The production rate, P₀, is the result of nuclide production by nucleons, stopped and fast muons (Lal, 1991). The denudation rate, ε, is the combined effect of physical erosion and chemical weathering. In the remaining part of the paper, we refer to the in-situ produced ¹⁰Be denudation rates as cosmogenic ¹⁰Be-based erosion rates.

Fluvial sediments for ¹⁰Be analyses were collected in the active part of the stream channel from the Pangor trunk stream and its main tributaries (Figure 1). Samples were dried and sieved, and pure quartz was extracted from the 500-3000 µm grain size fraction of the alluvial sediment. We refer to (von Blanckenburg et al., 1996) for the separation methods to extract ¹⁰Be from pure quartz. Approximately 300 µg of 9Be carrier was added to the ≈20 to 40 g pure quartz. Ratios of ¹⁰Be/⁹Be were measured in BeO targets with the accelerator mass spectrometry at ETH zürich (Kubik & Christl, 2010). The ratios were normalised to the ETH in-house secondary standard S2007N with a nominal value of ¹⁰Be/⁹Be of 28.1×10^–12 (Kubik & Christl, 2010) which is in agreement with a half-life of 1.387 Myr (Chmeleff et al., 2010). From each sample, a laboratory blank with a ¹⁰Be/⁹Be ratio of 4.65 ± 0.50×10^–15 (n=9) was subtracted. The analytical uncertainties on the ¹⁰Be/⁹Be ratios of sample and blank were then propagated into the 1 SD (standard deviation) analytical uncertainty for nuclide concentrations.

In this study, the relative spallogenic and muogenic production rates are based on the empirical muogenicto-spallogenic production ratios established by Braucher et al. (2011), using a fast muon relative production rate at SLHL of 0.87% and slow muon relative production rate at SLHL of 0.27%. Corrections for skyline shielding were made following Norton & Vanacker (2009). Production rates were scaled following Lal (1991), with a sea level high-latitude production rate of 4.37±0.18 atoms.g^–1.yr^–1 (Martin et al., 2017). The latter represents the regionally averaged SLHL production rate for the Tropical Andes. As the sampled basins were not covered by snow or ice during the denudation integration time scales, no correction was made for snow shielding. The CosmoCalc add-in for Excel (Vermeesch, 2007) was used to calculate the scaling factors and catchment-average denudation rates. The total denudation is the result of physical erosion and chemical weathering. In the remaining part of the text, we refer to the denudation rates derived from in-situ ¹⁰Be nuclide as cosmogenic ¹⁰Be-based erosion rates.

3.1. Landslide-related sediment flux

The multi-temporal landslide inventory contains 283 landslides: 123 in the Virgen Yacu and 160 in the Panza catchment (Table 1). About 0.42% of the total surface area of 41 km² is affected by landslides. Most of the observed landslides can be classified as earth slides (translational slides), earth slumps (rotational slides) or rock falls according to the Cruden & Varnes (1996) classification scheme. Figure 2 shows that the landslide size, characterised by the surface area, is statistically similar in the two catchments (Mann Whitney test, p-value: 0.19). The largest slide is observed in the VirgenYacu catchment, with a surface area of 8522m².

Based on the field measurements on 32 accessible slides, the following empirical landslide volume-area relationship was established (Figure 3):

Figure 3: Empirical landslide volume-area relationship, established for the Ecuadorian Andes based on field measurements of 32 representative landslides. The 95% confidence and prediction intervals are plotted on the graph. Figura 3: Relación empírica entre el área de deslizamientos (A) y el volumen (V) para los Andes ecuatorianos, basados en mediciones en campo de 32 deslizamientos representativos. Los intervalos de confianza y de predicción del 95% están dibujados en el gráfico.

Where V is the volume of material (m²) displaced by the landslide, and A is the surface area of the landslide (m²). Although the landslide volume-area relationship is based on a limited number of field measurements, the parameter estimates for α and ϒ are in good agreement with the values that were reported by (Guzzetti et al., 2009; Larsen et al., 2010). A comparison of the landslide scaling relationship that we established for the Ecuadorian Andes with those proposed earlier in the literature is shown in Table 2.

The total landslide volume, obtained by integrating the individual landslide volumes estimated using Eq. 6 over all landslides observed, equals 3.25 x 10⁵ m³ for the Virgen Yacu and 3.39 x10⁵ m³ for the Panza catchment (Table 2). These estimates are within 1 SD error of the landslide volume estimates using the empirical landslide scaling relationship established by Larsen et al. (2010) for the himalaya.

When considering the overall time span of observations (1963-2010), the landslide mobilisation rates (derived using Eq. 4) equal 0.63 mm.yr^–1 for the Virgen Yacu and 0.24 mm.yr^–1 for Panza catchment. These rates are underestimating the real landslide mobilisation rates because the time series of remote sensing data is erratic: over the 48-year period, landslides were only detected at 4 moments of time (1963, 1977, 1989, and 2010) and correspond to morphologically fresh features. Field observations during the period 2008-2011 revealed that landslide scars are recolonised by vegetation within 2 to 3 years, making them indiscernible on optical remote sensing data after 3 to 4 years. As such, the inventory based on four discrete moments of time is representative for a 12-year time interval, instead of the full 48 years spanning the observation period. To correct for this sampling bias, we multiplied the landslide mobilisation rates by four assuming that landslide processes are homo-geneous through time. Table 3 shows the corrected landslide mobilisation rates with their uncertainties. The corrected values equal and mm.yr^–1 for the Virgen Yacu and Panza catchments respectively, and are of the same magnitude than those reported by Muenchow et al. (2012) for the southern Ecuadorian Andes. To facilitate comparison with sediment yields derived from sediment gauging and in-situ ¹⁰Be cosmogenic nuclides, the landslide mobilisation rate was converted to t.km^–2.yr^–1. here, a bulk density of 2 g.cm^–3 is used assuming that the mobilised material is a mixture of soil (with a mean density of 1.3 g.cm^–3) and weathered bedrock (with density of 2.7 g.cm^–3). After correcting for landslide runout, landslide-related sediment fluxes equal and t.km^–2.yr^–1 for the Virgen Yacu and Panza catchments respectively.

3.2. Sediment-gauging-based average sediment yield

The sediment rating curve is based on 109 measurements of suspended sediment concentration at the Pangor AJ Chimbo gauging station: 16 data points collected by INAMHI over the period 1979-2002, and 93 data points by our team over the period July 2011 and August 2012 (Figure 4). Figure 5 shows the best fit (R²=0.30) obtained for the log-log plot of suspended sediment concentration, SSC (kg.m^–3) against water discharge, Q (m³.s^–1):

Figure 4: Summary of the daily discharge data collected in the period (1974-2009) and the discharges at which 109 sediment samples were collected. The red dots correspond to the average values. Figura 4: Resumen de los datos de caudales diarios recolectados en el periodo 1974-2009, y los caudales durante las cuales se recolectó 109 muestras de agua. Los puntos rojos se refieren a los valores promedios.

Figure 5: Empirical relationship linking the suspended sediment concentration and water discharge (n=109). The 95% confidence and prediction intervals are plotted on the graph. Figura 5: Relación empírica entre el caudal y la concentración de sedimentos en suspensión en las muestras de agua del río (n=109). Los intervalos de confianza y de predicción del 95% están dibujados en el gráfico.

The relatively low R² statistic points to significant scatter in the data, which is not atypical for small to medium-sized mountain rivers (Latrubesse & Restrepo, 2014; Vanacker et al., 2007a; Rosas et al., 2020). This is probably a reflection of the flashy regime of the Pangor River, causing large variability in discharge and suspended sediment concentration through time.

The 12784 daily discharge data observed during the time period 1974-2009 were used to calculate the flow frequency distribution (in d.yr^–1). Using a class size of 0.2 log₁₀(Q), twenty-one classes cover the range of observed log₁₀(Q) (Figure 6). The best statistical adjustment corresponds to the combination of two normal distributions having approximately the same mean but strongly different dispersion. The empirical relationship provides a good representation of the reality:

where d (d.yr^–1) is the number of days per year when a discharge Q (m³.s^–1) occurred.

Figure 6: Gauging-based sediment load (kg/year, blue line) in function of log10(discharge). The mean sediment load per discharge class is obtained from the multiplication of the sediment rating curve (orange) and the flow frequency (grey bars). The dots represent the values for the midpoints of each discharge interval, and the curves are the distributions adjusted to the data. Figura 6: Carga de sedimentos en la estación de aforo (kg/año, linea azul) en función de log (caudal). La carga promedio de sedimentos por categoría de caudal ha sido obtenida de la multiplicación de la curva de descarga de sedimento (color naranja) y la curva de frecuencia de caudales (barres gris). Los puntos representan los valores para los puntos medios de cada intervalo de caudales, y las curvas son ajustadas a los datos.

To account for rare but extreme discharge values in the derivation of the average sediment yield, we extended the range of recurrence intervals to 1000 years (d = 0.001). Using Eq. 8, the discharge with a 1000-year recurrence interval was estimated at 398 m³.s^–1. At the river section of the Pangor AJ Chimbo gauging station, the 1000-year flood corresponds to a water height of 3.73 m. The high flood marks that we observed along the river channel in 2012 indicate that this flood peak discharge is entirely reasonable.

Figure 6 shows the yearly sediment loads per discharge class, and the associated 1 SD errors on the estimates. The values are obtained from multiplying the daily-suspended sediment loads estimated with Eq.7 with the flow frequency estimated with Eq.8. The statistical fit shows a humped curve, with a sediment load peak around Q = 25 m³.s^–1 (Figure 6), corresponding to the effective discharge. The average transported load decreases at higher discharge, because of rapidly decreasing flow frequency. The integration of the yearly sediment loads over the discharge classes gives the sediment-gauging average sediment yield that equals t.km^–2.yr^–1, with the uncertainties calculated based on the standard errors on the scaling exponents in Eq.7 and 8.

3.3. Cosmogenic ¹⁰Be-based erosion rates

The measured ¹⁰Be/⁹Be ratio, the ¹⁰Be concentrations, the corresponding analytical errors and uncertainties are listed in Table 4. The ¹⁰Be/⁹Be ratios of the six alluvial samples are all below 1×10^–134. There is no measurable grain size effect on the ¹⁰Be concentrations, based on the measurements of ¹⁰Be in the sand (0.5 – 2.0 mm) and fine gravel (2.0 – 3.0 mm) fractions. For two of the six samples (PG-HY¹, PG-MA), the ¹⁰Be concentrations are below the detection limit (LOD). The long-term catchment-average erosion rates range between 400 and 5357 mm. kyr^–1, and integrate over 160 to 2000 years. For the samples with ¹⁰Be/⁹Be ratios below 1×10^–14, the background of 4.65 ± 0.50×10^–15 limits the accuracy of the ¹⁰Be analyses and derived erosion rates. The lowest rates correspond to the tributary catchments of Panza (400 ± 37 mm.kyr^–1) and Virgen Yacu (758 ± 154 mm.kyr^–1), on the east of the trunk stream. The erosion rates for the trunk stream are systematically higher, with rates of 1404 ± 401 mm.kyr^–1 at the middle segment of the trunk stream and 5357 ± 5997 mm.kyr^–1 at the outlet of the Pangor basin.

4.1. Proxies for sediment yield at catchment scale

Sediment yields based on sediment gauging are first-order proxies for the suspended sediment yield. By incorporating information on the flow frequency distribution, it is theoretically possible to estimate the average sediment yield at annual and centennial time scale. however, the derivation of centennial erosion rates strongly hinges on (i) the assumption of stationarity in the flood records, and (ii) the reliability of the sediment rating curve for extreme floods. First, it is likely that decadal climate variability due to the influence of ENSO in the Pangor basin will affect the assumption of stationary in our flood frequency analysis (Machado et al., 2015; Molina et al., 2015). Second, the sediment rating curve is based on 109 measurements of total suspended solids. Although the flow rates at the 109 sampling dates are not significantly different (at 95% confidence interval) from the daily flow rates observed between 1974 and 2009, it is clear from Figure 4 that they underrepresent the high and very high flow rates. The maximum discharge observed during the time period 1974 to 2009 is 66.1 m³.s^–1 (during the 1998 El Niño), while the water sample at highest flow rate was taken at a discharge of 17 m³.s^–1 only (March 1983). Also, sediment sampling was very limited during the most influential recorded El Niño-Southern Oscillation Event of 1997-98 (Tobar & Wyseure, 2018), when the annual precipitations of 1996 and 1997 were three times higher than the mean annual precipitation in Pangor (INAMHI, 2009). We, therefore, pose that the methods based on sediment gauging data strongly underestimate the suspended sediment yield. As such, the sediment-gauging average sediment yield that equals t.km^–2.yr^–1 is a minimum estimate of the average suspended sediment yield.

Besides, the suspended sediment represents only a fraction of the river’s total sediment load that comprises the dissolved, suspended load and bed load. Two recent studies in Ecuador (Armijos et al., 2013; Tenorio et al., 2018) highlighted the importance of dissolved load. They showed that for one ton of sediment transported in suspension, 0.5 to 2 ton of solutes are transported as dissolved load. As for the dissolved load, information on the material transported as bedload is unavailable for most gauging stations in the Andes (Latrubesse & Restrepo, 2014). Bedload transport in steep Andean rivers along the western margin of the Andes can be very high, as shown by Link et al. (2002). Based on the work pursued by Link et al. (2002) and Pepin et al. (2010) in the BioBio River, Latrubesse and Restrepo (2014) suggested that the bed load sediments might account for about 75% of the total sediment load in Andean rivers. By assuming a ratio of the dissolved-to-suspended yield of 2 (Tenorio et al., 2018) and bedload-to-suspended yield of 3 (Latrubesse & Restrepo, 2014), we obtain a preliminary estimate of the gauging-based total sediment yield of t.km^–2.yr^–1.

Erosion rates based on ¹⁰Be cosmogenic nuclides range between 400 ± 37 and 5357 ± 5997 mm.kyr^–1 or 800 ± 74 and 10714 ± 11996 t.km^–2.yr^–1 (Table 4). These rates are notably higher than the ¹⁰Be-based erosion rates that were published for river sediments in the southern Ecuadorian Andes (Vanacker et al., 2007b; 2015), northern Peru (Abbühl et al., 2011) and northern Chile (Starke et al., 2017). however, they are within the range of ¹⁰Be-based erosion rates of active mountain ranges. For example Lupker et al. (2012) reported rates of 500 to 2400 mm.kyr^–1 for the himalaya, and Norton et al. (2011) rates of 170 to 1400 mm.kyr^–1 for the Alps.

The derivation of erosion rates from ¹⁰Be concentrations in river sediments is conditional on several assumptions described in Dosseto & Schaller (2016). In the Pangor catchment, the assumptions of minimal sediment storage, uniform quartz-bearing lithology, and absence of quartz enrichment during erosion processes are all validated. Yet, the sediment cosmogenic nuclide budget might be affected by the occurrence of landslides that supply material that is shielded at depth and poor in cosmogenic nuclides to the river. Model simulations (Niemi et al., 2005; Yanites et al., 2009) showed that this landslide bias in the isotopic system can lead to an overestimation of the ¹⁰Be-derived erosion rates, particularly when (1) deep landslides (depth > 5m) dominate the landslide signal and (2) catchments drain areas smaller than 50 to 75 km². Figure 2 illustrates that the majority of landslides are smaller than 10³ m², and remove material over a depth of less than 0.8 m. As the landscape is dominated by shallow landslides, the estimated ¹⁰Be-derived erosion rates are reliable proxies of the long-term catchment-average erosion.

When comparing the stream gauge-derived and ¹⁰Be-derived erosion rates, we observe significant discrepancies between them (Figure 7). The ¹⁰Be-derived erosion rates are 5 to > 50 times larger than the gauging-based total erosion rates (Table 4). As they integrate over 160 to 2000 years, they include the impact of extreme phases of the El Niño-Southern Oscillation that were reported to have occurred in 1790-91, 1877-78 (Arteaga et al., 2006), 1924-25 and 1997-98 (Tobar & Wyseure, 2018). During extreme ENSO events, the annual suspended sediment yield can increase by 3 to 60 times, as shown by Morera et al. (2017) in the northwestern part of Peru. Contrary to the ¹⁰Be-derived erosion rates, the stream-gauge derived sediment yield does not adequately account for sediment transport during high and very high floods. Hence, we consider that the ¹⁰Be-derived erosion rates provide more realistic proxies of the geomorphic process rates.

Figure 7: Comparison of landslide-related sediment flux with 10Be-derived erosion rates for four locations (PG-YC, PG-PZ, PG-PU and PG-HY). The sediment-gauging average sediment yield (± 1 SD) is shown by the blue line. Figura 7: Comparación entre las tasas de erosión causada por deslizamientos y las tazas de erosión derivadas de isotopos cosmogénicos de 10Be de cuatro puntos de muestreo (PG-YC, PG-PZ, PG-PU and PG-HY). La carga promedio de sedimentos de la estación de aforo (± 1 desviación estándar) está indicada por una línea azul.

4.2. Contribution of landslide-related erosion to catchment-average erosion rate

From the landslide inventory, the landslide mobilisation rate was estimated at and t.km^–2. yr^-1 for the Virgen Yacu and Panza catchment, respectively. The landslide mobilisation rate is systematically higher than the landslide-related erosion, as only a fraction of the displaced material is removed by fluvial processes and exported from the catchment. In the study sites, the landslide run-out deposits are located at less than 75m from the river network for > 50% of the inventoried slides (Figure 2). hence, we can reasonably assume that the landslide deposits are evacuated within 2 to 3 years after landslide occurrence, in accordance with previous work on storm-triggered landslides in the Peruvian Andes by Clark et al. (2016). After correction, we obtain landslide-related erosion rates of and t.km^–2.yr^–1 for the Virgen Yacu and Panza catchments respectively (Table 3). The difference in landslide-related sediment flux between the two tributary catchments is a result of the higher landslide density in VirgenYacu compared to Panza catchment (Table 1).

For the two subcatchments with detailed landslide inventories, the landslide-related erosion rates are similar to the CRN-derived erosion rates (Figure 7). It shows that landslides are a major source of sediments in these upland catchments and a main, if not the main, hillslope process in this high-relief tropical environment. This result is in agreement with data from the Luquillo experimental catchment in Puerto Rico, where Larsen (1997) showed that episodic landslides are the dominant source of sediment and they contribute to more than 90% of the total sediment yield at the catchment scale.

The difference between the average landslide-derived sediment flux and the ¹⁰Be-derived erosion rates seems to increase with increasing catchment area. Figure 7 shows that the difference becomes particularly large at the outlet of the Pangor catchment, where the ¹⁰Be-derived erosion rate is about 9 times greater than the average landslide-derived erosion. Given the active tectonic setting of the area, with the Pangor River draining the western Andes mountain front, the increase in CRN-derived erosion rate downstream is not surprising as was found by Vanacker et al. (2015). Downstream of sampling point PG-PU, the Pangor River rapidly steepens in a broad knickzone (Figure 1) where a variety of knickpoints are present in the trunk stream and tributaries. The enhanced fluvial erosion in the knickzone requires the adjacent hillslopes to adjust by frequent landsliding.