Climate and Land Use Change Effects on Sediment Production in a Dry Tropical Forest Catchment

Abstract

Understanding the natural and anthropogenic drivers that influence erosion and sediment transport is a key prerequisite for adequate management of river basins, where, especially in tropical catchments, there are few direct measurements or modeling studies. Therefore, this study analyzed the effect of human-induced land-use changes and natural ENSO (El Niño-Southern Oscillation) related changes in rainfall patterns on soil erosion and catchment-scale sediment dynamics with the SEDD (Sediment Delivery Distributed) model. In the 393 km² Tonusco river basin, representative of tropical, mountainous conditions, daily rainfall data were used to quantify changes in rainfall erosivity and satellite images for the evaluation of cover factor changes between 1977 and 2015. The final model combined soil loss, calculated by RUSLE, with a sediment routing-based delivery ratio, that was calibrated and validated with data from the sediment load recorded at the basin outlet. The results detected a great reduction of the vegetation cover in the catchment during the last decade of from 79.5 to 29.5%, and the influence of important runoff and erosion events linked to La Niña episodes. Soil erosion rates were locally very high, of over 120 Mg ha⁻¹yr⁻¹, and sediment yields were estimated at the range of 6.17–8.23 Mg ha⁻¹yr⁻¹.

1. Introduction

Soil erosion is one of the greatest threats to the future of our agriculture and environment, but also to our society. As Sposito [1] indicated, soil is a natural resource that cannot be taken for granted. Several recent studies warn about the risks of erosion for global soil resources, but also water resources and ecosystems [2,3,4,5]. To tackle on and off-site damage due to erosion, adequate land and watershed management is necessary [6]. This requires information on the spatial and temporal variability of soil erosion and sediment yield, in order to tailor soil conservation measures and downstream sediment management.

For watershed management assessment, runoff and sediment yield can be estimated with hydrologic and erosion parametric models, either based on physical principles, as reviewed by Pandey et al. [7], or on simple schemes, such as the revised universal soil loss equation (RUSLE) or any of its modifications, as used by Zi et al. [8]. The advantage of the latter is the possibility of its application to areas with reduced information. Recent studies have extended the potential of the RUSLE model, representing diffuse water erosion, by combining it with other models. For example, Gudino-Elizondo et al. [9] applied the AGNPS (Agricultural Non-Point Source) model for gully erosion combined with the RUSLE model for diffuse erosion. An interesting approach was designed by Chiu et al. [10] using RUSLE and a previous version including runoff characteristics (MUSLE) coupled with a landslide sediment production function linked to a sediment transport one. Other studies [11] used RUSLE together with a distributed sediment delivery model (SEDD) of Ferro and Porto [12], or SEDD combined with some locally adapted models such as the Chinese soil loss equation, (CSLE) by Wu et al. [13]. A review by Alewell et al. [14] revealed that many or most of the erosion and sediment delivery models incorporate USLE- type modeling.

The Magdalena River is the main drainageway of the Northern bifurcation of the Andean cordillera in Colombia. The Tonusco River is a tributary of the Magdalena River left margin affluent, the Cauca River. The Magdalena River discharged an average of 228 km³yr⁻¹ in the period 1975–1995, with a mass of sediment of 144 Tg yr⁻¹ [15] which implied a sediment yield of 0.560 kg m⁻²yr⁻¹ in the whole watershed. A posterior estimate [16], for the period 2000–2010, increased the sediment yield to 0.710 kg m⁻²yr⁻¹, although the runoff yield for the period 1958–2000, was not much different, i.e., 205 km³yr⁻¹, according to Restrepo et al. [17]. That sediment yield rate is very high, even for the South American continent, and the perspectives of possible climate fluctuations are of concern because of their potential effect on runoff and sediment yield. This climate effect is well known, and several studies have documented how El Niño and La Niña events have caused extreme weather and an increase in runoff and sediment yield all along the Andes, although their impact varies significantly from a spatial perspective [18]. Regional climatic conditions of Central and South America, East Asia, and Oceania are influenced by the trade wind circulation along the equator area of the Pacific Ocean. Trade winds transport moist air from the ocean, regulating the rainfall regime in these zones. The extreme change in sea surface temperature alters the normal wind circulation, causing weather anomalies called ENSO phenomena. These events have two opposing phases known as El Niña and La Niña. During El Niño, the average temperature of the surface sea increases causing trade winds circulation to decrease. La Niña has the opposite effect of El Niño because the sea temperature decreases, and trade winds are even stronger than usual [19]. In the Andes region of Colombia, these events can change trade winds speed affecting the moisture transport dynamics in the Pacific Jet stream [20]. These effects can either reduce or increase the frequency and volume of rainfalls [21]. During El Niño and La Niña events, the annual rainfall level decreases between 200 mm and 400 mm and increases between 200 mm and 1000 mm compared to the average, respectively. In northwestern Colombia, there are few studies on it but Restrepo and Kjerfve [15] found that extreme El Niño events caused lower discharge and sediment yield in the Magdalena River in Colombia; whereas in western Ecuador and northwestern Peru, such events produce an increase in discharge and suspended sediment yield of 5.4–11 times [22,23]. Restrepo et al. [17] detected the runoff enhancement effects of La Niña years on the continental rivers. The reasons for such high sediment yield rates were, basically, the intense rainfall events, and the steep slopes of most of the contributing watersheds.

