**Susceptibility Mapping of Soil Water Erosion Using Machine Learning Models**

#### **Amirhosein Mosavi 1,2, \* , Farzaneh Sajedi-Hosseini 3, \*, Bahram Choubin 4 , Fereshteh Taromideh 5 , Gholamreza Rahi <sup>6</sup> and Adrienn A. Dineva 7**


#### Received: 3 May 2020; Accepted: 10 July 2020; Published: 14 July 2020

**Abstract:** Soil erosion is a serious threat to sustainable agriculture, food production, and environmental security. The advancement of accurate models for soil erosion susceptibility and hazard assessment is of utmost importance for enhancing mitigation policies and laws. This paper proposes novel machine learning (ML) models for the susceptibility mapping of the water erosion of soil. The weighted subspace random forest (WSRF), Gaussian process with a radial basis function kernel (Gaussprradial), and naive Bayes (NB) ML methods were used in the prediction of the soil erosion susceptibility. Data included 227 samples of erosion and non-erosion locations through field surveys to advance models of the spatial distribution using predictive factors. In this study, 19 effective factors of soil erosion were considered. The critical factors were selected using simulated annealing feature selection (SAFS). The critical factors included aspect, curvature, slope length, flow accumulation, rainfall erosivity factor, distance from the stream, drainage density, fault density, normalized difference vegetation index (NDVI), hydrologic soil group, soil texture, and lithology. The dataset cells of samples (70% for training and 30% for testing) were randomly prepared to assess the robustness of the different models. The functional relevance between soil erosion and effective factors was computed using the ML models. The ML models were evaluated using different metrics, including accuracy, the kappa coefficient, and the probability of detection (POD). The accuracies of the WSRF, Gaussprradial, and NB methods were 0.91, 0.88, and 0.85, respectively, for the testing data; 0.82, 0.76, and 0.71, respectively, for the kappa coefficient; and 0.94, 0.94, and 0.94, respectively, for POD. However, the ML models, especially the WSRF, had an acceptable performance regarding producing soil erosion susceptibility maps. Maps produced with the most robust models can be a useful tool for sustainable management, watershed conservation, and the reduction of soil and water loss.

**Keywords:** water erosion; susceptibility; Gaussian process; climate change; radial basis function kernel; weighted subspace random forest; extreme events; extreme weather; naive Bayes; feature selection; machine learning; hydrologic model; simulated annealing; earth system science

#### **1. Introduction**

Soil conservation is of utmost importance for sustainable development, food security, and environmental protection [1]. Understanding soil erosion is considered to be an essential practice for soil conservation programs around the world [2]. Currently, soil erosion has increasingly become known as a severe concern for sustainable agriculture, water resource management, and modern civilization [3]. Soil erosion is a significant menace for soil, ecology, and for humanity since the long-term production of soil productive capacity is profoundly affected by the destruction and leaching of soil's organic and topsoil matters [4]. Soil erosion is an intricate process that depends on the plant cover and land use, watershed topography, soil properties, climate, and land management practices. In the last century, soil erosion has intensified due to human activity and is an environmental problem [5]. Primary soil segregates when the rainfall or water flow power is greater than the soil's resistance to corrosion [6]. Generally, there are different types of water erosion, such as sheet, gully, landslide, debris flow, streambank, etc. [7].

In semiarid regions, such as Iran, soil erosion is a significant crisis [8] and can be considered to be one of the critical problems concerning agricultural development, natural resources, and the environment [9]. In such regions, water is limited, and there are many sources of sediment [10]. The high input of sediment in upstream rivers increases the water turbidity, reduces the lifespan of dams owing to reservoir siltation, and negatively affects water quality and biological activity [8]. According to scholars, the mean annual rate of soil erosion in Iran is about 25 tons/ha/year, which is four times more than the mean yearly rate around the world [11,12]. Therefore, the susceptibility mapping of soil erosion is necessary for controlling this critical problem.

Rather than using traditional and experimental models, such as the universal soil loss equation (USLE) [13] and multi-criteria decision-making methods [14], that have been used in water erosion assessments, machine learning (ML) models are known to be successful methods [15,16]. Different ML methods, such as support vector machine (SVM), boosted regression trees (BRT), random forest (RF), naive Bayes (NB), and artificial neural network (ANN), have been used for landslides [17–23], debris flows [24–26], and gully erosion [27–30]. For instance, Angileri et al. [15] used the stochastic gradient tree boost (SGT) for water erosion susceptibility mapping in central-northern Sicily, Italy. The results indicated that the applied model had excellent reliability (accuracy from 0.87 to 0.92). Recently, Garosi et al. [31] applied the RF, SVM, and NB models, along with the generalized additive model (GAM), to predict the gully erosion susceptibility in the Ekbatan Dam drainage basin, Iran. The results indicated that the RF model had the highest performance (accuracy = 92.4%) among the models tested. Svoray et al. [16] used different ML models, namely, SVM, ANN, and decision trees (DT), for predicting the gully erosion in a watershed scale in Israel and compared them with the results from topographic threshold (TT) and analytic hierarchy process (AHP) methods. The results indicated that the ML models produced better performances than the AHP and TT methods. Mao et al. [32] evaluated the soil erosion in the Shiqiaopu catchment, Hubei province, China, using SVM and ANN models. They optimized the parameters of the SVM using the particle swarm optimization (PSO) algorithm. The results indicated that the SVM had higher accuracy in comparison with the ANN model. Rahmati et al. [28] compared the ML models of SVM, ANN, RF, and BRT when predicting the gully erosion susceptibility in the Kashkan watershed, Iran. The results indicated that the performance of the RF and SVM models for predicting the gully occurrences in the watershed were better than the other models.

Due to the advancement of ML models, applying and evaluating novel methods in water erosion studies can help to accurately predict hazardous areas, especially in developing countries where soil erosion data are incomplete. The current study tried to predict water erosion susceptibility using two novel ML models, namely, a weighted subspace random forest (WSRF) and a Gaussian process with a radial basis function kernel (Gaussprradial), for the first time and compared their results with the NB model. Therefore, the primary purposes of this study were: (i) to identify the more significant factors regarding soil erosion through feature selection, (ii) to compare the performance of the novel predictive models (i.e., WSRF and Gaussprradial) with a model previously used for this application (i.e., NB), and (iii) the prediction of the spatial susceptibility of soil erosion induced by water.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The Nur-Rood watershed is located in the southwest of the Haraz watershed, in the north of Iran. The watershed lies within 51◦26′–52◦19′ E and 36◦01′–36◦16′ N (Figure 1). The elevation of the watershed ranges from 732 to 4333 m. There are six rain gauge stations in the region provided the long-term mean annual data from 1976 to 2016 that were were used in this study. The study area is about 1297 km<sup>2</sup> and is located upstream of the Haraz dam. The main application of this dam is to provide drinking and agriculture water for five cities (i.e., Amol, Babol, Babolsar, Nur, and Mahmoodabad) in the Mazandaran province. According to the literature, the watershed generates water with a high sediment load such that it causes a reduction in the dam's capacity [8,33]. Therefore, identifying the hazardous areas can help to control the upstream erosion and aid with providing sustainable watershed management in the Nur-Rood watershed. ′ ′ ′ ′

