Since the GCPs placed at the bottom of the gullies caused a degradation in the quality of their point network, it is necessary to ensure that the error in the SfM-MVS-4D results is acceptable.

<table><caption>**Table 3.** Errors in X, Y and Z of the three sites.</caption><thead><tr><th rowspan="2">Sites</th><th colspan="4">B.G-1</th><th colspan="3">B.G-2</th><th colspan="3">B.G-3</th></tr><tr><th>M-I</th><th>M-II</th><th>M-III</th><th>M-I</th><th>M-II</th><th>M-III</th><th>M-I</th><th>M-II</th><th>M-III</th></tr></thead><tbody><tr><td>Missions</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>X (cm)</td><td>1.93</td><td>2.96</td><td>2.17</td><td>3.26</td><td>4.94</td><td>4.47</td><td>1.46</td><td>1.01</td><td>1.18</td></tr><tr><td>Y (cm)</td><td>2.19</td><td>3.13</td><td>2.37</td><td>3.36</td><td>4.53</td><td>4.46</td><td>1.27</td><td>1.60</td><td>1.82</td></tr><tr><td>Z (cm)</td><td>5.32</td><td>5.13</td><td>4.87</td><td>5.3</td><td>5.79</td><td>5.67</td><td>1.36</td><td>1.44</td><td>1.98</td></tr></tbody></table>

Accuracy analysis of the control points showed a strong degradation in the accuracy of the control points depending on the depths of the gullies. This is due to the poor reception of the GNSS signal at the bottom of the gullies, which led us not to use GCPs for point cloud validation. In this study we chose to evaluate the quality of the data by three methods.

## The Uncertainty of Coregistration

We presented in the methodology the procedure used for the coregistration of dense point clouds using the ICP algorithm which allows to considerably increase the accuracy of georeferencing of SfM-MVS-4D data in an external georeferencing system. One of the main sources of error in detecting changes between two point clouds is uncertainty in coregistration. Even with low coregistration uncertainties between each pair of dense point clouds, the risk of the uncertainty amplifying increases with the temporal resolution of the point cloud collection. To reduce the risk, we proceeded by coregistration in pairs of point clouds. Here we seek to assess the quality of dense point clouds from a dual coregistration (GCP and ICP). Data uncertainty is modeled using three spatially variable models: FIS and cross-validation.

## Cross Validation

The training and validation data are constructed by resampling without replacement, which allows us to calculate the uncertainty on the different cloud points. This method is more robust, more stable and more thorough than Bootstrapping (resampling with discount). To our knowledge, this method has never been used to quantify the uncertainty of point clouds. The quantification of the uncertainty of the point clouds with the cross validation tools was successfully carried out for the three bank gullies. Analysis of the uncertainty results obtained by the cross-validation method shows that the absolute difference can go up to 0.5 m, but the frequency distribution shows an dissymmetry with an absolute average for the three gullies not exceeding 0.08 m and a standard deviation of 0.09. The repetitive analysis of uncertainty on the three gullies has brought to light aspects of the spatial structure of uncertainty. The visual interpretation of the results (Figure 9) aimed that the high uncertainty values tend to be concentrated on the talwegs (main and secondary), and on some parts of the slopes. On the other hand, the low values are generally found in the flat areas and at the tops of the gullies. The uncertainty values are plotted against the slope and the roughness. On the Figure 9, it appears that all the high values of uncertainties are strongly correlated with the steep slopes and the zones of strong roughness what met in relation to the spatial distribution of the uncertainties with the topographic complexity and the roughness of ground. It must be concluded that roughness and slope are the main sources of errors. In several studies carried out, in order to determine the accuracy of the SfM-MVS data, the DGPS surveys of the GCPs are compared with the final product (DEM-SFM-MVS). This confusion is deliberately maintained, as if the uncertainty of the SfM-MVS data was only related to the quality of the DGPS surveys, and that the geometry of the photos taken could not improve the accuracy of the data. In our case of river gullies, the possibility of using a high number of photos, of high overlap, of including oblique photos as well as the double coregistration of point clouds with the ICP algorithm allowed us to compensate and to improve the georeferencing weaknesses by the GCPs.
<figure><img src="image_85.png" /><figcaption>**Figure 9.** Example of point cloud uncertainty estimated by cross-validation on the gully B.G-1.</figcaption></figure>

FIS

The results of the detection of geomorphological changes by the DoD technique are influenced by the type of error model used (spatially uniform or variable). Spatially uniform models can lead to overestimation or underestimates of topographic changes.

The FIS model, on the other hand, allows the changes to be reflected thanks to its better resolution, which leads to a better interpretation of the results. This is due to its ability to exclude noise, which is more important than for other models, also it leaves a morphological signature. Finally the results are transcribed in the form of a histogram, often bimodal. This tool facilitates the quantification of real changes in topographic models.

For the current study, on the three bank gullies, we detected geomorphological changes in the gullies with spatially variable FIS models. As the resolution of the data influences the detection of changes, we chose to use this model with a resolution of 5 cm. This technique allowed us to detect all possible changes, even the smallest ones, for our digital terrain models.

## 3.2. *Detection of Geomorphological Changes from UAV Data*

All change detection measurements over a period of time represent the result of the aggregation of all geomorphological processes. Positive values are assigned to all sedimentation processes, while negative values are appropriate for all erosion processes. According to the change map established by the DoD technique and the FIS propagated error, on both sides of the B.G-1 gully, erosion occurs with a significant difference between the two. This difference is mainly due to the exposure of the two slopes. The results of the two change detection techniques are similar but differ slightly in the quantification of sediment budgets. This small difference corresponds to the error model used by each technique. The two change detections are implemented by spatially variable error models. M3C2 uses roughness and coregistration error as error model variables. The FIS combines all the potential sources of errors (slope, point density, interpolation error, DTM-SfM-MVS error) to obtain a single error model. However, it is possible that the M3C2 technique underestimates the SfM-MVS error model.
### 3.2.1. DOD in Site 1

In general, monitoring of the sediment balance on the bank gully (B.G-1) allows us to notice significant changes during the year 2018. The bank gully (B.G-1) is located 250 m from the reservoir with a V shape.

Upstream the B.G-1 gully is covered by tree strata forming a protective layer. In this area of the gully, the changes are not significant, but photo-interpretation of the PNOA orthophotographs shows a strong degradation of the vegetation cover due to the digging of the head of the gully. The orthophotographs generated by the SfM-MVS data showed significant immobilization of the wood planks on the gully bed (Figure 10). On the south-western south facing slope, the annual rate of erosion between the 1st and 3rd missions varied between 30 and 60 cm/year with vertical erosion. Conversely, a low erosion rate of no more than 10 cm/year was recorded on the north-northeast facing slope. However, the rate of erosion is not constant on the same slope. That said, the rate of erosion or sedimentation varies according to the slope, roughness, and exposure of each slope (Figure 10).

It should be noted that the sediment balance on the gully (B.G-1) during the wet season shows a dominance of erosion activity over the entire gully. The south-western south slope shows the highest erosion rates of the gully. On the other hand, there are spatial variations in the erosion rate within the same slope, highlighting the contrast in the existence of land use variation upstream of the gully.

On the gully bed this dynamic is reversed. The transfer of sediment to the bed has generated a fattening of the bed. Between the first and the third mission, the calculated sedimentation rate is of the order of 353.85 T/ha/year (442.31 m³/ha/year) (Figure 10) which does not reach the reservoir.