An additional factor contributing to the high soil losses and posterior sediment yield is the diminishing soil cover in northern Colombia. The natural vegetation of the region, in particular in the mountain watersheds, is the tropical forest, or, more precisely, the seasonally dry tropical forest [24], consisting of many plant species which shed their leaves during part of the year [25,26]. Therefore, there is an increased erosion risk during the period of the year when the vegetation cover is sparser. The progressive human occupation of the watershed aggravates this erosion risk, given that much of the natural vegetation is being cleared. Etter et al. [27] indicated the accelerated land-use change in Colombia in the last forty years and quantified the forest area reduction as being 60% in the Andean Region. Restrepo and Escobar [28], in an analysis of sediment load trends in 21 tributaries of the Magdalena River basin during the last 30 years, reported an important increase in sediment yield between 1990–2000 and 2005–2010, that they were able to link to rapid deforestation phases. These authors concluded that the increasing trends were mainly linked to human-induced forest clearance and should not be attributed to climate change. In a follow-up study, they quantified the deforestation in the Magdalena basin from satellite imagery and found a 241% increase in deforestation between 1990 and 2010, which led to a 33% increase in erosion rates.

The Magdalena River is the principal river of Colombia and its extremely high contribution of sediments to the Atlantic Ocean is a problem that needs to be tackled. This requires a good understanding of the spatial and temporal variations in erosion and sediment yield, and their response to land use and climate change. Currently, no robust or complete understanding exists for northwestern Colombia at the scale of small river catchments that could be used for detailed environmental planning.

The general objective of this study was to analyze the soil erosion and sediment yield dynamics in a representative watershed of northwestern Colombia, the Tonusco basin, and to determine the effect of human-induced land-use changes and ENSO-related climatic changes. This will eventually enable resource managers to identify priority pollution sources and to plan appropriate conservation actions. The framework set out in this study can then be applied to other, similar catchments. The specific objectives are:

(i)

To analyze the changes in rainfall erosivity, vegetation cover, and soil loss rates as a function of land-use, and ENSO-induced seasonal variations.

(ii)

To calibrate the spatially distributed SEDD model, by contrasting it with sediment yield measured at the catchment outlet.

(iii)

To apply the SEDD model to estimate the spatial and temporal distribution of sediment yield.

2. Materials and Methods

2.1. Study Area

The Tonusco River basin’s catchment area is 393.4 km² and it is located within the western Andes Mountains in Antioquia-Colombia (Figure 1). This river starts at the “Tonusco Arriba” Mountains, (6°34′21″ N, 76°0′40″ W) and flows into the Cauca river (6°30′49″ N, 75°49′11″ W). It has a high altitudinal gradient and descends from 3425 to 450 m amsl with an average slope of 39.5% The watershed lies on the accretion of oceanic rocks of the Cretaceous age [29], with some outcrops of the Paleozoic-Cretaceous basement [30]. The soils are very recent, of the Andept and Orthent suborder [31], in the upper basin, and formed on sediments, suborder Ochrept, in the lower part [32]. Currently, the land use in the basin corresponds mainly to forests (54.5%), followed by pastures and grasses (28.1%), crops (14.2%), bare soils (2.9%), and urban area (0.3%), with important land-use changes in recent decades. Restrepo et al. [16] have reported a significant increase in deforestation in the greater Magdalena basin, and the same is true for this sub-basin, according to the Alexander von Humboldt Institute for research [33], although this will be analyzed in detail below. The predominant climate can be classified as tropical savannah, Aw type in the Köppen-Geiger scheme [34]. The annual rainfall varies from 1000 to 1900 mm, with two rainy periods in April–May and October–November, and a marked dry interval from December to February. The mean annual reference evapotranspiration is 819 mm. During the drier interval, the scant rainfall and high evapotranspiration rates reduce the vegetation cover, especially in the lower part of the watershed. The vegetation is typical of the seasonally dry tropical forest, with species like the Bursera simaruba and Cnidosculus urens.

