Exotic tree plantations in the Chilean Coastal Range: balancing the effects of discrete disturbances, connectivity, and a persistent drought on catchment erosion

The Chilean Coastal Range, located in the Mediterranean segment of Chile, is a soil-mantled landscape with the potential to store valuable freshwater supplies and support a biodiverse native forest. Nevertheless, human intervention has been increasing soil erosion for ∼ 200 years, culminating in the intense management of exotic tree plantations throughout the last ∼ 45 years. At the same time, this landscape has been severely affected by a prolonged megadrought. As a result, this combination of stressors complicates disentangling the effects of anthropogenic disturbances and hydroclimatic trends on sediment fluxes at the catchment scale.

In this study, we calculate decennial catchment erosion rates from suspended-sediment loads and compare them with a millennial catchment denudation rate estimated from detrital ¹⁰Be. We then contrast both of these rates with the effects of discrete anthropogenic-disturbance events and hydroclimatic trends. Erosion and denudation rates are similar in magnitude on decennial and millennial timescales, i.e., 0.018 ± 0.005 and 0.024 ± 0.004 mm yr⁻¹, respectively. Recent human-made disturbances include logging operations throughout all seasons and a dense network of forestry roads, thereby increasing structural sediment connectivity. Further disturbances include two widespread wildfires (2015 and 2017) and an earthquake with an M_w value of 8.8 in 2010.

We observe decreased suspended-sediment loads during the wet seasons for the period 1986–2018, coinciding with declining streamflow, baseflow, and rainfall. The low millennial denudation rate aligns with a landscape dominated by slow diffusive soil creep. However, the low decennial erosion rate and the decrease in suspended sediment disagree with the expected effect of intense anthropogenic disturbances and increased structural (sediment) connectivity. Such a paradox suggests that suspended-sediment loads, and thus respective catchment erosion, are underestimated and that decennial sediment detachment and transport have been masked by decreasing rainfall and streamflow (i.e., weakened hydroclimatic drivers). Our findings indicate that human-made disturbances and hydrologic trends may result in opposite, partially offsetting effects on recent erosion, yet both contribute to landscape degradation.

1 Introduction

Over 75 % of the Earth's ice-free land has been altered by humans (Ellis and Ramankutty, 2008), resulting in severe consequences for sediment transport during the Anthropocene (Syvitski et al., 2022). Land-use and land-cover change (LULCC) is an important factor in increasing soil erosion (Borrelli et al., 2020). Human-made forests – or better, tree plantations (DellaSala, 2020) – are frequently disturbed by logging and the implementation of forestry roads. Such disturbances may intensify soil erosion (e.g., Schuller et al., 2013; Sidle and Ziegler, 2012), as can heavy-machinery traffic (e.g., Malmer and Grip, 1990), wildfires, and terracing (e.g., Martins et al., 2013). Short rotational cycles, i.e., the period between planting, harvesting, and replanting tree plantations, also change hillslope stability through cycles of root strength decay and recovery, which, in turn, promote landsliding and debris flows (Imaizumi et al., 2008; Montgomery et al., 2000). Ultimately, all these processes may modify sediment trajectories and storage on hillslopes and along rivers (Wainwright et al., 2011), with long-lasting impacts on sediment yields over periods of 10–100 years (Moody and Martin, 2009; Bladon et al., 2014).

Figure 1 Study region. (a) Land cover in the Chilean Coastal Range (Zhao et al., 2016) and catchments with published detrital ¹⁰Be denudation rates, outlined in black (Carretier et al., 2018). The Purapel catchment, whose denudation rate is presented in this work, is represented in red. (b) Purapel catchment. All the detected forestry roads intersecting tree plantations and the locations of where the photos in panels (c) and (d) were taken are shown. Elevation data come from a lidar digital terrain model (DTM) with a 5 m resolution that was obtained in 2009. (c) Photos captured on hillslopes in the Purapel catchment after the fire. (d) Photos taken after the fire showing bare topsoil on hillslopes with (4) granitic and (5) metamorphic lithologies.

The segment of the Chilean Coastal Range (CCR) that is located between 35–37.5° S is characterized a landscape of gentle and largely convex hillslopes (Fig. 1) under Mediterranean conditions. Its morphology results from relatively slow denudation rates caused by soil creep on regolith-mantled landscapes (Roering et al., 2007), yet it is modified by the underlying bedrock (Gabet et al., 2021). Here, forests, soils, and water are closely coupled (Galleguillos et al., 2021). Currently, the remnants of secondary native forests (i.e., successional forests growing in areas where forest cover was removed at some point in the past) stand on soils that are up to 2 m in depth (Soto et al., 2019), suggesting that soils under natural vegetation cover could have an even greater thickness. In the absence of snow storage, these soils form a major freshwater supply along the Mediterranean CCR, which many rural communities rely on. Thus, decision-making regarding land management is critical for the resilience of these communities (e.g., Gimeno et al., 2022), especially given the recent (Garreaud et al., 2020) and projected (IPCC, 2021) conditions of water scarcity.

The CCR has experienced deforestation for more than 200 years (Armesto et al., 2010), intensifying soil erosion, as recognized by Bianchi-Gundian (1947) and Chilean governments in the middle of the 20th century (IREN, 1965). From the beginning of the 20th century, governments blamed environmental issues caused by deforestation to promote the expansion of tree plantations (e.g., CONAF and MINAGRI, 2016; Pizarro et al., 2020). The most relevant transformation of land cover began with a decree law (DL 701 in 1974), which aimed to subsidize the forestry sector (Manuschevich, 2020). This law and subsequent political action accelerated LULCC, which, in practice, transformed degraded lands, shrublands, and native forests into industrially managed tree plantations (Heilmayr et al., 2016). From ∼ 450 000 ha of tree plantations in 1974 (Barros, 2018), their spatial extent increased to at least some 2.8 ± 0.2 ×10⁶ ha in 2011 (Heilmayr et al., 2016), mostly within the Mediterranean CCR (Fig. 1).

Figure 2 Forest roads in the Purapel catchment at different stages of the tree plantation rotation cycle and their connections to streams. (a–d) Google images (©Google Maps 2016, ©Google Maps 2021, ©Google Maps 2022, and ©Google Earth 2016, respectively). (e, f) Photographs of a gravel road and its connection to a stream taken during a storm in July 2017.

In Chile, tree plantations are managed mostly as monocultures of fast-growing Eucalyptus spp. or Pinus radiata stands. For these species, the rotation cycles are as short as 9–12 and 18–25 years, respectively (INFOR, 2004; Gerding, 1991). Harvesting commonly occurs by means of clear-cutting, and the extent usually expands across entire hillslopes (Fig. 2; Video supplement S1, (Tolorza, 2023)). This practice is permitted by current Chilean law, which mandates environmental impact assessments for clear-cutting only when harvest areas are ≥ 500 or ≥ 1000 ha yr⁻¹ in Mediterranean and temperate regions, respectively (Artículo Primero, Título I, Artículo 3, m.1 in accordance with Chilean Law 19.300, 2013). As a consequence, the CCR ranks among regions with the highest forest loss and gain worldwide (Hansen et al., 2013).