**Figure 1.** Location of the Nur-Rood watershed, Mazandaran province, Iran.

#### *2.2. Methodology*

The methodology consisted of several fundamental building blocks to ensure the accuracy of the susceptibility prediction. Figure 2 presents the schematic of the methodology workflow from data sampling to the susceptibility prediction. The method consisted of five sections: (i) preparation and collection of the relevant factors for soil erosion modeling; (ii) extraction of the erosion and non-erosion locations by the field observations; (iii) selection of the essential factors using the simulated annealing feature selection (SAFS) algorithm; (iv) water erosion modeling using the Gaussprradial, NB, and WSRF models in the Nur-Rood watershed; and (v) evaluating the models' performance.

**Figure 2.** Schematic representation of the proposed method.

#### 2.2.1. Field Data

It was necessary to know the locations of eroded and non-eroded areas for susceptibility mapping of the Nur-Rood watershed. Therefore, the locations (i.e., x and y coordinates) of 227 area (116 erosion locations and 111 non-erosion locations) were sampled through field surveys to model the water erosion susceptibility based on a binary scale (occurrence/non-occurrence). According to Sajedi-Hosseini et al. [8], the recorded soil erosion areas include different kinds of water erosions (such as sheet, rill, gully, and mass movements).

#### 2.2.2. Predictive Variables

In this study, according to the literature review, 19 relevant factors regarding soil erosion were collected and prepared, including the topographic, hydro-climatic, geological, and land cover factors. The attributes of the factors are presented in Table 1. A brief description of each of the predictive factors is presented afterward.


**Table 1.** Characteristics of the considered factors for susceptibility mapping of water erosion.



<sup>1</sup> Silt loam types of soils with a moderate infiltration rate. <sup>2</sup> Sandy clay loam types of soils with low infiltration rates. <sup>3</sup> Clay loam, silty clay loam, sandy clay, silty clay, or clay with the highest runoff potential. <sup>4</sup> Definition of the geological factors include; TRJs: Dark grey shale and sandstone; Pr: Dark grey medium - bedded to massive limestone; Mm.s.l: Marl, calcareous sandstone, sandy limestone and minor conglomerate; Pd: Red sandstone and shale with subordinate sandy limestone: Odi: Diorite; Tre: Thick bedded grey o'olitic limestone; PZ2bvt: Basaltic volcanic tuff; Tre1: Thin bedded, yellow to pinkish argillaceous limestone with worm tracks; Qs.D: Unconsolidated wind-blown sand deposit including sand dunes; Ebv: Basaltic volcanic rocks; Tra.bv: Triassic, andesitic and basaltic volcanics; Jl: Light grey, thin-bedded to massive limestone; Ek: Well bedded green tuff and tuffaceous shale; K1bvt: Basaltic volcanic tuff; Ktzl: Thick bedded to massive, white to pinkish orbitolina bearing limestone; Pldv: Rhyolitic to rhyodacitic volcanics; Jk: Conglomerate, sandstone and shale with plantremains and coal seams; K2l2: Thick-bedded to massive limestone; Eksh: Greenish-black shale, partly tuffaceous with intercalations of tuff.

#### Topographic Parameters

The topographic parameters included the elevation, slope, aspect, slope length (SL), and curvature (Figure 3). These factors are influential regarding soil erosion velocity [34]. The different elevations (Figure 3a), aspects (Figure 3c), and curvature (Figure 3e) cause different conditions of evaporation, soil temperature, soil moisture, and solar radiation, which have different effects on the soil erosion. Furthermore, slope (Figure 3b) and SL (Figure 3d) affect the runoff velocity and volume, where a steeper slope or a longer SL can increase the soil erosion by water [8].

The topographic factors were produced using a digital elevation model (DEM) with a cell size of 30 m in the ArcGIS 10 software (Environmental Systems Research Institute, Redlands, CA, USA).

#### Hydro-Climate Factors

The hydro-climate factors included the drainage density (DD), distance from the stream (DFS), topographic wetness index (TWI), stream power index (SPI), flow accumulation (FA), precipitation (PCP), rainfall erosivity factor (R), and hydrologic soil group (HSG) (Figure 4). The DD (Figure 4a) is calculated from the sum of the length of all streams in the watershed area. The DD values depend on the permeability and resistance of the surface and deeper soil layers that affect water erosion [8]. Regarding the DFS (Figure 4b), the regions near streams are more susceptible to soil erosion [35]. The DD and DFS layers were created using line density and Euclidian distance tools, respectively, in geographic information system (GIS). TWI (Figure 4c) shows the soil moisture and water-saturated area of the watershed. SPI (Figure 4d) indicates the potential for erosion due to the water flow, in which higher values indicate a higher potential. TWI and SPI were produced using the SAGA GIS 2.0.7 software (SAGA User Group Association, Hamburg, Germany). The flow accumulation (FA) function (Figure 4e) computes the sum of the weight of all accumulated pixels upstream [36], which is most important for showing the water-accumulated pixels that affect the water erosion. The PCP (Figure 4f) and R (Figure 4g) were the climate factors considered to affect soil erosion. Their effects depend on soil attributes such as the soil texture, soil organic matter, and soil structure. The PCP map is produced by the mean annual precipitation of the gauge stations in the study area. The R factor is directly related to the soil erodibility. The best method for calculating it is a direct measurement of soil erosion in plots [37]. However, in this study, according to Takal et al. [38], an empirical equation was used to calculate this factor, as follows.

$$R = 0.0483 P^{1.61} \text{ \AA} \tag{1}$$

where *<sup>R</sup>* is the precipitation erosivity index (MJ·mm·ha−<sup>1</sup> ·hr−<sup>1</sup> ) and *P* is the mean annual precipitation (mm).