<figure><img src="image_86.png" /><figcaption>**Figure 10.** Bank gully (B.G-1)(A) DoD Mission I and II. (B) DoD Mission II and III. (C) DoD Mission I and III. (D) Graph of DoD volume. (E) Bar chart of DoD volume change distribution (erosion volume in red, deposition volume in blue.</figcaption></figure>

### 3.2.2. DOD in Site 2

The U-shaped B.G-2 gully, located near the reservoir, is highly exposed to erosion. The results measured by the DoD method have further highlighted the variation in the erosion rate between the two sides of the gully. The results show a strong activity of the annual sedimentation rate on the
north-facing slope as opposed to the south-facing slope where the annual erosion rate exceeds 90 cm of digging. Similarly, on the map, we observe significant changes at the head of the gully. These changes are reflected by an erosion rate of up to 2 m/year (Figure 11).

On the gully bed we observe the same sedimentation phenomenon as in Gully B.G-1 since this sector represents a sedimentary basin of the elements eroded by the banks and the head of the gully. The results of the DoD method show that the annual sedimentation rate at the bottom of the gully is of the order of 4264.1 T/ha/year (5330.1 m³/ha/year) (Figure 11D). At the time of maximum filling of the lake the bottom of the gully corresponds to the bottom of the lake.

<figure><img src="image_87.png" /><figcaption>**Figure 11.** Bank gully (B.G-2)(A) DoD Mission I and II. (B) DoD Mission II and III. (C) DoD Mission I and III. (D) Graph of DoD volume. (E) Bar chart of DoD volume change distribution (erosion volume in red, deposition volume in blue.</figcaption></figure>

### 3.2.3. DOD in Site 3

For the gully (B.G-3), the result of the sediment balance reveals that the rate of erosion is slightly lower than the rate of accumulation during the wet season. This is because the water level decreased by 1 m between the first and second field missions, which limits the sediment transport capacity by littoral drift and waves, but also that, in its bed, the gully experienced a strong sedimentation activity after the mixing of the lake water (regression deposits). Thus, during this season, the lakeshores were less watered. The sedimentation rate is the highest on all the slopes. The volume of sediment decreases from upstream to downstream of the gully bed. During the dry season, the banks were watered more, resulting in a high erosive potential (Figure 4) and a higher rate of slope erosion, estimated at 1772 m³. Thus, the dry season of 2018 contributed to the triggering of the gullying process after episodes of heavy rainfall. The particles eroded on the slopes are transported by the water to the gully bed and the lake. The rate of sediment accumulated in the gully bed is 1993 m³. The sedimentation rate of the gully (B.G-3) is estimated at 138.3 T/ha/year (172.7 m³/ha/year) (Figure 12).
<figure><img src="image_88.png" /><figcaption>**Figure 12.** Bank gully (B.G-3)(A) DoD Mission I and II. (B) DoD Mission II and III. (C) DoD Mission I and III. (D) Graph of DoD volume. (E) Bar chart of DoD volume change distribution (erosion volume in red, deposition volume in blue.</figcaption></figure>

### 3.2.4. M3C2 in Site 3

The Gully (B.G-3) is formed by overhanging slopes. To quantify the sedimentation rate of this particular morphology we used a 3D comparison of the point clouds using the M3C2 algorithm, in order to take into account the share of sediments coming from the overhanging areas. This method seems to us more complete than the DoD method, which allowed us to highlight the difference between 2.5 D digital terrain models. This gully is typical of heterogeneous soil formations rich in clay. It is developed on the concave bank of the Rambla de Algeciras where the reservoir constitutes a sedimentary basin of products from the eroded slopes of the gully. The gully (B.G-3) is the result of a combined erosion of the slopes and tunnels. This complex shape is representative of the fragile geomorphological bank gullies. The tunnels formed on Gully B.G-3 are related to the high clay contents in the upper soil horizon (river deposits). In this part of the gully the water carries away the clays and forms tunnels that develop until the ground collapses. From this we can deduce several results about the erosion rate illustrated in particular by distance measurements. Concerning the tunnel zones, the changes are almost zero. It appears that they remain stable, but these areas are very dangerous because their evolution is unpredictable. At the bottom of the gully, the sedimentation rate is very high, exceeding 250 m³ over an area of 2.28 ha. We notice that the erosion rate on both sides of the gully is similar, which seems logical for an area with steep slopes and no vegetation cover. It is clear from the distance measurements that the sedimentation rate at the bed of gully B.G-3 is higher than for the other gullies.

### 3.3. Detection of Geomorphological Changes from LiDAR-PNOA Data

The Figure 13 shows DoD values for all bank gullies along the lakeshore between 2009 and 2016. The value (minLOD = 0.28 m) was determined according to the equation (Equation 2). According to the gross DoD, about 45% of the lakeshore area has been gully eroded and the rest of the area has been
sedimented. The histogram (Figure 14B) of the distribution of elevation changes in the figure reflects
the positive bias of the eroded areas.

<figure><img src="image_89.png" /><figcaption>Figure 13. DoD in the lakeshore area for the survey PNOA 2009 and PNOA 2016.</figcaption></figure>

Interpretations of the signature of topographic changes can be justified, however, taking into
account the volume and area of the areas that have changed. The histogram (Figure 14A) shows
that the total surface of the gullies undergoes changes in altitude due to erosion. On the histogram,
we can visualize an imbalance between the rate of sedimentation and erosion on the lakeshore.
This unbalanced sediment balance is represented by the asymmetric distribution of the elevation
values around the value 0 with an erosion bias. Visualization of the histogram of the distribution of
eroded or sedimented volumes on the lakeshore provides a more detailed view of the quantity of
sediment deposited in the lake: between 2009 and 2016, 177,641 m³ (39 T/ha/year) of sediment was
deposited (Figure 14A).

<figure><img src="image_90.png" /><figcaption>**Figure 14.** (A) Bar chart of DoD volume change distribution. (B) Bar chart of DoD elevation change distribution.</figcaption></figure>
From the histogram (Figure 15), it can be seen that the highest morphological changes are in the eroded part with moderate concentrations but with a large amplitude reaching −0.98 m over a surface area of 900,000 m². The peak of the changes is not quite in the same area. This peak corresponds to a dominance of low amplitude accumulation processes of 0.7 m over an area of 1,300,000 m².

<figure><img src="image_91.png" /><figcaption>Figure 15. Histogram of elevation distribution.</figcaption></figure>

From the multimodal distribution of the histogram (Figure 15), it is easy to conclude the complexity of the morphodynamic events of the embankment gullies, especially since the significant temporal variability of the R-factor agrees on this point.

With the LiDAR PNOA data two sediment budgets were established. The first calculated balance is really the result of the mean RMSEZ error of 0.2 m. The second is inherent to the actual cubic content of the balance with statistical thresholds that vary little. The results of these two balances vary slightly at the gully scale, but on all lakeshores the difference between the results of these two balances is significant.

# 4. Discussion

## 4.1. Data Acquisition

The acquisitions of SfM-MVS-4D data by UAV met our needs in terms of accuracy, and frequency of monitoring the gullies along the lakeshore. Throughout the shooting, the flight altitude remained at a quasi-constant level at each site during the three missions, allowing us to keep the same pixel resolution for each gully. The UAV flew over the gullies at an altitude varying between 28 and 37 m, which allowed us to guarantee adequate resolution for monitoring topographic changes in the gullies. We also noted that the choice of these altitudes allowed us to ensure a small B/h ratio of the order of 0.6, in order to minimize as much as possible the surface of the hidden areas.

The results of the dense point clouds obtained by a vertical shot provide significant information on the topography of the gullies (B.G-1 & 2) compared to the gully (B.G-3). However, it should be noted that the photos taken on the gully (B.G-3) with a vertical camera do not cover the entire studied area. The overhang, having a vertical geometry, could not be taken into consideration. Moreover, the photos taken with an oblique camera were very useful to link the base and the top of the gully overhang slopes. Several previous studies have shown that the systematic errors of the dome effect significantly decrease or even eliminate after the use of an oblique and nadir image acquisition configuration, regardless of the image overlap rate, the camera angle and the configuration of the oblique and nadir images [46,47]. On the gullies (B.G-1 and B.G-2) the aerial photos are taken at nadir. The application of the integrated approach of the two flight configurations, i.e. nadir (90°) and oblique (30°) acquisition, reinforces the geometry of the photo network on the gully (B.G-3). This technique of oblique photo
acquisition increases the overlap rate between the photos, making it possible on the one hand to densify the density of the point clouds and on the other hand to fill the void created by the topographical complexity of the gullies (overhang). This therefore reduces georeferencing errors linked mainly to the low precision of the GCPs. We noted that the integrated approach has considerably reduced the RMSE in XYZ of the GCPs on the gully (B.G-3) with an average that does not exceed 1.4 cm in X, 1.2 cm in Y and 1.3 cm in Z during the three missions, bearing in mind that the GCPs went from 16 points to 12 points during the last mission. For the gullies (B.G-1; 2) it is the opposite, since the RMSE errors are high and vary between 2 and 5 cm.

## 4.2. Data Processing

Vertical and oblique photos are somewhat contrasting since the angle between the two shooting modes fluctuates very strongly along the same swath. Alignment in the PhotoScan software of the photos from two angles failed. The alignment results of the two shooting modes contain only the sparse point cloud from the nadir photos. To solve this problem we realigned all the photos using the common GCPs on both photos and the geolocation data of the photos before applying the bundle adjustment model.

Then we included the following camera parameters: Fit F, Cx, Cy, K1, K2, K3, P1, P2 (Table 4) which improved the camera model and reduced the RMSE errors of the photos as well as the deviations at the GCP point networks. The parameters were selected based on correlation analysis where high correlation parameters were eliminated to avoid over-configuration of the camera model. In this case study Mean precision across the camera is 0.10, this precision value will allow us to evaluate the correlation between the camera parameters. The low correlation between the parameters indicates a good image network except for the radial parameters, however when the correlation between F and Cy is very high this is explained by the absence of large camera rotation and a slight weakness of the image network. In our case study, the value is very close to the accuracy which illustrates the good image network taken in the field.

<table><caption>**Table 4.** Correlation between the camera parameters.</caption><thead><tr><th></th><th>F</th><th>Cx</th><th>Cy</th><th>K1</th><th>K2</th><th>K3</th><th>P1</th><th>P2</th></tr></thead><tbody><tr><th>F</th><td>1</td><td><strong>-0.61</strong></td><td><strong>-0.13</strong></td><td>-0.07</td><td>0.06</td><td>-0.08</td><td><strong>0.51</strong></td><td><strong>0.22</strong></td></tr><tr><th>Cx</th><td></td><td>1</td><td>-0.04</td><td>-0.05</td><td>0.8</td><td>-0.06</td><td><strong>0.33</strong></td><td>-0.03</td></tr><tr><th>Cy</th><td></td><td></td><td>1</td><td>0.01</td><td>0.02</td><td>-0.02</td><td><strong>0.17</strong></td><td><strong>0.76</strong></td></tr><tr><th>K1</th><td></td><td></td><td></td><td>1</td><td><strong>-0.80</strong></td><td><strong>0.90</strong></td><td>-0.03</td><td>0.02</td></tr><tr><th>K2</th><td></td><td></td><td></td><td></td><td>1</td><td><strong>-0.89</strong></td><td>0.02</td><td>0.00</td></tr><tr><th>K3</th><td></td><td></td><td></td><td></td><td></td><td>1</td><td>0.00</td><td>0.01</td></tr><tr><th>P1</th><td></td><td></td><td></td><td></td><td></td><td></td><td>1</td><td><strong>0.17</strong></td></tr><tr><th>P2</th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>1</td></tr></tbody></table>

## 4.3. Error Modeling

We took advantage of the Lidar PNOA data collected between 2009 and 2016 to continue the analysis of the morpho-sedimentary dynamics of the three gullies. Similarly, the work was enriched by taking into account the spatio-temporal variability of the dynamics of all the bank gullies as a function of topographical influence and climate action. We present and analyze here the results of the raw DoD method obtained from the difference between the two DEMs Lidar 2009 and 2016 of the PNOA project. however, the error of the two Lidar missions PNOA 2009 and 2016 is of the order of 0.28 m this value is considered spatially homogeneous and uniform for two main reasons:

(1) The low density of the Lidar point clouds (between 0.5 and $1/m^2$): the denser the point cloud is, the more robust it is to accurately determine terrain parameters such as slope or roughness.

(2) The absence of GCPs: the georeferencing errors of the GCPs are essential to estimate an error model of the Lidar PNOA data.
To quantify the uncertainty of the SfM-MVS data accurately, it is necessary to validate the SfM-MVS data by a comparison with previous SfM-MVS data or with a very high precision field truth (DGPS, TS). The method of validation is often dictated by the extent, scale, field survey method, data format, and morphological characteristics of the study site.

In this paper we first chose to quantify spatially variable SfM DEM uncertainty, using the cross-validation method instead of Bootstrapping which has already been used by Wheaton and [35] and Rossi [48]. Then the uncertainty values will be used to build an error model that takes into account the different sources of uncertainty using the FIS model.

Several works have shown that the vertical accuracy on the GCPs control points is about two to three times higher than the accuracy obtained by the products derived from the SfM-MVS data (DTM and orthophotography) [49,50]. The table (Table 5) shows that the FIS error propagation model estimates higher vertical errors on GCPs points than those obtained after georeferencing (Table 3). The accuracy of the SfM-MVS data obtained during this study is comparable to that of the GNSS/TS data [35,51].

<table><caption>Table 5. FIS errors in X, Y and Z of the GCPs.</caption><thead><tr><th rowspan="2">Sites</th><th colspan="3">B.G-1</th><th colspan="3">B.G-2</th><th colspan="3">B.G-3</th></tr><tr><th>M-I</th><th>M-II</th><th>M-III</th><th>M-I</th><th>M-II</th><th>M-III</th><th>M-I</th><th>M-II</th><th>M-III</th></tr></thead><tbody><tr><td>Missions</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>X (cm)</td><td>3.62</td><td>3.63</td><td>4.02</td><td>6.21</td><td>6.33</td><td>6.52</td><td>2.63</td><td>2.21</td><td>2.54</td></tr><tr><td>Y (cm)</td><td>4.45</td><td>6.12</td><td>3.24</td><td>4.06</td><td>6.87</td><td>6.77</td><td>2.31</td><td>3.12</td><td>3.05</td></tr><tr><td>Z (cm)</td><td>10.53</td><td>12.69</td><td>9.35</td><td>11.55</td><td>12.81</td><td>11.31</td><td>2.80</td><td>3.32</td><td>3.54</td></tr></tbody></table>

### 4.4. Detection of Geomorphological Changes

From the Figure 16 we observe a ramified gully system on the entire lakeshore, with a gully density varying between 18Km/Km² and 59 Km/Km² with an annual erosion rate of 39 T/ha/year these values are extreme are of the same order as some areas that reach values varying between 15% and 60% [52–55].

<figure><img src="image_92.png" /><figcaption>Figure 16. Bank gullies density.</figcaption></figure>

Globally, the volume eroded by gullies varies between 0.002 and 47,430 m³/year with an average of 358.6 m³/year. In the Murcia region, research on the gully erosion is limited. As a result, there is very little data on the gully erosion budget. Indeed, in the literature, we find two references on the sediment balance of gullies and bank gullies in the Guadalentin watershed downstream of the
Rambla de Algeciras watershed, namely the work of L. Vandekerckhove [56] and Marzolff [57]. Before comparing the results of this study with those of Vandekerckhove and Marzolff, it should be pointed out that the very complex morphology, the very friable pedology and the strong anthropic activity in the Rambla de Algeciras catchment area allowed us to consider that the Guadalentin basin is less vulnerable to gully erosion than the Rambla de Algeciras basin. Despite this observation, the volumes calculated in this work are close to the results estimated in this study. The annual volume calculated by Vandekerckhove over some time varies between 40 and 43 years based on aerial photos and field measurements vary between 13 and 3467 m³/year. Similarly, Marzolff's work in 2011 using aerial photos taken by the UAV showed that the rate varies between 0.5 and 100 m³/year. On the gully (B.G-3) we have estimated the volume eroded from the calculation of the distance between the 3D SfM-MVS point clouds according to the direction of the surface normals, rather than the difference between two 2.5 D rasters. This technique takes into account the volume of sediments eroded on the gully overhang area. The M3C2 algorithm calculates the distance between the multi-temporal point clouds, which allows to estimate the eroded volume with high accuracy. The results of the M3C2 technique thus testify on the good choice of an oblique and nadir acquisition configuration on the gully (B.G-3) to generate accurate and complete SfM-MVS data.

The analysis of the geomorphological changes estimated by the UAV and LiDAR are very different, with the results of the LiDAR data the percentage of the eroded area over the whole period 2009–2016 is not more than (45%) however the deposit areas are more important with a percentage of (55%). On the other hand, the estimates of the sedimentary balance during the year 2018 with the UAV data allow us to observe an erosion activity that affected on average more than 60% of the surface of the studied gullies. The results of the LiDAR data appeared aberrant but the monitoring of the fluctuation of the lake level explains the strong sedimentation activity on the lake banks between 2009 and 2016. The rise of lake level is related to the high trend in rainfall inputs between 2010 and 2015 (Figure 17). It should also be noted that the lake level decreased from 260 m in 2014 to 240 m in 2018. This is related to the decrease in rainfall inputs between 2015 and 2017. This water level fluctuation has caused the settling of regression deposits. After the period of maximum lake filling at 260.12 m (4 May 2014), the lakeshore has experienced silt sedimentation in flat areas and gully beds. These sediments have four origins:

* Gullying (Figure 18)

* Bluff slumping (Figure 19)

* Biogenic (biological activity)

* Detritic ( torrential inputs of Rambla)

<figure><img src="image_93.png" /><figcaption>Figure 17. Lake level fluctuation between 2006 and 2019.</figcaption></figure>

The photo-interpretation of PNOA orthophotography illustrates an imbalance in the
morphodynamics of the lakeshore.
When the lake level reached the rills and small gullies, the vertical movement of the lake water overwhelmed them and caused them to lose their erosion ability (Figure 18). This has led to the temporal disappearance of these small forms of erosion. They were then replaced by the sedimentary activity of waves loaded with suspended matter (a process linked to longshore drift).

<figure><img src="image_94.png" /><figcaption>Figure 18. Temporal disappearance of the small bank gullies.</figcaption></figure>

The missing data on the solid transport at gauging stations tributaries of the Rambla, prevent from quantifying the sediment inputs into the Rambla de Algeciras, as it has been done in Morocco by Ezzaouini et al [58], and unfortunately also, we have not been able to install sedimentometres [59] in the lake to collect and analyze the nature and percentage of detrital and biogenic sediments, mainly due to time and money issues.

The bluff slumping of the lakeshore, observed through photo-interpretation of the PNOA aerial orthophotographs between 2009 and 2016, is a phenomenon triggered during the maximum filling of the lake. The former overhanging areas are destabilized by moisture and weight. This phenomenon is localized in particular on the meanders and concave banks of the Rambla (Figure 19).

<figure><img src="image_95.png" /><figcaption>Figure 19. Bluff slumping.</figcaption></figure>
**5. Conclusions**

The friable nature of the soil, the fluctuation of the lake level, anthropic activity and climate
changes marked by increasingly rare but very violent rainfall episodes have made it possible to carve
out several forms of gullies on the lakeshores of Rambla de Algeciras. The main goal of this paper is to
demonstrate whether the SfM-MVS-4D model generated from photos taken by a UAV and the LiDAR
PNOA data, are able to accurately quantify the sediment budget on the banks gullies of the lakeshores.
The answer is yes. Indeed, for bank gully monitoring on lakeshores, we detected geomorphological
changes on three types of banks gullies. The models for the detection of geomorphological changes
were constructed using multi-temporal numerical terrain models coupled with spatially variable
uncertainty models. The rate of geomorphological changes detected on the three gullies in 2018 is
mainly related to orientation, land use and the shape of the gullies.

The aerial photos taken by the UAV are capable of generating SfM-MVS-4D models with very
high spatial and temporal resolution, which allows to highlight the movement of sediments along
the slopes of the gullies. The use of the nadir and oblique configuration for taking pictures has
improved the accuracy of the SfM-MVS models taken on very steep slopes or with overhangs. This
accuracy was measured using several techniques, namely georeferencing error estimation and the
statistical cross-validation approach and finally the FIS technique. This was useful for detecting
geomorphological changes on banks gullies. The accuracy of the SfM-MVS data varies between 1 and
5 cm.

The data from the PNOA project (orthophotographs and LiDAR) provided reliable and accurate
data to study the effects of lake level fluctuation, climate change and land use change on the initiation
of lake bank gullying. The DOD technique applied on the LiDAR data from the PNOA allowed the
estimation of an annual sedimentation rate of 39 T/ha/year (49 m³/ha/year) from the gullies of the
Rambla de Algeciras lakeshores.

We also note that the approach proposed in this paper will be deployed in Morocco on the
Sidi Mohamed Ben Abdellah dam, which shares the same bank gully on lakeshores issues and
geo-environmental characteristics as our study site.

**Author Contributions:** Conceptualization, V.M. and R.H.; Data curation, R.C.; Funding acquisition, R.C. and F.C.; Methodology, V.M., R.H. and F.C.; Project administration, G.M. and E.R.; Resources, V.M, E.R. and R.C.; Software, F.C.; Validation, R.H.; Writing—original draft, R.H.; Writing—review and editing, V.M, G.M. and R.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** This work has been financed with the following funds: GEOLAB laboratory and CASBA project and is gratefully acknowledged. The authors also want to thank Confederación Hidrográfica del Segura in Murcia, Mr Pablo Albaladejo and all colleagues, engineer and technicians participating in the field campaigns supporting successful data acquisition. Constructive and valuable feedback given from the anonymous reviewers helped to increase the quality of the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

**References**

1. Hentati, A.; Kawamura, A.; Amaguchi, H.; Nakagawa, N. Erosion assessment at small hillside river basins in semiarid region of Tunisia. *Annu. J. Hydraul. Eng. Isce* **2010**, accepted under publication

2. Antipolis, S. Les menaces sur les sols dans les pays méditerranéens (Etude bibliographique). *Les Cahiers du 690 Plan Bleu* **2003**, 2, 44–48.

3. Inoubli, N. Ruissellement et Éronsion Hydrique en Milieu Méditerranéen Vertique: Approche Expérimentale et Modélisation. Ph.D. Thesis, Montpellier, SupAgro, Montpellier, France, 2016.

4. Castillo, C.; Gómez, J. A century of gully erosion research: Urgency, complexity and study approaches. *Earth-Sci. Rev.* **2016**, *160*, 300–319. [CrossRef]

5. Poesen, J.; Nachtergaele, J.; Verstraeten, G.; Valentin, C. Gully erosion and environmental change: importance and research needs. *Catena* **2003**, *50*, 91–133. [CrossRef]
6. Bull, L.J.; Kirkby, M.J. *Dryland Rivers: Hydrology and Geomorphology of Semi-Arid Channels*. JohnWiley & Sons: Chichester, UK, 2002.

7. Vanmaercke, M.; Poesen, J.; Van Mele, B.; Demuzere, M.; Bruynseels, A.; Golosov, V.; Bezerra, J.F.R.; Bolysov, S.; Dvinskih, A.; Frankl, A.; et al. How fast do gully headcuts retreat? *Earth-Sci. Rev.* **2016**, *154*, 336–355. [CrossRef]

8. Nazari Samani, A.; Tavakoli Rad, F.; Azarakhshi, M.; Reza Rahdari, M.; Rodrigo-Comino, J. Assessment of the sustainability of the territories affected by gully head advancements through aerial photography and modeling estimations: A case study on Samal Watershed, Iran. *Sustainability* **2018**, *10*, 2909. [CrossRef]

9. Valentin, C. *Les Sols au Coeur de la Zone Critique 5: DÉgradation et Réhabilitation*; ISTE Group: Santa Clara, CA, USA, 2018; Volume 5.

10. Wu, Y.; Cheng, H. Monitoring of gully erosion on the Loess Plateau of China using a global positioning system. *Catena* **2005**, *63*, 154–166. [CrossRef]

11. Mockus, V. *National Engineering Handbook*; U.S. Soil Conservation Service: Washington, DC, USA, 1964; Volume 4.

12. Thompson, J.R. Quantitative effect of watershed variables on rate of gully-head advancement. *Trans. Asae* **1964**, *7*, 54-0055. [CrossRef]

13. Seginer, I. Gully development and sediment yield. *J. Hydrol.* **1966**, *4*, 236–253. [CrossRef]

14. Yuan, Y.; Bingner, R.; Rebich, R. Evaluation of AnnAGNPS on Mississippi delta MSEA watersheds. *Trans. Asae* **2001**, *44*, 1183. [CrossRef]

15. Nachtergaele, J.; Poesen, J.; Vandekerckhove, L.; Oostwoud Wijdenes, D.; Roxo, M. Testing the ephemeral gully erosion model (EGEM) for two Mediterranean environments. *Earth Surf. Process. Landf. J. Br. Geomorphol. Res. Group* **2001**, *26*, 17–30. [CrossRef]

16. Sidorchuk, A. Dynamic and static models of gully erosion. *Catena* **1999**, *37*, 401–414. [CrossRef]

17. Tucker, G.; Lancaster, S.; Gasparini, N.; Bras, R. The channel-hillslope integrated landscape development model (CHILD). In *Landscape Erosion and Evolution Modeling*; Springer: Cham, Switzerland, 2001; pp. 349–388.

18. Evans, M.; Lindsay, J. High resolution quantification of gully erosion in upland peatlands at the landscape scale. *Earth Surf. Process. Landf.* **2010**, *35*, 876–886. [CrossRef]

19. Maleval, V. *Les Apports Sédimentaires Directs au Barrage SMBA*; Workshop SIGMED. 2012. Available online: https://hal-bnf.archives-ouvertes.fr/GEOLAB/hal-00924798v1/ (accessed on 29 September 2020).

20. D'Oleire Oltmanns, S.; Marzolff, I.; Peter, K.D.; Ries, J.B. Unmanned aerial vehicle (UAV) for monitoring soil erosion in Morocco. *Remote Sens.* **2012**, *4*, 3390–3416. [CrossRef]

21. Castillo, C.; James, M.; Redel-Macías, M.; Pérez, R.; Gómez, J. SF3M software: 3-D photo-reconstruction for non-expert users and its application to a gully network. *Soil* **2015**, *1*, 583. [CrossRef]

22. Wang, R.; Zhang, S.; Pu, L.; Yang, J.; Yang, C.; Chen, J.; Guan, C.; Wang, Q.; Chen, D.; Fu, B.; et al. Gully erosion mapping and monitoring at multiple scales based on multi-source remote sensing data of the Sancha River Catchment, Northeast China. *Isprb Int. J. Geo-Inf.* **2016**, *5*, 200. [CrossRef]

23. Di Stefano, C.; Ferro, V.; Palmeri, V.; Pampalone, V.; Agnello, F. Testing the use of an image-based technique to measure gully erosion at Sparacia experimental area. *Hydrol. Process.* **2017**, *31*, 573–585. [CrossRef]

24. Frankl, A.; Stal, C.; Abraha, A.; Nyssen, J.; Rieke-Zapp, D.; De Wulf, A.; Poesen, J. Detailed recording of gully morphology in 3D through image-based modelling. *Catena* **2015**, *127*, 92–101. [CrossRef]

25. Christian, P.; Davis, J. Hillslope gully photogeomorphology using structure-from-motion. *Z. Geomorphol. Suppl. Issues* **2016**, *60*, 59–78. [CrossRef]

26. Koci, J.; Jarihani, B.; Leon, J.X.; Sidle, R.C.; Wilkinson, S.N.; Bartley, R. Assessment of UAV and ground-based structure from motion with multi-view stereo photogrammetry in a gullied savanna catchment. *Isprsb Int. J. Geo-Inf.* **2017**, *6*, 328. [CrossRef]

27. Stocker, T.F.; Qin, D.; Plattner, G.K.; Tignor, M.; Allen, S.K.; Boschung, J.; Nauels, A.; Xia, Y.; Bex, V.; Midgley, P.M.; et al. Climate change 2013: The physical science basis. *Contrib. Work. Group Fifth Assess. Rep. Intergov. Panel Clim. Chang.* **2013**, 1535.

28. Zhang, Z. Iterative point matching for registration of free-form curves and surfaces. *Int. J. Comput. Vis.* **1994**, *13*, 119–152. [CrossRef]

29. Stöcker, C.; Eltner, A.; Karrasch, P. Measuring gullies by synergetic application of UAV and close range photogrammetry—A case study from Andalusia, Spain. *Catena* **2015**, *132*, 1–11. [CrossRef]
30. Micheletti, N.; Chandler, J.H.; Lane, S.N. Investigating the geomorphological potential of freely available and accessible structure-from-motion photogrammetry using a smartphone. *Earth Surf. Process. Landf.* **2015**, *40*, 473–486. [CrossRef]

31. Lague, D.; Brodu, N.; Leroux, J. Accurate 3D comparison of complex topography with terrestrial laser scanner: Application to the Rangitikei canyon (NZ). *Isprs J. Photogramm. Remote. Sens.* **2013**, *82*, 10–26. [CrossRef]

32. Mosbrucker, A.R.; Major, J.J.; Spicer, K.R.; Pitlick, J. Camera system considerations for geomorphic applications of SfM photogrammetry. *Earth Surf. Process. Landf.* **2017**, *42*, 969–986. [CrossRef]

33. Cucchiaro, S.; Maset, E.; Fusiello, A.; Cazorzi, F. 4D-SFM photogrammetry for monitoring sediment dynamics in a debris-flow catchment: Software testing and results comparison. In Proceedings of the 2018 ISPRS TC II Mid-term Symposium “Towards Photogrammetry 2020”. International Society for Photogrammetry and Remote Sensing, Riva Del Garda, Italy, 4–7 June 2018; Volume 42, pp. 281–288.

34. Vericat, D.; Smith, M.; Brasington, J. Patterns of topographic change in sub-humid badlands determined by high resolution multi-temporal topographic surveys. *Catena* **2014**, *120*, 164–176. [CrossRef]

35. Wheaton, J.M.; Brasington, J.; Darby, S.E.; Sear, D.A. Accounting for uncertainty in DEMs from repeat topographic surveys: improved sediment budgets. *Earth Surf. Process. Landf. J. Br. Geomorphol. Res. Group* **2010**, *35*, 136–156. [CrossRef]

36. Brasington, J.; Vericat, D.; Rychkov, I. Modeling river bed morphology, roughness, and surface sedimentology using high resolution terrestrial laser scanning. *Water Resour. Res.* **2012**, *48*. [CrossRef]

37. Bangen, S.; Hensleigh, J.; McHugh, P.; Wheaton, J. Error modeling of DEMs from topographic surveys of rivers using fuzzy inference systems. *Water Resour. Res.* **2016**, *52*, 1176–1193. [CrossRef]

38. Lane, S.N.; Chandler, J.H. The generation of high quality topographic data for hydrology and geomorphology: new data sources, new applications and new problems. *Earth Surf. Process. Landf. J. Br. Geomorphol. Res. Group* **2003**, *28*, 229–230. [CrossRef]

39. Milan, D.J.; Heritage, G.L.; Large, A.R.; Fuller, I.C. Filtering spatial error from DEMs: Implications for morphological change estimation. *Geomorphology* **2011**, *125*, 160–171. [CrossRef]

40. Passalacqua, P.; Belmont, P.; Staley, D.M.; Simley, J.D.; Arrowsmith, J.R.; Bode, C.A.; Crosby, C.; DeLong, S.B.; Glenn, N.F.; Kelly, S.A.; et al. Analyzing high resolution topography for advancing the understanding of mass and energy transfer through landscapes: A review. *Earth-Sci. Rev.* **2015**, *148*, 174–193. [CrossRef]

41. Victoriano, A.; Brasington, J.; Guinau, M.; Furdada, G.; Cabré, M.; Moysset, M. Geomorphic impact and assessment of flexible barriers using multi-temporal LiDAR data: The Portainé mountain catchment (Pyrenees). *Eng. Geol.* **2018**, *237*, 168–180. [CrossRef]

42. Brasington, J.; Rumsby, B.; McVey, R. Monitoring and modelling morphological change in a braided gravel-bed river using high resolution GPS-based survey. *Earth Surf. Process. Landf. J. Br. Geomorphol. Res. Group* **2000**, *25*, 973–990. [CrossRef]

43. Brasington, J.; Langham, J.; Rumsby, B. Methodological sensitivity of morphometric estimates of coarse fluvial sediment transport. *Geomorphology* **2003**, *53*, 299–316. [CrossRef]

44. Lane, S.N.; Westaway, R.M.; Murray Hicks, D. Estimation of erosion and deposition volumes in a large, gravel-bed, braided river using synoptic remote sensing. *Earth Surf. Process. Landf. J. Br. Geomorphol. Res. Group* **2003**, *28*, 249–271. [CrossRef]

45. James, M.R.; Robson, S.; Smith, M.W. 3-D uncertainty-based topographic change detection with structure-from-motion photogrammetry: precision maps for ground control and directly georeferenced surveys. *Earth Surf. Process. Landf.* **2017**, *42*, 1769–1788. [CrossRef]

46. Nesbit, P.R.; Hugenholtz, C.H. Enhancing UAV-SFM 3D model accuracy in high-relief landscapes by incorporating oblique images. *Remote. Sens.* **2019**, *11*, 239. [CrossRef]

47. James, M.R.; Robson, S. Mitigating systematic error in topographic models derived from UAV and ground-based image networks. *Earth Surf. Process. Landf.* **2014**, *39*, 1413–1420. [CrossRef]

48. Rossi, R.K. Evaluation of Structure-from-Motion from a Pole-Mounted Camera for Monitoring Geomorphic Change. Ph.D. Thesis, Utah State University, Logan, UT, USA, 2018.

49. Fernández, T.; Pérez, J.L.; Cardenal, J.; Gómez, J.M.; Colomo, C.; Delgado, J. Analysis of landslide evolution affecting olive groves using UAV and photogrammetric techniques. *Remote Sens.* **2016**, *8*, 837. [CrossRef]
50. Fernández, T.; Pérez, J.L.; Colomo, C.; Cardenal, J.; Delgado, J.; Palenzuela, J.A.; Irigaray, C.; Chacón, J. Assessment of the evolution of a landslide using digital photogrammetry and LiDAR techniques in the Alpujarras region (Granada, southeastern Spain). *Geosciences* **2017**, *7*, 32. [CrossRef]

51. Rumsby, B.; Brasington, J.; Langham, J.; McLelland, S.; Middleton, R.; Rollinson, G. Monitoring and modelling particle and reach-scale morphological change in gravel-bed rivers: Applications and challenges. *Geomorphology* **2008**, *93*, 40–54. [CrossRef]

52. Dewitte, O.; Daoudi, M.; Bosco, C.; Van Den Eeckhaut, M. Predicting the susceptibility to gully initiation in data-poor regions. *Geomorphology* **2015**, *228*, 101–115. [CrossRef]

53. Martínez-Casasnovas, J. A spatial information technology approach for the mapping and quantification of gully erosion. *Catena* **2003**, *50*, 293–308. [CrossRef]

54. Perroy, R.L.; Bookhagen, B.; Asner, G.P.; Chadwick, O.A. Comparison of gully erosion estimates using airborne and ground-based LiDAR on Santa Cruz Island, California. *Geomorphology* **2010**, *118*, 288–300. [CrossRef]

55. Nadal-Romero, E.; Revuelto, J.; Errea, P.; López-Moreno, J. The application of terrestrial laser scanner and SfM photogrammetry in measuring erosion and deposition processes in two opposite slopes in a humid badlands area (central Spanish Pyrenees). *Soil* **2015**, *1*, 561. [CrossRef]

56. Vandekerckhove, L.; Poesen, J.; Govers, G. Medium-term gully headcut retreat rates in Southeast Spain determined from aerial photographs and ground measurements. *Catena* **2003**, *50*, 329–352. [CrossRef]

57. Marzolff, I.; Ries, J.B.; Poesen, J. Short-term versus medium-term monitoring for detecting gully-erosion variability in a Mediterranean environment. *Earth Surf. Process. Landf.* **2011**, *36*, 1604–1623. [CrossRef]

58. Ezzaouini, M.A.; Mahé, G.; Kacimi, I.; Zerouali, A. Comparison of the MUSLE Model and Two Years of Solid Transport Measurement, in the Bouregreg Basin, and Impact on the Sedimentation in the Sidi Mohamed Ben Abdellah Reservoir, Morocco. *Water* **2020**, *12*, 1882. [CrossRef]

59. Maleval, V.; Pitois, F. Fonctionnement hydro-sédimentaire et bilan sédimentaire du lac de Saint-Germain-de-Confolens (Charente, France): paramètres géomorphologiques à prendre en considération dans la gestion du lac. *Phys.-Géo. Géo. Phys. Environ.* **2017**, 197–227. [CrossRef]

<figure><img src="image_96.png" /></figure>

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license ([http://creativecommons.org/licenses/by/4.0/](http://creativecommons.org/licenses/by/4.0/)).
Article

# Understory Limits Surface Runoff and Soil Loss in Teak Tree Plantations of Northern Lao PDR

Layheang Song ¹,², Laurie Boithias ¹,*, Oloth Sengtaheuanghoung ³, Chantha Oeurng ², Christian Valentin ⁴, Bounthan Souksavath ⁵, Phabvilay Sounyafong ⁶, Anneke de Rouw ⁴, Bounsamay Soulileuth ⁶, Norbert Silvera ⁴, Bounchanh Lattanavongkot ⁵, Alain Pierret ⁶ and Olivier Ribolzi ¹

1. GET, Université de Toulouse, CNRS, IRD, UPS, 31400 Toulouse, France; layheangsong@itc.edu.kh (L.S.); olivier.ribolzi@ird.fr (O.R.)

2. Research and Innovation Center (RIC), Faculty of Hydrology and Water Resources Engineering, Institute of Technology of Cambodia, Phnom Penh 12156, Cambodia; chantha@itc.edu.kh

3. Department of Agricultural Land Management (DALaM), Ministry of Agriculture and Forestry, Vientiane 01000, Laos; oloth_s@hotmail.com

4. Institute of Ecology and Environmental Sciences of Paris (iEES-Paris), Sorbonne Université, Univ Paris Est Creteil, IRD, CNRS, INRA, 75005 Paris, France; christian.valentin@ird.fr (C.V.); anneke_de_rouw@ird.fr (A.d.R.); norbert.silvera@ird.fr (N.S.)

5. Luang Prabang Teak Program (LPTP-TFT), Luang Prabang Provincial Forestry Section, Luang Prabang 06000, Laos; b.souksavath@earthworm.org (B.S.); chanh.lptp@gmail.com (B.L.)

6. IRD, IEES-Paris UMR 242, c/o National Agriculture and Forestry Research Institute, Vientiane 01000, Laos; phabvilay_laopdr@yahoo.com (P.S.); sbounsamay@gmail.com (B.S.); alain.pierret@ird.fr (A.P.)

* Correspondence: laurie.boithias@get.omp.eu; Tel.: +33-561-332-712

Received: 14 July 2020; Accepted: 15 August 2020; Published: 19 August 2020

**Abstract:** Many mountainous regions of the humid tropics experience serious soil erosion following rapid changes in land use. In northern Lao People's Democratic Republic (PDR), the replacement of traditional crops by tree plantations, such as teak trees, has led to a dramatic increase in floods and soil loss and to the degradation of basic soil ecosystem services such as water filtration by soil, fertility maintenance, etc. In this study, we hypothesized that conserving understory under teak trees would protect soil, limit surface runoff, and help reduce soil erosion. Using 1 m² microplots installed in four teak tree plantations in northern Lao PDR over the rainy season of 2017, this study aimed to: (1) assess the effects on surface runoff and soil loss of four understory management practices, namely teak with no understory (TNU; control treatment), teak with low density of understory (TLU), teak with high density of understory (THU), and teak with broom grass, *Thysanolaena latifolia* (TBG); (2) suggest soil erosion mitigation management practices; and (3) identify a field visual indicator allowing a rapid appraisal of soil erosion intensity. We monitored surface runoff and soil loss, and measured teak tree and understory characteristics (height and percentage of cover) and soil surface features. We estimated the relationships among these variables through statistics and regression analyses. THU and TBG had the smallest runoff coefficient (23% for both) and soil loss (465 and 381 g·m⁻², respectively). The runoff coefficient and soil loss in TLU were 35% and 1115 g·m⁻², respectively. TNU had the highest runoff coefficient and soil loss (60%, 5455 g·m⁻²) associated to the highest crusting rate (82%). Hence, the soil loss in TBG was 14-times less than in TNU and teak tree plantation owners could divide soil loss by 14 by keeping understory, such as broom grass, within teak tree plantations. Indeed, a high runoff coefficient and soil loss in TNU was explained by the kinetic energy of rain drops falling from the broad leaves of the tall teak trees down to bare soil, devoid of plant residues, thus leading to severe soil surface crusting and soil detachment. The areal percentage of pedestal features was a reliable indicator of soil erosion intensity. Overall, promoting understory, such as broom grass, in teak tree plantations would: (1) limit surface runoff and improve soil infiltrability, thus increase soil water stock available for both root absorption and groundwater recharge; and (2) mitigate soil loss while favoring soil fertility conservation.
**Keywords:** overland flow; inter-rill erosion; teak tree plantation; understory; broom grass; South-East Asia; land management; soil erosion

# 1. Introduction

Mountain areas of the humid tropics are characterized by steep slopes and heavy rains [1]. Hence these regions are prone to high surface runoff and soil erosion [2]. Conversion of natural forest to e.g., agricultural land exacerbates runoff production and soil erosion, leading to physicochemical and biological changes of the altered ecosystems [3,4].

On-site effects of increased soil erosion include the reduction of soil quality impacting the sustainability of agricultural production [5], and economics, due to the loss of ecosystem services [6]. The tremendously higher rate of soil loss compared to its formation rate threatens food production and environmental quality (water, soil, and air) [7–9].

Off-site effects comprise floods, depletion of groundwater recharge, degradation of stream environments, and downstream sedimentation [5,10]. In addition, the adsorption of organic and inorganic matter on soil particles and suspended sediments plays a leading role in the transport of nutrients [11], radionuclides [12], metals [13], pesticides [14], and bacterial pathogens [15,16]. Higher contamination of rivers by fecal bacteria is often correlated to higher in-stream suspended sediment concentration [17,18]. Sediment loads cause massive accumulation of fine sediments to river beds and cause the siltation of irrigation canals and dam reservoirs [19–21], thus reducing their life spans. All of these off-site effects lead to high economic and ecological costs (i.e., sedimentation, flooding, landslides, water eutrophication, biodiversity loss, land abandonment, destruction of infrastructures) [6].

Natural forests are known for their protective effect against erosion [22]. Forests have favorable hydrodynamic properties for surface infiltration [23], subsurface drainage [24] due to the biological activity, and the development of macropores in the litter and underneath. A fraction of rainfall, known as throughfall, is intercepted by the canopy and the underlying vegetative strata [25]. The amount of water and the kinetic energy of drops corresponding to throughfall is lower than that of rainfall, thus reducing soil erosion. Kinetic energy is the main factor initiating soil erosion; raindrops hit the soil surface and disaggregate the soil structure [26] resulting in redistribution of soil material by splash effect [5]. Kinetic energy is controlled by rainfall (amount, drop size, fall velocity), vegetation characteristics (height, cover, residues) [27], and slope [26].

In the mountainous region of northern Lao People's Democratic Republic (PDR), widespread agricultural practices are known to increase surface runoff, soil erosion, and in-stream suspended sediment concentrations compared to natural forest [3,23]. In particular, commercial perennial monocultures, such as teak tree plantations, lead to a drastic surface runoff increase [28] and exacerbate soil loss along hillslopes [29,30]. In teak tree plantations, such as the old teak trees (high timbers with broad leaves) growing in the Houay Pano catchment, surface runoff and soil erosion are likely related to recurrent understory and leaf litter suppression by burning at the end of the dry season, leaving the soil bare and exposed to the kinetic energy of the raindrops [31].

Land use, understory cover, and soil surface condition, are the main drivers influencing the infiltration of water into the soil [32,33]. Living plant roots modify both mechanical and hydrological characteristics of the soil matrix and negatively influence the soil erodibility [34,35]. Understory and plant residues on soil surface are known to efficiently attenuate the effect of splash and thus soil detachment [36–40]. Soil surface features have a strong impact on surface runoff and soil erosion: soil aggregation limits surface runoff generation and soil loss [41,42], while crusted soils can be self-protecting from erosion because of their high surface shear stress resistance, but they also promote surface runoff, hence downstream erosion, especially on non-crusted areas [26,30]. For example, in northern Lao PDR, splash in teak tree plantations increases soil erosion, and the absence of
understory enhances the effect of splash [29]. By adapting the model developed by [43], Patin et al. [31] showed that plant residues and weed cover at soil surface were the main attenuation factors of soil erosion. In northern Thailand, soil erosion in rubber tree plantations decreased when understory was grown [44]. Indeed, in northern Lao PDR headwater catchments, the runoff coefficient is approximately 55% under teak trees at the plot scale [30], whereas it nearly doubled from 16% to 31% within 13 years at the catchment scale [29], resulting from the absence of understory in the teak tree plantation area. In particular, mature teak trees with limited understory were shown to export respectively 5.5- and 31-times more water and soil than broom grass (*Thysanolaena latifolia*) at the plot scale [30].

Hence, in this paper, we hypothesized that conserving understory such as broom grass, which provides income to farmers through broom making and selling [45], protects soil, limits surface runoff, and helps to reduce soil erosion in teak tree plantations. The overarching goal of this study, conducted in the teak tree plantations of Ban Kokngew, northern Lao PDR, was thus to find the best strategy to contribute to sustainable agricultural production. In this context, the three objectives of our 1 m² microplot experiment performed during the June to October monsoon period of 2017 were to: (1) assess the effects on surface runoff and soil loss of four understory management practices, namely teak with no understory (TNU; control treatment), teak with low density of understory (TLU), teak with high density of understory (THU), and teak with broom grass (TBG); (2) suggest soil erosion mitigation management practices; and (3) identify a field visual indicator allowing a rapid appraisal of soil erosion intensity.

## 2. Materials and Methods

### 2.1. *Study Area and Experimental Plots*

We conducted the experiment in 2017 in teak tree plantations surrounding Ban Kokngew, a village located in Luang Prabang Province, northern Lao PDR, and predominantly situated over Acrisol soil and Carboniferous and Permian limestones (Figure 1). The climate is sub-tropical humid and is characterized by a monsoon regime with a dry season from November to May, and a rainy from June to October. Mean annual rainfall recorded at Luang Prabang from 1960 to 2006 was 1268 mm, about 76% of which falls during the rainy season. The mean annual temperature is 25.3 °C. Mean annual reference evapotranspiration is 1116 mm. The study area belongs to the mountainous region of northern Laos PDR. More specifically, the area is located within the "Luang Prabang mountain rain forest" ecoregion [46]. The area has been experiencing dramatic land use changes in the last decade with the introduction of the teak tree plantations [29,47].

We selected this area because it presents, over short distances, diversely managed teak tree plantations. This area used to be shifting cultivation land. This last decade, teak has gradually replaced most of the fields and spontaneous forest regrowth because of land degradation, lack of labor and the expectation of profitability. Teak timber is valued for its durability and water resistance; it is used for furniture and construction in the rapidly developing city of Luang Prabang. Farmers spontaneously adopted different management practices within their plantations, and we selected four sites corresponding to actual contrasted situations that we intended to test and compare.

Aside from the most common situation which is teak with no understory (TNU), i.e., teak tree plantations where soil is kept bare, often by burning the leaf litter and understory, and which represented our control situation, we considered the three following alternative treatments: teak trees grown with high density understory (THU), teak trees grown with low density and/or periodically pruned understory (TLU), and teak trees grown with broom grass, *Thysanolaena latifolia* (TBG). In each treatment, we set up 1 m² microplots (Figure 2) with six replicates per treatment.
<figure><img src="image_97.png" /><figcaption>**Figure 1.** Study site in Ban Kokngew, Luang Prabang Province, Lao People's Democratic Republic, with location of experimental microplots, treatments, and soil types.</figcaption></figure>

The age of the teak trees in the four treatments varied between 12 years in TBG and 18 years in TLU (Table 1). Elevation above sea level ranged between 316 m in TLU and 358 m in TBG, whereas TNU and THU were both at 325 m. The slopes of the microplots ranged between 39% in TNU and 46% in TBG. The slope difference between TNU and TBG was not considered a limitation for the comparison of treatments in this study since the effect of slope is known to be imperceptible for these slope ranges [31]. The size of the teak tree plantations was 1.87 ha, 3.78 ha, 1.21 ha, and 1.32 ha in TNU, TLU, THU, and TBG, respectively. We installed a rain gauge near the treatments TNU, TLU, and THU, and approximately 500 m from TBG (Figure 1).

<table><caption>**Table 1.** Characteristics of the four experimental sites in 2017, measured on 14 December 2017, in Ban Kokngew, Luang Prabang Province, Lao People's Democratic Republic. TNU: teak with no understory; TLU: teak with low density of understory; THU: teak with high density of understory; TBG: teak with broom grass.</caption><thead><tr><th>Treatment</th><th>Teak Height (m)</th><th>Teak Stem Diameter (cm)</th><th>Teak Cover (%)</th><th>Tree Density (tree·ha<sup>-1</sup>)</th><th>Teak Age (Years)</th><th>Altitude (m)</th><th>Slope (%)</th><th>Latitude (°)</th><th>Longitude (°)</th></tr></thead><tbody><tr><td>TNU</td><td>22</td><td>15.3</td><td>35</td><td>1200</td><td>15</td><td>325</td><td>39.2</td><td>19.8577</td><td>102.21269</td></tr><tr><td>TLU</td><td>22</td><td>17.0</td><td>30</td><td>1000</td><td>18</td><td>316</td><td>40.7</td><td>19.85641</td><td>102.21202</td></tr><tr><td>THU</td><td>20</td><td>15.4</td><td>70</td><td>800</td><td>15</td><td>325</td><td>41.5</td><td>19.85735</td><td>102.21237</td></tr><tr><td>TBG</td><td>20</td><td>16.4</td><td>60</td><td>1000</td><td>12</td><td>358</td><td>45.5</td><td>19.85475</td><td>102.21666</td></tr></tbody></table>
<figure><img src="image_98.png" /><figcaption>**Figure 2.** (a) Sketch of microplot of 1 × 1 m metal frame connected to a bucket through a pipe for surface runoff and sediment collection (source: [31]). (b) Microplot of TNU: teak with no understory. (c) Microplot of TLU: teak with low density of understory. (d) Microplot of THU: teak with high density of understory. (e) Microplot of TBG: teak with broom grass. Microplots were installed on 9th and 10th of May 2017, in Ban Kokngew, Luang Prabang Province, Lao People's Democratic Republic.</figcaption></figure>
## 2.2. Teak Trees and Understory Structure Assessment

We measured teak tree height, canopy cover, and stem diameter at 1.6 m height in 10 × 10 m plots enclosing the microplots. At the time of measurement, on 14 December 2017, i.e., during the dry season, some teak trees were already shedding their leaves (30% and 35% cover in TLU and TNU, respectively) but in other plantations, many teak leaves were still attached (70% and 60% cover in THU and TBG, respectively). Average teak heights ranged from 20 m to 22 m. Average diameter at 1.6 m height was 16 cm. The original planting density had been 2000 tree·ha⁻¹ in TNU and TBG, from which the actual density was obtained by thinning 40% and 50%, respectively. The original planting density had been 1600 tree·ha⁻¹ in THU and only 1200 tree·ha⁻¹ in TLU; here the actual densities were obtained by thinning 50% and 16%, respectively (Table 1).

We considered two categories of understory: one consisting of weed and low vegetation, hereafter referred to as “understory” and another consisting of purposely planted broom grass (Thysanolaena latifolia) monocrop, hereafter referred to as “broom grass”. Both kinds of understories are spontaneous vegetation, but farmers enhance broom grass propagation by cutting. The TBG treatment in this study was made of replanted broom grass. Farmers cut the inflorescences of broom grass to make brooms and regularly prune the grass. We described the structure of the understory in each of the four selected treatments during the 2017 rainy season by combining the use of visual inspections (percentage of cover assessment) and measuring tapes (girth, height). We estimated the mean understory height and cover in a representative area of 18 m² encompassing the three microplot replicates on 4 June and 27 October 2017.

## 2.3. Rainfall Measurements

We measured rainfall by using a rain gauge (Campbell BWS200 equipped with ARG100, 0.2 mm capacity tipping-bucket; Figure 1).

## 2.4. Surface Runoff, Soil Loss, and Soil Surface Features Assessment

We collected surface runoff through a small channel at the lower side of each microplot, connected to a large plastic bucket by a plastic pipe (Figure 2a). The height of the buckets was 0.45 m while their diameter at the bottom was 0.32 m and the diameter at the top was 0.38 m. We emptied the buckets after every major rainfall event, or after a series of two to ten smaller rainfall events. We calculated the runoff coefficient as the ratio between the total surface runoff depth and the total rainfall depth, expressed in percentage. We also computed a cumulated runoff coefficient over the rainy season.

Soil loss was the total weight of sediment collected each time the buckets were emptied. It was measured after flocculation, filtration, and oven dehydration. We calculated suspended sediment concentration by dividing sediment mass by surface runoff depth and we cumulated soil loss over the rainy season.

Soil surface features were assessed at the beginning (4 June) and at the end (27 October) of the rainy season using the method proposed by [48] and extensively used by [23,26,28–33,44,49–54]. We calculated the average of each soil surface feature for each treatment by first calculating the averages of the two measurements dates per treatment, and then by calculating the averages among the six replicates per treatment. Surface features include understory, residues (leaves, branches, seeds), constructions by soil macro-organisms like earthworms and termites, moss and algae, charcoals, free aggregates, free gravel, and three types of crust: structural, erosion, and gravel crusts. In these soils, structural crusts result from the packing of highly stable micro-aggregates [32,50]. Compacted by raindrops and smoothed by surface runoff, this structural crust gradually transforms into an erosion crust characterized by a thin and very compacted smooth plasmic layer [55]. When they include gravels, structural or erosion crusts become a gravel crust [51]. Additionally, we assessed the percentage areas of soil corresponding to pedestal features [5].
## 2.5. Statistical Analysis and Modelling

We performed non-parametric Wilcoxon tests (R version 3.5.3, The R Foundation for Statistical Computing) to compare the distribution between paired groups of the four treatments for three variables (surface runoff depth, suspended sediment concentration, and runoff coefficient).

We conducted a correlation analysis between the measured variables, namely Rc: seasonal runoff coefficient; Sl: seasonal soil loss; and soil surface features expressed in areal percentage: Fa: free aggregates; Fg: free gravel; Tc: total crust; Sc: structural crust; Ec: erosion crust; Gc: gravel crust; Cha: charcoals; Res: residues; Wor: worm casts; Alg: algae; Mos: mosses; Ped: pedestals; Und: understory. To meet the distributional and variance assumptions required for linear statistical models [56], we transformed the variables prior to the analysis; variables expressed in percent were normalized using the arcsine of the square root, which is a classical transformation for percentages, while Sl was scaled with logarithm transformation [31]. The main objective of these transformations was to make each distribution symmetrical. After transformation, we calculated Pearson correlation coefficients and significance levels (XLSTAT Premium version 20.1.1., Addinsoft, Paris, France) in order to test the correlation between variables.

We also calculated partial least squares regression (PLSR) analysis [57] on the measured variables (XLSTAT Premium version 20.1.1.) in order to model Rc and Sl depending on soil surface features: Fa: free aggregates; Fg: free gravel; Tc: total crust; Cha: charcoals; Res: residues; Wor: worm casts; Alg: algae; Mos: mosses; Ped: pedestals; Und: understory. PLSR has the advantage to be little sensitive to multi-collinearity and can be used with datasets where the number of observations is close to the number of variables, or even smaller. The importance of each projected variable is estimated by the variable importance in the projection number (VIP). In order to limit the uncertainty related to the variables that bring little information to the model, and consequently to limit the distortion of the results, we discarded the VIP values below 0.8 [58,59].

# 3. Results

## 3.1. Rainfall

Accumulated rainy season rainfall was 1133 mm from 4 June to 15 October 2017. A total of 22 major rainfall events occurred during the same period (Figure S1). Minimum rainfall depth was 17 mm whereas maximum was 93 mm. Average rainfall depth was 52 mm. About 36% of the major rainfall events occurred in July, which represents 33% of the accumulated rainy season rainfall depth.

## 3.2. Height and Cover of Teak Trees and of Understory

The average height of the teak trees ranged from 20 m in THU and TBG to 22 m in TNU and TLU. Teak cover was 30% and 35% in TLU and TNU, respectively, and was 60% and 70% in TBG and THU, respectively. The average height of understory varied between 0.6 m in TNU and 4 m in TLU. Understory cover varied between 30% in TNU and 90% in THU (Figure 3).

## 3.3. Soil Surface Features and Pedestal Features

The average height of pedestal features ranged from 1.2 cm in TBG to 2.1 cm in TNU (Figure 4a,b). The pedestal features cover was high in TNU (50%) compared to the other treatments (3–6%; Figure 4b).

Figure 4c illustrates the percentage area of total crust (erosion crust, structural crust, and gravel crust), free aggregates, free gravel, charcoals, and residues in each treatment. Soil surface across the treatments excluding TNU shared similar conditions, such as total crust (8.5%), free aggregates (32.5%), and residues (58.5%). On the contrary, TNU exhibited high total crust (82.5%), free gravel (9.3%), and little residues (2.91%). Structural crust accounts for more than 90% of total crust in each treatment and charcoals are negligible (less than 0.03%). Nevertheless, erosion crust was 0.54% in TNU, and negligible in the other three treatments (less than 0.05%).
Figure 5 shows the relationship between the percentage of cover of pedestal features and both surface runoff (logarithm model, R² = 0.82) and soil loss (linear model, R² = 0.91). The percentage of cover of pedestal features exhibits a consistent increase with increasing surface runoff and soil loss.

<figure><img src="image_99.png" /></figure>

**Figure 3.** Percentage of cover (%) by teak trees and understory, and mean height (m) of teak trees and understory in each treatment measured on 14 December 2017, in Ban Kokngew, Luang Prabang Province, Lao People's Democratic Republic. TNU: teak with no understory; TLU: teak with low density of understory; THU: teak with high density of understory; TBG: teak with broom grass.

<figure><img src="image_100.png" /></figure>

**Figure 4.** *Cont.*
<figure><img src="image_101.png" /><figcaption>**Figure 4.** (a) Example of pedestal features circled by red lines. (b) Percentage of cover of pedestal features (%) and pedestal features' height (m; logarithmic scale). (c) Cumulative percentage areas (%) of the soil surface features in 2017 in Ban Kokngew, Luang Prabang Province, Lao People's Democratic Republic. TNU: teak with no understory; TLU: teak with low density of understory; THU: teak with high density of understory; TBG: teak with broom grass.</figcaption></figure>

<figure><img src="image_102.png" /><figcaption>**Figure 5.** (a) Relationship between runoff coefficient (%) and percentage of cover of pedestal features (%). (b) Relationship between soil loss (g·m⁻²) and percentage of cover of pedestal features (%) measured in 2017 in Ban Kokngew, Luang Prabang Province, Lao People's Democratic Republic. TNU: teak with no understory; TLU: teak with low density of understory; THU: teak with high density of understory; TBG: teak with broom grass.</figcaption></figure>

### 3.4. Relationship between Surface Runoff and Soil Loss across Four Treatments

Figure 6 depicts the linear correlation between accumulated surface runoff and accumulated soil loss on log scale during the rainy season in 2017 for all the replicates. Total surface runoff averaged over the rainy season ranged from approximately 170 mm in THU and TBG to 480 mm in TLU, while it reached up to 680 mm in TNU. The average value in the four treatments was 370 mm. TBG and THU
provided the same amount of surface runoff. Soil loss varied between 249 g·m⁻² in TBG to 6012 g·m⁻² in TNU, with an average of 1848 g·m⁻² in all treatments. The highest sediment concentration was 9.60 g·L⁻¹ in TNU, while the other three treatments had an average of 2.34 g·L⁻¹.

<figure><img src="image_103.png" /><figcaption>**Figure 6.** Surface runoff (mm) and soil loss (g·m⁻², logarithmic scale) measured from 4 June to 15 October 2017, in the four treatments with six replicates in Ban Kokngew, Luang Prabang Province, Lao People's Democratic Republic. TNU: teak with no understory; TLU: teak with low density of understory; THU: teak with high density of understory; TBG: teak with broom grass.</figcaption></figure>

### 3.5. Effect of Understory on Soil Loss and Surface Runoff Generation

The Wilcoxon test applied to the cumulative surface runoff, the average suspended sediment concentration, and the runoff coefficient for the four treatments highlighted three significantly (*p*-value > 0.05) different categories of treatments: from the little erosive treatments (THU and TBG) to highly erosive treatments (TNU; Figure 7).

The surface runoff pattern matched the rainfall pattern, while the soil detachment pattern did not perfectly match rainfall pattern (Figure S1). The highest median values of surface runoff and runoff coefficient (44 mm mostly in TNU and about 110% in TBG and TNU, respectively) were observed from 15 July to 15 August and this pattern was less clear for soil loss. We found the maximum median value of soil loss in TNU (551 g·m⁻²). For one replicate of TNU, the soil loss reached 1054 g·m⁻² (8 July).

Figure S2 shows the cumulative surface runoff and cumulative soil loss in relation to cumulative rainfall over the 2017 rainy season for the different treatments. Surface runoff and soil loss in TBG were 242 mm and 381 g·m⁻², respectively. Surface runoff and soil loss in THU were 242 mm and 465 g·m⁻², respectively. Surface runoff and soil loss in TLU were 358 mm and 1115 g·m⁻², respectively. TNU produced more surface runoff (612 mm) and much more soil loss (5455 g·m⁻²) than the other treatments. Hence, the surface runoff in TNU was approximately 2.5-times higher than in THU and TBG, and 1.7-times higher than in TLU. The soil loss in TNU was 13- and 14-times higher than in THU and TBG, respectively. TNU had the sharpest rise of soil loss among all the treatments. Median runoff coefficients for TBG, THU, TLU, and TNU were 23%, 23%, 25%, and 60%, respectively.

### 3.6. Runoff Coefficient and Soil Loss in Relation to Soil Surface Features and Understory Cover

Table 2 shows the relationship between the surface runoff coefficient, soil loss, and soil surface features areal percentages. Erosion crust, charcoals, worm casts, algae, and mosses show weak correlations with the surface runoff coefficient while the other soil surface features and understory cover provided Pearson correlation coefficient (*r*) above 0.7. The features with lower *r* had lower percentages, which may cause poor relation with the surface runoff coefficient. Free aggregates, residues and understory had a strong and negative correlation with surface runoff with *r* between -0.72 and -0.91. Free gravel, total crust, structural crust, gravel crust, and pedestal features cover exhibited a strong correlation with surface runoff (*r* = 0.84, 0.89, 0.89, 0.98, and 0.89, respectively).
We also found a strong correlation between soil loss and surface runoff (r = 0.97). Soil surface features and understory cover exhibited significant inter-correlation except for charcoals, worm casts, algae, and mosses, which were weakly related to other types of soil surface features or understory cover.

<figure><img src="image_104.png" /></figure>

**Figure 7.** Boxplots of (a) cumulative surface runoff (mm), (b) average suspended sediment concentration (g·L⁻¹), and (c) runoff coefficient (%) in each treatment measured from 4 June to 15 October 2017, in Ban Kokngew, Luang Prabang Province, Lao People's Democratic Republic. Each rainfall bar represents the accumulated rainfall over the period prior to the sampling. Each boxplot contains the extreme of the lower whisker (vertical line), the lower hinge (thin line), the median (bold line), the upper hinge (thin line), the extreme of the upper whisker (vertical line), and the outliers (black dots) with *p*-values from Wilcoxon tests between two groups of treatments. The whiskers extend to the most extreme data point, which is no more than 1.5-times the interquartile range from the box. TNU: teak with no understory; TLU: teak with low density of understory; THU: teak with high density of understory; TBG: teak with broom grass. The runoff coefficient is the ratio in percentage between total surface runoff depth and total rainfall depth.
<table><caption>**Table 2.** Correlation (Pearson coefficient) between seasonal runoff coefficient (Rc), seasonal soil loss (Sl), and soil surface features: Fa: free aggregates; Fg: free gravel; Tc: total crust; Sc: structural crust; Ec: erosion crust; Gc: gravel crust; Cha: charcoal; Res: residues; Wor: worm casts; Alg: algae; Mos: mosses; Ped: pedestals; Und: understory. All variables were measured in 2017 in Ban Kokngew, Luang Prabang Province, Lao People's Democratic Republic.</caption><thead><tr><th>Variables</th><th>Rc</th><th>Sl</th><th>Fa</th><th>Fg</th><th>Tc</th><th>Sc</th><th>Ec</th><th>Gc</th><th>Cha</th><th>Res</th><th>Wor</th><th>Alg</th><th>Mos</th><th>Ped</th></tr></thead><tbody><tr><td>Sl</td><td>0.97 ***</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Fa</td><td>-0.72 ***</td><td>-0.77 ***</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Fg</td><td>0.84 ***</td><td>0.86 ***</td><td>-0.72 ***</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Tc</td><td>0.89 ***</td><td>0.90 ***</td><td>-0.87 ***</td><td>0.84 ***</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Sc</td><td>0.89 ***</td><td>0.90 ***</td><td>-0.87 ***</td><td>0.83 ***</td><td>0.99 ***</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Ec</td><td>0.62 **</td><td>0.64 ***</td><td>-0.58 **</td><td>0.44 *</td><td>0.67 ***</td><td>0.66 ***</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Gc</td><td>0.87 ***</td><td>0.90 ***</td><td>-0.78 ***</td><td>0.90 ***</td><td>0.91 ***</td><td>0.90 ***</td><td>0.70 ***</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Cha</td><td>0.48 *</td><td>0.44 *</td><td>-0.39</td><td>0.32</td><td>-0.97 ***</td><td>0.52 **</td><td>0.49 *</td><td>0.43 *</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Res</td><td>-0.87 ***</td><td>-0.89 ***</td><td>0.77 ***</td><td>-0.88 ***</td><td>-0.97 ***</td><td>0.52 **</td><td>-0.66 ***</td><td>-0.94 ***</td><td>-0.53 **</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Wor</td><td>-0.47 *</td><td>-0.42 *</td><td>-0.07</td><td>-0.46 *</td><td>-0.27</td><td>-0.27</td><td>-0.15</td><td>-0.37</td><td>-0.13</td><td>0.39</td><td></td><td></td><td></td><td></td></tr><tr><td>Alg</td><td>0.55 **</td><td>0.53 **</td><td>-0.31</td><td>0.49 *</td><td>0.45 *</td><td>0.45 *</td><td>0.16</td><td>0.44 *</td><td>0.28</td><td>-0.47 *</td><td>-0.23</td><td></td><td></td><td></td></tr><tr><td>Mos</td><td>-0.26</td><td>-0.34</td><td>0.27</td><td>-0.25</td><td>-0.39</td><td>-0.38</td><td>-0.32</td><td>-0.37</td><td>-0.10</td><td>0.42 *</td><td>0.26</td><td>-0.13</td><td></td><td></td></tr><tr><td>Ped</td><td>0.89 ***</td><td>0.91 ***</td><td>-0.80 ***</td><td>0.82 ***</td><td>0.94 ***</td><td>0.94 ***</td><td>0.73 ***</td><td>0.91 ***</td><td>0.54 **</td><td>-0.92 ***</td><td>-0.30</td><td>0.38</td><td>-0.31</td><td></td></tr><tr><td>Und</td><td>-0.91 ***</td><td>-0.93 ***</td><td>0.89 ***</td><td>-0.86 ***</td><td>-0.98 ***</td><td>-0.98 ***</td><td>-0.63 **</td><td>-0.91 ***</td><td>-0.49 *</td><td>0.95 ***</td><td>0.26</td><td>-0.47 *</td><td>0.36</td><td>-0.92 ***</td></tr></tbody></table>

Significance level: *** *p* < 0.0001; ** *p* < 0.001; * *p* < 0.05.
Figure 8a shows PLSR biplot of inter-correlation between runoff coefficient, soil loss, understory, and soil surface features. The runoff coefficient and soil loss were positively correlated with total crust, pedestal features cover, free gravel, charcoals, and algae, but were negatively correlated with understory, residues, and free gravel. Worm casts and mosses had no significant relation with the other variables.

<figure><img src="image_105.png" /><figcaption>**Figure 8.** (a) Partial least squares regression (PLSR) biplot for seasonal surface runoff coefficient (Rc) and seasonal soil loss (Sl) in relation with soil surface features and understory. (b) Variable importance for the projection (VIP) score plot of each variable contributing the most to the models of Rc and Sl. Und: understory; Tc: total crust; Ped: pedestals; Res: residues; Fg: free gravel; Fa: free aggregates; Alg: algae; Cha: charcoals; Wor: worm casts; Mos: mosses. All variables were measured in 2017 in Ban Kokngew, Luang Prabang Province, Lao People's Democratic Republic.</figcaption></figure>

Figure 8b displays the score plot of variable importance in projection (VIP) for each soil surface feature and understory treatment. This plot allows the rapid identification of the variables that
contribute the most to the models of runoff coefficient and soil loss (Figure S3). Understory, total crust, pedestal features cover, and residues, were considered the most important variables and therefore the best predictors in the model. The coefficients of each variable contributing to the models' equation are listed in Table S1.

Figure 9 shows the observed and the modelled soil loss proposed by [31]. The soil loss is a function of the runoff coefficient with the equation $\ln(\text{SI}) = -1.30 + 2.36 \ln(\text{Rc})$, where SI is the soil loss and Rc is the runoff coefficient. The model provided significant statistics ($R^2 = 0.91$, $p < 0.0001$) and is a promising framework for the prediction of soil loss based on the runoff coefficient.

<figure><img src="image_106.png" /><figcaption>**Figure 9.** Observed and modelled seasonal soil loss. Each point represents the seasonal soil loss of the 24 microplots. Values (in g·m⁻²) were predicted using the simple formulation $\ln(\text{SI}) = -1.30 + 2.36 \ln(\text{Rc})$, where Rc is the seasonal runoff coefficient and SI is the seasonal soil loss, as proposed by [31]. Observed Rc and SI were measured from 4 June to 15 October 2017, in Ban Kokngew, Luang Prabang Province, Lao People's Democratic Republic.</figcaption></figure>

# 4. Discussion

## 4.1. Understory Limits Surface Runoff and Soil Erosion

Soil loss is strongly related to surface runoff (Figures 6 and 8, and Table 2). Soil loss increased with increasing surface runoff, especially in TNU, which produced the most surface runoff, over the rainy season (Figure S1). Patin et al. [31] successfully applied the Terrace Erosion and Sediment Transport (TEST) model developed by [43] to estimate soil detachment observed for various land covers in the Houay Pano catchment, northern Lao PDR. They also proposed a simple soil loss model describing the soil loss as a function of the runoff coefficient. We applied this simplified model and also found that seasonal surface runoff coefficient is a reliable predictor of the seasonal soil loss (Figure 9). This result shows that, at the plot scale, surface runoff is statistically the main process responsible for the transfer of solid particles. Our results show that under the studied conditions, sediment production, which is high, is not the limiting factor to soil loss. We found that soil loss had a strong correlation with structural crust and a strong negative correlation with free aggregates and residues. Overall, these findings are in accordance with the findings reported by [30], namely that soil loss increases with an increasing surface runoff coefficient proportionally predicted by structural crust. The percentage of erosion crust is one of the main drivers of soil erosion [30]. At the plot scale, soil detachment firstly depends on
splash [5]. The packing of the soil particles, which is due to the compaction by large drops, leads to crust forming [5]. Soil detachment by splash eventually leads to pedestal features [5]. We found that the structural crust represents in average 93–100% of the total crust, similar to the 91% reported by [30] for mature teak trees. Consistently, the percentage of pedestal features cover (Figure 4a) and total crust (82% of the percentage area of soil surface feature; Figure 4b) was relatively high in TNU. The highest surface runoff and soil loss were observed in TNU, due to the highly crusted area together resulting from less residues cover and the absence of understory cover. Hence, runoff generation occurs on the crusted soil area which restricts infiltration [60,61], depending on the presence of plant residues at soil surface and understory [62,63]. Residues do not only prevent crust formation [62], but also favor the biological activity of organisms which decompose organic matter from the residues, especially macro-faunal communities (also referred to as soil engineers) [64]. This biological process contributes to increase the infiltration rate through the development of soil biological porosity and breakdown of soil crusts resulting from previous inappropriate land use [63,65]. The surface runoff coefficient is a function of structural crust [30] and our results confirmed this finding (Table 2). A similar result was obtained by [61,66,67] who found that a soil crusted area enhances surface runoff and soil erosion. Crusted soil has less influence on soil detachment due to its physical protection [33]. This implies that structural crust is not the main cause of soil erosion but the main driver of surface runoff generation [30]. Once generated, surface runoff carries towards downstream the soil particles detached by splash [5].

Our models based on PLSR analysis confirmed the importance of each soil surface feature, including understory cover, in predicting seasonal surface runoff coefficient and seasonal soil loss at the plot scale (Figure S3 and Table S1). High soil detachment and surface runoff could be alleviated by the physical protection of both understory covers, including broom grass, and plant residues (originating from teak leaves and understory). Soil splash depends on understory cover and on the percentage of cover by residues [5], while surface runoff decreases when understory cover increases [68]. The understory and a thick residue layer reduces the rain splash effect and allows a high rate of infiltration, which then reduces surface runoff and soil erosion [69]. The soil is protected by understory, which dissipates the high raindrops' kinetic energy under teak trees [5]. In the absence of understory, the height of trees plays a decisive role with respect to soil detachment (Figure 3). Tall tree canopy in dense tree plantations intercepts the raindrops. Since teak tree leaves are broad, leaves dramatically increase the diameter of raindrops with throughflow, hence their kinetic energy and erosivity [27,30,70]. Plant residue attenuates the effect of rain splash and reduces the percentage area of total crust [30]. A number of authors have shown that the removal of understory by slashing and burning induces increased surface runoff and soil erosion [52,65,71,72].

Runoff coefficients in the four treatments varied in time and increased in the middle of the rainy season (July 2017, Figure S1), which is explained by antecedent rainfall, i.e., the rain falling before any rainfall event of interest [73]. Rainfall events happen more frequently in the middle of the rainy season. Since there is less time between two subsequent rainfall events to dry the soil, soil moisture increases, and soils get saturated or nearly saturated by water. Consequently, surface runoff generation is enhanced, and the runoff coefficients increase in the middle of the rainy season. Surface runoff in TNU plots were much higher than in all other treatments over the rainy season (Figure S1). The runoff coefficient in TNU was higher than in the other treatments since the very beginning of the season despite the low soil moisture and the low hydraulic conductivity of soil. Overall, the runoff coefficient in TNU was the highest with 59.7% (Figure S2), which is comparable to the values found by [30] in mature teak tree plantations with similar conditions of understory. On the contrary, the runoff coefficients of THU and TBG were around 23% (Figure S2), which is comparable to the overall runoff coefficient calculated for the treatments associating young teak trees with a variety of understory in [30].

The least erosive treatment in this study, i.e., TBG (381 $g·m^{-2}$), produced soil loss about 7.5-times higher than the broom grass treatment in Houay Dou catchment in northern Lao PDR [30]. The difference between our finding and broom grass treatment in [30] can be explained by the lower residue cover
and the greater height of teak trees in our study, and by the fact that broom grass canopy was less dense because it was grown under the shadow of the teak trees and frequently harvested. Our values in TBG treatment were also greater than the values reported by [66,74] in northern Vietnam, with similar understory, topography, and climate, but with different vegetation types, namely mixed tree plantation and fallow, respectively. Similarly, the most erosive treatment in this study, i.e., TNU (5455 g·m⁻²), produced soil loss about 3.5-times higher than the mature teak tree treatment in [30]. The values in TNU were also higher than the values reported by [30,33,66] in the same region with similar understory conditions. Such discrepancies may be related to soil surface conditions or plant and understory characteristics (height, plant density, percentage of cover), which were usually discussed in this study. However values 5.4-times higher were measured by [75] in Costa Rica, where antecedent land use as rangeland most likely greatly and durably reduced soil hydraulic conductivity [72]. Overall, the soil loss in TBG was 14-times less than in TNU. Hence, teak tree plantation owners could divide soil loss by 14 by keeping understory, such as broom grass, within teak tree plantations.

Considering the finding of our study at the plot scale, the farmers would also have saved about 40–50 ton·ha⁻¹ of soil over the 2017 rainy season if they had grown understory under teak trees. However, this result must be qualified since the value cannot be extrapolated at scales larger than the plot scale, such as the catchment scale, because erosion processes at catchment scale are different (suspended sediments deposition or resuspension, gully formation on steep slope). At the catchment scale, a supplementary mitigation measure to trap water and eroded soil particles before entering the river network would be riparian buffers and vegetation filter strips [54,76].

## 4.2. Broom Grass Grown in Teak Tree Plantations: Agronomic Aspects and Ecosystem Services

Tree spacing is a determinant factor of teak tree productivity [65]. In our study, densities of 800–1200 tree·ha⁻¹ [47] (Table 1) were relatively high compared to densities reported for other plantation areas (<800 tree·ha⁻¹) [65,77,78]. Planting trees at lower densities in our study site should be seriously considered as not only would it warrant higher productivity [65], possibly higher than in natural forest (40 tree·ha⁻¹) [77], but it would also allow growing intercropped understory adding economic value to the overall plantation yield.

The introduction of agroforestry into agricultural practices ensures an increased food security by a restored soil fertility for food crops [79] and a sustainable production of wood [64]. Broom grass is neither a food nor a feed crop but similar to agroforestry systems, it is grown under trees. Thus, the management practice may possibly increase the overall productivity (biomass, economic yield [45,47]) of the plot, by making an optimal use of the resources (water and nutrients), that would otherwise not be utilized by a single crop, consumed by plant species at different soil layers depending on their root length [80]. Although it seems to be more productive when grown in full light [30], broom grass can be grown under teak trees, similar to other shade-tolerant crops, such as patchouli [78].

Growing broom grass in teak tree plantations may supply a range of ecosystem services that is not limited to the supply of raw material for brooms and to the prevention of high surface runoff and soil loss in teak tree plantations. At the plot scale, the litter from broom grass leaves may increase top soil organic matter content [64], thus improving soil structure, soil nutrient availability [65], soil carbon sequestration [81], and increasing water infiltration and soil moisture retention [44,63]. Diversifying vegetation strata may diversify the habitats for bird species and other forest-dependent species [77,82,83], thus increasing predator biodiversity [84,85] and reducing the need for chemical inputs (insecticides, herbicides, etc.) in e.g., surrounding annual crop plots. At the catchment scale, favoring water infiltration and reducing surface runoff may mitigate natural disasters (floods and droughts) [65,81] and increase the transfer time of contaminants deposited at the soil surface [86]. Decreasing the sediment supply to the stream network would increase the life span of dam reservoirs [87] and avoid dredging costs [88–90], an issue of particular concern along the Mekong River where the number of dams is dramatically increasing [91,92].
4.3. *Percentage of Cover of Pedestal Features: An Indicator of Soil Erosion*

Both the PLSR model and the simplified model proposed by [31] that we presented in this study
provided promising results to predict soil loss. Nevertheless, such models require input data from
experimental microplots which are rather difficult and time consuming to be implemented, particularly
within the framework of a large-scale approach. It would be more practical to use simple indicators
that are known to be proxies of soil loss, such as soil surface features and/or understory characteristics
(Table 2). On one hand, understory cover varies throughout the year, which causes difficulties in
visual observation and inconsistent estimation. On the other hand, observing soil surface features
(such as total crust, free aggregates, free gravel, residues, etc.) requires some expertise. However,
among soils surface features, pedestal feature cover is a reliable proxy of soil erosion intensity in the
field (Figure 4a) that is easy to assess visually. Hence, lay people could use the proxy of the percentage
of cover of pedestal features to identify the impact of the agricultural management of their land on soil
degradation through surface runoff and soil loss.

**5. Conclusions**

We investigated the impact of different types of understory on surface runoff and soil loss in a teak
tree plantation of the mountainous region of northern Lao PDR. We analyzed the relationship between
understory management, soil surface features (including pedestal features), and surface runoff and soil
loss at the plot scale. This paper clearly demonstrates that teak tree plantations, especially in steep
sloping lands, can be prone to considerable soil loss if not properly managed. Our main findings,
which are graphically synthesized in Figure 10, are that:

* Understory cover acts as an umbrella that protects soil surface from rain splash despite the height of teak trees and the large size of their leaves which contribute to produce raindrops of high kinetic energy. Teak tree plantation owners could divide soil loss by 14 by keeping understory, such as broom grass, within teak tree plantations. Hence, growing understory under teak trees is a mitigation management practice that can be reliably promoted to limit surface runoff and soil erosion.

* Residues from both teak tree leaves and understory not only protect the soil but also enhance the infiltrability of water into the soil. In contrast, the main driver of surface runoff and soil erosion is the percentage of crusted area.

* Understory such as broom grass provides several benefits to the relevant stakeholders in the area, in terms of incomes and ecosystem services. For example, the farmers can sell the brooms made from broom grass.

* The percentage of cover of pedestal features appears as a good indicator of soil erosion that farmers and teak tree plantations owners could easily use to assess the degradation of their land.

In such a context, to minimize surface runoff and soil erosion in steep slope areas such as the
montane regions of south-east Asia, decision makers should, if not legally enforce the maintenance
of understory strata in teak tree plantation, at least recommend the plantations owners to maintain
understory and avoid understory and plant residue layers burning.

At the plot scale, our findings are relevant to farmers concerned about soil loss and soil fertility
in their teak tree plots. At catchment scale, they are relevant to decision makers concerned with
the management of costs (such as the cost of water treatment, and/or infrastructure rehabilitation
such as dam reservoir dredging) resulting from soil loss induced by upslope activities such as tree
plantations and improper agricultural land management. In addition to maintaining the understory
strata, encouraging the use of e.g., riparian zone buffers along the streams [93,94] could also be
recommended to trap soil particles from the cultivated hillslopes and favor runoff infiltration, and thus
ensure the sustainability of the system.
<figure><img src="image_107.png" /><figcaption>**Figure 10.** Graphical abstract. TNU: teak with no understory; TLU: teak with low density of understory; THU: teak with high density of understory; TBG: teak with broom grass.</figcaption></figure>

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/12/9/2327/s1, Figure S1: (a) Cumulated rainfall (mm). (b) boxplots of runoff coefficient (%). (c) boxplot of surface runoff (mm). (d) boxplot of soil loss (g·m⁻²) in each treatment measured from 4 June to 15 October 2017, in Ban Kokngew, Luang Prabang Province, Lao People's Democratic Republic. TNU: teak with no understory; TLU: teak with low density of understory; THU: teak with high density of understory; TBG: teak with broom grass. Each rainfall bar represents the accumulated rainfall over the period previous to the sampling. Each boxplot contains the extreme of the lower whisker (dashed line), the lower hinge (thin line), the median (bold line), the upper hinge (thin line), and the extreme of the upper whisker (dashed line). The whiskers extend to the most extreme data point, which is no more than 1.5-times the interquartile range from the box. Figure S2: (a) Cumulative surface runoff (mm) versus cumulative rainfall and runoff coefficient (%, total surface runoff divided by total rainfall), and (b) cumulative soil loss (g·m⁻²) versus cumulative rainfall (mm), measured from 4 June to 15 October 2017, in Ban Kokngew, Luang Prabang Province, Lao People's Democratic Republic. TNU: teak with no understory; TLU: teak with low density of understory; THU: teak with high density of understory; TBG: teak with broom grass. Figure S3: Observed and modelled seasonal runoff coefficient (Rc; in %) and seasonal soil loss (S1; in g·m⁻²) from partial least squares (PLS) regression method. Observed Rc and S1 were measured from 4 June to 15 October 2017, in Ban Kokngew, Luang Prabang Province, Lao People's Democratic Republic. Table S1: Coefficients of each variables in the runoff coefficient (%) and soil loss (g·m⁻²) partial least squares (PLS) models. SD: standard deviation; Fa: free aggregates; Fg: free gravel; Tc: total crust; Cha: charcoals; Res: residues; Wor: worm casts; Alg: algae; Mos: mosses; Ped: pedestals; Und: understory.

**Author Contributions:** Conceptualization, O.S., B.S. (Bounthan Souksavath), B.L., and O.R.; methodology, A.d.R., C.V., and O.R.; formal analysis, L.S., O.R., and L.B.; data curation, L.S., C.V., P.S., A.d.R., B.S. (Bounsamay Soulileuth), N.S., and O.R.; writing—original draft preparation, L.S.; writing—review and editing, L.B., O.R., A.P., C.V., and A.d.R.; visualization, L.S. and O.R.; supervision, L.B., O.R., and C.O.; funding acquisition, O.R. and A.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the French National Research Agency TecItEasy project, grant number ANR-13-AGRO-0007, and by a scholarship of French government and Ministry of Education, Youth and Sports of Cambodia (L.S.).

**Acknowledgments:** This study was supported by the Luang Prabang Teak Programme (LPTP), the French Institute for Sustainable Development (IRD) through the international joint laboratory LUSES (Dynamic of Land Use Changes and Soil Ecosystem Services), and the Faculty of Agriculture of Lao PDR and Lao students, namely VACHOIMA Sythong and TONGSENG Chong Moua. The authors sincerely thank the M-TROPICS Critical Zone Observatory (https://mtropics.obs-mip.fr/), which belongs to the French Research Infrastructure OZCAR (http://www.ozcar-ri.org/), and the Lao Department of Agricultural Land Management (DALaM) for their support, including granting the permission for field access.

**Conflicts of Interest:** The authors declare no conflict of interest.
References

1. Gerstengarbe, F.W.; Werner, P.C. Precipitation Pattern. In *Encyclopedia of Ecology*; Jørgensen, S.E., Fath, B.D., Eds.; Elsevier Science: Amsterdam, The Netherlands, 2008; pp. 2916–2923. [CrossRef]

2. Descroix, L.; González Barrios, J.L.; Viramontes, D.; Poulenard, J.; Anaya, E.; Esteves, M.; Estrada, J. Gully and sheet erosion on subtropical mountain slopes: Their respective roles and the scale effect. *Catena* **2008**, *72*, 325–339. [CrossRef]

3. Valentin, C.; Agus, F.; Alamban, R.; Boosaner, A.; Bricquet, J.-P.; Chaplot, V.; De Guzman, T.; De Rouw, A.; Janeau, J.-L.; Orange, D. Runoff and sediment losses from 27 upland catchments in Southeast Asia: Impact of rapid land use changes and conservation practices. *Agric. Ecosyst. Environ.* **2008**, *128*, 225–238. [CrossRef]

4. Nandi, A.; Luffman, I. Erosion related changes to physicochemical properties of Ultisols distributed on calcareous sedimentary rocks. *J. Sustain. Dev.* **2012**, *5*, 52. [CrossRef]

5. Valentin, C.; Rajot, J.L. Erosion and Principles of Soil Conservation. In *Soils as a Key Component of the Critical Zone 5: Degradation and Rehabilitation*; Johns Hopkins University Press: Baltimore, MD, USA, 2018; Volume 5, pp. 39–82.

6. Sartori, M.; Philippidis, G.; Ferrari, E.; Borrelli, P.; Lugato, E.; Montanarella, L.; Panagos, P. A linkage between the biophysical and the economic: Assessing the global market impacts of soil erosion. *Land Use Policy* **2019**, *86*, 299–312. [CrossRef]

7. Panagos, P.; Borrelli, P.; Poesen, J.; Ballabio, C.; Lugato, E.; Meusburger, K.; Montanarella, L.; Alewell, C. The new assessment of soil loss by water erosion in Europe. *Environ. Sci. Policy* **2015**, *54*, 438–447. [CrossRef]

8. Pimentel, D.; Burgess, M. Soil erosion threatens food production. *Agriculture* **2013**, *3*, 443–463. [CrossRef]

9. Pimentel, D. Soil erosion: A food and environmental threat. *Environ. Dev. Sustain.* **2006**, *8*, 119–137. [CrossRef]

10. Owens, P.; Batalla, R.; Collins, A.; Gomez, B.; Hicks, D.; Horowitz, A.; Kondolf, G.; Marden, M.; Page, M.; Peacock, D. Fine-grained sediment in river systems: Environmental significance and management issues. *River. Res. Appl.* **2005**, *21*, 693–717. [CrossRef]

11. Yuan, S.; Tang, H.; Xiao, Y.; Xia, Y.; Melching, C.; Li, Z. Phosphorus Contamination of the Surface Sediment at a River Confluence. *J. Hydrol.* **2019**, *573*, 568–580. [CrossRef]

12. Chartin, C.; Evrard, O.; Onda, Y.; Patin, J.; Lefèvre, I.; Ottlé, C.; Ayrault, S.; Lepage, H.; Bonté, P. Tracking the early dispersion of contaminated sediment along rivers draining the Fukushima radioactive pollution plume. *Anthropocene* **2013**, *1*, 23–34. [CrossRef]

13. Turner, J.N.; Brewer, P.A.; Macklin, M.G. Fluvial-controlled metal and As mobilisation, dispersal and storage in the Rio Guadiamar, SW Spain and its implications for long-term contaminant fluxes to the Donana wetlands. *Sci. Total. Environ.* **2008**, *394*, 144–161. [CrossRef] [PubMed]

14. Domagalski, J.L.; Kuivila, K.M. Distributions of pesticides and organic contaminants between water and suspended sediment, San Francisco Bay, California. *Estuaries* **1993**, *16*, 416–426. [CrossRef]

15. Gateuille, D.; Evrard, O.; Lefèvre, I.; Moreau-Guigon, E.; Alliot, F.; Chevreuil, M.; Mouchel, J.M. Mass balance and decontamination times of Polycyclic Aromatic Hydrocarbons in rural nested catchments of an early industrialized region (Seine River basin, France). *Sci. Total. Environ.* **2014**, *470*, 608–617. [CrossRef] [PubMed]

16. Ribolzi, O.; Cuny, J.; Sengsoulichanh, P.; Mousques, C.; Soulileuth, B.; Pierret, A.; Huon, S.; Sengtaheuanghoung, O. Land use and water quality along a Mekong tributary in northern Lao P.D.R. *Environ. Manag.* **2011**, *47*, 291–302. [CrossRef] [PubMed]

17. Boithias, L.; Choisy, M.; Souliyasen, N.; Jourdren, M.; Quet, F.; Buisson, Y.; Thammahacksa, C.; Silvera, N.; Latsachack, K.; Sengtaheuanghoung, O.; et al. Hydrological Regime and Water Shortage as Drivers of the Seasonal Incidence of Diarrheal Diseases in a Tropical Montane Environment. *PLoS. Negl. Trop. Dis.* **2016**, *10*, e0005195. [CrossRef] [PubMed]

18. Kim, M.; Boithias, L.; Cho, K.H.; Sengtaheuanghoung, O.; Ribolzi, O. Modeling the Impact of Land Use Change on Basin-scale Transfer of Fecal Indicator Bacteria: SWAT Model Performance. *J. Environ. Qual.* **2018**, *47*, 1115–1122. [CrossRef]

19. Thothong, W.; Huon, S.; Janeau, J.-L.; Boonsaner, A.; de Rouw, A.; Planchon, O.; Bardoux, G.; Parkpian, P. Impact of land use change and rainfall on sediment and carbon accumulation in a water reservoir of North Thailand. *Agric. Ecosyst. Environ.* **2011**, *140*, 521–533. [CrossRef]
20. Zarris, D.; Vlastara, M.; Panagoulia, D. Sediment Delivery Assessment for a Transboundary Mediterranean Catchment: The Example of Nestos River Catchment. *Water Resour. Manag.* **2011**, *25*, 3785. [CrossRef]

21. Efthimiou, N.; Lykoudi, E.; Panagoulia, D.; Karavitis, C. Assessment of soil susceptibility to erosion using the EPM and RUSLE Models: The case of Venetikos River Catchment. *Global NEST J.* **2016**, *18*, 164–179.

22. Sidle, R.C.; Ziegler, A.D.; Negishi, J.N.; Nik, A.R.; Siew, R.; Turkelboom, F. Erosion processes in steep terrain—Truths, myths, and uncertainties related to forest management in Southeast Asia. *For. Ecol. Manag.* **2006**, *224*, 199–225. [CrossRef]

23. Patin, J.; Mouche, E.; Ribolzi, O.; Chaplot, V.; Sengtahevanghoug, O.; Latsachak, K.O.; Soulileuth, B.; Valentin, C. Analysis of runoff production at the plot scale during a long-term survey of a small agricultural catchment in Lao PDR. *J. Hydrol.* **2012**, *426*, 79–92. [CrossRef]

24. Ziegler, A.D.; Giambelluca, T.W.; Tran, L.T.; Vana, T.T.; Nullet, M.A.; Fox, J.; Vien, T.D.; Pinthong, J.; Maxwell, J.F.; Evett, S. Hydrological consequences of landscape fragmentation in mountainous northern Vietnam: Evidence of accelerated overland flow generation. *J. Hydrol.* **2004**, *287*, 124–146. [CrossRef]

25. Ziegler, A.D.; Giambelluca, T.W.; Nullet, M.A.; Sutherland, R.A.; Tantasarin, C.; Vogler, J.B.; Negishi, J.N. Throughfall in an evergreen-dominated forest stand in northern Thailand: Comparison of mobile and stationary methods. *Agric. For. Meteorol.* **2009**, *149*, 373–384. [CrossRef]

26. Valentin, C. (Ed.) Soil Surface Crusting of Soil and Water Harvesting. In *Soils as a Key Component of the Critical Zone 5: Degradation and Rehabilitation*; John Wiley & Sons: New York, NY, USA, 2018; Volume 5, pp. 21–38.

27. Goebes, P.; Seitz, S.; Kühn, P.; Li, Y.; Niklaus, P.A.; Oheimb, G.V.; Scholten, T. Throughfall kinetic energy in young subtropical forests: Investigation on tree species richness effects and spatial variability. *Agric. For. Meteorol.* **2015**, *213*, 148–159. [CrossRef]

28. Lacombe, G.; Ribolzi, O.; de Rouw, A.; Pierret, A.; Latsachak, K.; Silvera, N.; Pham Dinh, R.; Orange, D.; Janeau, J.L.; Soulileuth, B.; et al. Contradictory hydrological impacts of afforestation in the humid tropics evidenced by long-term field monitoring and simulation modelling. *Hydrol. Earth Syst. Sci.* **2016**, *20*, 2691–2704. [CrossRef]

29. Ribolzi, O.; Evrard, O.; Huon, S.; de Rouw, A.; Silvera, N.; Latsachack, K.O.; Soulileuth, B.; Lefevre, I.; Pierret, A.; Lacombe, G.; et al. From shifting cultivation to teak plantation: Effect on overland flow and sediment yield in a montane tropical catchment. *Sci. Rep.* **2017**, *7*, 3987. [CrossRef]

30. Lacombe, G.; Valentin, C.; Sounyafong, P.; de Rouw, A.; Soulileuth, B.; Silvera, N.; Pierret, A.; Sengtaheuanghoug, O.; Ribolzi, O. Linking crop structure, throughfall, soil surface conditions, runoff and soil detachment: 10 land uses analyzed in Northern Laos. *Sci. Total. Environ.* **2018**, *616*, 1330–1338. [CrossRef]

31. Patin, J.; Mouche, E.; Ribolzi, O.; Sengtahevanghoug, O.; Latsachak, K.O.; Soulileuth, B.; Chaplot, V.; Valentin, C. Effect of land use on interrill erosion in a montane catchment of Northern Laos: An analysis based on a pluri-annual runoff and soil loss database. *J. Hydrol.* **2018**, *563*, 480–494. [CrossRef]

32. Ribolzi, O.; Patin, J.; Bresson, L.M.; Latsachack, K.O.; Mouche, E.; Sengtaheuanghoug, O.; Silvera, N.; Thiébaux, J.P.; Valentin, C. Impact of slope gradient on soil surface features and infiltration on steep slopes in northern Laos. *Geomorphology* **2011**, *127*, 53–63. [CrossRef]

33. Chaplot, V.; Khampaseuth, X.; Valentin, C.; Bissonnais, Y.L. Interrill erosion in the sloping lands of northern Laos subjected to shifting cultivation. *Earth Surf. Process. Landf.* **2007**, *32*, 415–428. [CrossRef]

34. Vannoppen, W.; Vanmaercke, M.; De Baets, S.; Poesen, J. A review of the mechanical effects of plant roots on concentrated flow erosion rates. *Earth Sci. Rev.* **2015**, *150*, 666–678. [CrossRef]

35. Shinohara, Y.; Otani, S.; Kubota, T.; Otsuki, K.; Nanko, K. Effects of plant roots on the soil erosion rate under simulated rainfall with high kinetic energy. *Hydrol. Sci. J.* **2016**, *61*, 2435–2442. [CrossRef]

36. Roosevelt, E. *Land Husbandry: Components and Strategy*; FAO: Rome, Italy, 1996.

37. Ehigiator, O.A.; Anyata, B.U. Effects of land clearing techniques and tillage systems on runoff and soil erosion in a tropical rain forest in Nigeria. *J. Environ. Manag.* **2011**, *92*, 2875–2880. [CrossRef] [PubMed]

38. Cerdà, A.; Rodrigo-Comino, J.; Giménez-Morera, A.; Novara, A.; Pulido, M.; Kapović-Solomun, M.; Keesstra, S.D. Policies can help to apply successful strategies to control soil and water losses. The case of chipped pruned branches (CPB) in Mediterranean citrus plantations. *Land Use Policy* **2018**, *75*, 734–745. [CrossRef]

39. Li, X.; Niu, J.; Xie, B. The effect of leaf litter cover on surface runoff and soil erosion in Northern China. *PLoS ONE* **2014**, *9*, e107789. [CrossRef]
40. Durán Zuazo, V.H.; Martinez, J.R.F.; Pleguezuelo, C.R.R.; Martinez Raya, A.; Rodriguez, B.C. Soil-erosion and runoff prevention by plant covers in a mountainous area (se spain): Implications for sustainable agriculture. *Environmentalist* **2006**, *26*, 309–319. [CrossRef]

41. Poesen, J. Conditions for gully formation in the Belgian loam belt and some ways to control them. In *Proceedings of the of Soil Erosion Protection Measures in Europe*. Proc. EC Workshop, Freising, Germany, 24–26 May 1988; pp. 39–52.

42. Valentin, C.; Ruiz Figueroa, J. Effects of kinetic energy and water application rate on the development of crusts in a fine sandy loam soil using sprinkling irrigation and rainfall simulation. *Micromorphol. Sols* **1987**, *401–408*.

43. Van Dijk, A.; Bruijnzeel, L.; Eisma, E. A methodology to study rain splash and wash processes under natural rainfall. *Hydrol. Process.* **2003**, *17*, 153–167. [CrossRef]

44. Neyret, M.; Robain, H.; de Rouw, A.; Janeau, J.-L.; Durand, T.; Kaewthip, J.; Trisophon, K.; Valentin, C. Higher runoff and soil detachment in rubber tree plantations compared to annual cultivation is mitigated by ground cover in steep mountainous Thailand. *Catena* **2020**, *189*, 104472. [CrossRef]

45. Pachas, A.N.A.; Newby, J.C.; Siphommachan, P.; Sakanphet, S.; Dieters, M.J. Broom grass in Lao PDR: A market chain analysis in Luang Prabang Province. *For. Trees Livelihoods* **2020**, 1–18. [CrossRef]

46. Olson, D.M.; Dinerstein, E.; Wikramanayake, E.D.; Burgess, N.D.; Powell, G.V.N.; Underwood, E.C.; D'amico, J.A.; Itoua, I.; Strand, H.E.; Morrison, J.C.; et al. Terrestrial Ecoregions of the World: A New Map of Life on Earth: A new global map of terrestrial ecoregions provides an innovative tool for conserving biodiversity. *Bioscience* **2001**, *51*, 933–938. [CrossRef]

47. Pachas, A.; Sakanphet, S.; Midgley, S.; Dieters, M. Teak (Tectona grandis) silviculture and research: Applications for smallholders in Lao PDR. *Aust. For.* **2019**, *82*, 94–105. [CrossRef]

48. Casenave, A.; Valentin, C. A runoff capability classification system based on surface features criteria in semi-arid areas of West Africa. *J. Hydrol.* **1992**, *130*, 231–249. [CrossRef]

49. Ribolzi, O.; Lacombe, G.; Pierret, A.; Robain, H.; Sounyafong, P.; de Rouw, A.; Soulileuth, B.; Mouche, E.; Huon, S.; Silvera, N. Interacting land use and soil surface dynamics control groundwater outflow in a montane catchment of the lower Mekong basin. *Agric. Ecosyst. Environ.* **2018**, *268*, 90–102. [CrossRef]

50. Janeau, J.L.; Bricquet, J.P.; Planchon, O.; Valentin, C. Soil crusting and infiltration on steep slopes in northern Thailand. *Eur. J. Soil. Sci.* **2003**, *54*, 543–554. [CrossRef]

51. Valentin, C.; Casenave, A. Infiltration into sealed soils as influenced by gravel cover. *Soil Sci. Soc. Am. J.* **1992**, *56*, 1667–1673. [CrossRef]

52. Lacombe, G.; Ribolzi, O.; de Rouw, A.; Pierret, A.; Latsachak, K.; Silvera, N.; Pham Dinh, R.; Orange, D.; Janeau, J.L.; Soulileuth, B.; et al. Afforestation by natural regeneration or by tree planting: Examples of opposite hydrological impacts evidenced by long-term field monitoring in the humid tropics. *Hydrol. Earth Syst. Sci.* **2015**, *12*, 12615–12648. [CrossRef]

53. Janeau, J.-L.; Gillard, L.-C.; Grellier, S.; Jouquet, P.; Le, T.P.Q.; Luu, T.N.M.; Ngo, Q.A.; Orange, D.; Pham, D.R.; Tran, D.T.; et al. Soil erosion, dissolved organic carbon and nutrient losses under different land use systems in a small catchment in northern Vietnam. *Agric. Water Manag.* **2014**, *146*, 314–323. [CrossRef]

54. Vigiak, O.; Ribolzi, O.; Pierret, A.; Sengtaheuanghoung, O.; Valentin, C. Trapping efficiencies of cultivated and natural riparian vegetation of northern Laos. *J. Environ. Qual.* **2008**, *37*, 889–897. [CrossRef]

55. Valentin, C.; Bresson, L.-M. Morphology, genesis and classification of surface crusts in loamy and sandy soils. *Geoderma* **1992**, *55*, 225–245. [CrossRef]

56. Quinn, G.; Keough, M. *Experimental Design and Data Analysis for Biologists*; Cambridge University Press: Cambridge, UK, 2002.

57. Abdi, H. Partial least squares regression and projection on latent structure regression (PLS Regression). *Wiley Interdiscip. Rev. Comput. Stat.* **2010**, *2*, 97–106. [CrossRef]

58. Zaldívar Santamaría, E.; Molina Dagá, D.; Palacios García, A.T. Statistical Modelization of the Descriptor "Minerality" Based on the Sensory Properties and Chemical Composition of Wine. *Beverages* **2019**, *5*, 66.

59. Wold, S. PLS for multivariate linear modeling. *Chemom. Methods Mol. Des.* **1995**, 195–218.

60. Bu, C.-F.; Wu, S.-F.; Yang, K.-B. Effects of physical soil crusts on infiltration and splash erosion in three typical Chinese soils. *Int. J. Sediment Res.* **2014**, *29*, 491–501. [CrossRef]

61. Fox, D.; Bryan, R.; Fox, C. Changes in pore characteristics with depth for structural crusts. *Geoderma* **2004**, *120*, 109–120. [CrossRef]
62. Neave, M.; Rayburg, S. A field investigation into the effects of progressive rainfall-induced soil seal and crust development on runoff and erosion rates: The impact of surface cover. *Geomorphology* **2007**, *87*, 378–390.
[CrossRef]

63. Rey, F.; Ballais, J.-L.; Marre, A.; Rovéra, G. Rôle de la végétation dans la protection contre l'érosion hydrique de surface. *C. R. Geosci.* **2004**, *336*, 991–998. [CrossRef]

64. Angelsen, A.; Kaimowitz, D. Is agroforestry likely to reduce deforestation. In *Agroforestry and Biodiversity Conservation in Tropical Landscapes*; Schroth, G., Fonseca, G.A.d., Harvey, C.A., Gascon, C., Vasconcelos, H.L., Izac, A.-M.N., Eds.; Island Press: Washington, DC, USA, 2004; pp. 87–106.

65. Blanco, J.A.; Lo, Y.-H. *Forest Ecosystems: More Than Just Trees*; BoD-Books on Demand: Norderstedt, Germany, 2012.

66. Podwojewski, P.; Orange, D.; Jouquet, P.; Valentin, C.; Janeau, J.; Tran, D.T. Land-use impacts on surface runoff and soil detachment within agricultural sloping lands in Northern Vietnam. *Catena* **2008**, *74*, 109–118.
[CrossRef]

67. Kosmas, C.; Danalatos, N.; Cammeraat, L.H.; Chabart, M.; Diamantopoulos, J.; Farand, R.; Gutierrez, L.; Jacob, A.; Marques, H.; Martinez-Fernandez, J. The effect of land use on runoff and soil erosion rates under Mediterranean conditions. *Catena* **1997**, *29*, 45–59. [CrossRef]

68. Nouwakpo, S.K.; Weltz, M.A.; Green, C.H.M.; Arslan, A. Combining 3D data and traditional soil erosion assessment techniques to study the effect of a vegetation cover gradient on hillslope runoff and soil erosion in a semi-arid catchment. *Catena* **2018**, *170*, 129–140. [CrossRef]

69. Calder, I.R.; Hall, R.L.; Prasanna, K. Hydrological impact of Eucalyptus plantation in India. *J. Hydrol.* **1993**, *150*, 635–648. [CrossRef]

70. Geißler, C.; Kühn, P.; Böhnke, M.; Bruelheide, H.; Shi, X.; Scholten, T. Splash erosion potential under tree canopies in subtropical SE China. *Catena* **2012**, *91*, 85–93. [CrossRef]

71. Le, H.T.; Rochelle-Newall, E.; Ribolzi, O.; Janeau, J.L.; Huon, S.; Latsachack, K.; Pommier, T. Land use strongly influences soil organic carbon and bacterial community export in runoff in tropical uplands-. *Land. Degrad. Dev.* **2020**, *31*, 118–132. [CrossRef]

72. Fernández-Moya, J.; Alvarado, A.; Forsythe, W.; Ramírez, L.; Algeet-Abarquero, N.; Marchamalo-Sacristán, M. Soil erosion under teak (*Tectona grandis* Lf) plantations: General patterns, assumptions and controversies. *Catena* **2014**, *123*, 236–242. [CrossRef]

73. Jadidoleslam, N.; Mantilla, R.; Krajewski, W.F.; Goska, R. Investigating the role of antecedent SMAP satellite soil moisture, radar rainfall and MODIS vegetation on runoff production in an agricultural region. *J. Hydrol.* **2019**, *579*, 124210. [CrossRef]

74. Phan Ha, H.A.; Huon, S.; Henry des Tureaux, T.; Orange, D.; Jouquet, P.; Valentin, C.; De Rouw, A.; Tran Duc, T. Impact of fodder cover on runoff and soil erosion at plot scale in a cultivated catchment of North Vietnam. *Geoderma* **2012**, *177*, 8–17. [CrossRef]

75. Santamaría Leandro, F. Evaluación de la Pérdida de Suelo en Plantaciones de Teca, Bajo la Aplicación de Sistemas de Conservación de Suelos en Nicoya, Guanacaste. 1992. Available online: http://www.sidalc.net/cgi-bin/wxis.exe/?IsisScript=rednia.xis&method=post&formato=2&cantidad=1&expresion=mfn=005930 (accessed on 2 February 2019).

76. Cao, X.; Song, C.; Xiao, J.; Zhou, Y. The optimal width and mechanism of riparian buffers for storm water nutrient removal in the Chinese eutrophic Lake Chaohu watershed. *Water* **2018**, *10*, 1489. [CrossRef]

77. Koonkhunthod, N.; Sakurai, K.; Tanaka, S. Composition and diversity of woody regeneration in a 37-year-old teak (*Tectona grandis* L.) plantation in Northern Thailand. *For. Ecol. Manag.* **2007**, *247*, 246–254. [CrossRef]

78. Kumar, D.; Bijalwan, A.; Kalra, A.; Dobriyal, M.J. Effect of shade and organic manure on growth and yield of patchouli [Pogostemon cablin (blanco) benth.] under teak (*Tectona grandis* Lf) based agroforestry system. *Indian For. *2016*, *142*, 1121–1129.

79. Schroeder, P. Agroforestry systems: Integrated land use to store and conserve carbon. *Clim. Res.* **1993**, *3*, 53–60. [CrossRef]

80. Chitra-Tarak, R.; Ruiz, L.; Dattaraja, H.S.; Kumar, M.M.; Riotte, J.; Suresh, H.S.; McMahon, S.M.; Sukumar, R. The roots of the drought: Hydrology and water uptake strategies mediate forest-wide demographic response to precipitation. *J. Ecol.* **2018**, *106*, 1495–1507. [CrossRef]

81. Balmford, A.; Whitten, T. Who should pay for tropical conservation, and how could the costs be met? *Oryx* **2003**, *37*, 238–250. [CrossRef]
82. Harvey, C.A.; Villalobos, J.A.G. Agroforestry systems conserve species-rich but modified assemblages of tropical birds and bats. *Biodivers. Conserv.* **2007**, *16*, 2257–2292. [CrossRef]

83. Wunderle, J.M., Jr. The role of animal seed dispersal in accelerating native forest regeneration on degraded tropical lands. *For. Ecol. Manag.* **1997**, *99*, 223–235. [CrossRef]

84. Imron, M.; Tantaryzard, M.; Satria, R.; Maulana, I.; Pudyatmoko, S. Understory Avian Community In A Teak Forest Of Cepu, Central Java. *J. Trop. For. Sci.* **2018**, *30*, 509–518. [CrossRef]

85. Perrin, R. Pest management in multiple cropping systems. *Agro-Ecosystems* **1976**, *3*, 93–118. [CrossRef]

86. Monaghan, R.; Smith, L.; Muirhead, R. Pathways of contaminant transfers to water from an artificially-drained soil under intensive grazing by dairy cows. *Agric. Ecosyst. Environ.* **2016**, *220*, 76–88. [CrossRef]

87. Annandale, G.W. Reservoir sedimentation. In *Encyclopedia of Hydrological Sciences*; Wiley: New York, NY, USA, 2006. [CrossRef]

88. Boithias, L.; Terrado, M.; Corominas, L.; Ziv, G.; Kumar, V.; Marques, M.; Schuhmacher, M.; Acuna, V. Analysis of the uncertainty in the monetary valuation of ecosystem services—A case study at the river basin scale. *Sci. Total. Environ.* **2016**, *543*, 683–690. [CrossRef]

89. Crowder, B.M. Economic costs of reservoir sedimentation: A regional approach to estimating cropland erosion damage. *J. Soil Water Conserv.* **1987**, *42*, 194–197.

90. McHenry, J.R. RESERVOIR SEDIMENTATION 1. *J. Am. Water Resour. Assoc.* **1974**, *10*, 329–337. [CrossRef]

91. Dang, T.D.; Cochrane, T.A.; Arias, M.E. Quantifying suspended sediment dynamics in mega deltas using remote sensing data: A case study of the Mekong floodplains. *Int. J. Appl. Earth Obs. Geoinf.* **2018**, *68*, 105–115. [CrossRef]

92. Arias, M.E.; Cochrane, T.A.; Kummu, M.; Lauri, H.; Holtgrieve, G.W.; Koponen, J.; Piman, T. Impacts of hydropower and climate change on drivers of ecological productivity of Southeast Asia’s most important wetland. *Ecol. Model.* **2014**, *272*, 252–263. [CrossRef]

93. Bhat, S.A.; Dar, M.U.D.; Meena, R.S. Soil erosion and management strategies. In *Sustainable Management of Soil and Environment*; Springer: New York, NY, USA, 2019; pp. 73–122.

94. Ahmad, N.S.B.N.; Mustafa, F.B.; Gideon, D. A systematic review of soil erosion control practices on the agricultural land in Asia. *Int. Soil Water Conserv. Res.* **2020**. [CrossRef]

<figure><img src="image_108.png" /></figure>

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
<figure><img src="image_109.png" /></figure>
<figure><img src="image_110.png" /></figure>

Article

# Comparison of the MUSLE Model and Two Years of Solid Transport Measurement, in the Bouregreg Basin, and Impact on the Sedimentation in the Sidi Mohamed Ben Abdellah Reservoir, Morocco

Mohamed Abdellah Ezzaouini ¹,²,*, Gil Mahé ³, Ilias Kacimi ¹ and Abdelaziz Zerouali ²

¹ Geoscience, Water and Environment Laboratory, Faculty of Sciences, Mohammed V University, Avenue Ibn Batouta, Rabat 10100, Morocco; iliaskacimi@yahoo.fr

² Basin of the Bouregreg and Chaouia Agency, Benslimane BP 262, Morocco; azlso@yahoo.fr

³ Laboratoire HydroSciences, Université Montpellier, IRD, 34090 Montpellier, France; gilmahe@hotmail.com

* Correspondence: ezzaouini70@yahoo.fr; Tel.: +212-661-110-897

Received: 18 May 2020; Accepted: 15 June 2020; Published: 1 July 2020

**Abstract:** The evaluation and quantification of solids transport in Morocco often uses the Universal Soil Loss Model (USLE) and the revised version RUSLE, which presents a calibration difficulty. In this study, we apply the MUSLE model to predict solid transport, for the first time on a large river basin in the Kingdom, calibrated by two years of solid transport measurements on four main gauging stations at the entrance of the Sidi Mohamed Ben Abdellah dam. The application of the MUSLE on the basin demonstrated relatively small differences between the measured values and those expected for the calibrated version, these differences are, for the non-calibrated version, +5% and +102% for the years 2016/2017 and 2017/2018 respectively, and between -33% and +34% for the calibrated version. Besides, the measured and modeled volumes that do not exceed $1.78 \times 10^6 \text{ m}^3/\text{year}$ remain well below the dam's siltation rate of $9.49 \times 10^6 \text{ m}^3/\text{year}$, which means that only 18% of the dam's sediment comes from upstream. This seems very low because it is calculated from only two years. The main hypothesis that we can formulate is that the sediments of the dam most probably comes from the erosion of its banks.

**Keywords:** modeling; MUSLE; erosion; solid transport; dam; Bouregreg; Morocco

## 1. Introduction

Erosion is a natural phenomenon that reduces the capacity of dam reservoirs around the world. The natural erosive process is aggravated by anthropogenic activities including pastoral activity [1], deforestation [2–4], and climate change [5] with the advent of periods of heavy rainfall and increasingly frequent dry periods. This phenomenon constitutes a major challenge for water resource management at the scale of the Bouregreg basin [6,7] in northern Morocco.

The Sidi Mohamed Ben Abdellah (SMBA) dam, commissioned in 1974 and raised in 2007, is intended solely to supply drinking water to the coastal area between Rabat and Casablanca, which represents nearly eight million inhabitants. It has a relatively low silting rate compared to other dams in the Kingdom [8]. However, it has lost $132 \text{ Mm}^3$ since its commissioning of which 58% of this volume was lost before rising the dam height, with this loss constituting a very significant reduction in its storage capacity. Given the magnitude of this situation, the modeling of soil losses in the basin aims at achieving the following objectives:

* analyzing the biophysical environment;
* describing and evaluating the erosive processes affecting the Bouregreg Basin, and the solid transport by the main tributaries to the SMBA dam;

* and identifying the priority areas contributing to siltation, in order to better guide spatial planning actions.

Universal Soil Loss Equation (USLE) was first established in USA to model erosion in small agricultural catchment [9]. It is based on several parameters linked to climate, soil cover and properties, topography, and human activities. This equation has been modified and adapted several times. The MUSLE model includes the use of water flow rates [10]. Despite the difficulties encountered in calibrating and adapting the Universal Land Loss Model (USLE) to conditions in Morocco [11], most studies on watershed management in Morocco [12–15] and bordering Mediterranean regions, in particular in Algeria and Tunisia [15–17], continue to use USLE, and the revised version RUSLE [18–21], or an event normalized plot soil loss estimated by a modified USLE model—USLE-MM—as in Italy [22], most often for small basins of much less than 5000 km². As USLE [23] and RUSLE [24] were developed for the rough assessment of annual land loss at the scale of small plots, their application to large areas leads to rather large errors [25,26]. However, their accuracy increases when coupled with hydrological models [27]. To overcome the difficulties in assessing the accuracy of using a simple erosion equation like USLE, Alewell et al. [28] recommend to strengthen and extend measurement and monitoring programs to build up validation data sets.

Thus, Williams [10] developed a modified version of the USLE (MUSLE) that takes into account the flow load at the outlet by taking into account the biophysical characteristics of the watershed. This model has already been applied to micro-watersheds [29,30] and gave very reliable results compared to measurements. Indeed, in the Sidi Sbaa basin in Morocco, the deviation of the results compared to the MUSLE model was –4% by underestimating the solid inputs [31]. Samaras and Koutitas et al. [32] use MUSLE with SWAT to simulate the potential impact of land cover change on sediment yields to the sea in Greece, but with no observed validation data; while other authors like Fang [33] use the WaTEM/SEDEM model, which includes the RUSLE formula, to estimate erosion, but always without observations to compare. Only a few studies compare the erosion rates with the three USLE formulas. In Maghreb, only one study by Djoukbala et al. [34] compared them on the small basin of 384 km² in the north of Algeria, with erosion rates quite similar between the three methods, but they were slightly superior in the case MUSLE. Unfortunately, they could not compare their results with observation data, which does not allow an assessment of the validity of the erosion rates produced in regard to real natural processes.

## 2. Materials and Methods

### 2.1. Study Area

The Bouregreg Basin, located in central western Morocco, covers a total area of approximately 10,130 km² at its mouth in the north of the city of Rabat (Figure 1). The main rivers of the basin are the Bouregreg, the Grou, and the Korifla and its tributary, the Machraa.
<figure><img src="image_111.png" /><figcaption>Figure 1. Situation of the Bouregreg Basin on the north Morocco Atlantic coast.</figcaption></figure>
The climate in the study area is Mediterranean with oceanic influence, with an average annual rainfall over the basin varying from 450 mm in Rabat in the north-west, to nearly 750 mm in the mountainous area in the south-east.

Rain events generate water volume of 680 × 10⁶ m³/year, i.e., an annual mean of 22 m³/s. It is regulated by the capacity of the SMBA dam, located downstream of the confluence between the Bouregreg, Korifla/Machraa, and Grou rivers. Table 1 summarizes the characteristics of the SMBA dam.

<table><caption>Table 1. Sidi Mohamed Ben Abdellah (SMBA)dam data.</caption><thead><tr><th></th><th>Watershed Area (km²)</th><th>Initial Dam Capacity (10⁶ m³)</th><th>Inter-Annual Water Resources (10⁶ m³)</th><th>Opening Date</th><th>Number of Bathymetric Measurements</th></tr></thead><tbody><tr><td>SMBA at its construction</td><td>9800</td><td>508.6</td><td>680</td><td>1974</td><td>5</td></tr><tr><td>SMBA after raising its dike</td><td></td><td>974.0</td><td></td><td>2007</td><td>3</td></tr></tbody></table>

## 2.2. Basic Data

### 2.2.1. Discharges Measurements

The hydrometric network of the Bouregreg basin upstream of the SMBA dam is composed of three major tributaries: Bouregreg, Grou, and Korifla (including its tributary the Machraa). Four hydrological stations located on these tributaries at the entrance of the dam's reservoir were chosen to carry out the measurements of solid transport. These hydrological stations control a basin of about 8521 km², i.e., 87% of the catchment area of the SMBA dam. They also have human and material means to ensure the measurement of flows, rainfall, and the concentration of suspended solids. Table 2 gives the characteristics of the four hydrological stations.

<table><caption>Table 2. Characteristics of the four hydrological stations studied.</caption><thead><tr><th rowspan="2">Hydrological Station</th><th rowspan="2">Name of River</th><th rowspan="2">Date of Commissioning</th><th rowspan="2">Averaging Period</th><th rowspan="2">N°IRE ABHBC Code</th><th rowspan="2">Watershed Area (km²)</th><th colspan="3">Lambert Coordinates</th></tr><tr><th>X</th><th>Y</th><th>Z</th></tr></thead><tbody><tr><td>Aguibat Zear</td><td>Bouregreg</td><td>1975</td><td>1975/2018</td><td>3118/13</td><td>3681</td><td>394.500</td><td>368.150</td><td>90</td></tr><tr><td>Ras Fathia</td><td>Grou</td><td>1975</td><td>1975/2018</td><td>989/20</td><td>3485</td><td>394.250</td><td>351.800</td><td>100</td></tr><tr><td>Ain Loudah Sidi Mohamed Cherif</td><td>Korifla</td><td>1971</td><td>1971/2018</td><td>2673/20</td><td>699</td><td>373.750</td><td>329.150</td><td>175</td></tr><tr><td></td><td>Machraa</td><td>1971</td><td>1971/2018</td><td>2674/21</td><td>656</td><td>385.850</td><td>328.200</td><td>270</td></tr></tbody></table>

Figure 2 shows that the two years, 2016/2017 and 2017/2018, are dry years which did not record significant water inflows.

Examination of the historical measurement data from the four observation stations shows that 80% of the inflows to the dam are recorded between October and May (Figure 3). The Bouregreg and Grou rivers contribute more than 90% of the inflows generated in the basin.

Table 3 shows the maximum flows, volumes, and number of flood events recorded during the 2016/2017 and 2017/2018 hydrological years at the four hydrological stations studied. The number of recorded events is between six and seventeen, depending on the observation station and the year.
<figure><img src="image_112.png" /><figcaption>Figure 2. Annual discharges at the 4 hydrological stations studied.</figcaption></figure>

<figure><img src="image_113.png" /><figcaption>Figure 3. Mean monthly discharges at the hydrological stations studied.</figcaption></figure>

<table><caption>Table 3. Flood events for the four hydrological stations studied.</caption><thead><tr><th>Hydrological Station</th><th>Hydrological Year</th><th>Maximum Discharge Recorded Per Year (m³/s)</th><th>Date of Flood Event</th><th>Volume (10⁶ m³)</th><th>Number of Flood Events</th></tr></thead><tbody><tr><td rowspan="2">Aguibat Zear</td><td>2016/2017</td><td>263.4</td><td>24/02/2017 at 23h00</td><td>52</td><td>13</td></tr><tr><td>2017/2018</td><td>265.9</td><td>07/03/2018 at 14h00</td><td>77</td><td>14</td></tr><tr><td rowspan="2">Ras Fathia</td><td>2016/2017</td><td>256.4</td><td>25/02/2017 at 07h30</td><td>33</td><td>11</td></tr><tr><td>2017/2018</td><td>423.8</td><td>07/03/2018 at 09h00</td><td>59</td><td>10</td></tr><tr><td rowspan="2">Ain Loudah</td><td>2016/2017</td><td>24.5</td><td>13/10/2016 at 11h00</td><td>0.7</td><td>10</td></tr><tr><td>2017/2018</td><td>198.7</td><td>24/04/2018 at 23h00</td><td>7.6</td><td>17</td></tr><tr><td>Sidi Mohamed</td><td>2016/2017</td><td>62.6</td><td>13/10/2016 at 07h00</td><td>1.8</td><td>6</td></tr><tr><td>Cherif</td><td>2017/2018</td><td>48.0</td><td>11/12/2017 at 17h00</td><td>1.28</td><td>6</td></tr></tbody></table>

### 2.2.2. Concentration of Suspended Solids (CSS)

Observers at the four hydrological stations studied, who are contracted by the Hydraulic Basin of the Bouregreg and Chaouia Agency (ABHBC), carry out daily sampling during low-flow periods and hourly sampling during periods of flooding. At each instantaneous sampling, the date, time, and scale
rating are noted on the specimen bottles catalogued. The samples are then analyzed in the laboratory and filtered under a vacuum using filtering membranes (0.45µm). The available data for measuring the concentration of suspended solids cover two hydrological years, 2016/2017 and 2017/2018. The number of measurement samples is presented in Table 4.

<table><caption>Table 4. Number of concentration of suspended solids (CSS) samples and floods events per station.</caption><thead><tr><th rowspan="2">Name of Hydrological Station</th><th rowspan="2">River</th><th rowspan="2">Watershed Area (km²)</th><th colspan="2">2016/2017 and 2016/2018</th></tr><tr><th>Number of CSS Samples</th><th>Number of Flood Events</th></tr></thead><tbody><tr><td>AguibatZear</td><td>Bouregreg</td><td>3681</td><td>727</td><td>27</td></tr><tr><td>Ras Fathia</td><td>Grou</td><td>3485</td><td>636</td><td>21</td></tr><tr><td>Ain Loudah</td><td>Korifla</td><td>699</td><td>233</td><td>27</td></tr><tr><td>Sidi Mohamed Cherif</td><td>Machraa</td><td>656</td><td>789</td><td>12</td></tr></tbody></table>

### 2.2.3. Bathymetric Data

In order to draw up an inventory of the silting of the SMBA dam reservoir, the bathymetric surveys carried out by the Moroccan Directorate of Water Research and Planning were collected and analyzed. A total of seven bathymetric surveys were collected, the oldest dating back to 1974, while the most recent were carried out in 2013. Analysis of the bathymetric data shows that the silting up of the SMBA dam reservoir was of the order of $2.65 \times 10^6 \text{ m}^3/\text{year}$, i.e., a specific degradation of $270.4 \text{ m}^3/\text{km}^2/\text{year}$ before raising. After raising, the silting increased to $9.49 \times 10^6 \text{ m}^3/\text{year}$, i.e. a specific degradation of $968.37 \text{ m}^3/\text{km}^2/\text{year}$, which increased the silting rate by 400% [4]. However, the data is subject to doubt for values in the early 2000s [35]. The difference between the beginning and the end of the chronicle is of reliable quality. Therefore, the total silting of the SMBA dam reached $132 \times 10^6 \text{ m}^3$, i.e., an average silting rate of $3.7 \times 10^6 \text{ m}^3/\text{year}$ since its commissioning. Subsequently, the dam has exceeded its dead unit, which is sized for $100 \times 10^6 \text{ m}^3$ and has lost $32 \times 10^6 \text{ m}^3$ of its useful reserve since its commissioning. Table 5 summarizes the silting status of the dam.

<table><caption>Table 5. Summary of SMBA dam siltation calculation results since its impoundment.</caption><thead><tr><th>Name of the Dam</th><th>Initial Dam Capacity (10<sup>6</sup> m<sup>3</sup>)</th><th>Total Siltation (10<sup>6</sup> m<sup>3</sup>)</th><th>Number of Years</th><th>Rate of Siltation (10<sup>6</sup> m<sup>3</sup>)</th><th>Dead-unit Volume at Dam Building (10<sup>6</sup>m<sup>3</sup>)</th><th>Lost Volume (%)</th><th>Current Capacity (10<sup>6</sup> m<sup>3</sup>)</th></tr></thead><tbody><tr><td>SMBA before raising</td><td>508.60</td><td>76.88</td><td>29</td><td>2.65</td><td rowspan="2">100.00</td><td>15</td><td>431.72</td></tr><tr><td>SMBA after raising</td><td>974.79</td><td>55.28</td><td>6</td><td>9.49</td><td>6</td><td>919.51</td></tr></tbody></table>

The silting of the dam reservoir represents a real threat to the sustainability of the mobilization of surface water resources in the Bouregreg basin to satisfy required needs. The regulation of the SMBA dam, prior to its raising, was done on a seasonal basis (capacity lower than the annual inflow). In other words, the dam had to discharge the surplus inflows most often from flood spillway or bottom discharge. This technique of management favored the elimination of solid deposits and thus a reduction in the siltation rate of the dam. After the dam was raised and consequently the water capacity of the reservoir increased, the SMBA dam moved to multi-year regulation (capacity greater than annual inflows). By using the restitution devices of the dam, this method of regulation favors storage to the detriment of evacuation, which accentuates the silting rate. Figure 4 illustrates the evolution over time of the normal capacity of the SMBA dam.
<figure><img src="image_114.png" /><figcaption>Figure 4. Evolution of the normal capacity of the SMBA dam reservoir, in millions of m³.</figcaption></figure>

## 2.3. Methods

### 2.3.1. Model Selection

In this study we apply the Williams [10] model based on the Modified Universal Soil Loss Equation (MUSLE), integrated in the ArcGIS Geographic Information System for the determination of the soil loss potential at the level of the Bouregreg watershed to the four hydrological stations located immediately upstream of the SMBA dam, as previously done by Khali Issa et al. [36] in another region of Morocco in the North of the country, over a very small basin of 38 km².

This model evaluates the average annual rate of erosion at the outlet of the basin. It uses hydrological parameters, measured at the four hydrological stations, taking into account biophysical characteristics. Thus, the model equation is as follows:

$$A = a(Q_{max} \times V_t)^b K \times LS \times C \times P, \quad (1)$$

where A: amount of sediment produced at the outlet in tons, a and b: in this study, we used the scale factor values of the Sidi Sbaa micro-basin [30] (a = 11.8 and b = 0.56), Qmax: maximum flow rate in m³/s, Vt: total volume of runoff water in m³, K: average soil erodibility (mg MJ-1mm-1), LS: average topographic factor, C: average vegetation cover factor, P: average cultural practices and amenities factor.

### 2.3.2. Methods Selection

We adopted the suite of methods and operations explained in the flowchart (Figure 5) below to assess siltation rates at hydrological stations upstream of the SMBA dam [37].
<figure><img src="image_115.png" /><figcaption>Figure 5. Methodology adopted for the assessment of erosion in the Bouregreg basin.</figcaption></figure>

## 3. Results

### 3.1. *Rainfall and Hydrometric Analysis*

Solid transports are calculated at the four hydrological stations concerned in the hydrological years 2016/2017 and 2017/2018. These two years were respectively dry and wet in terms of rainfall. The average rainfall recorded at all the rainfall stations in the basin reached 354 mm in 2016/2017 and 496 mm in 2017/2018. Thus, the rainfall differences recorded in relation to the arithmetic mean of the data from the rainfall stations located in the Bouregreg basin varied respectively by −10% and +26% since the commissioning dates of the stations, based on data available to study from the ABHBC.

Table 6 summarizes the rainfall variations recorded in relation to the average rainfall of the stations in the Bouregreg basin.

<table><caption>Table 6. Rainfall for the 2016/2017 and 2017/2018 water years.</caption><thead><tr><th rowspan="2">Rain Gauge Station</th><th rowspan="2">Basin</th><th colspan="2">Cumulative Rainfall (mm)</th><th rowspan="2">Averaging Period</th><th rowspan="2">Mean Year (mm)</th><th colspan="2">Deviation from the Mean (%)</th></tr><tr><th>2016/2017</th><th>2017/2018</th><th>2016/2017</th><th>2017/2018</th></tr></thead><tbody><tr><td>LallaChafia</td><td>Bouregreg</td><td>409</td><td>569</td><td>1972/2019</td><td>486</td><td>-16</td><td>17</td></tr><tr><td>AguibatZear</td><td>Bouregreg</td><td>482</td><td>578</td><td>1975/2019</td><td>429</td><td>12</td><td>35</td></tr><tr><td>Sidi Jabeur</td><td>Grou</td><td>260</td><td>432</td><td>1971/2019</td><td>317</td><td>-18</td><td>36</td></tr><tr><td>Tsalat</td><td>Bouregreg</td><td>454</td><td>622</td><td>1977/2019</td><td>469</td><td>-3</td><td>33</td></tr><tr><td>Roumani</td><td>Machraa</td><td>261</td><td>431</td><td>1933/2019</td><td>342</td><td>-24</td><td>26</td></tr><tr><td>Ras Fathia</td><td>Grou</td><td>290</td><td>505</td><td>1975/2019</td><td>387</td><td>-25</td><td>30</td></tr><tr><td>S. M Cherif</td><td>Machraa</td><td>290</td><td>495</td><td>1971/2019</td><td>367</td><td>-21</td><td>35</td></tr><tr><td>Barrage SMBBA</td><td>Bouregreg</td><td>518</td><td>479</td><td>1984/2019</td><td>469</td><td>10</td><td>2</td></tr><tr><td>OuljatHaboub</td><td>Grou</td><td>321</td><td>350</td><td>1972/2019</td><td>311</td><td>3</td><td>12</td></tr><tr><td>Ain Loudah</td><td>Krofla</td><td>258</td><td>498</td><td>1971/2019</td><td>347</td><td>-26</td><td>43</td></tr><tr><td>Total</td><td></td><td>354</td><td>496</td><td></td><td>392</td><td>-10</td><td>26</td></tr></tbody></table>
With the exception of the Korifla sub-basin at the right of the Ain Loudah station which recorded a surplus of nearly 24% during 2017/2018, the contributions were in deficit on the rest of the sub-basins during the two hydrological years with more pronounced deficits during the dry year 2016/2017. Table 7 summarizes the inflows recorded at the four stations of the study.

<table><caption>Table 7. Water supplies to the four hydrological stations studied.</caption><thead><tr><th rowspan="2">Hydrological Station</th><th rowspan="2">River</th><th colspan="3">Water Supply (Mm³)</th><th colspan="2">Deviation from the Mean (%)</th></tr><tr><th>2016/2017</th><th>2017/2018</th><th>Averaging Period</th><th>Annual Average (mm)</th><th>2016/2017</th><th>2017/2018</th></tr></thead><tbody><tr><td>AguibatZear</td><td>Bouregreg</td><td>122.2</td><td>221.5</td><td>1975/2019</td><td>227.6</td><td>-46</td><td>-3</td></tr><tr><td>Ras Fathia</td><td>Grou</td><td>95.0</td><td>162.8</td><td>1975/2019</td><td>194.2</td><td>-51</td><td>-16</td></tr><tr><td>Ain Loudah Sidi</td><td>Korifla</td><td>13.3</td><td>23.5</td><td>1971/2019</td><td>18.9</td><td>-29</td><td>+24</td></tr><tr><td>Mohammed Cherif</td><td>Machraa</td><td>3.7</td><td>12.9</td><td>1971/2019</td><td>16.7</td><td>-77</td><td>-23</td></tr></tbody></table>

## 3.2. Analysis of the Biophysical Environment

The analysis of the biophysical environment consists of determining the average factors used by the model for each sub-basin. These factors, that affect soil erosion, are: soil type, topography, land use and cropping practices, and erosion control facilities [38].

### 3.2.1. Soil Erodibility Factor (K)

Erodibility is defined as the degree to which soils are resistant to erosion. The factors that have a major influence on the response of soils to erosion, namely the detachment and transport of particles by rain and runoff, are texture, structure, organic matter, and permeability. The methodology for estimating RUSLE K has been applied, and is written as follows [39]:

$$K = \frac{2.1 \times 10^{-4}(12 - \text{MO})M^{1.14} + 3.25(S - 2) + 2.5(P - 3)}{10}, \quad (2)$$

with K: soil erodibility expressed in t.ha.h/ha.MJ.mm (tonne. hectare. hour/hectare. mega joule. millimeter); MO: percentage of organic matter; M: textural term % fine sand + % silt; S: structure class code 1 to 4, with 1 fragmented structure and 4 coarse structure, soil structure affects both landslide susceptibility and infiltration, the profiles described on the Bouregreg have a subangular polyhedral structure and fall under class (3); and P: the permeability code (1 to 6), its value can be inferred indirectly from the organic matter content by calculating the infiltration given by the Equation [40]:

$$Y = 3.53 \times X + 2.08, \quad (3)$$

with: Y = infiltration in cm/h, X = organic matter in %.

Thereafter, the soil erodibility factor K will be calculated using the Harmonized World Soil Database (HWSD), developed by the Food and Agriculture Organization of the United Nations (FAO) (http://www.fao.org/soils-portal/fr/). This database gives the distribution of silt, sand, and clay soil compositions by soil type [41].

Soil is composed of organic and mineral matter. Its texture is determined by the size of the soil particles and their respective quantities. There are three categories of particles that determine soil texture: sand, silt, and clay. They are distinguished by particle diameter: sand 0.05 mm to 2 mm, silt 0.002 mm to 0.05 mm, clay 0.002 mm and less. The modified illustrated triangular graph [42] was used to determine soil texture classification. Soil texture is classified according to the percentage of silt, sand, and clay.
Once the textures have been determined, it is possible to establish the correspondence between the standard texture and the K-factor [43]. These values are given in tons/ha and ton/acre (US system). Although this methodology provides an approximation in the calculation of the K-factor, it has the advantage of lending itself to the constraints imposed by the study area. The rate of organic matter at the level of each watershed for each texture value is calculated by converting organic carbon into organic matter. The conversion factor of 1.724 is commonly used to convert the organic carbon content of a soil sample to organic matter. The conversion factor is old and has survived the test of time and modern analytical methods. According to authors [44], this conventional factor is attributed to the 19th century authors Van Bemmelen [45], Wolff [46], or even Sprengel [47]. It is based on well-established and very old studies showing that soil organic matter contains 58% carbon. Since the C/MO ratio would be equal to 0.58, the MO/C ratio would be equal to 1.724.

Thus, the calculation of the organic matter is carried out by the following formula, which uses the value 1.724 and is widely used in Morocco:

$$MO = CO \times 1.724, \qquad (4)$$

the average erodibility factor K per basin is calculated by the formula:

$$K_{\text{average}} = \frac{\sum (K \times \text{number})}{\sum \text{number}}, \qquad (5)$$

with K: K-factor per value, number: counting of pixels with the same K-value in the Arcmap allocation table.

The results of the calculation of the average K-factor for each watershed are in Table 8.

<table><caption>**Table 8.** Average erodibility factor for each watershed.</caption><thead><tr><td>Name Watershed</td><td>Erodibility Factor K<sub>average</sub></td></tr></thead><tbody><tr><td>Bouregreg</td><td>0.352</td></tr><tr><td>Grou</td><td>0.350</td></tr><tr><td>Ain Loudah</td><td>0.348</td></tr><tr><td>Sidi Mohamed Cherif</td><td>0.348</td></tr></tbody></table>

The Kmoy factor used is 0.35. This value corresponds to the silt, with the latter being a sedimentary formation whose grains are of intermediate size between clays and sands. The loss of silt leads to a decrease in the water retention capacity of the soil. The result is an increased erodibility and an increased risk of erosion. Because silt is often suspended in water, it is easily transported by floods and can contribute to the siltation of dam reservoirs. Figure 6 shows the distribution of the K factor over the Bouregreg basin for the four hydrological stations located immediately upstream of the SMBA dam.
<figure><img src="image_116.png" /><figcaption>Figure 6. K-factor of the Bouregreg basin at the four hydrological stations.</figcaption></figure>

### 3.2.2. Topographic Factor (LS)

The topographic factor LS is an essential parameter of the model. It expresses the result of erosion due to the combined effect of the degree of slope and its length. Most recent studies for the determination of the soil loss potential at the watershed level by the MUSLE model utilize the equation established by Wischmeier and Smith [9,23], which is expressed as:

$$LS = \left( \frac{\kappa}{22.13} \right)^m \times 65.4 \sin^2 \beta + 4.56 \sin \beta + 0.0654, \quad (6)$$

with: $κ$: length of the slope in m, $β$: slope in degrees, $m = 0.5$ if $β \ge 5\%$, $m = 0.3$ if $1 < β < 5\%$, $m = 0.2$ if $β \le 1$.

First, the slope map was established by the spatial analyst tool of the ArcGis software using a 90 m resolution DTM. Then, the LS map (Figure 7) by basin was generated using the Wishmeierand Smith formula [9,23]. The mean slopes and mean LS values are given in Table 9.

<table><caption>Table 9. Mean slope and mean topographic factor (LS).</caption><thead><tr><th></th><th>Mean Slope<br>(%)</th><th>Mean Topographic Factor (Ls)</th></tr></thead><tbody><tr><td>Bouregreg à Lala Chafia</td><td>16.6</td><td>0.48</td></tr><tr><td>Grou à Ras Fathia</td><td>14.6</td><td>0.49</td></tr><tr><td>Korifla à Ain Loudah</td><td>9.8</td><td>0.27</td></tr><tr><td>Machraa à Sidi Mohamed Cherif</td><td>12.2</td><td>0.58</td></tr></tbody></table>
<figure><img src="image_117.png" /><figcaption>Figure 7. LS factor of the Bouregreg basin at the four hydrological stations.</figcaption></figure>

### 3.2.3. Land Cover Factor (C)

The soil protection factor (C) indicates the degree of soil protection by vegetation cover. This factor has undergone several changes since the establishment of the universal soil loss equation. The vegetation cover is—after topography—the second most important factor controlling the risk of soil erosion. The value of C depends mainly on the percentage of vegetation cover and the growth phase. The Cfactor map for the Bouregreg catchment area was derived from the land use maps. These were determined from the use of remote sensing data and field observations [48]. The land cover map was extracted from SPOT satellite images at 20 m resolution combined with recent Landsat ETM+ images (2011/2012) using the supervised classification method [41]. Another approach could be to derive the C factor from NDVI maps as practiced on the Wadi Mina in Algeria by Toumi et al. [16].

The distribution of surfaces according to the nature of the vegetation cover is carried out by the ArcGis tool, which subdivides the catchment area into several polygons. Each polygon corresponds to a specific type of vegetation cover. Each vegetation cover corresponds to a factor given by Wischmeir and Smith according to the theme [49]. A color is assigned to each polygon that represents a land cover type. For the calculation of the Cfactor, an independent calculation is required for each sub-basin. Indeed, the sensitivity to erosion of the different classes is determined from the main land use themes, i.e., forest formations, rangelands, agricultural land, arboriculture, water, and bare soil, whose values
vary between 0 and 1. The calculation of the average vegetation cover factor depends on the factor given for each land use and surface area by the following formula:

$$C_{\text{average}} = \frac{\sum S_i \cdot C_i}{\sum S_i}, \qquad (7)$$

with $S_i$: partial polygon area, $C_i$: value of the C factor according to the theme. Figure 8 shows the C factors for each sub-basin of the Bouregreg River.

<figure><img src="image_118.png" /><figcaption>Figure 8. Factor C of the Bouregreg basin at the four hydrological stations.</figcaption></figure>

The application of the formula leads to the results presented in Table 10.
<table><caption>**Table 10.** Average vegetation cover factor C for each watershed.</caption><thead><tr><th rowspan="2">Watershed</th><th>Land Cover Factor (C)</th></tr><tr><th>C<sub>average</sub></th></tr></thead><tbody><tr><td>Bouregreg</td><td>0.26</td></tr><tr><td>Grou</td><td>0.45</td></tr><tr><td>Korifla</td><td>0.41</td></tr><tr><td>Machraa</td><td>0.33</td></tr></tbody></table>

The vegetation cover factor (C) can vary from close to 0 for well-protected soils, to 1 for striated surfaces that are very sensitive to gully erosion. The determination of this factor for the Bouregreg watershed and sub-basins is based on the density of vegetation and the height of the vegetation strata. These data were deduced from the field updating of the land use map by the Royal Centre for Remote Sensing in Space of Morocco (CRTS) available at the ABHBC. The values assigned to the different land use patterns are based on Wischmeier and Smith's tables for forests, matorrals, and pastures. The results show that the values of the C factor range from 0.26 for the sub-basin located in the northeast region, the Bouregreg, the wettest, to 0.45 for the sub-basin that extends furthest south, which is more arid [50,51]. The spatial distribution of the vegetation cover index by class for the Bouregreg watershed shows on 50.45% of the total surface area, i.e., nearly 181,556 ha, an index lower than or equal to 0.2, indicating good protection, while 49.5% of the surface area shows very low protection against erosion, i.e., about 177,986 ha. It can also be noted that the upper northeastern basin, which is both the most humid and the least covered with vegetation, is not, however, subject to high erodibility [52].

### 3.2.4. Average Cultivation Practices and Amenities Factor P

The P factor is a dimensionless factor expressing soil protection through agricultural practices (P). This factor takes into account purely anti-erosion practices. Contour, strip or terrace cultivation, bench planting, and ridging are the most effective soil conservation practices. These practices proportionally affect erosion by altering the flow pattern or direction of surface runoff and by reducing the amount and speed of runoff. P values are less than or equal to 1. A value of 1 is assigned to lands where none of the practices listed are used. P values vary between 0 and 1 depending on the practice adopted and also on the slope. In the Bouregreg watershed, there are few anti-erosion facilities and farmers do not use anti-erosion cultivation practices. These actions are small-scale and do not have a major impact on reducing erosion due to the size of the basin. As a result, a P value equal to 1 has been assigned to the entire area of the basin.

### 3.3. Application and Calibration of the MUSLE Model

#### 3.3.1. Application of the MUSLE Model

The parameters required for the application of the MUSLE model for the evaluation of solid inputs by the MUSLE model are summarized in Table 11.

<table><caption>**Table 11.** Averaging factors for the universal soil loss equation modified per sub-basin.</caption><thead><tr><th></th><th>Mean Topographic Factor (LS)</th><th>Mean Erodibility Factor (K)</th><th>Mean Vegetation Cover Factor (C)</th><th>Mean Development Factor (P)</th></tr></thead><tbody><tr><td>Bouregreg</td><td>0.48</td><td>0.32</td><td>0.26</td><td>1</td></tr><tr><td>Grou</td><td>0.49</td><td>0.32</td><td>0.45</td><td>1</td></tr><tr><td>Korifla</td><td>0.27</td><td>0.34</td><td>0.33</td><td>1</td></tr><tr><td>Machraa</td><td>0.58</td><td>0.31</td><td>0.41</td><td>1</td></tr></tbody></table>
By introducing the various factors of the MUSLE model, the amount of solid inputs produced at the outlet of each watershed is calculated using the flow and volume measurement data for each flood, according to the model equation. The solids transport is then calculated using the measured suspended solids concentrations for the same floods. The concentrations are linearly interpolated over the range of flows in each hydrograph. The bottom solid transport is taken as 10% of the total volume of suspended solids [53], which remains an empirical but widely used value.

### 3.3.2. Calibration of the MUSLE Model

Calibration of the model at each hydrological station is carried out using the percentage bias method (PBIAS), which consists of minimizing the difference between the observed mean and the mean predicted by the model by acting on the model's scaling parameters a and b:

$$ \text{PBIAS} = \left[ \frac{\sum_{i=1}^{n} (Y_{i}^{\text{obs}} - Y_{i}^{\text{sim}}) \times (100)}{\sum_{i=1}^{n} (Y_{i}^{\text{obs}})} \right], \quad (8) $$

with $Y_i^{\text{obs}}$: observed solid transport, $Y_i^{\text{sim}}$: solid transport simulated by MUSLE, n: number of observations.

In our case, a and b are respectively equal to 11.80 and 0.56. The calibration period is spread over two years of observations available. The optimal value of PBIAS is 0, where PBIAS has the ability to clearly indicate poor model performance [54]. To ensure model reliability, the model is also calibrated using the KGE (Kling-Gupta Efficiency) Method [55]:

$$ \text{KGE} = 1 - \sqrt{(r-1)^2 + \left(\frac{\sigma_{\text{sim}}}{\sigma_{\text{Obs}}} - 1\right)^2 + \left(\frac{\mu_{\text{sim}}}{\mu_{\text{Obs}}} - 1\right)^2}, \quad (9) $$

with r: simple correlation coefficient between observed and simulated values, $\sigma_{\text{Obs}}$: standard deviation of observed solid transport values, $\sigma_{\text{sim}}$: standard deviation of simulated solid transport values, $\mu_{\text{Obs}}$: mean of observed solid transport values, $\mu_{\text{sim}}$: mean of simulated solid transport values. The optimal value of the KGE factor is 1.

The PBAIS calibration of the model shows that the scaling factor (a) does not change while (b) is more sensitive, but remains close to the initial value. Table 12 shows the results of the PBIAS model calibration.

<table><caption>Table 12. Model calibration results at the four stations studied by percentage bias method (PBIAS).</caption><thead><tr><th rowspan="2">Name of Basin</th><th colspan="2">Scaling Parameters after Calibration</th><th colspan="3">Calibration Criteria</th></tr><tr><th>A</th><th>B</th><th>Bias</th><th>KGE</th><th>r</th></tr></thead><tbody><tr><td>Bouregreg</td><td>11.799</td><td>0.548</td><td>0.00</td><td>0.48</td><td>0.58</td></tr><tr><td>Grou</td><td>11.799</td><td>0.535</td><td>0.00</td><td>0.11</td><td>0.26</td></tr><tr><td>Korifla</td><td>11.800</td><td>0.541</td><td>0.00</td><td>0.49</td><td>0.57</td></tr><tr><td>Machraa</td><td>11.799</td><td>0.502</td><td>0.00</td><td>0.65</td><td>0.70</td></tr></tbody></table>

Calibration of the model by the KGE method shows that the scaling factors a and b are very sensitive, with the terms being very far from the initial values, and are quite or widely and varyingly different from one basin to another. Table 13 shows the results of model calibration using the KGE method.
<table><caption>**Table 13.** Model calibration results for the four stations studied by KGE.</caption><thead><tr><th rowspan="2">Name of Basin</th><th colspan="2">Scaling Parameters after Calibration</th><th colspan="3">Calibration Criteria</th></tr><tr><th>A</th><th>b</th><th>Bias</th><th>KGE</th><th>r</th></tr></thead><tbody><tr><td>Bouregreg</td><td>0.039</td><td>0.807</td><td>0.024</td><td>0.601</td><td>0.603</td></tr><tr><td>Grou</td><td>11.799</td><td>0.546</td><td>-0.236</td><td>0.149</td><td>0.272</td></tr><tr><td>Korifla</td><td>1.727</td><td>0.643</td><td>-0.055</td><td>0.537</td><td>0.545</td></tr><tr><td>Machraaf</td><td>1.421</td><td>0.620</td><td>-0.017</td><td>0.689</td><td>0.690</td></tr></tbody></table>

### 3.4. Comparison of Observations with MUSLE Model Results

#### 3.4.1. Aguibat Zear Hydrological Station on the Bouregreg River

During the two hydrological years studied, twenty-seven floods were recorded. The total solid volumes observed were 742,334 tons in 2016/2017 and 439,735 tons in 2017/2018. Table 14 summarizes the results from applying the model to the Aguibat Zear hydrological station.

<table><caption>**Table 14.** Results of the MUSLE model at Aguibat Zear on the Bouregreg and a comparison with observations.</caption><thead><tr><th rowspan="2">Hydrologic Year</th><th colspan="2">Observed Solid Transport (tons)</th><th colspan="3">Solid Transport (MUSLE) (tons)</th><th colspan="3">Difference Observation-MUSLE (%)</th></tr><tr><th></th><th>Uncalibrated</th><th>Calibration PBAIS</th><th>Calibration KGE</th><th>Uncalibrated</th><th>Calibration PBAIS</th><th>Calibration KGE</th></tr></thead><tbody><tr><td>2016/2017</td><td>742 334</td><td>619 023</td><td>487 494</td><td>436 922</td><td>-17</td><td>-34</td><td>-41</td></tr><tr><td>2017/2018</td><td>439 735</td><td>888 011</td><td>694 575</td><td>716 809</td><td>+102</td><td>+58</td><td>+63</td></tr></tbody></table>

Examination of the results of the MUSLE model calibrated by the two approaches—PBAIS and KGE—shows that calibration by the PBAIS method better simulates solid transport. Figure 9 compares the simulated and observed solid transport. The latter events are less well simulated by MUSLE than the former.

<figure><img src="image_119.png" /><figcaption>**Figure 9.** Comparison between MUSLE and observations at Aguibat Zear on the Bouregreg.</figcaption></figure>
### 3.4.2. Ras Fathia Hydrological Station on the Grou River

During the two hydrological years studied, twenty-one floods were recorded. The total observed solid volumes amount to 517,194 tons in 2016/2017 and 896,438 tons in 2017/2018. Table 15 summarizes the results of applying the model to the Ras Fathia hydrological station for the two versions of the model calibration.

<table><caption>Table 15. Results of the MUSLE model at Ras Fathia on the Grou and comparison with observations.</caption><thead><tr><th rowspan="2">Hydrologic Year</th><th colspan="2">Observed Solid Transport (tons)</th><th colspan="2">Solid Transport (MUSLE) (tons)</th><th colspan="2">Difference Observation-MUSLE (%)</th></tr><tr><th></th><th>Uncalibrated</th><th>Calibration PBAIS</th><th>Calibration KGE</th><th>Uncalibrated</th><th>Calibration PBAIS</th><th>Calibration KGE</th></tr></thead><tbody><tr><td>2016/2017</td><td>673 551</td><td>875 366</td><td>517 194</td><td>649 389</td><td>+30</td><td>-23</td><td>-4</td></tr><tr><td>2017/2018</td><td>740 083</td><td>1 446 831</td><td>896 438</td><td>1 098 213</td><td>+95</td><td>+21</td><td>+48</td></tr></tbody></table>

Examination of the results of the MUSLE model calibrated by the two approaches PBAIS and KGE, shows that calibration by the PBAIS method better simulates solid transport. Figure 10 shows the comparison of simulated and observed solid transport. As for the Bouregreg, the events of the second year are less well represented by MUSLE.

<figure><img src="image_120.png" /><figcaption>Figure 10. Comparison between MUSLE and observations at Ras Fathia on the Grou.</figcaption></figure>

### 3.4.3. Sidi Mohamed Cherif Hydrological Station on the Machraa River

During the two hydrological years studied, twelve floods were recorded. The total solid volumes observed were 28,155 tons in 2016/2017 and 48,422 tons in 2017/2018. Table 16 summarizes the results of the application of the model to the hydrological station Sidi Mohamed Cherif.
<table><caption>**Table 16.** Results of the MUSLE model at Sidi Mohamed Cherif on the Machraa and comparison with observations.</caption><thead><tr><th rowspan="2">Hydrologic Year</th><th colspan="2">Observed Solid Transport (tons)</th><th colspan="3">Solid Transport (MUSLE) (tons)</th><th colspan="2">Difference Observation-MUSLE (%)</th></tr><tr><th></th><th>Uncalibrated</th><th>Calibration PBAIS</th><th>Calibration KGE</th><th>Uncalibrated</th><th>Calibration PBAIS</th><th>Calibration KGE</th></tr></thead><tbody><tr><td>2016/2017</td><td>28,155</td><td>75,194</td><td>27,431</td><td>26,225</td><td>+167</td><td>-3</td><td>-7</td></tr><tr><td>2017/2018</td><td>48,422</td><td>141,249</td><td>49,146</td><td>51,657</td><td>+192</td><td>+1</td><td>+7</td></tr></tbody></table>

Examination of the results of the MUSLE model calibrated by the two approaches—PBAIS and KGE—shows that calibration by the PBAIS method better simulates solid transport. Figure 11 shows the comparison of simulated and observed solid transport. Again, the MUSLE simulations are significantly too high for the events in the recording portion.

<figure><img src="image_121.png" /><figcaption>**Figure 11.** Comparison between MUSLE and observations at Sidi Mohamed Cherif on the Machraa.</figcaption></figure>

### 3.4.4. Ain Loudah Hydrological Station on the Korifla River

During the two hydrological years studied, twenty-seven floods were recorded. The total solid volumes observed were 68,759 tons in 2016/2017 and 92,324 tons in 2017/2018. Table 17 summarizes the results of the application of the model to the Ain Loudah hydrological station.

<table><caption>**Table 17.** Results of the MUSLE model at Ain Loudah on the Korifla and comparison with observations.</caption><thead><tr><th rowspan="2">Hydrologic Year</th><th colspan="2">Observed Solid Transport (tons)</th><th colspan="3">Solid Transport (MUSLE) (tons)</th><th colspan="2">Difference Observation-MUSLE (%)</th></tr><tr><th></th><th>Uncalibrated</th><th>Calibration PBAIS</th><th>Calibration KGE</th><th>Uncalibrated</th><th>Calibration PBAIS</th><th>Calibration KGE</th></tr></thead><tbody><tr><td>2016/2017</td><td>68,759</td><td>41,083</td><td>28,881</td><td>27,739</td><td>-40</td><td>-58</td><td>-60</td></tr><tr><td>2017/2018</td><td>92,324</td><td>191,859</td><td>132,202</td><td>142,138</td><td>+108</td><td>+43</td><td>+54</td></tr></tbody></table>

Examination of the results of the MUSLE model calibrated by the two approaches PBAIS and KGE, shows that calibration by the PBAIS method better simulates solid transports. Figure 12 shows the comparison of simulated and observed solid transports. Out of the four stations studied, it is on the Korifla at Ain Loudah that the simulations seem to be the most efficient.
<figure><img src="image_122.png" /><figcaption>Figure 12. Comparison between MUSLE and observations at Ain Loudah on the Korifla.</figcaption></figure>

From the deviations from observations at the four hydrological stations, it can be deduced that in general, the model calibrated by the PBAIS method gives better results than the model calibrated by the KGE method. Thus, from the observation we note that the differences varied between −58% and +58% for the model calibrated by PBAIS, between −60% and +63% for the model calibrated by the KGE method, and between −40% and +192% for the model not calibrated.

### 3.5. Comparison of Observations and Results of the MUSLE Model at the Four Hydrological Stations Upstream from the SMBA Dam, with Bathymetric Data

To validate the calibration method, the Nash and Sutcliffe NSE index was used to measure the performance of the model. According to Nash–Sutcliffe [56], NSE is defined as:

$$NSE = 1 - \frac{\sum_{i=1}^{n} (Y_{i}^{\text{obs}} - Y_{i}^{\text{sim}})^{2}}{\sum_{i=1}^{n} (Y_{i}^{\text{obs}} - \mu_{\text{Obs}})^{2}}, \quad (10)$$

With NSE: Nash coefficient, $Y_{i}^{\text{obs}}$: observed solid transport, $Y_{i}^{\text{sim}}$: solid transport simulated by MUSLE, $\mu_{\text{Obs}}$: average of the observed solid transport values.

The NSE index varies from −∞ to 1, such that if NSE = 1, then the modeled values match the observations perfectly, while a value above 0 shows a relationship between simulation and reality, and a value below 0 shows that there is no relationship between the two. In other words, the closer the efficiency is to 1, the more the model is observed to be accurate.

Table 18 shows the different calculated values of the Nash–Sutcliffe index, which compares observations to the values calculated by the uncalibrated MUSLE model and observations to the values of the calibrated MUSLE model.
<table><caption>**Table 18.** Nash–Sutcliffe index (NSE) at study stations.</caption><thead><tr><th rowspan="2">Name of Basin</th><th colspan="2">Scaling Parameters after Calibration</th><th colspan="2">Nash–Sutcliffe Index (NSE)</th></tr><tr><th>a</th><th>B</th><th>MUSLE</th><th>MUSLE calibrated</th></tr></thead><tbody><tr><td>Bouregreg à Aguibat Zear</td><td>11.799</td><td>0.548</td><td>0.38</td><td>0.46</td></tr><tr><td>Grou à Ras Fathia</td><td>11.799</td><td>0.535</td><td>-0.07</td><td>0.02</td></tr><tr><td>Korifla à Ain Loudah</td><td>11.800</td><td>0.541</td><td>-0.35</td><td>0.30</td></tr><tr><td>Machraa à Sidi</td><td>11.799</td><td>0.502</td><td>-9.8</td><td>0.47</td></tr></tbody></table>

In general, the uncalibrated MUSLE model has negative NSE values except for those in the Bouregreg basin, which shows that the model values are far from reality, whereas the NSE values of the calibrated model are positive, varying between 0.02 and 0.47, which shows that the calibrated model is superior in reliably representing reality.

Though, for the continuity of the experiment we are continuing the calculation on the four basins in order to complete the comparisons with the volumes of sediment deduced from the bathymetric measurements carried out in the SMBA dam reservoir. For the calculation of the volume of suspended matter, an earthy density of 1.5 t/m³ is adopted [57]. The siltation of the dam is considered to be equal to 9.49 × 10⁶ m³/year. Table 19 summarizes the results of the calculations for all the hydrological stations.

<table><caption>**Table 19.** Summary table of observation results, application of the MUSLE, calibrated MUSLE\_model and bathymetry for the four sub-basins of the Bouregreg.</caption><thead><tr><th rowspan="2">Name of Basin</th><th colspan="3">Solid Transport 2016/2017</th><th colspan="3">Solid Transport 2017/2018</th></tr><tr><th>Observation</th><th>MUSLE</th><th>MUSLE calibrated</th><th>Observation</th><th>MUSLE</th><th>MUSLE calibrated</th></tr></thead><tbody><tr><td>Bouregreg</td><td>742,334</td><td>619,023</td><td>487,494</td><td>439,735</td><td>888,011</td><td>694,575</td></tr><tr><td>Grou</td><td>673,551</td><td>875,366</td><td>517,194</td><td>740,083</td><td>1,446,831</td><td>896,438</td></tr><tr><td>Korifla</td><td>68,759</td><td>41,083</td><td>28,881</td><td>92,324</td><td>191,859</td><td>132,202</td></tr><tr><td>Machraa</td><td>28,155</td><td>75,194</td><td>27,431</td><td>48,422</td><td>141,249</td><td>49,146</td></tr><tr><td>Total (tons)</td><td>1,512,799</td><td>1,610,666</td><td>1,061,000</td><td>1,320,564</td><td>2,667,950</td><td>1,772,361</td></tr><tr><td>Total (10<sup>6</sup> m<sup>3</sup>)</td><td>1.01</td><td>1.07</td><td>0.71</td><td>0.88</td><td>1.78</td><td>1.18</td></tr></tbody></table>

Thus, the results listed in Tables 19 and 20 show that the calibrated model is closer to the measurement than the non-calibrated model, and that the observed values of sediment transport at the four stations, which vary between 1.01 × 10⁶ m³ and 0.88 × 10⁶ m³ for the two hydrological years 2016/2017 and 2017/2018 respectively, and remain largely lower than the quantity of sediments obtained by bathymetric measurements in the SMBA dam, which is approximately 9.49 × 10⁶ m³.

<table><caption>**Table 20.** Comparison of observation results, application of the calibrated MUSLE, MUSLE\_model, and bathymetry for the four sub-basins of the Bouregreg during 2016/2017 and 2017/2018.</caption><thead><tr><th>Hydrologic Year</th><th>Bathymetry (10<sup>6</sup> m<sup>3</sup>)</th><th>Difference MUSLE/ Observation</th><th>Difference MUSLE Calibrated/ Observation</th><th>Difference Observation/ Bathymetry</th><th>Difference MUSLE/ Bathymetry</th><th>Difference MUSLE Calibrated/ Bathymetry</th></tr></thead><tbody><tr><td>2016/2017</td><td rowspan="2">9.25</td><td>+5%</td><td>-30%</td><td>-89%</td><td>-89%</td><td>-93%</td></tr><tr><td>2016/2018</td><td>+102%</td><td>+34%</td><td>-91%</td><td>-81%</td><td>-88%</td></tr></tbody></table>

## 4. Discussion and Conclusions

The relationships between man, environment, and solid transport are well understood at the spatial scale of catchment areas, and modeling tools using satellite data are particularly suitable [58]. The results obtained from the application of the model (MUSLE), show that the soils of the Bouregreg catchment area are affected by several factors engendering erosion, i.e., steep slopes, low vegetation
cover and soil erodibility. Thus, the processing of spatial remote sensing images to determine the biophysical characteristics of the watershed has made it possible to identify areas of high erodibility, low vegetation cover, and steep slopes, which can be used for potential development purposes.

The analysis of solid transports observed during two years at the four hydrological stations that control the main tributaries of the SMBA dam reservoir showed that the Bouregreg and the Grou generated 1,415,885 tons of suspended matter, i.e., 94% and 956,929 tons of suspended matter, i.e., 73% of the total suspended matter generated by the four basins during the 2016/2017 and 2017/2018 hydrological years. Compared to bathymetry, the four sub-basins generate only 12% of the solid transport to the dam in the best case. Hence the relevance of concentrating anti-erosion actions on the dam banks and the intermediate basin which contribute to more than 88% of the solid transport. These actions will have to be carried out in areas of high erodibility and slope to ensure the effectiveness of the intervention.

On the other hand, examination of the figures for MUSLE and calibrated MUSLE values shows that the model in these two versions does not predict peaks in solid transport well. It overestimates low values and underestimates high solid transport values most of the time. Indeed, the deviations from observations range from +5% for the dry year, to +102% for the wet year for the non-calibrated model. These results are similar to several applications of the model in different areas around the world [40]. The deviations recorded for the calibrated model are significantly improved. They vary respectively between -30% and +34% [59]. Thus, the calibrated version of the model gives better results with less than +55% error in sediment prediction [60]. These results are interesting on an interannual average, but do not allow the prediction of quantities transported to the stations during unit floods due to the very high dispersion of the relationship between observation and modeling, which results in average-to-low Nash index values. The large positive differences are due to the underestimation of solid transport peaks by the model. Indeed, the model involves the entire basin without concern for the spatial and temporal distribution of rainfall.

Although the model overestimates values relative to observations, comparison of the model to bathymetric data shows that the differences between mean siltation and model results range from -89% to -81% for the non-calibrated version of the model, and from -93% to -88% for the calibrated version. This underestimation of the inflows to the dam in relation to the mean silting is explained by the fact that the upstream basins do not control the whole SMBA dam basin and that the basins immediately surrounding the banks of the dam reservoir bring significant quantities of sediment to the dam, as assumed by the first results of Maleval [61] on the SMBA dam reservoir, and this is mainly after its elevation, for reasons that are perhaps not only related to the physics of the environment, but also perhaps to the evolution of agro-sylvo-pastoral practices in an environment undergoing rapid socio-economic change due to the proximity of the capital and the development of the "Grand Rabat" and new traffic routes [62,63]. The gaps between the observed solid transport and the silting rate, which are respectively -89% and -91%, confirm this observation, which is also the situation observed by Hallouz et al. [64] on the Wadi Mina basin and the Sidi Mohamed Ben Aouda dam, of a size comparable to the Bouregreg in Western Algeria.

In order to improve the model's output, its scaling parameters should be reviewed and its sensitivity to different physical factors including basin size by adopting a distributed form of sub-basin modeling should be tested [65]. This approach requires a long series of observations that are currently unavailable, but which are being acquired since the ABHBC, which manages the water resource on the Bouregreg catchment area and the SMBA dam, has set up permanent concentration of suspended solids observation on the stations of the basin, and soon a long time series will be available to refine these first results. Michalec et al. [66] have pointed out the necessity to calibrate the MUSLE model with long time series of observations in order to reduce the uncertainty of simulated erosion provided, even if in some cases very good correlations have been obtained from a calibration with only 10 events like in Nigeria [30]. These measurements will allow a monitoring of the contribution of each sub-basin in order to better direct intervention measurements towards the most productive areas on the one
hand, thus reducing the overall cost of the measures to reduce erosion, and to improve modeling and subsequent solid transport forecasting on the other hand.

Eventually, the results of the MUSLE model confirm the relevance of its application in watershed management studies once calibrated, since it integrates all forms of erosion observed at the watershed scale, coupled with hydrological data, which greatly improves accuracy compared to USLE [67] and to RUSLE [68]. Indeed, the average specific degradation at the level of the Bouregreg basin alone at the Aguibat Zear station is of the order of 13.81 t/ha/year according to Moussebbih et al. [69] when calculated by the RUSLE method, whereas those of calibrated MUSLE and observation are of the order of 1.6 t/ha/year. The studies carried out on the wadi Mina basin in Algeria, a basin similar to the Bouregreg basin [64], confirm this finding with deviations of the RUSLE model results from the measurement of suspended solids concentrations, which can reach 79%. In Spain, Ramos-Diez et al. [70] have used the USLE formula to assess the amount of sediment inputs to 25 very small check dams used in an important restoration project over a 9 km² area in the north of Spain, showing interesting but still mitigated results and concluding that the size and shape of the dams had an impact on the quality of the USLE assessment.

Kronvang et al. [71] showed that bank erosion was the dominant sediment source (90–94%) in the River Odense catchment in Denmark during three study years. They add that in-channel and overbank sediment sinks and storage dominated the sediment budget, as 79–94% of the sediment input from all sources was not exported from the catchment during the three study years. This is in accordance with our results, as the regular bathymetric survey to monitor the silting up rate of the dam have shown that half of the sediment input to the reservoir comes from the dam banks. Our results are very important for the forthcoming works for erosion mitigation on the Bouregreg catchment; they will enable more efficient orientation of the soil conservation and restoration works that may be carried out in the future, with priority being given to the environment close to the dam, in order to preserve the water capacity of the SMBA dam reservoir while optimizing the financial resources mobilized. Palazon and Navas [72] have shown the same interest to monitor the sediment sources draining to a large reservoir in the Esera River, Ebro basin in Spain. The small basin size (1500 km²) allowed them to use the SWAT model as an alternative to MUSLE with two gauging stations, each one being covered by different types of soil. This might make it possible to implement the SWAT model on the Bouregreg basin, when the ABHBC will have started the monitoring of sediment transport at all the stations of the basin.

To conclude, it should be pointed out that in the context of climate change, which predicts, in Morocco, an increase in temperatures and a drop in rainfall [73–75], the projections of erosion and sediment transport evolution that are possible from climate model outputs, and which are based on USLE or SWAT type approximation models [76], will need observed data to calibrate and validate their results.

**Author Contributions:** Conceptualization, M.A.E. and G.M.; methodology, M.A.E.; software, M.A.E.; validation, G.M., I.K. and A.Z.; formal analysis, M.A.E.; investigation, M.A.E.; resources, A.Z.; data curation, A.Z.; writing—original draft preparation, M.A.E. and G.M.; writing—review and editing, M.A.E.; visualization, M.A.E.; supervision, G.M., I.K. and A.Z.; project administration, G.M., I.K. and A.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** Thanks to ABHBC for the hydrological and rainfall data.

**Conflicts of Interest:** The authors declare no conflict of interest.

**References**

1. Mahe, G.; Emran, A.; Brou, Y.T.; Tra Bi, A.Z. Analyse statistique de l'évolution de la couverture végétale à partir d'images MODIS et NOAA sur le bassin versant du Bouregreg (Maroc). *Géo Obs.* **2012**, *20*, 33–44.
2. Laouina, A.; Aderghal, M.; Al Karkouri, J.; Antari, M.; Chaker, M.; Laghazi, Y.; Machmachi, I.; Machouri, N.; Nafaa, R.; Naïmi, K. The efforts for cork oak forest management and their effects on soil conservation. *For. Syst.* **2010**, *19*, 263–277. [CrossRef]

3. Schmidt, S.; Alewell, C.; Meusburger, K. Monthly RUSLE soil erosion risk of Swiss grasslands. *J. Maps* **2019**, *15*, 247–256. [CrossRef]

4. Yahiaoui, S.; Zerouali, A. Etude de l'évolution de l'occupation du sol sur deux grands bassins d'Algérie et du Maroc, et relation avec la sédimentation dans les barrages. In *Considering Hydrological Change in Reservoir Planning and Management*; IAHS Publ 362; Schumann, A., Belyaev, V.B., Gargouri, E., Kucera, G., Mahe, G., Eds.; IAHS Press: Wallingford, UK, 2013; pp. 115–124.

5. Khomsi, K.; Mahe, G.; Sinan, M.; Snoussi, M. Hydro-climatic variability in two Moroccan watersheds: A comparative analysis of temperature, rain and flow regimes. In *Climate and Land Surface Changes in Hydrology*, IAHS Publ. 359; Boegh, E., Blyth, E., Hannah, D.M., Hisdal, H., Kunstmann, H., Su, B., Yilmaz, K.K., Eds.; IAHS Press: Wallingford, UK, 2013; pp. 183–190.

6. Khomsi, K.; Mahe, G.; Tramblay, Y.; Sinan, M.; Snoussi, M. Regional impacts of global change: Seasonal trends in extreme rainfall, runoff and temperature in two contrasted regions of Morocco. *Nat. Hazards Earth Syst.* **2016**, *16*, 1079–1090. [CrossRef]

7. Mahe, G.; Emran, A.; Brou, Y.T.; Tra Bi, A.Z. Impact de la variabilité climatique sur l'état de surface du bassin versant du Bouregreg (Maroc). *Eur. J. Sci. Res.* **2012**, *84*, 417–425.

8. ABV. *Plan National d'Aménagement des Bassins Versant*; Haut-Commissariat aux Eaux et Forêts et à la Lutte Contre la Désertification: Rabat, Morocco, 1996.

9. Wischmeier, W.H.; Smith, D.D. Predicting Rainfall-Erosion Losses from Cropland East of the Rocky Mountains: Guide for Selection of Practices for Soil and Water Conservation; US Department of Agriculture: Washington, DC, USA, 1965.

10. Williams, J.R. Sediment-Yield Prediction with Universal Equation Using Runoff Energy Factor; ARS-S. Southern Region, Agricultural Research Service, US Department of Agriculture: Washington, DC, USA, 1975; Volume 40, p. 244.

11. El Bahi, S.; El Wartiti, M.; Yassin, M.; Calle, A.; Casanova, J.L. Applying revised universal loss equation model to forest lands in Central Plateau of Morocco. *Rev. Teledetección* **2005**, *23*, 89–97.

12. Chadli, K. Estimation of soil loss using RUSLE model for Sebou watershed (Morocco). *Model. Earth Syst. Environ.* **2016**, *2*, 51. [CrossRef]

13. Elaloui, A.; Marrakchi, C.; Fekri, A.; Maimouni, S.; Aradi, M. USLE-based assessment of soil erosion by water in the watershed upstream Tessaoute (Central High Atlas, Morocco). *Modeling Earth Syst. Environ.* **2017**, *3*, 873–885. [CrossRef]

14. Hara, F.; Achab, M.; Emran, A.; Mahe, G.; Fhel, B.E. Estimate the Risk of Soil Erosion Using USLE through the Development of an Open Source Desktop Application: DUSLE (Desktop Universal Soil Loss Equation). In Proceedings of the 3rd International Conference on African Large River Basin Hydrology, Alger, Algeria, May 2018; Available online: https://hal.archives-ouvertes.fr/hal-02397686/ (accessed on 17 June 2020).

15. Boufala, M.; El Hmaidf, A.; Chadli, K.; Essahlaoui, A.; El Ouali, A.; Lahjouj, A. Assessment of the Risk of Soil Erosion Using RUSLE Method and SWAT Model at the M'dez Watershed, Middle Atlas, Morocco. In Proceedings of the E3S Web Conf., The Seventh International Congress "Water, Waste and Environment" (EDE7-2019), Salé, Morocco, 20–22 November 2019; Volume 150. [CrossRef]

16. Toumi, S.; Meddi, M.; Mahe, G.; Brou, Y.T. Cartographie de l'érosion dans le bassin versant de l'Oued Mina en Algérie par télédétection et SIG. *Hydrol. Sci. J.* **2013**, *58*, 1–17. [CrossRef]

17. Mihi, A.; Benarfa, N.; Arar, A. Assessing and mapping water erosion-prone areas in northeastern Algeria using analytic hierarchy process, USLE/RUSLE equation, GIS, and remote sensing. *Appl. Geomat.* **2019**. [CrossRef]

18. Gaubi, I.; Chaabani, A.; Ben Mammou, A.; Hamza, M.H. A GIS-based soil erosion prediction using the Revised Universal Soil Loss Equation (RUSLE) (Lebna watershed, Cap Bon, Tunisia). *Nat. Hazards* **2017**, *86*, 219–239. [CrossRef]

19. Toubal, A.K.; Achite, M.; Ouillon, S.; Dehni, A. Soil erodibility mapping using the RUSLE model to prioritize erosion control in the Wadi Sahouat basin. North-West Alger. *Environ. Monit. Assess.* **2018**, *190*, 4. [CrossRef] [PubMed]
20. El Aroussi, O.; Mesrar, L.; El Garouani, A.; Lahrach, A.; Beaabidaate, L.; Akdi, B.; Jabrane, R. Predicting the potential annual soil loss using the revised universal soil loss equation (RUSLE) in the oued El Malleh catchment (Prerif, Morocco). *Present. Environ. Sustain. Dev.* **2011**, *5*, 5–16.

21. Lamyaa, K.; M'bark, A.; Brahim, I.; Hicham, A.; Soraya, M. Mapping Soil Erosion Risk Using RUSLE, GIS, Remote Sensing Methods: A Case of Mountainous Sub-watershed, Ifni Lake and High Valley of Tifnoute (High Moroccan Atlas). *J. Geogr. Environ. Earth Sci. Int.* **2018**, *14*, 1–11. [CrossRef]

22. Bagarello, V.; Ferro, V.; Giordano, G.; Mannocchi, F.; Todisco, F.; Vergni, L. Statistical check of USLE-M and USLE-MM to predict bare plot soil loss in two Italian environments. *Land Degrad. Dev.* **2018**, *29*, 2614–2628. [CrossRef]

23. Wischmeier, W.H.; Smith, D.D. Predicting Rainfall Erosion Losses: A Guide to Conservation Planning; Department of Agriculture, Science and Education Administration: Washington, DC, USA, 1978.

24. Renard, K.G. Predicting Soil Erosion by Water: A Guide to Conservation Planning with the Revised Universal Soil Loss Equation (RUSLE); United States Government Printing: Washington, DC, USA, 1997.

25. Kinnell, P.I.A. Why the Universal Soil Loss Equation and the Revised Version of It Do Not Predict Event Erosion Well. *Hydrol. Process.* **2005**, *19*, 851–854. [CrossRef]

26. Haan, C.T.; Barfield, B.J.; Hayes, J.C. Design Hydrology and Sedimentology for Small Catchments. *J. Hydrol.* - *Amst. 1963*, *176*, 296–297.

27. Novotny, V. Water Quality: Prevention, Identification and Management of Diffuse Pollution; Van Nostrand-Reinhold Publishers: New York, NY, USA, 1994; ISBN 978-0442005597.

28. Alewell, C.; Borrelli, P.; Meusburger, K.; Panagos, P. Using the USLE: Chances, challenges and limitations of soil erosion modelling. *Int. Soil Water Conserv. Res.* **2019**, *7*, 203–225. [CrossRef]

29. Odongo, V.O.; Onyando, J.O.; Mutua, B.M.; van Oel, P.R.; Becht, R. Sensitivity analysis and calibration of the Modified Universal Soil Loss Equation (MUSLE) for the upper Malewa Catchment, Kenya. *Int. J. Sediment Res.* **2013**, *28*, 368–383. [CrossRef]

30. Adegede, A.P.; Mbajiorgu, C.C. Event-based sediment yield modeling for small watersheds using MUSLE in north-central Nigeria. *Agric. Eng. Int.*: CIGR J. **2019**, *21*, 7–17.

31. Yassin, M.; Pépin, Y.; El Bahi, S.; Zante, P. Evaluation de l'érosion au micro-bassin de Sidi Sbaa: Application du modèle MUSLE. *Ann. Rech. For. Maroc* **2013**, *42*, 171–181.

32. Samaras, A.G.; Koutitas, C.G. Modeling the impact of climate change on sediment transport and morphology in coupled watershed-coast systems: A case study using an integrated approach. *Int. J. Sediment Res.* **2014**, *29*, 304–315. [CrossRef]

33. Fang, H. Impact of land use changes on catchment soil erosion and sediment yield in the northeastern China: A panel data model application. *Int. J. Sediment Res.* **2020**. [CrossRef]

34. Djoukbala, O.; Hasbaia, M.; Benselama, O.; Mazour, M. Comparison of the erosion prediction models from USLE, MUSLE and RUSLE in a Mediterranean watershed, case of Wadi Gazouana (N-W of Algeria). *Modeling Earth Syst. Environ.* **2018**. [CrossRef]

35. Mahe, G.; Benabdelfadel, H.; Dieulin, C.; Elbaraka, M.; Ezzaouini, M.; Khomsi, K.; Rouche, N.; Sinan, M.; Snoussi, M.; Tra Bi, A.; et al. Evolution des débits liquides et solides du Bouregreg. In *Gestion Durable des Terres. Proceedings de la Réunion Multi-Acteurs Sur le Bassin du Bouregreg. CERGéo, Faculté des Lettres et Sciences Humaines, Université Mohammed V-Agdal, Rabat, 28 mai 2013*; Laouina, A., Mahe, G., Eds.; Edité par ARGDT: Rabat, Morocco, 2014; pp. 21–36, ISBN 978-9954-33-482-9.

36. Khali Issa, L.; Ben Hamman Lech-Hab, K.; Raissouni, A.; El Arrim, A. Cartographie quantitative du risque d'érosion des sols par approche SIG/USLE au niveau du bassin versant Kalaya (Maroc Nord Occidental). *J. Mater. Environ. Sci.* **2016**, *7*, 2778–2795.

37. El Jazouli, A.; Barakat, A.; Ghafiri, A.; El Moutaki, S.; Ettaqy, A.; Khellouk, R. Soil erosion modeled with USLE, GIS, and remote sensing: A case study of Ikkour watershed in Middle Atlas (Morocco). *Geosci. Lett.* **2017**, *4*, 1. [CrossRef]

38. Sabri, E. Impact des Processus Erosifs sur les Ressources Naturelles des Bassins Hydrographiques et sur L'état D'envasement des Retenues des Barrages: Application au Bassin Versant de L'oued El Abid en Amont du Barrage de Bin El Ouidane- Maroc. Ph.D. Thesis, Univ. Cadi Ayyad, Marrakech, Morocco, 2018.

39. Sadeghi, S.H.R.; Gholami, L.; Khaledi Darvishan, A.; Saeidi, P. A review of the application of the MUSLE model worldwide. *Hydrol. Sci. J.* **2014**, *59*, 365–375. [CrossRef]
40. Aboulabbes, O. Infiltration Characteristics on a Small Watershed. Master's Thesis, Utah State University, Logan Utah, UT, USA, 1984.

41. El Gaatib, R.; Erraji, A.; Larabi, A. Impact des processus érosifs sur les ressources naturelles des bassins hydrographiques et sur l'état d'envasement des retenues de barrages: Application au Bassin Versant de l'Oued Beht en amont du barrage El Kansra, (Maroc). *Géo Obs.* **2014**, *21*, 4.

42. Jamagne, M. *Bases et Techniques d'une Cartographie des Sols*'; Institut National de la Recherche Agronomique: Paris, France, 1967.

43. Stone, R.P.; Hilborn, D. *Universal Soil Loss Equation, USLE*; Factsheet; Government of Ontario, Ministry of Agriculture, Food and Rural Affairs: Canada, 2000. Available online: http://www.omafra.gov.on.ca/english/engineer/facts/12-051.htm (accessed on 17 June 2020).

44. Pribyl, D.W. A critical review of the conventional SOC to SOM conversion factor. *Geoderma* **2010**, *156*, 75–83.
[CrossRef]

45. Van Bemmelen, J.M. "Über die Bestimmung des Wassers, des Humus, des Schwefels, der in den colloïdalen Silikaten gebundenen Kieselsäure, des Mangansuswim Ackerboden". *Die Landwirtschaftlichen Vers.-Stn.* **1890**, *37*, 279–290.

46. Wolff, E. Entwurf zur boden analyse. Z. Anal. Chem. **1864**, *3*, 85–115. [CrossRef]

47. Sprengel, C. Ueber Pflanzenhumus, Humussäure und humussaure Salze. Arch. für die Gesammtze **1826**, *8*, 145–220. [CrossRef]

48. Tra Bi, Z.A.; Emran, A.; Brou, Y.T.; Mahé, G. Cartographie par arbre de décision de la dynamique de l'occupation du sol du bassin versant du Bouregreg, en région semi-aride au centre Nord-Ouest du Maroc, Abidjan, Côte d'Ivoire. *Rev. Sci. Int. Géomatique* **2014**, *1*, 43–54.

49. Sadiki, A.; Bouhajlassa, S.; Auajjar, J. Utilisation d'un SIG pour l'évaluation et la cartographie des risques d'érosion par utilisation d'un SIG pour l'évaluation et la cartographie des risques d'érosion par l'équation universelle des pertes en sol dans le Rif oriental (Maroc). *Bull. de l'Inst. Sci.* **2004**, *26*, 69–79.

50. Goussot, E.; Brou, Y.T.; Laouina, A.; Chaker, M.; Emran, A.; Machouri, N.; Mahe, G.; Sfa, M.; Tra Bi, A. Land covers dynamic and agricultural statistics on the Bouregreg watershed in Morocco. *Eur. J. Sci. Res.* **2014**, *126*, 191–205.

51. Tra Bi, A.Z.; Brou, T.; Emran, A.; Mahe, G. Remote sensing and GIS analysis of the dynamic of vegetation in climate variability context on the Bouregreg watershed. In *Climate and Land Surface Changes in Hydrology*, IAHS Publ. 359; Boegh, E., Blyth, D.M., Hannah, H., Hisdal, H., Kunstmann, B., Su, K.K., Yilmaz, Eds.; IAHS Press: Wallingford, UK, 2013; pp. 403–410.

52. El Hadraoui, Y. *Étude Diachronique de L'occupation du sol et de Modélisation des Processus Erosifs du Bassin Versant du Bouregreg (Maroc) à Partir des Données de L'observation de la Terre*; Mémoire d'ingénieur du CNAM, Spécialité Géomètre et Topographe: Le Mans, France, 2013; Available online: https://dumas.ccsd.cnrs.fr/dumas-00920465 (accessed on 17 June 2020).

53. Lahlou, A. *La Dégradation Spécifique des Bassins Versants et son Impact sur L'envasement des Barrages*; IAHS Publ. 137; IAHS Press: Wallingford, UK, 1982; pp. 163–169. Available online: https://pdfs.semanticscholar.org/cddc/4f4b21f52e34043e1b0547d00d0e9dd389b9.pdf (accessed on 17 June 2020).

54. Gupta, H.V.; Sorooshian, S.; Yapo, P.O. Status of automatic calibration for hydrologic models: Comparison with multilevel expert calibration. *J. Hydrol. Eng.* **1999**, *4*, 135–143. [CrossRef]

55. Gupta, H.V.; Kling, H.; Yilmaz, K.K.; Martinez, G.F. Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling. *J. Hydrol.* **2009**, *377*, 80–91.
[CrossRef]

56. Nash, J.E.; Sutcliffe, J.V. River flow forecasting through conceptual models part I—A discussion of principles. *J. Hydrol.* **1970**, *10*, 282–290. [CrossRef]

57. Lahlou, A. *Etude du transport solide à la station Dar Es Soltane sur l'oued Bouregreg*; Rapport de travaux du service GDE; Gestion des eaux, Ministère des travaux publics et des communications, Direction de l'hydraulique: Rabat, Morocco, 1971.

58. Mahe, G.; Aksoy, H.; Brou, Y.T.; Meddi, M.; Roosevelt, E. Relationships among man, environment and sediment transport: A spatial approach. *Rev. Sci. Eau* **2013**, *26*, 235–244. [CrossRef]

59. De Vente, J.; Verduyn, R.; Verstraeten, G.; Vanmaercke, M.; Poesen, J. Factors controlling sediment yield at the catchment scale in NW Mediterranean geoecosystems. *J. Soils Sediments* **2011**, *11*, 690–707. [CrossRef]