El Niño” or “La Niña” Southern Oscillation (ENSO) phenomenon exerts strong impacts on the interannual Tonusco basin weather variability. Restrepo and Kjerfve [15] computed a reduced mean discharge rate of 76.5% of the average 1975–1995 period during El Niño years and 121% of it as a mean value during La Niña years at the Magdalena River. In the same Magdalena River Basin, Poveda et al. [21] detected Normalized Difference Vegetation Index (NDVI) reductions of up to 30% of the average values during El Niño events, indicating vegetation stress and a sparser vegetation cover due to the absence of rain.

2.2. Rainfall, Flow Rate, and Sediment Concentration Measurements

Rainfall data were taken from six weather stations close to the watershed (

Table 1. Location details and daily records available in the period 1977–2015 of weather (W) and stream gauge (S) stations.

Stations	Type	Coordinates (UTM)		Records
		X	Y
Cotové	W	806,138.1	1,215,249.5	12,691
Guasabra		792,864.5	1,209,840.7	13,171
Giraldo		792,898.0	1,230,809.9	13,473
Buriticá		797,313.3	1,235,292.9	13,240
Abriaquí		779,613.9	1,225,348.3	13,366
Caicedo		789,472.9	1,201,028.1	13,250
La Galera	S	802,304.7	1,216,969.4	11,286

). Water discharge rates and sediment concentrations were measured at the gauging station of La Galera, located in the lower part of the Tonusco River, close to its confluence with the Cauca River (Table 1). The five rain gauges and the stream gauge stations are run by the Meteorological and Environmental Institute of Colombia (IDEAM) The location of all the gauges is shown in Figure 1.

The data series were recorded between 1977 and 2015, with the first year being 1977 because it was the first one with complete records in all the stations.

The rainfall, discharges, and sediment load raw data series were incomplete, or abnormal in some periods, which should be corrected. The data correction method is based on three criteria [35]:

(a)

Accept the years with more than 95% of daily records (>347 data).

(b)

Reject the years with less than 70% of daily records (<256 data).

(c)

Accept the years with a percentage of records between 70 and 90%, if their average value falls into the confidence interval (>95%) established using the standard deviation of these.

2.3. Potential Soil Erosion Rate Modelling

The revised universal soil loss equation or RUSLE [36], was selected for the evaluation of the potential erosion at the watershed scale, at a spatial resolution of 25 m. Annual soil loss per unit area and year, A, is expressed by the product of six factors representing rainfall-runoff erosivity factor, R; soil erodibility, K; slope length and steepness, L and S; cover-management factor, C; and the adoption of any conservation practice, P. Since the value of the latter factor is usually taken as being 1 in the absence of conservation practices, which is the case here, it will be ignored in the following discussion.

$$A=R·K·LS·C·P$$

(1)

All the factors were taken as being constant over time, except the rainfall erosivity factor R and the cover-management factor C, which permitted the evaluation of the effect of both the climate variations on R and C, and the human-induced land-use changes on C.

The rainfall-runoff erosivity factor, R, is computed by integrating the unit-specific energy of the different pulses of a rainfall event with their respective duration and multiplying it later by the maximum intensity in a 30-min period, which is known as EI₃₀ (MJ-mm-ha⁻¹h⁻¹yr⁻¹), for each rainstorm. Summing EI₃₀ for all events during a particular year yields the annual R (MJ-mm-ha⁻¹h⁻¹). In the absence of high-resolution rainfall records, as was the case here, an alternative procedure can be used, in which the R-value is estimated by a potential relationship between the EI₃₀, and the daily rainfall depth measured, P_d (mm). In this study, the equation proposed by Hoyos et al. [37] for a neighboring catchment was applied:

$$E{I}_{30}=0.422{P_d}^{1.769}$$

(2)

The annual rainfall erosivity factor R (MJ-mm-ha⁻¹yr⁻¹) is equal to the cumulative annual sum of EI₃₀ values.

$$R=\frac{{{\displaystyle \sum }}_{n=1}^N{\left(E{I}_{30}\right)}_n}{Y}$$

(3)

where N is the number of rainfall events in Y years.

The soil erodibility factor, K, was estimated from a multiple regression equation with the concentration of several soil particle sizes and organic matter. Since no detailed information on the particle size distribution of the watershed was available, the proposed values for different soil textural classes provided by Mitchell and Bubenzer [38], and the Colombian soil map system data generated by the Agustin Codazzi Geographic Institute of Colombia [39] were used. The erodibility values selected were similar to those obtained by Bonilla and Johnson [40] in Central Chile and vary between 0.0135 to 0.0544 Mg-ha-MJ⁻¹mm⁻¹.

The topographic factor LS was estimated following the method indicated by Fernández et al. [41].

The soil cover management factor, C, was evaluated by resorting to its relationship with the NDVI, measured from Landsat satellite imagery, as outlined by Van Der Knijff et al. [42] and Carvalho et al. [43]. This index is the ratio between the difference and the sum of the spectral reflectance measurements acquired in near-infrared and red (visible) regions, respectively. The satellite images were taken by Landsat 4–5, Landsat 7, and Landsat 8, corresponding to 1986, 1991, 1997, 2001, 2009, and 2015, were downloaded from the EarthExplorer web portal [44] and processed with the software ENVI [45]. The results of these specific years were interpolated in order to obtain continuous C values during the study period 1977 to 2015. The interpolation was made using C factor values estimated from NDVI for 1986, 1991, 1997, 2001, 2009, and 2015, assuming linear relationships between every year.

The vegetation cover and the NDVI can be expected to vary as a function of (i) overall land-use changes; and (ii) inter- and intra-annual climate variability resulting in soil moisture differences. Therefore, the long-term evolution of NDVI was always determined at the transition between the wet and dry seasons.

Next, a relationship between C and NDVI was established in order to assess the temporal evolution of C-values. Therefore, first, 9 different land use classes were identified, based on their NDVI values, ranging from bare soil to forest, according to Van der Knijff [42] and Dalezios et al. [46] articles. Second, a different C-value was assigned to each one of the land cover classes, based on the suggestions of Benavidez et al. [47] and RUSLE literature [36].

Table 2. Land cover classifications with corresponding NDVI and C factor values, used in this study.

Classification	NDVI	C Factor
Water	<0	-
Bare Soil	0.00–0.10	0.25–1.0
Light vegetation	0.11–0.20	0.20–0.24
Moderate vegetation	0.21–0.30	0.08–0.19
Medium vegetation	0.31–0.40	0.030–0.07
High vegetation	0.41–0.50	0.026–0.029
Dense vegetation	0.51–0.60	0.017–0.025
Forest	0.61–0.75	0.0042–0.016
Dense Forest	0.76–1.0	<0.0041

shows the NDVI and C factor interval classifications that were used. Figure 2 shows the NDVI-C factor relation obtained, and in the rest of the analysis the following equation was employed to relate remote-sensing-derived NDVI values to C-values directly (R² = 0.99):

$$C=0.997e^{-6.27NDVI}$$

(4)

The sediment production in landslide and gully erosion is a non-relevant process in SEDD model analyses because the latter evaluates the water erosion process. Sediment contribution from a landslide is a complex process that needs to be investigated in specific studies.

2.4. Distributed Sediment Delivery Modeling