Tree plantations are frequently intersected by dense networks of logging roads. These roads facilitate access to heavy forest machinery, storage, and transport of timber, as well as subsequent (re)plantation. Logged hillslopes, such as logging roads, are important sediment sources and routes during storms and after wet-season clear-cutting (Schuller et al., 2013, 2021; Aburto et al., 2020). This is not surprising since they remain bare and prone to compaction by heavy-machinery transit. Post-harvest erosion is mainly triggered by rainfall (Aburto et al., 2020; Schuller et al., 2013) after exceeding soil-hydrology-specific rainfall intensity thresholds (Mohr et al., 2013). The erosional work efficacy depends on the logging season and is higher for wet-season logging (Mohr et al., 2014). Nevertheless, post-harvest erosion may last more than one season. Indeed, Aburto et al. (2020) reported the highest post-harvest soil loss in a catchment sustaining a 1-year-old plantation.

During the plantation cycle, roads are the prime sources of and pathways for sediments in catchments covered by tree plantations (Schuller et al., 2013). These roads often intersect streams, which form bypasses that preferentially transport sediment (Fig. 2), increasing the efficacy of mass transfer within a geomorphic system. This is also called sediment connectivity (Wohl et al., 2019). In this case, road networks modify the pathways of runoff and sediments and may also modify the thresholds of rainfall needed to trigger sediment detachment and transport (for example, due to soil compaction), potentially affecting both the structural and functional components of sediment connectivity, as defined by Wainwright et al. (2011). This shift is also relevant to constraining off-site impacts of soil erosion (Boardman et al., 2019).

Despite the increase in structural connectivity, sediment mobilization depends mostly on specific thresholds of rainfall. For example, to initiate runoff in recently logged areas, (functional) hydrologic connectivity required a threshold of 20 mm h⁻¹ in rainfall simulations conducted on tree plantations near Nacimiento (Mohr et al., 2013). In the absence of long-term records of rainfall intensities, hydroclimatic trends in rainfall and streamflow are relevant to the context of catchment erosion. In central Chile (30–39° S), rainfall decreased at a rate of approximately 4 % per decade between 1960 and 2016 (Boisier et al., 2018), culminating in an unprecedented megadrought, which started in 2010 (Garreaud et al., 2020).

While the erosional impacts of logging are largely indisputable, hydrologic responses to tree harvesting are ambiguous. On the one hand, logging may increase streamflow discharge, particularly for peak flows (Iroumé et al., 2006). On the other hand, logging may decrease streamflow discharge due to enhanced groundwater recharge immediately after logging (Mohr, 2013). The specific response most likely varies with tree species and age; harvest size; forestry treatment (thinning, clear-cutting, and replanting); riparian-buffer width; and, in particular, soil water storage decrease caused by recent drought conditions, which exacerbated declines in runoff (Iroumé et al., 2021).

In the context of the megadrought, recent increases in both the magnitude and frequency of wildfires have affected more tree plantations than alternative types of land cover (Bowman et al., 2019). This is likely because fuel is more abundant under dense plantation cover, linking extensive and uninterrupted stretches of the landscape. In contrast, native species distribute more patchiness (Gómez-González et al., 2017, 2018).

While the observed disturbances to vegetation cover commonly increase soil loss (Aburto et al., 2020) and fluvial-sediment yields (e.g., Reneau et al., 2007; Brown and Krygier, 1971), the long and persistent decline in rainfall (Méndez-Freire et al., 2022; Tolorza et al., 2019), together with the high water demands of tree plantations, is expected to reduce water availability for both sediment detachment and subsequent mobilization. In most fluvial catchments, long-term (10³–10⁴ years) denudation rates exceed short-term rates (Covault et al., 2013). This situation, however, may flip if anthropogenic soil erosion is high (Hewawasam et al., 2003; Vanacker et al., 2007). To evaluate this conundrum, we explore the catchment-scale erosion and denudation of the Purapel river. To this end, we establish a long-term benchmark based on detrital ¹⁰Be denudation rates for comparison with recent sediment yields. We also include discrete disturbance events (i.e., wildfires that occurred in 2017 and 2015), wet-season logging, and (dynamic) sediment connectivity associated with forestry roads. In a recent analysis of this river, Pizarro et al. (2023) concluded that the concomitant afforestation and the decrease in sediment export out of the catchment may form a causal relationship. Here, we analyze a wide range of meteorological and vegetational indexes over a 42-year period. We show that causes other than afforestation may explain the observed decline in sediment discharge, particularly the drought that this catchment has been experiencing since 2010.

2 Materials and methods

2.1 The Purapel catchment

The Purapel river drains the eastern flank of the CCR. The climate exhibits a Mediterranean type. The mean annual rainfall is 845 mm, and the mean minimum and maximum air temperatures are 7.2 and 20.3 °C, respectively. The hydrological regime is exclusively pluvial (Álvarez-Garreton et al., 2018). The catchment is 406 km² and dominated by metamorphic (48 %) and granitic (44 %) lithologies. Elevation ranges from 164 to 747 m above sea level (a.s.l). Most hillslopes are gentle (with gradients around 16 %), largely convex, and incised by gullies that converted this landscape into badlands (Fig. 1c). CIREN (2021) classified most of these hillslopes as severely affected by soil erosion. The dominant soil types are Inceptisols and Alfisols (Bonilla and Johnson, 2012). Soil properties are highly variable across space. Soils under tree plantations are generally thinner, are more depleted in soil organic matter, and exhibit lower invertebrate diversity than soils under native forests (Cifuentes-Croquevielle et al., 2020; Soto et al., 2019). Throughout the entire soil column, the soil bulk densities of Eucalyptus spp. stands (1.38 ± 0.08 to 1.58 ± 0.12 g cm⁻³) and Pinus radiata stands (1.28 ± 0.18 to 1.53 ± 0.13 g cm⁻³) are 40 %–80 % higher than in native forests (0.89 ± 0.27 to 1.25 ± 0.24 g cm⁻³) (Soto et al., 2019). After the mulch layer burned away in the 2017 wildfire, the bared surfaces exposed coarse topsoils (Fig. 1d).

2.2 Hydrometeorological data and analysis