The HSG (Figure 4h) indicates the infiltration and runoff generation rates that affect soil erosion. This layer is extracted from the digital soil map of the world [39] and it includes three groups: B, C, and D. Group B has moderately low runoff potential when completely humid. Soils in this group have 50 to 90% sand, 10 to 20% clay, and have sandy loam or loamy sand textures. Water transition across the soil is unrestricted. Group C soils have moderately high runoff potential when completely humid. They have less than 50% sand, 20 to 40% clay, and include sandy clay loam, silty clay loam, loam, silt loam, and clay loam textures. Group D soils have high runoff potential and the infiltration across the soil is very limited [40].

#### Geological Factors

Geological factors include the fault density (FD), lithology, and soil texture (Figure 5). The FD affects infiltration and runoff, which can affect soil erosion. Furthermore, the existence of a fault can accelerate the mass movements [41]. The layer of FD (Figure 5a) was produced in the ArcGIS environment by using the line density tool on the fault layer. The lithology has the greatest effect on erosion control. Erosion depends on the exposed material weathering attributes or the lithology [42,43]. The lithology map (Figure 5b) was taken from a geological survey done by the Iranian department of environment and had a scale of 1:100,000. The other important factor is soil texture (Figure 5c). Porosity and soil texture, along with the soil profile and surface, are the dominant soil attributes that influence soil erosion. An increase in the clay value of the soil causes a decrease in soil erosion [44]. The soil textures of the study area were clay, clay loam, loam, loamy sand, sandy clay loam, and sandy loam (Figure 5c).

#### Land Cover Factors

The land cover factors considered were the normalized difference vegetation index (NDVI), land-use, and distance from road (DFR) (Figure 6). The NDVI (Figure 6a) was extracted from Landsat satellite images for June 2018. The NDVI values range from −1 to 1 [45]. A watershed with a higher NDVI provides higher resistance against soil erosion [9,46]. The land uses of the study area included rangeland, residential, forest, agriculture, and rock (Figure 6b). The land-use map was received from the Iranian Water Resources Management Company (IWRMC). Roads are one of the man-made features that increase the availability of materials for transformation and increase the sediment yield in the watershed. Moreover, roads increase the runoff speed through collecting and concentrating the surface runoff in the given areas (such as near bridges); therefore, faster flows increase the erosion. The DFR layer (Figure 6c) was calculated using the line density tool within the ArcGIS environment.

**Figure 3.** Topographic factors: (**a**) elevation, (**b**) slope, (**c**) aspect, (**d**) slope length (SL), and (**e**) curvature.

**Figure 4.** Hydro-climate factors: (**a**) drainage density (DD), (**b**) distance from the stream (DFS), (**c**) topographic wetness index (TWI), (**d**) stream power index (SPI), (**e**) flow accumulation (FA), (**f**) precipitation (PCP), (**g**) rainfall erosivity factor (R), and (**h**) hydrologic soil group (HSG).

**Figure 5.** Geological factors: (**a**) fault density (FD), (**b**) lithology, and (**c**) soil texture.

**Figure 6.** Land cover factors: (**a**) normalized difference vegetation index (NDVI), (**b**) land use, and (**c**) distance from road (DFR).

#### 2.2.3. Feature Selection

To select the most important factors in the water erosion of soil based on parsimonious objectives from the large number of factors considered, the simulated annealing feature selection (SAFS) model was used. The SAFS method is based on the minimum energy configuration theory, whereby a solid is gradually cooled such that its structure is frozen [47]. Many studies have used this method for feature selection in environmental fields, such as flash-flood hazard assessment [48], dust and air quality evaluation [49], and earth fissure hazard prediction [50]; see Bertsimas and Tsitsiklis [47] for more details of the SAFS method.

In the current research, the SAFS was conducted using the k-fold (k = 10) cross-validation methodology and it was implemented in the Caret package [51] of the R software (4.0.2, R Core Team, Vienna, Austria).

#### 2.2.4. Weighted Subspace Random Forest (WSRF)

Xu et al. [52] suggested a new random forest, namely, the WSRF model, which involves weighting the input variables and afterward opting for the variables that ensure each subspace always includes informative attributes. The WSRF model is implemented as multi-thread processes. This algorithm categorizes very high-dimensional data and sparse data with random forests made using small subspaces. A new variable weighting manner is applied for the variable subspace choice rather than the traditional random variable sampling in the random forest model [53]. More details of the WSRF model are presented in Xu et al. [52] and Zhao et al. [53]. The WSRF model was implemented using the "wsrf" package [53] in the R software using the k-fold (k = 10) cross-validation procedure.

#### 2.2.5. Naive Bayes (NB)

The NB classifiers are a set of assortment algorithms that use Bayes' Theorem. This is a family of algorithms, where every pair of features being categorized is independent of each other; on the other hand, all of them share a common principle. The dataset is categorized into two sections:


According to the primary naive Bayes hypothesis, each element must be independent and equal [54,55]; see Webb et al. [56] for more details of the NB model. The NB model was done using the k-fold (k = 10) cross-validation method in the "klar" package [57,58] within the R software.

#### 2.2.6. The Gaussian Process with a Radial Basis Function Kernel (Gaussprradial)

Gaussian process regression is a vigorous, non-parametric Bayesian method used for solving regression problems and modeling unknown functions [59,60]. It can capture the different relationships between inputs and output variables by applying a hypothetically infinite number of parameters and allowing the dataset to determine the level of complexity via Bayesian inference [61]. The Gaussian process is parametrized using a kernel. One of the benefits of Gaussian process regression is the flexibility in choosing the kernel; furthermore, the different kernels can be combined to perform the regression [59]. In this study, the radial basis function network (RBF) was used to perform the Gaussian process. The Gaussprradial was performed in the R software using the "kernlab" package [62] using the k-fold (k = 10) cross-validation approach.

#### 2.2.7. Model Calibration and Validation

The database, including the predictand and predictors, was randomly divided into the training (70%) and testing (30%) sets. A k-fold (k = 10) cross-validation methodology was used to calibrate the models. The models were assessed using testing datasets after the calibration using the features selected by the SAFS. Here, for the assessment of the models' performances, three classification evaluation metrics were used: accuracy, kappa, and the probability of detection (POD). The models' performances were represented as accuracy percentages. Kappa indicates the probability of agreement by chance using the likelihood of the model classification [63]. The metrics are computed as follows:

$$\text{Accuracy} = \frac{H + \text{CN}}{H + \text{FA} + \text{M} + \text{CN}'} \tag{2}$$

where *H* (the number of hits), *FA* (the number of false alarms), *M* (the number of misses), and *CN* (the number of correct negatives) were computed from a contingency table.

$$\text{Kappa} = \frac{\text{Accuracy} - \text{Pe}}{1 - \text{Pe}} \,\tag{3}$$

where *Pe* is the expected probability of chance agreement [64] that is computed using Equation (4):

$$Pe = \frac{(H + FA)(H + M) + (M + CN)(FA + CN)}{(H + FA + M + CN)^2} \tag{4}$$

The POD is a metric used to quantify the possibility of finding a specific detect. The POD is significantly linked to the subject of risk evaluation and probabilistic analyses of the components' integrity. The POD is the ratio of the correct predicted data to the total number observed occurrences. It ranges from 0 to 1, where 1 indicates a perfect score [49,50]. The metric is calculated using Equation (5):

$$\text{POD} = \frac{H}{H+M}.\tag{5}$$

#### **3. Results and Discussion**

#### *3.1. Feature Selection Results*

A relatively large number of factors, such as elevation, slope, aspect, SL, curvature, DD, DFS, TWI, SPI, FA, PCP, R, HSG, FD, lithology, soil texture, NDVI, land use, and DFR, were used in the current study to predict water erosion. The results of the feature selection using the SAFS algorithm are shown in Table 2. As can be seen, the minimum and maximum selected features were 8 and 14 variables, respectively, in the folds number of 8 (accuracy = 0.84, Kappa = 0.67) and 3 (accuracy = 0.92, Kappa = 0.83). The fold number 6 provided the worst performance (accuracy = 0.74, Kappa = 0.48), whereas the fold number 10 provided the best performance (accuracy = 0.93, Kappa = 0.87).


**Table 2.** Selected factors in each fold using the simulated annealing feature selection (SAFS) method.

According to the 10-fold results, the selected factors should be between the minimum and maximum selected features and should be mostly equal to the mean selected factors across all folds. However, the percentage of selected factors in all folds can be a good criterion for selecting the final variables [48–50]. Figure 7 shows the percentage of selected factors in all folds. Twelve variables had a frequency of at least 50% across all folds. As can be seen, the NDVI with a 100% frequency (F) was selected in all folds. However, the important role of the vegetation and NDVI is obvious and shown in previous studies [8,9,46]. Followed the NDVI (F = 100%), the variables of lithology (F = 80%), R (F = 80%), SL (F = 80%), FD (F = 70%), FA (F = 70%), soil texture (F = 60%), DFS (F = 60%), aspect (F = 50%), curvature (F = 50%), HSG (F = 50%), and DD (F = 50%) were selected.

Although the feature selection has largely not been done in studies on the water erosion of soil, the importance of these selected variables in the water erosion of soil is demonstrated by previous studies, such as those of De Baets et al. [46], Md. Rejaur et al. [9], Sajedi-Hosseini et al. [8], Di Stefano et al. [42], Arabameri et al. [43], Lin et al. [37], Choubin et al. [41], Auzet et al. [44], and Nekhay et al. [35].

**Figure 7.** The percentage of selected factors in all folds.

#### *3.2. Results of Water Erosion Modeling*

The calibration of models was conducted using the "tunelength" function in the Caret R package [51]. The performance results of the three models (WSRF, Gaussprradial, and NB) were evaluated using the three statistics of accuracy, kappa, and the probability of detection (POD), which are presented in Table 3.


**Table 3.** The performances of the models using the testing dataset.

WSRF: Weighted Subspace Random Forest, NB: Naive Bayes.

As can be seen from Table 3, the evaluation of the models' performance indicated that the WSRF model had a higher accuracy (accuracy = 0.91), followed by the Gaussprradial (accuracy = 0.88) and NB (accuracy = 0.85) models. According to Monserud and Leemans [65], the kappa values indicated that all three models were in the "very good" degree of agreement (i.e., 0.70 < Kappa < 0.85) (Table 3). However, like the accuracy, the kappa statistic for the WSRF model (Kappa = 0.82) was more than the Gaussprradial (Kappa = 0.76) and NB (Kappa = 0.71) models. Regarding the POD, the WSRF, Gaussprradial, and NB models showed an equal performance (POD = 0.94) (Table 3).

Generally, the evaluation of the applied machine learning (ML) models in this study indicated an acceptable performance for all the ML models. However, regarding the accuracy and kappa values, the models' performances were ranked as follows: WSRF > Gaussprradial > NB. A direct comparison between the results of this study and previous ones is not possible because the application of the WSRF and Gaussprradial models was undertaken for the water erosion of soil for the first time. However, two novel ML models (WSRF and Gaussprradial) applied in this study indicated a better performance than the NB model that has previously been used in this field. Previous studies have indicated the accurate performance of the NB model in the assessment of soil erosion, such as Weihua et al. [66] and Nhu et al. [67].

#### *3.3. Spatial Prediction of Water Erosion Susceptibility*

After the calibration and validation of the models, the maps of the soil erosion probability were predicted using the values of the pixels throughout the study area. Then, the probability maps were classified into five susceptibility classes of very low, low, medium, high, and very high based on the classification method of natural breaks through the ArcGIS software (Figure 8).

**Figure 8.** Spatial prediction of water erosion using various methods: (**A**) Gaussprradial, (**B**) NB, and (**C**) WSRF.

The area of the susceptibility classes found using each model is presented in Table 4. As can be seen, the Gaussprradial model predicted most of the area in the moderate class (about 450 km<sup>2</sup> , 34.69% of the study area). The sum of the areas for low and very low classes was less than 7% (about 85 km<sup>2</sup> ) (Figure 8A). According to the NB model, more than 65% of the study area (about 850 km<sup>2</sup> ) was located in very high susceptibility zones (Figure 8B). Results of the WSRF model indicated that the classes of the very low, low, moderate, high, and very high susceptibilities covered 11.91% (154.52 km<sup>2</sup> ), 9.76% (126.6 km<sup>2</sup> ), 20.25% (262.64 km<sup>2</sup> ), 28.66% (371.87 km<sup>2</sup> ), and 29.42% (381.7 km<sup>2</sup> ) of the study area, respectively (Figure 8C).


**Table 4.** Area of the susceptibility classes found using each model.

Although the predicted models indicated different areas for each class, there was something in common for all predicted maps. By comparing the predicated maps (Figure 8) with the NDVI map (Figure 5a), it was clear that the susceptibility maps were approximately matched with the NDVI and lithology maps. For example, the green areas in the east of the region on the NDVI map (Figure 6a) had higher values of NDVI, which correspond to lower values on the water erosion susceptibility maps (Figure 7). Furthermore, the higher susceptibly values (Figure 8) corresponded to the TRJ lithology (Figure 5b). TRJs include dark grey shale, claystone, siltstone, and sandstone of the Shemshak formation. In this formation, various kinds of water erosion, such as rill, riverbank, gully, and badland erosions can be seen. This agrees with the SAFS results, which indicated that the NDVI and lithology were the most important variables during the feature selection.

#### **4. Conclusions**

This study focused on the probability of water erosion occurring in the Nur-Rood watershed. Using the SAFS model, the most important factors were selected among nineteen parameters, namely, NDVI, lithology, R, SL, FD, FA, soil texture, DFS, aspect, curvature, HSG, and DD. Based on the performance analysis of the machine learning (ML) models, the two novel applied ML models of WSRF (accuracy = 0.91, Kappa = 0.82) and Gaussprradial (accuracy = 0.88, Kappa = 0.76) displayed better performances than the NB (accuracy = 0.85, Kappa= 0.71) model that has previously been used in this field. The predicted maps created using the ML models indicated the different areas for each susceptibility class but it was obvious that the susceptibility maps were approximately matched with the NDVI and lithology maps (which were identified as the most important variables). One of the main limitations in this study that also occurs in other spatial modeling studies is that different scales are used for the input variables all over the world. Although all of the input variables were resampled into the same spatial resolution, the data collection and sampling of them were not on the same scale; this is an inevitable limitation for the time being. It may be the case that the data availability of the NDVI (30 m resolution) helped this variable to be the most important variable during the soil erosion modeling compared with the other variables (such as the soil dataset and lithology with scales of 1:100,000 or more). Despite these limitations, producing the water erosion susceptibility maps in developing countries can be a useful tool for sustainable management, the conservation of watersheds, the reduction of soil degradation, and alleviating water quality decline.

**Author Contributions:** Conceptualization, A.M.; data curation, F.S.-H. and B.C.; formal analysis, F.S.-H., B.C., and G.R.; investigation, F.S.-H., B.C., F.T., and G.R.; methodology, F.S.-H. and B.C.; project administration, F.S.-H. and F.T.; resources, F.T. and G.R.; software, B.C. and G.R.; supervision, A.M. and A.A.D.; validation, F.T.; visualization, G.R. and A.A.D.; writing—original draft, A.M.; writing—review and editing, F.T. All authors have read and agreed to the published version of the manuscript.

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

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

### **References**


© 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/).