Once detached, the fragments of soil aggregates travel downslope carried by water until they settle at the bottom of hillslopes or within channels. Consequently, only a fraction of the original soil loss leaves the watershed. The concept of sediment delivery ratio (SDR), or the ratio of the mass of sediment leaving the watershed over the averaged in-situ soil loss, has been widely employed to express sediment connectivity within watersheds. In the early studies of sediment transport, SDR was estimated by some empirical potential relationships with the watershed area [48]. Williams [49] suggested a negative exponential of the product of the sediments travel time along the watershed and the square root of its effective diameter as a representation of the probability of its exit at the outlet. More recently, Ferro and Porto [12] formulated a simple model based on the Williams suggestion, but removing the sediment effective diameter from the exponent and estimating the travel time by adopting the inverse of the square root of the slope as a proxy of the water velocity in the watershed. The sediment delivery distributed model (SEDD), contains a parameter β to match the missing factors in the sediment delivery ratio from the i-th geomorphological unit to the lower end of the watershed.

$$SD{R}_i=e^{-\beta t_{p,i}}$$

(5)

The particles’ travel time, t_p,i could be calculated using the length (l_p,i) and slope (s_p,i) of each morphologic unit i of the basin.

Although the SEDD model was originally developed in a small agricultural catchment, other studies have successfully applied it to larger basins. For example, Fernández et al. [41] and Boakye et al. [50] estimated sediment yield using the SEDD model in Idaho basins (USA) and the Pra river (Guinea), with areas of 114 km² and 23,200 km², respectively. Di Stefano et al. [51], tested it in 7 Sicilian basins, up to 213 km². Diwediga et al. [52] used SEDD to model erosion dynamics in the 1485 km² large Mo river in the central part of Togo. Similarly, Liu et al. [53] applied it to study the sediment contribution process in a basin with an area of 130 km² located in Germany.

From the estimated potential soil loss of the basin, A_total (Equation (1)), and the measured sediment yield, Y_total, the sediment delivery ratio, SDR, was calculated as

$$SDR=Y_{total}/A_{total}$$

(6)

This enabled the calibration of the β coefficient by using the measured sediment load and estimated soil loss rate, for a part of the measurement period, leaving the other part for model validation. The Nash and Sutcliffe [54] efficiency coefficient was used to evaluate the model’s performance.

Sub-watersheds do not have sediment load measurements, so that the β coefficients calculated for the Tonusco Basin were employed to establish its SDRs, by using the length-l_p,i and slope-s_p,i of each subwatershed, following Equation (5). l_p,i and s_p,i was established for each pixel of the sub-watershed with digital cartography of the zone contained in the Geographic Information System. Originally, subwatershed sediment loads were established by multiplying SDRs by the soil loss potentials A_p of each one.

3. Results and Discussion

3.1. Observed Runoff and Sediment Loads

The annual rainfall, discharge, and sediment load between 1977 and 2015, during La Niña, El Niño (ENSO years), and Normal periods are shown in Figure 3. ENSO years were detected quantifying Oceanic indexes based on the surface temperature of the Pacific Ocean southern region, such as Oceanic Niño Index (ONI) and Cold Tongue Index (Nct) and to confirm weather features related to these periods [55]. The Niño and La Niña periods are based on a threshold of +/−0.5 values of these indexes [19], and color codes of light blue, red, and dark blue were used in Figure 3 to indicate whether the precipitation (P), discharge (Q) and sediment load (W) of a particular year belong to a normal, El Niño and la Niña years.

The variation in rainfall reflects the climate variability expected as a result of ENSO phenomena very well Maximum values of rainfall were recorded in 1999, 2010, and 2011. These years are all La Niña years, marked in dark blue. It can be seen how these La Niña years are generally wetter compared to normal years, marked in light blue, while El Niño years, marked in red, are drier compared to normal years. The trends in discharge and sediment load are more complex. In the La Niña years 1999, 2010, and 2011, the maximum flow rates measured correspond well to the higher rainfall in those years. Also, the sediment load then was greater compared to the years before and after. However, during 1983, 1984, and 1987, high discharges were measured, although these were considered El Niño years. The same occurs with the sediment load data, which were unusually high in the entire period from 1983–1987. It is difficult to find a good explanation for these data, which could suggest the mobilization of sediments that had accumulated in the river valley during previous dry years. Morera et al. [18], proposed a similar mechanism in the western Peruvian Andes, where they observed an accumulation of sediment during dry, normal years and a rapid mobilization during El Niño events, which, in their study area, correspond to the wetter years. However, there were high loads during four consecutive years, which seems to indicate the presence of an additional source of sediments, which could come from a large-scale disturbance in the catchment, such as road construction, or a large landslide, or a combination of both. In the middle and upper part of the Tonusco Basin, there is, in fact, evidence of construction activity in that period. In the early 1980s, the “Carretera al Mar” highway was constructed, which is likely to have also caused additional landslides and increased soil erosion.