We downloaded daily rainfall and potential-evapotranspiration (PET) data (1979–2020) for the Purapel catchment from the CAMELS-CL data set (https://camels.cr2.cl/, last access: 30 March 2022). Daily rainfall data are derived from the CR2MET v2 precipitation product, which combines ERA5 reanalysis with local topographic data calibrated by an updated national rain-gauge network (DGA, 2017). In addition, we downloaded daily rainfall time series (1956–2019) for the Nirivilo rain gauge from Mawün (https://mawun.cr2.cl/, last access: 30 March 2022). We analyzed trends in annual and seasonal rainfall for several periods using both data sets. We calculated annual aridity trends using an aridity index (AI), computed as the annual rainfall divided by the annual potential evapotranspiration ($\mathrm{AI}=P/\mathrm{PET}$), calculated using the CAMELS-CL data set. For the annual analysis, we considered years with more than 330 daily data points and months with more than 27 daily data points. For the seasonal analysis (∼ 90 d), we included data for seasons with more than 60 daily data points – that is, data for more than ∼ 66.6 % of the days for an assumed 3-month season. We defined fall as March–May (MAM), winter as June–August (JJA), spring as September–November (SON), and summer as December–February (DJF). For the pre- and post-fire hydrologic years (March–February), we analyzed daily rainfall as well. To this end, we used the gap-free gridded CR2METv2 product.

Streamflow and suspended-sediment data are available on the Chilean General Water Directorate (DGA) site (https://mapas2.mop.gob.cl/, last access: 1 July 2024). The DGA estimated the daily streamflow discharge (Q) based on single gauge-stage readings using calibrated rating curves. Roughly once a month, the rating curves are updated with current manual meter measurements. The suspended-sediment concentrations (SSCs) were also sampled on a daily basis. All samples were obtained close to the water surface in the vicinity of the water stage. The samples were filtered using cotton linter cellulose filter papers with an 80 % collection efficiency for particles larger than 0.3 µm (Advantec Qualitative Filter Papers (Grade 2); written communication from the DGA operator, 2021). They were then dried and burned for 2 h at 550–600 °C in DGA laboratories (Solar, 1999).

We calculated the daily suspended-sediment discharge (SSD; t d⁻¹) as the product of Q (m³ s⁻¹) and SSCs (mg L⁻¹), assuming these measurements to be representative of the entire day (Pepin et al., 2010). This assumption is based on the catchment size and measurements of streamflow velocities taken in this landscape (Sect. S1 in the Supplement). In addition, we calculated the number of data points, the percentiles, and the mean SSD value for individual hydrologic years. We calculated the daily baseflow at the Purapel en Sauzal station using a Lyne and Hollick filter (Ladson et al., 2013), which is a standard approach used in several studies (e.g., Li et al., 2022; Huang et al., 2021; Teutschbein et al., 2015; Zhang et al., 2017). For baseflow separation, we used several α values between 0.5 and 0.95 and “n.reflected” = 30 d as parameters.

The Q and SSC time series contain gaps that are not seasonally clustered. In particular, gaps during the dry season are mostly related to ceased Q (personal communication from the DGA operator, 2021). We calculated the mean value and the number of daily Q, baseflow, and SSC data points for each month and season. Based on the number of data (Sect. S2), we discarded the monthly analysis of trends. We then calculated seasonal trends in streamflows and suspended sediment for seasons with more than 60 data points. All trends were computed with the Mann–Kendall test (Helsel et al., 2020) and plotted using additional locally weighted scatterplot smoothing (LOWESS) (Cleveland, 1981). In addition, we computed the quantiles and averages of SSD for hydrologic years with more than 185 data points.

We also examined the existence of differences between the daily data points for the hydrologic years before and after the fire, i.e., 2014 and 2016 versus 2015 and 2017, respectively (n > 262 for each group of data). Given the non-normal distribution of all hydrological data sets, we used the Wilcoxon–Mann–Whitney test instead of the traditional t test. This test requires the two compared populations to have the same shape and variance (homoscedasticity). We used the non-parametric Fligner–Killeen test to determine whether our samples were homoscedastic (p value larger than 0.05) or not (p value lower than 0.05). If the compared populations exhibited unequal shape or variance (heteroscedasticity), we used Welch's t tests (Skovlund and Fenstad, 2001).

2.3 Catchment-wide erosion and denudation rates

We obtained catchment-wide erosion rates for the Purapel river at the Purapel en Sauzal gauge using two approaches for different timescales – short-term erosion (decennial timescale) from suspended sediment and long-term erosion (10³ to 10⁴ years) from detrital ¹⁰Be. A limitation of our approach is the fact that detrital ¹⁰Be rates include physical-erosion and chemical-weathering rates (von Blanckenburg and Willenbring, 2014), while suspended-sediment yields account only for the physical erosion of very fine sediment (Summerfield and Hulton, 1994), which excludes bed loads and dissolved loads. Thus, we regard our short-term erosion rates as the minimum rates for landscape lowering.

For short-term erosion, we calculated the mean of the specific SSD ($\mathrm{t}{\mathrm{km}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) as the average of all records (June 1985 to November 2018) on a yearly scale, normalized by catchment area (Pepin et al., 2010). We estimated the resulting erosion rate (mm yr⁻¹), assuming a mean soil bulk density of 2.6 g cm⁻³ (Carretier et al., 2018).

For long-term erosion, we assumed that the ¹⁰Be concentrations within fluvial sands are proportional to the catchment-wide-averaged denudation rate (von Blanckenburg, 2005; Granger and Schaller, 2014). This rate integrates over a characteristic timescale that is inversely proportional to the denudation rate. These timescales are commonly longer than 1000 years (Covault et al., 2013). We therefore regard the ¹⁰Be-derived rates as a reference that largely excludes recent human disturbances but includes low-frequency and high-magnitude erosion events (Kirchner et al., 2001; Carretier et al., 2013). We obtained a bulk sample of fluvial sands from the active riverbed along a cross section close to the “Purapel en Sauzal” water stage, collecting sands from the surface at three locations within ∼ 10 m of distance. We mixed all samples and sieved them to a grain size fraction of 0.5–1 mm.

The mixed-sand sample was processed at the French accelerator mass spectrometry (AMS) ASTER facility at the Centre Européen de Recherche et d'Enseignement des Géosciences de l'Environnement (CEREGE). In order to convert ¹⁰Be concentration, C, into a mean denudation rate for the catchment, we neglected radioactive decay and assumed a steady state of ¹⁰Be concentration, leading to the following classical equation:

$$\begin{array}{}\text{(1)}& \mathit{ϵ}={\displaystyle \frac{\mathrm{1}}{\mathit{\rho }C}}{P}_{\mathrm{SLHL}}(f_{\text{sp}}{S}_{\text{sp}}{\mathrm{\Lambda }}_{\text{sp}}+f_{\text{sm}}{S}_{\text{sm}}{\mathrm{\Lambda }}_{\text{sm}}+f_{\text{fm}}{S}_{\text{fm}}{\mathrm{\Lambda }}_{\text{fm}}),\end{array}$$

where P_SLHL = 4 $\mathrm{at}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ represents the sea-level high-latitude total production rate of the considered nuclide (Martin et al., 2017). Moreover, f_sp, f_sm, and f_fm are the fractions of this production rate due to spallation, slow-muon capture, and fast-muon capture averaged over the catchment area, respectively (Braucher et al., 2011). Furthermore, S_sp, S_sm, and S_fm are scaling factors that depend on latitude and elevation averaged over the catchment area (Stone, 2000), and ρ = 2.6 g cm⁻³. No geometric-shielding correction for topography was applied (horizon < 20° in all directions). The uncertainty in the denudation rate accounts for the propagation of the analytical uncertainty and an assumed uncertainty of 15 % in the production rate.

2.4 Land-cover changes

The Purapel catchment has experienced high rates of LULCC since the 19th century. This has been largely due to the extensive increase in wheat production caused by the gold rushes in California and Australia (Cortés et al., 2022). Later on, between 1955 and 2014, tree plantations increased from a minimum of 10.27 (Hermosilla-Palma et al., 2021) to 203.5 km² (Zhao et al., 2016). Recently, two large wildfires burned the catchment. In 2015, 14 % of the catchment area burned, and in 2017, almost the entire catchment burned (95 %) (Tolorza et al., 2022a).

To describe recent LULCC in this catchment, we use land-cover maps from compiled sources (1955, 1975, and 2017) and from our own sources (1986, 2000, 2005, 2010, and 2015):

• The 1955 and 1975 land-cover maps of Hermosilla-Palma et al. (2021) cover the headwaters of the Purapel catchment (157 km²). The 1955 map was created by interpreting land cover from a 1 : 70 000 aerial photograph (Hycon flight), and the 1975 map was created by interpreting land cover from the 60 m resolution Landsat-2 multispectral scanner (MMS) and 1:30 000 aerial photographs from 1978 (CH-30 flight).

• We used Landsat surface reflectance products to identify land-cover classes during the dry seasons of 1986, 2000, 2005, 2010, and 2015. We classified unburned land cover using a maximum-likelihood classifier (Chuvieco, 2008) which we trained and validated with 20 polygons for one class and 10 polygons for the other. We validated the results with field observations from 2014–2015. We subclassified the burned surfaces into low, moderate, and severe fire categories based on differences in the normalized burn ratio (NBR) index between pre- and post-fire images (thresholds: 0.1–0.269, 0.27–0.659, and 0.66–1.3) (Key and Benson, 2006).

• The 2017 land-cover map was made by Tolorza et al. (2022a) using pre-fire Sentinel and lidar data. Here, this classification was resampled to a 30 m resolution to be compatible with Landsat classifications.

2.5 Logging roads and sediment connectivity

To identify changes in structural connectivity, we applied the connectivity index (IC; dimensionless) using the weighting factor (W; dimensionless) from Cavalli et al. (2013). The IC is a semi-quantitative approach to describe the degree of coupling between hillslopes and a target (for example, the stream network) and is expressed as

$$\begin{array}{}\text{(2)}& \mathrm{IC}={\log }_{\mathrm{10}}\left({\displaystyle \frac{\stackrel{\mathrm{‾}}{WS}\sqrt{A}}{{\displaystyle }\sum _i{\displaystyle \frac{d_i}{W_i{S}_i}}}}\right),\end{array}$$

where $\stackrel{\mathrm{‾}}{W}$ and $\stackrel{\mathrm{‾}}{S}$ (m m⁻¹) are the average weighting factor and slope gradient of the upslope contributing area (A; m²), respectively. Moreover, d_i (m), W_i (dimensionless), and S_i (m m⁻¹) are the path length, the weighting factor, and the slope gradient for the ith cell downslope towards a target.

W is calculated from a DTM to account for the effect of topographic roughness. The roughness index (RI; m) is the standard deviation of the residual topography. The residual topography refers to the difference between the original DTM and a smoothed version obtained by averaging DTM values within a moving window with 5 × 5 (25) cells. The RI is given as

$$\begin{array}{}\text{(3)}& \mathrm{RI}=\sqrt{{\displaystyle \frac{{\displaystyle }\sum _i^{\mathrm{25}}(x_i-x_{\mathrm{m}}{)}^{\mathrm{2}}}{\mathrm{25}}}},\end{array}$$

where x_i (m) is the value of one specific cell of the residual topography within the moving window and x_m (m) is the mean of all 25 window cells. The weighting factor is calculated as

$$\begin{array}{}\text{(4)}& W=\mathrm{1}-{\displaystyle \frac{\mathrm{RI}}{{\mathrm{RI}}_{\text{max}}}},\end{array}$$

where RI_max is the maximum value of RI in the study area.

The routing of sediments on hillslopes is likely to change where a new temporary sink is closer downslope. Both streams and forestry roads may behave as temporary sinks of sediments after detachment from hillslopes (Schuller et al., 2013, 2021), and this behavior can be addressed using these elements (streams and forestry roads) as targets for the IC. We propose here an approach to estimate the related changes in sediment connectivity by considering the difference between two different targets, which is expressed as

$$\begin{array}{}\text{(5)}& \mathrm{RC}={\mathrm{IC}}_{\text{rs}}-{\mathrm{IC}}_{\mathrm{s}},\end{array}$$

where the subscripts of IC refer to the targets for its computation – rs refers to the network formed by streams and roads, and s refers to the network formed only by streams. We supplied the model with a mapped forestry road network obtained from images available from the OpenLayers plug-in of QGIS and post-2017-fire Sentinel compositions.

This approach only accounts for changes in sediment connectivity when forestry roads behave as sinks of sediments, omitting their role as sources.

2.6 Disturbances in vegetation

We used the Breaks For Additive Season and Trend (BFAST) algorithm (Verbesselt et al., 2010) on a Landsat collection to detect disturbances in vegetation at the pixel scale, i.e., ≥ 30 m. In the Purapel catchment, disturbances > 30 m are mostly due to wildfires and/or clear-cuts. Such disturbances rely on the seasonal behavior of the normalized difference vegetation index (NDVI) in a time series of Landsat surface reflectance (Level 2, Collection 2, Tier 1) for the period from September 1999 to October 2021. Clouds were filtered using the QA band, which uses the C Function of Mask (CFMask) algorithm (Foga et al., 2017). We used the same parameter as Cabezas and Fassnacht (2018), setting the threshold value for disturbances based on 93 manually labeled reference polygons with fire events, clear-cuts, and constant tree cover. It is worth mentioning that we applied a sieve filter to the results. Hence, only disturbances greater than 1 ha were considered. We trained the algorithm using the Landsat time series from 1999 to 2001. Given the disturbance regime of the Purapel catchment (two large wildfires and potential loggings every 9 to 25 years), we ran the BFAST algorithm, anticipating three possible breaks for the period 2002–2021. The accuracy assessment was performed on 35 manually drawn polygons that were randomly distributed across the catchment. Since we were dealing with imbalanced data (looking for disturbances assumed to be anomalies in the time series), we used balanced accuracy and F1 scores as the metrics for evaluating our data (Brodersen et al., 2010). Balanced accuracy was calculated as the arithmetic mean of sensitivity and specificity, while the F1 score was calculated as

$$\begin{array}{}\text{(6)}& \mathrm{F}\mathrm{1}\mathrm{score}=\mathrm{2}\cdot \left[\right(\text{precision}\cdot \text{recall})/(\text{precision}+\text{recall}\left)\right],\end{array}$$

where “precision” represents the amount of correct positive predictions ($\text{true\hspace{0.17em}\hspace{0.17em}positives}/(\text{true\hspace{0.17em}\hspace{0.17em}positives }+\text{ false\hspace{0.17em}\hspace{0.17em}positives})$) and “recall” represents how many positive predictions the model made divided by all positive cases ($\text{true\hspace{0.17em}\hspace{0.17em}positives}/(\text{true\hspace{0.17em}\hspace{0.17em}positives }+\text{ false\hspace{0.17em}\hspace{0.17em}negatives})$).

Figure 3 The annual and seasonal rainfall and the annual aridity index (AI = $P/\mathrm{PET}$) at the Purapel catchment. The main monotonic trends are tested using Mann–Kendall tests and LOWESS for 1962–2015 (blue), 1979–2019 (red), 1986–2018 (purple), and 2000–2018 (green). The unfilled circles represent discarded data. (a) The annual rainfall and AI time series. Note that “P” and “St” stand for Purapel and station, respectively. (b) The seasonal time series for the Nirivilo station.

3 Results

3.1 Hydroclimatic trends

At the annual scale at Nirivilo rainfall station, most data for the period 1962–2015 passed our completeness assessment criteria (53 out of 54 years). In the case of CR2MET, the longest period analyzed here is 1979–2019. Using Mann–Kendall tests and LOWESS, we did not find a single trend within the longest interval of records (1962–2015). Nevertheless, for the period after 1979, we obtained decreasing non-monotonic trends in rainfall at both the Nirivilo station and for the CR2MET product. This decrease is steeper for 2000–2019 but less pronounced for intermediate intervals, such as 1986–2018. During 1986–2018, however, a decrease in seasonal rainfall is observed for fall, which is at the beginning of the hydrologic year (Fig. 3). Generally, the AI follows similar decreasing trends to those pertaining to rainfall, indicating persistently dry conditions across this catchment. In only 2 separate years (1982 and 2002), the AI was higher than 1. During all other years, potential evapotranspiration was greater than rainfall.

Figure 4 The mean seasonal streamflow, baseflow, and suspended-sediment concentrations at the Purapel en Sauzal station on an annual basis. The main monotonic trends are tested using Mann–Kendall tests and LOWESS for 1986–2018 (purple) and 2000–2018 (green). The unfilled circles represent discarded data. (a) The mean seasonal streamflow at the Purapel en Sauzal station. (b) The mean seasonal baseflow at the Purapel en Sauzal station. (c) The mean seasonal suspended-sediment concentration at the Purapel en Sauzal station.

Streamflow data are available at the Purapel en Nirivilo station between 1979 and 2019 and at the Purapel en Sauzal station between 1981–2019. For the Purapel en Sauzal streamflow data, the results of baseflow separation are in Sect. S3 in the Supplement. We selected the results obtained using α = 0.7 for further trend analysis, considering the observed magnitudes and shape of the baseflow time series. For the data on suspended-sediment concentration, the longest period is 1985–2018. The seasonal analysis for the Purapel en Sauzal station is shown in Fig. 4. Although none of these time series are monotonic, the sharp decrease in suspended-sediment concentrations is clear for the three wetter seasons (fall, winter, and spring).

Carretier et al. (2018)Carretier et al. (2018)Carretier et al. (2018)

Table 1 Published and new detrital ¹⁰Be denudation rates in the Mediterranean CCR. Denudation rates and their uncertainties were calculated using procedures described in Carretier et al. (2018). Characteristic (char.) time refers to the quartz residence time within the mean free path of particles in rocks measuring 60 cm and represent a timescale for steady erosion (Lal, 1991). Details of ¹⁰Be concentrations and blanks analyzed in the CEREGE laboratory are available in Table S1 in the Supplement.

3.2 Catchment-wide erosion and denudation rates

The ¹⁰Be denudation rate was determined to be 0.024 ± 0.004 mm yr⁻¹ (Table 1), assuming a soil particle density of 2.6 t m⁻³. This rate translates into a specific sediment yield of 62.4 ± 10.4 $\mathrm{t}{\mathrm{km}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. This rate integrates over a characteristic timescale of ∼ 25 kyr. Together with published data, detrital ¹⁰Be denudation rates in the CCR range between 0.02 and 0.05 mm yr⁻¹ (Table 1;  Carretier et al., 2018).

Figure 5 Denudation and specific SSD (Spec Susp Sed Disch) at the Purapel en Sauzal gauge. Distributions of SSD for individual hydrologic years (March to February). The purple circles represent the percentiles (0.25, 0.5, 0.75, and 0.95), and the purple triangles each represent the mean of the daily data from each hydrologic year. The filled symbols represent years with more than 185 daily data points. Catchment erosion and denudation rates are indicated in red. The solid line represents the sediment yield equivalent to the ¹⁰Be denudation rate, while the dashed line represents the average of all suspended-sediment records. Both the daily and annual scales are expressed to facilitate comparisons. MK tau: Mann–Kendall tau; pval: p value.

Following Pepin et al. (2010), we calculated the mean of the specific SSD for all records between 1985 and 2018 as well as a 30 % of error (Pepin et al., 2010). For the Purapel catchment, we estimated 47 ± 14.1 $\mathrm{t}{\mathrm{km}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, equivalent to 0.018 ± 0.005 mm yr⁻¹, assuming the same soil bulk density. Both rates do not statistically differ (Fig. 5).

3.3 Decennial trends in SSD and comparison of pre- and post-fire hydrometric data

On a decennial timescale, we observe decreasing trends in mean and high (95th and 75th percentile) annual SSD values, based on Mann–Kendall p values for α = 0.05 (Fig. 5). Since 2010, the high and medium percentiles of SSD have been lower than most of the previous hydrologic years.

Figure 6 Daily hydrometric data for pre- and post-fire hydrologic years for the two large wildfire events. Rainfall, Q, SSC, and SSD from pre- and post-fire hydrologic years were compared in terms of homoscedasticity using the Fligner–Killeen test with α = 0.05. Heteroscedastic distributions were compared using Welch's t tests. Homoscedastic distributions were compared using Wilcoxon–Mann–Whitney tests. Eff. size: effect size. OA: overall accuracy.

On a shorter timescale, in Fig. 6, we compare pre- and post-fire hydrometric data for the 2015 and 2017 wildfires. The compared series were found to be homoscedastic for rainfall in both comparisons and for SSC in the comparison before and after the 2015 wildfire. No changes in median rainfall values can be inferred from the large p value of the Wilcoxon–Mann–Whitney tests. Before and after the 2015 wildfire, considering p values when α = 0.05 in Welch's t tests or Wilcoxon–Mann–Whitney tests, mean streamflow diminished, but the median SSC and mean SSD increased. Before and after the 2017 wildfire, mean streamflow increased, as did mean SSC and mean SSD.

3.4 Recent land-cover changes

We developed five land-cover maps for the period 1986–2015. The overall classification accuracy ranged between 83 % and 92 %. We distinguished between tree plantations, native forests, shrublands, and seasonal grasslands. Seasonal grasslands included bare surfaces, seasonal pasture, and sparse vegetation. We also classified seasonal grasslands to separate recently logged areas (clear-cuts) from other poorly vegetated areas.

Figure 7 Land-cover classification and transitions. Maps for 1955 and 1975 are from Hermosilla-Palma et al. (2021), maps for 1986–2015 are from this work (Sotomayor et al., 2022), and the 2017 map is from Tolorza et al. (2022a).

Figure 7 shows that the upper catchment was covered by a minimum of 1000 ha of tree plantations and 5500 ha of shrublands in 1955 (Hermosilla-Palma et al., 2021). Between the 1980s and the beginning of the 21st century, the most prominent change comprised the transition from seasonal grasslands and shrublands to tree plantations. The first two classes covered a minimum of 23 550 ha in 1986 and 13 050 ha in 2005. During the same period, tree plantations expanded from 8090 to 20 980 ha. Between the wildfires of 2015 and 2017, seasonal grasslands and shrublands together expanded to ∼ 20 300 ha (Fig. 7).

Figure 8 Sediment connectivity index (Cavalli et al., 2013) calculated using (a) the streams and (b) the streams and forest roads as targets. (c) Road connectivity (RC) is the difference between both models. A histogram illustrating RC is displayed in the bottom-left corner. Insets show the locations of maps in panels (d–g). (d) Detailed view of RC as an example of RC values, using the same color table as in panel (c), in relation to forestry roads and streams. (e–g) Details of hilltops with the highest values of RC (©Google Maps 2022; related dataset can be found in Tolorza et al., 2022b).

3.5 Landscape disturbances

The result of our mapped road network is illustrated in Fig. 1. Using this road network on a 5 m resolution lidar, we estimate approximately 18 000 ha of increased sediment connectivity (RC > 0) and 150 ha of decreased sediment connectivity (RC < 0) (Fig. 8). RC values exceeding the 95th percentile (> 3.12) cover 1986 ha. This area of high RC is primarily located on hilltops, with 1966 ha having an RC > 3.12 (i.e., 99 %) in upstream contributing areas < 1 ha. In particular, these topographic settings exceed an empirical threshold between high and low connectivity for a mountain catchment, i.e., −2.32 (Martini et al., 2022). In the Purapel catchment, the area above this threshold increased from 1120 to as much as 6570 ha simply due to the dense road network. This quantification, however, was performed with a digital terrain model with a coarser resolution compared to that used in the original study by Martini et al. (2022) – 5 m versus 0.5 to 2.5 m. The few pixels resulting in negative RC are located in first-order channels omitted from the stream network, specifically where these channels intersect more than one road or stream (Fig. 8d).

Figure 9 Detected disturbances from our BFAST modeling. (a) Map showing the number of disturbances in vegetation detected for the period 2002–2019. (b) Seasonality of the disturbance area detected within the Purapel catchment.

Based on our BFAST modeling, we obtained monthly time series of disturbances for 2002–2019, which we aggregated to the seasonal scale (Poblete-Caballero et al., 2022; Poblete-Caballero and Tolorza, 2023). We achieved a confusion matrix with a balanced accuracy of 0.86 and an F1 score of 0.69. For the complete period (Fig. 9a), 13 640 ha (33.7 %) of the Purapel catchment experienced one break in the NDVI time series, 16 810 ha (41.5 %) showed two breaks, and 5010 ha (12 %) presented three breaks. The undisturbed 12.8 % included tree plantation stands that remained unlogged and seasonal grasslands that remained poorly vegetated. Taking the seasonality into account (Fig. 9b), the periods with the largest disturbed areas were the summers (dry seasons) of 2015 and 2017 due to the wildfires. These wildfires were detected in ∼ 5000 and 24 000 ha, respectively. The subsequent disturbed areas occurred during the fall (wet season) of 2002 (∼ 2000 ha) and 2007 (1910 ha). The largest surface disturbances during winter and spring both corresponded to 770 ha in 2006 and 2009, respectively.

In comparison to the classification based on the differential normalized burn ratio (dNBR) for 2017 (Tolorza et al., 2022a), the BFAST results detected fewer burned areas for the 2017 wildfire (33 618 versus 24 299 ha). This difference could be explained by the better capabilities of the dNBR index with regard to detecting burned areas, since it is a dedicated method for classifying burned surfaces based on the NBR index of pre- and post-fire images (Key and Benson, 2006). In contrast, the BFAST algorithm applied here uses the NDVI, which is a more suitable index for detecting the density of vegetation and is thus more sensitive to clear-cuts.

4 Discussion

4.1 Long- and short-term catchment erosion

Both the ¹⁰Be denudation rate and suspended-sediment erosion rate are surprisingly similar (Fig. 5). Thus, the decennial average of the suspended sediment may capture at least the effects of erosion events recorded with ¹⁰Be in the long term. Both rates are low for fluvial catchments between 100–1000 km² (Covault et al., 2013). Yet, these rates are similar to those observed in three tributaries of the Biobío river that drain the eastern CCR, which range between 0.037 ± 0.006 and 0.042 ± 0.008 mm yr⁻¹ (Carretier et al., 2018).

The low ¹⁰Be denudation rate aligns with a landscape dominated by slow soil creep, with occasional mass wasting triggered by intense rainfall events and earthquakes. This long-term rate, however, may be overestimated: after 200 years of intense soil erosion, deep saprolite layers with low concentrations of ¹⁰Be are widely exposed at the Earth's surface. Those hillslopes might be depleted in ¹⁰Be, which may lead to an overestimation of the total denudation (Schmidt et al., 2016).

Short-term erosion does not exceed long-term denudation, as observed in other highly human-disturbed catchments (Hewawasam et al., 2003; Vanacker et al., 2007). However, we regard the decennial suspended-sediment erosion rate as a very conservative estimate for recent catchment-scale erosion for three reasons.

First, we argue that this estimate does not account for the possible effects of sediment storage and other transient processes that affect soils and streams. Over a certain period of time, only a fraction of the sediment detached within a catchment is transported to the outlet and, thus, captured in the sediment yield due to processes such as deposition and storage (Walling, 1983). The sediment delivery ratio (SDR) describes this fraction and is expressed as $\mathrm{SDR}=\mathrm{SY}/E$, where SY is the average annual sediment yield per unit area and E is the average annual erosion over that same area. Considering empirical data (Walling, 1983; Ferro and Minacapilli, 1995; Lu et al., 2006), an SDR of 0.7 may be used as a conservative first approach for catchments larger than 100 km², although this value may be a couple of orders of magnitude lower. Using the mean of the specific SSD as SY in the previous expression, our estimation of decennial erosion in the Purapel catchment increases to 67 ± 20 $\mathrm{t}{\mathrm{km}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (or 0.026 ± 0.007 mm yr⁻¹), which is even closer to the long-term denudation rate. Decennial catchment erosion may be much higher if the SDR is a few orders of magnitude lower. More reliable estimates of the SDR can be obtained using spatially distributed models, such as the linear lumped series model for hillslope storage and channel storage from Lu et al. (2006). This model considers a theoretical relationship between sediment delivery and certain physical parameters, such as particle size and rainfall duration.

Second, previous studies on rivers in the western Andes indicate that catchment-scale erosion from gauging data may be underestimated due to undersampling of extreme events (Vanacker et al., 2020; Carretier et al., 2018). Given the numerous data gaps in both streamflow and suspended sediments, we cannot rule this out. Also, we do not have subdaily or depth-integrated measurements of sediment concentrations.

Third, suspended sediments do not capture the effects of chemical weathering on denudation rates. This process seems to be relevant in the CCR. In the absence of spatially resolved data on regolith thickness, isolated observations indicate thick saprolite layers (at least) locally (Vázquez et al., 2016; Mohr et al., 2012; Krone et al., 2021), while weathering rates at plot scales range between 0 and ∼ 240 $\mathrm{t}{\mathrm{km}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ across other latitudes in the CCR (Schaller and Ehlers, 2022), i.e., 26, 29, 32, and 37° S). Thus, depending on the magnitude of mass loss due to chemical weathering, total denudation in the short term can be equal to or even higher than the reported long-term denudation. Local quantitative estimates of chemical weathering and soil production rates would help to constrain the magnitude of the difference between long- and short-term denudation. Indeed, deep chemical weathering may also lead to underestimating long-term denudation rates with ¹⁰Be. Such underestimation may occur if mass loss is significant below the depth at which most cosmogenic nuclides are produced, as observed in tropical landscapes dominated by dissolution (Campbell et al., 2022).

4.2 Disturbances and trends in suspended sediments

The Purapel river catchment has served as a staging ground for the rapid expansion of tree plantations and has experienced a number of disturbances during the period of suspended-sediment monitoring. This landscape has been affected by clear-cuts, two widespread wildfires, and one earthquake with an M_w value of 8.8. The expansion of tree plantations has mostly occurred at the expense of poorly vegetated surfaces (Fig. 7). Yet, plantation management includes extensive logging operations (Video supplement) – mostly conducted during wet seasons (Fig. 9) – and the construction and maintenance of forestry roads used by heavy machinery.

The distribution and density of a forestry road network significantly affect the routing of sediments because the roads may act as sources or sinks of the particles detached from hillslopes. Our results show changes in structural sediment connectivity across a surface measuring 18 150 ha, where roads act as target (i.e., behave as sinks of sediments). Most of the surface has experienced an increase in structural sediment connectivity and thus an increase in sediment-routing capacity downslope (Fig. 8). Consequently, detached soil is now (relatively) well connected to local sinks, even from hilltops (surfaces with upstream contributing areas < 1 ha). In this specific position in the landscape, soil moisture is likely at its lowest, and, consequently, we can assume that the soil production rate is slower compared to that in the mid-slope or toe positions across the CCR (Schaller and Ehlers, 2022). Thus, hilltop soils may be more difficult to recover over human timescales.

The decennial trends in mean and high suspended-sediment discharge (Fig. 5) and mean SSC during wet seasons (Fig. 4) are decreasing. Despite the disturbances, which cluster towards the end of our time series, the changes in post-disturbance suspended sediments did not overshadow this decennial trend. For example, after the 2015 wildfire, we found an increase in sediment concentrations but a decrease in mean streamflow (Fig. 6). Thus, we argue that the elevated sediment concentrations reflect the erosional signal of the wildfire. Nevertheless, the described changes in pre- and post-fire suspended sediments are limited since the high and medium percentiles of suspended-sediment discharge following the fire are lower than during most of the hydrologic years prior to 2010, the beginning of the megadrought (Fig. 5). Assuming that the suspended-sediment record is representative of the sediment yields from the Purapel river, the disturbance regime contrasts with the expected sediment mobilization reported in other landscapes (e.g., Reneau et al., 2007; Brown and Krygier, 1971).

4.3 Water availability, suspended sediments, and sediment storage

The arid conditions, i.e., the ratio between annual precipitation and evapotranspiration, indicate increasing water scarcity. Such increasing water scarcity is consistent with the decline in streamflow and baseflow during the wet seasons, which is mostly pronounced in the fall. Apart from the unprecedented drought that started in 2010, the high root water uptake by fast-growing tree plantations is likely a key factor contributing to the severe reduction in water availability. Indeed, both tree plantations and drought conditions may reduce groundwater recharge (Iroumé et al., 2021; Huber et al., 2010). Also, soil loss due to erosion may further reduce the water storage capacity (Ratta and Lal, 1998). When considered together, enhanced water scarcity due to drought and tree plantations may have reduced groundwater storage. This is consistent with the observed negative trend in baseflow. Streamflow is largely composed of groundwater, and as groundwater declines, less streamflow is available to perform geomorphic work, which may reduce sediment transport.

Under the described conditions of drought and land use, the frequency of the minimum rainfall intensity required to trigger runoff and soil erosion on hillslopes may also be declining. At the same time, as rainfall and direct runoff exert control over sediment fluxes in streams (Andermann et al., 2012; Tolorza et al., 2014), sediment mobilization under the current hydrological regime may likely remain low despite landscape disturbances. Although we observed an increase in structural sediment connectivity, it is possible that sufficient water was unavailable to effectively establish sediment connectivity, leading to increased residence times of sediments stored within the catchment. Storage may occur in diverse sediment compartments, including hillslopes, forestry roads, tributary junctions, fans, and flood plains.

Indeed, a model of post-fire sediment cascades indicates that, while post-fire erosion may be severe in source areas, a substantial fraction of the detached sediments may intermittently remain stored within valleys, with only moderate delivery to the river network (Murphy et al., 2019). A sustained decline in available water for streamflow generation in such a system would likely increase the residence times of sediments.

4.4 Sediments on hillslopes and streams over time

The biogeochemical processes impacting vegetation and soils alike, the amount and grain size distribution of sediments supplied by hillslopes, and the dynamics of sediment mobilization in streams are all closely related (Terweh et al., 2021). More than 200 years ago, the primary native forests probably grew on deep and carbon-rich soils, with surfaces likely much more enriched in fine sediments than present topsoils. In contrast, the landscape-scale wildfire of 2017 revealed coarse saprolite that was widely exposed at the Earth's surface (Fig. 1d). Roots of fallen trees dig up coarse sediments spread out over most hillslopes. If we assume that landscape lowering has denuded 2 m of soil over 200 years, physical erosion should have experienced one or more stages of much higher rates in the past and after the beginning of intense deforestation. Rates of 10 mm yr⁻¹ are required to produce such levels of denudation over 200 years, which exceed both our short- and long-term estimates by 3 orders of magnitude. During intense soil erosion, the most likely evolution of grain size distribution in topsoils involves a quick depletion of fine sediments and a relative enrichment in coarse and more weather resistant sediments. Together with deforestation and wildfires, this time-dependent coarsening of the inorganic component of topsoils would be described as a reduction in the supply of fine sediments from hillslopes to streams, which may lead to a transition from transport-limited to supply-limited catchment erosion. Because vegetation and soils have been systematically disturbed, and given the high degree of coupling of these landscape elements with sediment mobilization in streams, it is unlikely that the grain size distribution in rivers remain constant. Since short-term and long-term erosion are comparable, we interpret that, after a rapid response of erosion to deforestation, denudation rates decreased due to exhausted sediments. Hence, it is plausible that part of the observed decline in suspended sediment is related to changes in the grain size distribution of sediments detached from soils. In this interpretation, both long- and short-term erosion can be considered supply-limited erosion. In this former, this is due to the scarcity of landslides, and in the latter, this is due to the slow rate of soil (and fine-sediment) production.

4.5 Tree plantations, native forests, and landscape degradation in the CCR

The expansion of tree plantations has been proposed as a tool to mitigate soil erosion (CONAF and MINAGRI, 2016). Recently, plantations have been favored as a better solution for mitigating soil erosion compared to native forests in the same Purapel catchment (Pizarro et al., 2020, 2023). In the Purapel catchment, direct comparison between native forests and plantations cannot be achieved for the period 1986–2018 because the major source of land cover transitioned from poorly vegetated surfaces to tree plantations (Fig. 7). Nevertheless, we can discuss whether the observed land management is a suitable solution for soil erosion mitigation in the CCR. There is abundant evidence of accelerated soil erosion in Chilean tree plantations, such as truncated soil profiles in Eucalyptus spp. stands at 36°37^′ S (Banfield et al., 2018); a 4-fold increase in net soil loss under Pinus radiata stands compared to native forests in Talcamávida (37°7'S) and Nacimiento (37°30'S) (Aburto et al., 2020); and changes in nutrient cycles and increased sedimentation rates in numerous coastal lakes, such as the Laguna Matanzas (33°45^′ S; Fuentealba et al., 2020), Lago Vichuquén (34° S; Fuentealba et al., 2021), Laguna Chica de San Pedro (36°51^′; Cisternas et al., 2001), and Lanalhue Lake (37° S; Alaniz et al., 2021). Based on such strong empirical evidence across the CCR as well as our own results (Figs. 2 and 7–9), we argue that the observed ongoing forest management of tree plantations promotes soil erosion and landscape degradation. In addition, soils in tree plantations are depleted in carbon and nutrients (Soto et al., 2019; Banfield et al., 2018) and inhibit lower-invertebrate diversity (Cifuentes-Croquevielle et al., 2020) compared to soils under native forests. As a result, C and N stocks are relatively lower in tree plantations, extending down to deep soil compartments (> 120 cm) (Crovo et al., 2021). Soil organic matter is a key component of soil formation (Bernhard et al., 2018). For this reason alone, native forests rather than exotic tree plantations are a more appropriate source of land cover for regenerating soils and reversing, or at least decelerating, 200 years of intense soil erosion. Indeed, the protection and conservation of natural vegetation have the strongest effect on improving soil quality after water erosion (Vanacker et al., 2022). Also, empirical restoration examples show that transitioning from former Eucalyptus spp. plantations to native forests is valuable in terms of improving water availability (Lara et al., 2021).

5 Conclusions

The Purapel catchment, like other similar catchments along the CCR, experiences slow denudation on millennial timescales. The decennial-averaged suspended-sediment discharge, while of similar magnitude, likely underestimates total denudation. Accordingly, depending on the magnitude of the unconstrained portion of the denudation, the decennial lowering of the Earth's surface may be equal to or even higher than the long-term average.

Suspended-sediment transport decreased during the wet seasons between 1986 and 2018, which, at first glance, conflicts with the disturbances observed in vegetation, especially the intense and widespread wildfires. The decrease in several hydroclimatic variables and forcings, including baseflow and aridity, coincides with lower suspended-sediment loads. We argue that the low range of recent suspended-sediment results from the observed decline in water availability, thus limiting the detachment and transport of sediments. Or, in other words, the drought offsets the effects of disturbances and higher connectivity. Without sufficient water, the residence times of sediments are long. The contribution of tree plantations to reducing erosion, if any, is more closely related to their impact on water availability rather than directly to soil protection. Complementing this explanation, the coarsening of topsoils due to persistent soil erosion would lead to a reduction in the supply of fine-sized grains to the streams.

For most of the landscape disturbances described in this work, we cannot unambiguously quantify the overall effect on sediment fluxes. Yet, sediment fluxes are more efficient during periods of high flow, which correspond to wetter conditions. Consequently, the sediment stored in the valleys, which is highly rich in nutrients and carbon, can be resuspended during higher-discharge events, causing temporarily delayed off-site problems for several decades to come.

In the studied landscape (the CCR in the Maule Region) the surface lowering in the last 3 decades is similar to or higher than the long-term benchmark. Thus, we argue it may be considered high for this specific system. This conclusion and the documented effects of local tree plantations on soil organic carbon, soil density, and soil biodiversity are clear indicators of a degrading landscape.