### *Article* **Climate Change Impacts on Soil Erosion and Sediment Yield in a Watershed**

#### **Ching-Nuo Chen 1, \* , Samkele S. Tfwala <sup>2</sup> and Chih-Heng Tsai 3**


Received: 2 July 2020; Accepted: 7 August 2020; Published: 10 August 2020

**Abstract:** This study analyzed the influence of climate change on sediment yield variation, sediment transport and erosion deposition distribution at the watershed scale. The study was based on Gaoping River basin, which is among the largest basins in southern Taiwan. To carry out this analysis, the Physiographic Soil Erosion Deposition (PSED) model was utilized. Model results showed a general increase in soil erosion and deposition volume under the A1B-S climate change scenario. The situation is even worsened with increasing return periods. Total erosion volume and total sediment yield in the watershed were increased by 4–25% and 8–65%, respectively, and deposition volumes increased by 2–23%. The study showed how climate change variability would influence the watershed through increased sediment yields, which might even worsen the impacts of natural disasters. It has further illustrated the importance of incorporating climate change into river management projects.

**Keywords:** climate change; soil erosion; sediment yield; PSED Model

#### **1. Introduction**

Due to the increasing severity of global warming and climate change effect in recent years, extreme hydrographic phenomena have frequently been observed. Climate change has increased precipitation concentration, volume and intensity, which has significantly impacted runoff and soil erosion in many watersheds [1–3]. The sediments generated from watershed erosion are transported to rivers via surface runoff, and they are the main composition of river sediments and a major source of reservoir or river dam sediment deposition [4]. The degree of soil erosion has a significant impact on the evolution of river channels, influencing river stability, flood prevention safety and river remediation planning. Hence, the control of sediment yield is crucial in watershed management, especially since it usually involves high costs. In a review of sediment management strategies in Taiwan and the barriers to their implementation, Wang et al. [5] highlighted how technical barriers are driven primarily by engineering and costs. This was in reference to methods such as construction of upstream sediment structures and hydraulic and mechanical dredging. Moreover, in some regions, climate change is projected to decrease the overall soil erosion potential due to decrease in rainfall [6].

Several studies have focused on the impacts of climate change on precipitation volume, runoff volume, erosion volume and sediment yield. General Circulation Models (GCMs) have been applied to analyze the impacts of climate change on precipitation characteristic [7–9] and river runoff [10,11]. The Soil and Water Assessment Tool (SWAT) model seems to be the most favored by researchers when evaluating the impacts of climate change on flow rate, soil erosion and sediment yield in a

watershed area [12,13]. Thodsen et al. [14] applied the High Resolution Limited Area Model (HIRHAM) regional climate model to investigate the impact of climate change on suspended load transport rate of Danish rivers. In modeling flow rate and sediment yield for high flow-rate rivers under the A2 scenario in the rain season of 2050, Phan et al. [15] showed an increase of 11.4% and 15.3%, respectively. Cousino et al. [16] utilized SWAT to provide hydrological insights for the Maumee River watershed, showing a reduction by 10% in flow, while sediment yield increased by 11%. Most recently, Zhou et al. [17] used the SWAT model to evaluate the impacts of climate change on flow and sediment yield in northeast China. In northern Iran, Azari et al. [18] reported an annual increase of 5% in annual streamflow and more than 35% in sediment yield. Zhang et al. [19] generated climatic conditions for future periods 2020–2039, 2050–2069 and 2080–2099, and their results demonstrated an increase of 13% in sediment yield.

The abovementioned methods do not couple the computations of slope erosion and river sediment transport. Instead, watershed erosion is calculated first based on slope information, followed by watershed sediment yield based on sediment delivery ratio, a river sediment transport model or runoff volume with a flow rate-sediment transport rating curve. The estimation of slope erosion, regardless of applying the Universal Soil Loss Equation (USLE), Revised Universal Soil Loss Equation (RUSLE) [20] or Modified Universal Soil Loss Equation (MUSLE), is carried out by empirical model. Empirical models are established based on erosion data of the initial location via inductive analysis. Therefore, they have application and location limitations. Moreover, the erosion volume estimated by empirical models is the total erosion volume rather than the time-dependent or spatial-dependent erosion volume. Furthermore, watersheds and river systems are complex and sediment boundary conditions for carrying out river sediment transport simulations by the abovementioned models cannot be obtained directly. Typically, sediment transport volume is estimated by the flow rate-sediment transport rating curves established through hydrological observation stations. Nonetheless, sediment transport rates estimated by rating curves have large uncertainties, especially with high discharge values [21].