As mentioned earlier, the Tonusco River is a tributary to the Cauca-Magdalena Basin. Restrepo and Kjerfve [15] and Restrepo [56] reported maxima and minima values in water discharge and sediment load during La Niña and El Niño years, respectively. Here, the trend is less clear over the entire study period, with the exception of the period 1983–1987. During the most recent one, from 1990 to 2015, discharge and sediment loads have been increasing and decreasing in the same periods, when the La Niña and Niño phenomena, respectively, occur.

3.2. Modeled Potential Soil Erosion Rates

With respect to the influence of climate on potential soil erosion rates, the variation in the calculated rainfall erosivity factor (R) is shown in Figure 4. The average R-value over the entire study period was 4977 MJ-mm/ha-yr, with minimum and maximum values ranging between 3303 and 14,431 MJ-mm/ha-yr. The maximum erosivity values occurred during 1999, 2010, and 2011, and the minimum ones during 1977, 1979, 1987, 1992, 1997, 2002, 2014, and 2015.

The highest R-values are concurrent with the La Niña event years, while the lowest values correspond to El Niño” years. The average R-value during the La Niña years is 52.4% higher compared to the average of normal years, and during the extreme La Niña year of 2011, the difference amounted to 190% compared to normal years. The higher R-value obtained in 2011 is related to 3200 mm of annual rainfall, which is 104% higher than the average rainfall in the basin area.

The vegetation cover changes were estimated through the evolution of the NDVI values which were used as a proxy of the C-factor. The results show a twofold change in vegetation patterns. The first change is the periodic variation due to the succession of dry and rainy seasons. Figure 5A,B shows the contrast of NDVI, expressing the vegetationgreenness, between wet and dry seasons. It can be seen how during the wet season, almost the whole catchment has a high vegetation cover. In contrast, during the dry season, many areas with medium-low NDVI values occur in yellow (0.40 or lower), especially around the eastern part of the catchment. According to an ANOVA test (Figure 5C), the differences in NDVI values between dry and wet periods were not statistically significant (p = 0.00).

The second type of vegetation cover change is due to land-use variation over successive years [57,58]. Figure 6 shows the vegetation cover evolution in the study period. Significant changes can be observed, with an important increase in low and medium NDVI values (red to yellow colors, 0–0.40), corresponding to areas with bare soil or light to medium vegetation. These areas increased from less than 1.52% of the total area in 1997 to almost 11.7% in 2015. This land-use change has been a very fast one, as can be seen from the comparison of 2009 with 2015. The high NDVI values (0.60–1.00), which correspond to the dense forest, dropped from 79.5 to 54.1% in 1997 and 2015, respectively. In summary, there is an important land use effect on the temporal variation of C values, clearing dense forest vegetation and making the catchment more susceptible to soil erosion, on top of a wet-dry seasonal effect, that can be aggravated by ENSO events. Figure 6 shows the resulting evolution of C values between 1977 and 2015. The vegetation cover increases between the satellite images of 1986 and 1997 led to a clear drop in the C factor. The vegetation cover increases could be related to the decreasing of urban and agricultural activities in the zone, caused by the reduction in the rural population, which, in turn, is a consequence of the forced internal displacement of people occurring during the Colombian armed conflict from 1990 to 2000 [59]. After this period, Figure 7 shows how the decrease in vegetation cover led to an increase in C-values that was intensified up until 2015.

Finally, the spatial distribution of the estimated overall potential soil loss for the Tonusco basin is depicted in Figure 8, compared between 1997 and 2015. In such a period the area with higher soil loss potential has increased caused, probably, by the reduction of forest cover. There are important spatial differences that occurred within the catchment in the two periods. The largest soil losses, well over 120 Mg ha⁻¹yr⁻¹, were found on the steep slopes near the watercourses, in the northern and central region of the basin. This area corresponds to areas heavily deforested by anthropic intervention related to the urbanization activities executed near Giraldo municipality and the “Carretera al Mar” highway. These areas have low NDVI and high C values. Soil losses in the upper and lower parts are very small, with a gradual variation between the two extreme values. In the upper part of the catchment, the forest cover is relatively intact, as can be seen in Figure 8B. In the lower part, near the outlet of the catchment, there is a sparse vegetation cover, but low slopes. The estimated mean annual soil loss in the basin was 20.2 Mg ha⁻¹yr⁻¹.

The estimated average yearly soil loss value changed significantly over time, as indicated in Figure 9. Important peaks were found in 2010 and 2011, also La Niña years. However, the latter is not always related to higher soil erosion rates. The years 1995, 1996, 1998, 1999, and 2000 were also La Niña years, but, despite higher rainfall erosivities (Figure 4), the erosion rates were lower than in more recent, “normal”, and even El Niño years, in the period 2008–2015. This can be attributed to the lower C factors in most of the La Niña years related to increases in plant cover [21]. This results in a great variability of erosion rates in El Niño years, so that the differences from normal, and even La Niña, years were not statistically significant (p = 0.06) (Figure 9B).

More importantly, the pattern of simulated potential soil loss strongly resembles the temporal pattern of the measured sediment load at La Galera, with the exception of the unusually high sediment loads in the eighties. As was explained earlier, this is due to additional disturbances in the catchment related to road construction. This suggests that there is a relatively direct connection between sediment production and sediment export in this catchment, and that there is a limited accumulation of sediment within the river valley during dry months or years, [60]. This will be explored further in the next paragraph by fitting a distributed sediment delivery model.

3.3. Modeling Sediment Delivery

Figure 10 compares the annual sediment delivery and β coefficient values, estimated with Equation (5). The results show an important year-to-year variation, with SDR and β values ranging from 0.05 to 1.00 and 0.00 to 0.161.

As mentioned earlier, the sediment loads estimated between 1984 and 1986 were discarded because they were considered to be anomalous data. During this period, sediment loads and SDR values were maximum, which does not follow the trend of flow rate and rainfall values because those data were probably influenced by road construction.

To explain this temporal variation, we established a relationship between runoff coefficient (Cr) and the β parameter (Figure 11), according to the Ferro and Porto method [12], which we calibrated by using the results for 1977, 1980, 1983, 1992, 1995, 2001, 2004, 2006, 2008 and 2012. Cr and β coefficients result for 1978, 1979, 1987, 1997, 1999, 2002, 2005, 2007, 2009, 2010 and 2015 were used to validate it.

The obtained relationship between the parameter β, and the runoff coefficient, Cr is:

$$\beta =-0.0839 \times Ln\left(Cr\right)-0.0736$$

(7)

Figure 11 shows that the proposed Equation (7) is very close to the calibration and validation data. The results of this study are, therefore, similar to those established by Ferro and Porto [12]. The comparison between measured and calculated annual sediment load in the Tonusco basin, shown in Figure 12, indicates a good agreement, with an r²-value greater than 0.60, and a value of the Nash and Sutcliffe efficiency index of e_NS= 0.572, which is almost acceptable for hydrological research [61].

Figure 12 shows that the relationships have a good fit between 0 and 500 Gg yr⁻¹. In the same way, discarded data are higher than 500 Gg yr⁻¹, therefore the relationship can be applied in this interval. The Nash-Sutcliffe coefficient was thus recalculated discarding sediment loads of over 500, obtaining a value of 0.794. This good relationship between the runoff coefficient and β parameter, used to calculate the sediment delivery ratio, indicates that the sediment connectivity within the catchment varies as a function of runoff connectivity. It is well known that the fraction of the network connected to its uplands controls runoff and sediment connectivity [62], although empirical evidence such as that presented here is scarce and valuable at the catchment scale. It allows for taking into account the climate influence of ENSO [63] related years on the sediment transfer in this catchment, and this approach could be extrapolated to other, similar ones.

3.4. Annual Sediment Delivery in Tonusco River

Initially, the annual sediment deliveries rates-SY were estimated for the Tonusco River Basin and compared to measured deliveries from the Cauca and Magdalena rivers, obtained according to Niño [64] and Restrepo [56] research (Figure 13).

For many years, the annual SY in the Tonusco River Basin has been higher than those obtained in the Cauca and Magdalena rivers. The mean SY in the Tonusco river was 1.13 Gg km⁻² yr⁻¹, whereas in the Cauca and Magdalena rivers they were 0.834 Gg km⁻² yr⁻¹ and 0.710 Gg km⁻² yr⁻¹, respectively. Also,