To understand the impacts of climate change on erosion volume, sediment yield and erosion distribution for a watershed, this study utilized a Physiographic Soil Erosion-Deposition (PSED) model proposed by Chen, Tsai and Tsai [4]. Unlike empirical models, the PSED model has physical mechanisms that enable simultaneous computation of slope erosion and river sediment transport. Additionally, this model does not require the flow rate-sediment rating curve as a boundary condition, thus eliminating uncertainties that come with their adoption. Findings from this study will serve as a useful reference for decision-making authorities in planning appropriate strategies and corresponding measures to prevent water- and soil-related disasters.

#### **2. Materials and Methods**

#### *2.1. Study Area*

Gaoping River, which is the largest river in Taiwan, was used for the analysis. It is located in the southwest of Taiwan and has a total length of 171 km, covering 3257 km<sup>2</sup> . It originates from west of the Central Mountain Range near Yushan. The upstream end is connected to Laonong River, and its main branches include Qishan, Ailiao and Zhuokou (Figure 1). The topographic height varies significantly, descending along the direction from northeast to southwest, with a maximum difference of nearly 4000 m. The area above 1000 m accounts for 47.45% of the total drainage basin, that between 100 and 1000 m accounts for 32.38% and finally, the area below 100 m accounts for 20.17%. The average slope of the river bed is approximately 1/150, with 1/15, 1/100 and 1/1000 for the upstream, midstream and downstream sections, respectively. Additionally, the sections' lengths are 37, 68 and 66 km for the upstream, midstream and downstream sections, respectively.

Time- and location-based precipitation distribution in the basin varies widely. Near the Central Mountain Range it is large (~3400 mm), whereas in the plain and coastal areas it is significantly smaller (~2000 mm). Precipitation is concentrated from May to October, which accounts for 90% of the annual precipitation. Average annual runoff volume is ~8.45 billion m<sup>3</sup> , of which 7.69 billion m<sup>3</sup> (91%) occur in the wet season. The average annual sediment transport volume is 35.61 million tons, with 10,934 tons of sediment transport per km<sup>2</sup> of drainage basin area.

**Figure 1.** The drainage network in the Gaoping River basin.

– The Gaoping River basin is among the worst basins in Taiwan in terms of sediment yield and is highly vulnerable to sediment deposition. Statistical data from the Water Resources Agency of Taiwan, MOEA, estimate a total dredging sediment volume of 94,780,000 m<sup>3</sup> for Gaoping River basin between 2009 to 2014 [22]. The major contributions to such high rates include the high-slope landform and concentrated precipitation. The geological map in Figure 2 shows that the watershed mainly constitutes sand gravel and sandstone, rendering it vulnerable to erosion. Stefanidis and Stathis [6], in their assessment of soil erosion in a catchment, showed that vulnerable geological subsoil and the steep slopes favor the development of erosion phenomena. At Gaoping, about 80% of the annual precipitation is from an average of 3–4 typhoons per year [21], which fall between June and September each year. Precipitation in these events is characterized by high intensity and short duration, leading to enormous volumes of sediment yield. The situation is even worsened by additional factors, such as climate change, which has caused drastic changes in precipitation [23].

**Figure 2.** Geological map of Gaoping watershed [24].

#### *2.2. Long Term Climate and Hydrological Changes*

According to Taiwan's hydrological data According to Taiwan's hydrological data from the past 55 years, there has been a gradual increase in precipitation and typhoon intensity. Both the frequencies of occurrence and the severity of flooding are showing an increasing trend. Figures 3 and 4 show the long-term climate and hydrological data from 1962 to 2016 recorded by the Kaohsiung weather station in the study area. From Figure 3, the annual average temperature in the study area exhibited an increasing trend. Beginning in 1997, the annual average temperature was higher than the average temperature in the past 55 years (24.8 ◦C) and after 2000, the increase of annual average temperature was more significant. Precipitation after 1996 was more significant and most of the annual precipitation was higher than the long-term annual average precipitation of 1770 mm. The deviation between high and low precipitation became larger and the duration of high and low precipitation became shorter (Figure 4).

– **Figure 3.** Changes in annual average temperature recorded by the Kaohsiung weather station, 1960–2016.

– **Figure 4.** Changes in annual precipitation recorded by the Kaohsiung Weather Station, 1960–2016.

–

–

#### *2.3. The Physiographic Soil Erosion-Deposition Model (PSED Model)*

The PSED Model [4,25] is a physical mechanism-based model that was developed by integrating the Geographic Information System (GIS) with a Physiographic Precipitation-Runoff model. It incorporates the effects of slope and river channel erosion, entrainment and deposition on river bed erosion-deposition for a watershed. Based on topography, landform, river system, land use and soil characteristics of the watershed, the PSED Model utilizes GIS to partition the watershed area into non-structural computational cells. The computed cells are then classified into slope cells, river cells and special cells. Esri ArcMap 10.7 was utilized to obtain hydrological and physiographical data within each cell. In addition, the extension modules of ArcMap (spatial analysis, hydrologic model, 3D Analyst and Network Analyst) and their object-oriented programming language were used.

The model consists mainly of two parts; a water flow simulation and a soil erosion-deposition simulation. Water flow simulation calculates the transport of precipitation runoff in the watershed area. The continuity equation of water flow is as follows:

$$A\_i \frac{\partial h\_i}{\partial t} = \sum\_k Q\_{i,k}(h\_i, h\_k) + P\_{ei}(t) \tag{1}$$

where: *t* is time; *A<sup>i</sup>* is area of the *i* cell; *h<sup>i</sup>* and *h<sup>k</sup>* represent the water stage of the *i* and *k* cell, respectively; *Qi,k* denotes the flow rate (discharge) from the k cell into its neighboring *i* cell; and *Pei* expresses the effective rainfall volume per second in the *i* cell, which is equal to the effective rainfall per second in the *i* cell multiplied by its area. Depending on the topography and landform information, the watershed area can be divided into several computed cells and the water level change of each cell should then satisfy the continuity equation of water flow and flow rate simulation, as expressed in Equation (1).

In the soil erosion-deposition simulation for the watershed, simulations for slope cells and river cells were calculated separately. The sediment transport rate and river bed erosion profile of each cell were simulated by using the suspended load equation (Equation (2)), the river bed variation continuity equations (Equation (3)) [26], and the river bed load transport equation.