Table 3. Average sediment yield values estimated in tributaries catchments of Cauca and Magdalena rivers.

River	Basin	Area (km2)	SY (Gg km−2yr−1)
Samaná [56]	Magdalena	1710	0.525
Nare [56]	Magdalena	5564	0.404
Cocorná [56]	Nare-Magdalena	790	0.620
Gualí [56]	Magdalena	458	0.568
Coello [56]	Magdalena	1041	1.11
Arma [64]	Cauca	1031	0.430
San Juan [64]	Cauca	1108	0.554
Aurra [64]	Cauca	244	0.072

shows the mean SY values in different basin tributaries of the Magdalena River.

According to Table 3, the mean sediment delivery of the Tonusco river is comparatively much higher than the mean corresponding to other rivers located in the region.

3.5. Sediment Delivery in Sub-Watershed Tributaries of Tonusco River

Sub-basin sediment loads were established using SDRs and soil loss potentials A_p, respectively.

The estimated sediment load for one sub-watershed, La Pená, has been compared to the measured load of the whole watershed in Figure 14. Although there is a certain proportionality between the sub-watershed and the whole watershed, there were appreciable differences in some years, such as in the eighties, due to the influence of the specific erosivity and surface cover factors at that time.

This model has allowed the localization of the subwatersheds and areas with the highest sediment production in the Tonusco Basin. Also, the model provides an information base for understanding sediment connectivity in the catchment, according to studies by Bracken et al. [65] and Fryirs [66].

These results are a first step to the understanding of the erosion processes and connectivity in this watershed, which is one very small part of the whole Magdalena River basin, but, at the same time, highly representative of the water and sediment transport in it. Kettner et al. [67] have explained the main traits of the Magdalena River, but there are many factors whose influence needs a closer inspection, like that of humans in the deforestation and forest recovery processes and their implication [27,68].

4. Conclusions

The Tonusco River watershed, in the most erosive part of the Cauca and Magdalena Rivers, draining the steep slopes of the Western and Central Cordillera of the Andes, is an important contributor of sediments. In addition to the steepness of the land and the intense rainfall events, the vegetation cover of the watershed is not permanent, since the seasonally dry tropical forest reduces its protection of the soil during the dry season.

The rainfall erosivity was found to vary during the ENSO El Niño-La Niña events, with large estimated soil losses confirmed by the sediment loads measured at the outlet of the Tonusco basin, from 20 to 1010 Gg yr⁻¹. The high occupation rate of the watershed reduced the average surface cover, increasing the risk of soil loss by water erosion by a factor from 2 to 4 between 2002 and 2015.

The analysis of the evolution of NDVI and the C-factor shows that the forest and dense vegetation cover have been diminishing during the last 20 years, from 79.5 to 54.1% of the basin area, while grasses, crops, and bare soil cover areas have been increasing. The possible cause of that is the practice of agricultural activities in the zone. These changes in the cover have increased soil loss and sediment load in the Tonusco basin in the last two decades.

The actual maximum estimated soil loss rate in some parts of the Tonusco watershed is as high as 120 Mg ha⁻¹yr⁻¹and the mean is 20.2 Mg ha⁻¹yr⁻¹. Fortunately, part of this material does not reach the Cauca River but the specific averaged sediment delivery ratio SDR of the watershed is still in the range of 1.94–49.5%, with the exception of the period from 1983–1987, (these data were discarded).

The results obtained using the SEDD model show a good agreement between the sediment transport and climate and geomorphologic conditions of a Tropical Basin. This can be seen in sediment dynamics during the wet and dry periods and comparing it in mountainous and flat areas of the basin. The inter-annual rainfall changes related to ENSO phenomena events were reflected in the sediment load and delivery results because most of the sediment production and transport values show their maxima and minima during La Niña and El Niño periods, increasing to 172% and dropping to 64.5% with respect to the average values, respectively.

Most of the sediment load contribution to the Tonusco River is generated in the northwestern and central part of the watershed, in small sub-catchments located near the main river channel, with values between 25.0 and 38.0 Gg ha⁻¹yr⁻¹. Therefore, erosion control techniques should be applied to soil conservation in these areas.