$$\frac{\partial V\_{si}}{\partial t} = \sum\_{k} Q\_{\text{SC}\_{ik}} + Q\_{\text{sei}} - Q\_{\text{sdi}} + R\_{DTi} \tag{2}$$

$$(1 - \lambda) \frac{\partial V\_{di}}{\partial t} = \sum\_{k} Q\_{SB\_{ik}} - Q\_{sei} + Q\_{sdi} - R\_{DTi} \tag{3}$$

where: *Vsi* is the soil volume of water body in *i* cell (= *A<sup>i</sup>* × *D<sup>i</sup>* × *C<sup>i</sup>* ); *D* is the cell water depth; *C* is the suspended load volume concentration; λ is the porosity; *Vdi* is the alluvium volume of cell *i*; *QSCi*,*<sup>k</sup>* and *QSBi*,*<sup>k</sup>* denote the suspended and the river bed load flow rates, respectively, from the *k* cell into its neighboring *i* cell; *Qsei* represents the entrainment rate of ground surface soil or river bed sediment of the *i*-th cell; *Qsdi* expresses the deposition rate of river bed sediment for *i* cell; and *RDTi* is the precipitation separation rate of the *i* cell.

#### *2.4. Computational Cells*

Figures 5–7 show the Digital Elevation Model (DEM), soil map and land use maps of Gaoping River Basin, respectively. The basin ranges from 0 m to more than 3000 m above sea level and is dominated mostly by forest land at elevations higher than 500 m. Based on the DEM, land use, soil maps, road system maps, and slope maps (Figure 8) the basin was divided by Esri ArcMap 10.7 into 17,635 computational cells for further analysis.

**Figure 5.** Digital elevation model (DEM) of the Gaoping River watershed.

**Figure 6.** Soil map of the Gaoping River watershed [27].

**Figure 7.** Digital land use map of the Gaoping River watershed [28].

**Figure 8.** Slope map of Gaoping River watershed.

#### *2.5. Input Data and Model Setup*

In the basin, there are 28 precipitation monitoring stations established by the Central Weather Bureau; however, they are not evenly distributed. In order to minimize computational errors, we applied the Thiessen polygons method to determine the controlling area of each precipitation monitoring station. The reader should note that the Thiessen polygons were not used to derive the weighted precipitation. Instead, precipitation data for each station were used as the precipitation volume of computed cells in the control area for the same precipitation monitoring station (Figure 9).

**Figure 9.** Effective area for precipitation gauging stations in the Gaoping River watershed.

PSED simulated sediment transport was validated by comparing the flow rate and sediment transport data collected by sediment monitoring stations of the 7th River Management Office of Water Resources Agency, MOEA, which are located at the downstream of the Gaoping River and its tributaries. They include the Shanlin Bridge of Qishan River, Tachin Bridge of Laonong River, Sandimen Bridge of Ailiao River, and Lilin Bridge of Gaoping River (Figure 10).

The future scenario in this study was set from 2020 to 2039 and the corresponding baseline was set from 1980 to 1999. The Taiwan Climate Change Projection and Information Platform Project (TCCIP) provided downscaled precipitation data at 5 km<sup>2</sup> resolution. To process the data, TCCIP relies on 24 GCMs as described in the Intergovernmental Panel on Climate Change (IPCC), Fourth Assessment Report (AR4) [29]. Additionally, IPCC identifies A2, A1B and B1 as the most probable scenarios; hence, this study adopted the A1B-S for analysis, which is regarded as a worse scenario and is similar to the A1B scenario. The worst-case scenario is primarily obtained through subtracting or adding one standard deviation between the estimated values of GCMs from the multi-model ensemble of all GCMs [30]. Monthly precipitation scenario information was further combined with a weather generator to evaluate the impact of climate change on daily precipitation volume.

The baseline was defined by precipitation data from 1980 to 1999. Historical daily precipitation data from the monitoring stations were used as input files for the climate derived models. These were then applied to generate the daily precipitation data representing the future climate change scenario. In addition, the daily precipitation data of the baseline scenario were combined with the precipitation distribution in the watershed area to translate into precipitation profiles to be used by the PSED model.

**Figure 10.** Sediment monitoring stations on the Gaoping River and its tributaries.

#### *2.6. Model Verification*

–

Typhoon Morakot (2009), the most disastrous storm to have hit Taiwan in the last century, was used to validate runoff and suspended load hydrographs from the PSED model. We further compared actual discharge and sediment transport from each hydrological station to simulated data. Simulated and observed flow hydrograph from Lilin Bridge is shown in Figure 11. The peak of the simulated hydrograph coincided with that from the observed data, and the hydrograph shapes are similar. This suggests that the model can be successfully applied.

Since the hydrological monitoring stations along the Gaoping River system do not have suspended load concentration data hydrographs, the historical discharge and suspended load data from sediment monitoring stations downstream of the main branch of the Gaoping River were used to establish the correlation between discharge and sediment transport volume of each hydrological monitoring station. These served as the basis for validating the suspended load concentration hydrograph obtained from the numerical model. The simulated discharge and suspended load transport rate under Typhoon Morakot in 2009 were plotted onto the correlation diagram between the observed discharge and sediment transport rate for the Sanlin, Dajin, Sandimen and Lilin bridge Stations (Figures 12–15). In these figures, points are actual historical measured data. The solid line represents the regression relation between discharge and sediment transport rate for each hydrological monitoring station. It is noted from the figures that the simulated correlation between discharge and sediment transport rate was consistent with the correlation between discharge and sediment transport rate of the actually

observed data for all monitoring stations except Lilin Bridge station, particularly under high discharge. The model slightly overestimated data at this station (Figure 15). Nonetheless, the model also indicated reasonable estimates on suspended load and suspended load transport.

**Figure 11.** Simulated and observed discharge during Typhoon Morakot in 2009 at Lilin Bridge station.

**Figure 12.** Simulated and observed correlation between flow discharge and sediment transport rate at Sanlin Bridge station.

**Figure 13.** Simulated and observed correlation between flow discharge and sediment transport rate at Dajin Bridge station.

**Figure 14.** Simulated and observed relationship between flow discharge and sediment transport rate at Sandimen Bridge station.

**Figure 15.** Simulated and observed relationship between flow discharge and sediment transport rate at Lilin Bridge station.

#### **3. Results**

#### *3.1. Impacts of Climate Change on Erosion Volume and Sediment Yield*

Different precipitation types, rainfall distribution, rainfall intensities and precipitation volume will result in different runoff processes, leading to different erosion volumes and sediment yields in a watershed. The precipitation volumes collected by each precipitation monitoring station in the watershed area for each return period (2, 5, 10, 25, 50, 100 and 200 year return period) under the baseline and A1B-S scenarios were used to calculate the average maximum rainfall intensity and average rainfall for each return period via the controlled area weighted method. The above return periods were selected as they are the standards used for most engineering designs in Taiwan. The precipitation results for the baseline were then compared with the results of the A1B-S scenario. Table 1 shows the comparison of average maximum rainfall intensity and average rainfall between the baseline and selected scenarios of each return period. Average maximum rainfall intensity increased more than average rainfall, suggesting that under the influence of climate change, not only did the precipitation volume increase but also the precipitation intensity.



Simulated erosion volume and sediment yield under the baseline and A1B-S scenarios for various return periods are shown in Table 2. Total erosion volume and sediment yield under the A1B-S scenario for various return periods are greater than under the baseline. The increase in the total sediment yield rate was higher than that of the total erosion volume. The total erosion volume and total sediment yield increases by 4–25% and 8–65%, respectively, when compared to the baseline. This implies that climate change contributed to 15% and 36% increases in soil erosion volume and sediment yield, respectively, when compared to the baseline average.

**Table 2.** Soil erosion and sediment yield increase rates under baseline and A1B-S scenarios for various return periods.


#### *3.2. Climate Change E*ff*ect on Erosion and Erosion Distribution*

Soil erosion simulation results for the studied area were compiled and summarized in Table 3. The least return period indicated a 0.54% increase in area, while for under 200 year return period there was a 2.3% increase. The impacts of climate change on erosion were found to be lower when compared to other areas within the Asian region. Pal and Chakrabortty [31] simulated the impacts of climate change on soil erosion in a sub-tropical monsoon dominated watershed based on a RUSLE model, and found erosion to increase by 33% under a 15 year return period.

**Table 3.** Comparison of the increase rate of soil erosion area increase under the baseline and A1B-S scenarios for various return periods.


Climate change was also shown to cause larger erosion depths (Figure 16). Except in a few selected areas, a greater percentage indicates an increasing trend, and erosion depth increased with increasing return period. This is in line with observations made by Giang et al. [32], who predicted that Asian countries would be among the hardest hit regions globally.

**Figure 16.** Spatial distribution of erosion depth increase under the influence of climate change under (**a**) 10 year, (**b**) 25 year, (**c**) 100 year and (**d**) 200 year return periods.

#### *3.3. Climate Change E*ff*ects on Deposition Volume and the Deposition Distribution*

Deposition distributions under baseline and climate change scenarios for 10 and 100 year return periods are shown in Figures 17 and 18, respectively. High deposition was observed mainly at the confluences; between Gaoping and Qishan River (zone A), Gaoping River and Laonong River (zone B), Ailiao River Gaoshu Bridge and Ailiao weir (zone C) and in the middle and downstream of Gaoping River (zone D). Zone areas are shown in Figure 19. Large deposition at these areas is attributed to widening of the cross sections, low river bed slopes and low flow rates. Deposition in each computational cell was calculated by multiplying the deposition height of a cell by its area, and the total deposition of all cells was simply the summation of the volumes in each cell. A summary of the deposition under the different return periods is shown in Table 4.

Figure 20 shows the spatial distribution of deposition depth increase under the baseline and climate change scenarios. Similar increase patterns were observed between the simulated cases, with a larger deposition depth increase located in the middle and downstream of river channels. Deposition depth increase was highest in the main channel as expected, and increased with increasing return periods indicated by dark green areas in Figure 20.

**Table 4.** Increase in volume of estimated sediment deposition under baseline (1980–1999) and A1B-S (2020–2039) scenarios for various return periods.


**Figure 17.** Deposition depth distribution under 10 year return period for (**a**) baseline and (**b**) climate change scenarios.

**Figure 18.** Deposition depth distribution under 100 year return period for (**a**) baseline and (**b**) climate change scenarios.

**Figure 19.** Location of various river sections in the Gaoping River basin that are vulnerable to soil deposition.

**Figure 20.** Spatial distribution of deposition depth increase under the influence of climate change for (**a**) 10 year, (**b**) 25 year, (**c**) 100 year and (**d**) 200 year return periods.

Gaoping River basin exhibited high erosion volume and sediment yield. According to statistical data reported by the Water Resources Agency [22], 76,870,000 m<sup>3</sup> have been dredged between 2010 and 2013. Table 5 shows the dredged volume of each year, while Figure 21 indicates the location of the dredged site. Dredged locations coincide with high deposition areas computed by the PSED model, as illustrated by Figure 21.


**Table 5.** The actual dredged sediment deposition amounts from 2010 to 2013 in each river [22].

**Figure 21.** Dredged locations in Gaoping River basin from 2010 to 2013.

– – – – In each year, a huge amount of money is needed to carry out the dredging works at Gaoping River. With climate change increasing the deposition rates, not only will there be more pressure on financial resources, but flooding risk is also expected to increase. Hence, appropriate structures and policies should be put in place in order to redress climate change impacts. Proposed strategies include identifying and mapping areas more prone to soil erosion and implementing river management and stability measures. Planning the overall river basin operations and management strategies can effectively control the sediment yield in a watershed, hence reducing sediment deposition in river channels due to soil erosion and eliminating flooding disasters due to limited water passage.

#### **4. Conclusions**

This study applied a numerical model to investigate the impacts of climate change on erosion volume, sediment yield and erosion deposition in a watershed. The results showed that precipitation under the A1B-S climate change scenario would significantly increase soil erosion volume, sediment yield and sediment transport rate. Total erosion volume and total sediment yield in the watershed under the A1B-S scenario for various return periods increased by 4–25% and 8–65%, respectively, from 2 year to 200 year return periods. Climate change further increased deposition volume by 2–23% relative to the baseline and by 13% relative to the baseline average. Deposition was found to mostly occur at the river confluences, river middle and at the downstream end of Gaoping River. The study clearly revealed the adverse impacts climate change is likely to bring to this basin; hence, appropriate conservation measures are suggested.

**Author Contributions:** Conceptualization, methodology, formal analysis, C.-N.C.; data curation, C.-H.T.; validation, writing—review and editing, S.S.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Ministry of Science and Technology (MOST) of Taiwan, grant number MOST-107-2221E-020-006.

**Acknowledgments:** The authors extend their gratitude toward the anonymous reviewers for their valuable comments to improve the quality of the manuscript.

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

#### **References**


© 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/).
