Modelling and Mapping of Soil Erosion Susceptibility of Murree, Sub-Himalayas Using GIS and RS-Based Models

Abstract

Soil erosion is one of Pakistan’s most serious environmental threats. This study used geospatial modelling to identify the distinct zones susceptible to soil erosion in Murree, Pakistan. Using a machine learning technique in the Google Earth engine (GEE) and Google Earth, we identified 1250 soil erosion events. The inventory (dependent variable) was separated into two datasets, one for training (70%) and one for testing (30%). Elevation, slope, aspect, curvature, stream, precipitation, LULC, lithology, soil, NDVI, and distance to road were prepared in ArcGIS and considered as independent variables in the current research. GIS and RS-based models such as WOE, FR, and IV were used to assess the relationship between both variables and produce soil erosion susceptibility maps. Finally, the Area Under Curve (AUC) approach was used to confirm the research results. According to the validation data, the SRC for WOE, FR, and IV were 88%, 91%, and 87%, respectively. The present study’s validation results show that the PRC for WOE, FR, and IV are 92%, 94%, and 90%, respectively. Based on the AUC validation approach, we determined that the FR model had the highest accuracy when compared to the other two techniques, the WOE and IV models. The current analysis and final susceptibility maps of soil erosion could be useful for decision-makers in the future to prevent soil erosion and its negative repercussions.

1. Introduction

Soils are very fragile and vulnerable to environmental impacts [1,2]. Soil erosion is one of the most severe global environmental challenges, with a considerable impact on soil and water quality [3,4]. Soil erosion is the removal and transportation of surface material caused by natural and man-made processes [5]. Soil erosion can be accelerated by human factors such as land use change, deforestation, overgrazing, and over-cropping [6,7]. These actions affect soil properties and decrease the land cover that provides stability to the land.

Soil erosion probability is highly impacted by climatic, topographic, LULC, plant, and other biophysical parameters, such as soil characteristics. These factors contribute to the uncertainty and complexity of the soil erosion phenomenom [8,9]. Because the behavior of biophysical factors in one location may vary from that of another, sustainable erosion management necessitates unique region studies [10,11].

Soil erosion is mostly caused by natural causes such as water and wind. Water affects soil erosion by precipitation intrusion in sediments and runoff on the earth’s surface, resulting in rills and gullies [12,13]. Soil erosion and its negative impacts on the environment and agricultural land have increased dramatically in recent decades in developing countries as a result of poor policy planning for farmland, population growth, and urbanization [14]. In addition to polluting air, water and destroying soil, erosion has a big effect on how the world’s population changes [15,16,17]. Because of this, it’s important to come up with different ways to reduce the effects of erosion.

Therefore, the assessment of soil erosion is a crucial issue [18,19]. Researchers employed many methodologies and models to develop the Soil Erosion Susceptibility Map (SESM) [20]. Previously, scientists employed qualitative techniques such as Weighted Overlay Analysis (WOA) to develop SESM for the region of interest in order to identify soil erosion high-risk areas. Researchers used the WOA approach to analyze several predisposing variables such as rainfall, topography, and soil erodibility index and assign a value to each factor based on its effect on soil erosion [21].

Aside from the WOA method, various models such as the Universal Soil Loss Equation (USLE), the Revised Universal Soil Loss Equation (RUSLE), the Modified Universal Soil Loss Equation (MUSLE), the Erosion Productivity Impact Calculator (EPIC), the European Soil Erosion Model (EUROSEM), and the Water Erosion Prediction Project (WEPP) have been widely used by researchers to produce soil erosion maps. These approaches are employed by researchers for a variety of reasons, including their low cost, simplicity of implementation, and compatibility with GIS [22]. These approaches include drawbacks such as uncertainty and inaccuracy in the model’s real and predicted values [23]. Later, as geospatial technology advanced, researchers used GIS-based quantitative approaches to simulate soil erosion [24]. GIS-based technologies have been approved for computing and mapping watershed zones susceptible to soil erosion, and can help in the management of soil erosion and related dangers [25]. The causal elements, namely the watershed features, should be thoroughly evaluated [26]. A previous study revealed that several erosion-causing elements and approaches have been applied in soil erosion susceptibility mapping. Some of the models used in previous studies include logistic regression [27,28], the Analytical Hierarchy Process (AHP) technique [26], sensitivity analysis approaches [29], weighting overlay [30,31,32], a soft computing method [33], and artificial neural-network evaluation methods [34,35,36].

In the present study, we integrated GIS-based statistical models such as WOE and FR to forecast possible soil erosion areas. For soil erosion mapping, landslide inventory (dependent variable) and elevation, slope, aspect, curvature, drainage, road, precipitation, LULC, and soil (independent variables) are considered influential characteristics [31].

This study is unique in that it investigates the current study area for soil erosion using GIS-based statistical models such as FR, WOE, and IV that have previously been overlooked by researchers. In this study, we also look at the relationship between soil erosion and elevation, slope, aspect, curvature, drainage, roads, precipitation, LULC, and soil. The study adds significantly to the insufficient literature on soil erosion-related risks in Pakistan’s Murree area. The primary goal of this research was to identify and analyse the effects of soil erosion-causing factors and probable soil erosion zones in the Murree region of the Sub-Himalayas.

2. Study Area

The present study was carried out in Pakistan’s Murree area. Murree is a well-known tourist destination and is an attractive hill station in Pakistan [37]. The Murree region is geographically located in the Rawalpindi district of Punjab at 33.9078 °N and 73.3915 °E at a height of roughly 3000 m above sea level, as shown in Figure 1. Geographically, the study area is located to the north and west of Haripur, KP. The Jhelum river is to the east, and the Kotli Satian is to the south of tehsil Murree [38]. It is a subtropical highland zone with annual precipitation of up to 1800 mm and a mean annual temperature of 12.7 °C [39]. The Murree area was created by the collision of two plates, the Indian and Eurasian plates, and is located on the southern edge of the Sub-Himalayan foothills. The research area is often regarded as the most elevated zone [40]. This region is largely inhabited by molasses deposits of sandstone and siltstone, which are impacted by tectonic activity of the Main Boundary Thrust (MBT) [41].

3. Materials and Methods

3.1. Datasets

The datasets used to generate SESM in this investigation included both ground and RS data. For this study, relevant national and international institutions provided ground and RS data on soil erosion factors.

Table 1. Detailed satellite and ground data availability, statement, source, and purpose.

Data	Scale/Resolution	Availability of Data	Data availability Statement/Source	Parameter Maps
Sentinel -2	10 m	Data openly available	The data supporting this study’s findings are openly available in [USGS] at https://earthexplorer.usgs.gov (accessed on 1 July 2022). The spatial resolution is 30 m.	Land Use Land Cover
SRTM DEM	30 m	Data openly available	The data that support the findings of this study are openly available in [USGS] at https://earthexplorer.usgs.gov (accessed on 5 July 2022) Spatial resolution is 30 m.	Elevation, Slope, Aspect, Curvature, and drainage
CHIRPS	0.05°	Data openly available	The data supporting this study’s findings are available in [UCSB] at https://www.chc.ucsb.edu (accessed on 30 July 2022). The spatial resolution of CHIPRS is 0.05° (5.54 Km) and daily gridded.	Rainfall maps
Soil	1:2,000,000	Soil Survey of Pakistan	The data available from FAO and Soil Survey of Pakistan	Soil Texture Map
Geology	1:650,000	(Searle et al.)	Digitized from Geological Map of North Pakistan	Lithological Map
Sentinel 2	10 m	Data openly available	Data import to GEE using machine learning algorithm.	NDVI Map and soil erosion inventory

contains information about the data and sources.

3.2. Methods

In the present study, bivariate statistical models such as WOE, FR, and IV were used in a GIS context to produce a soil erosion map of the study region. Figure 2 depicts the entire technique for the present investigation.

3.2.1. Soil Erosion Inventory Map

The first and most important stage in producing SESM for the research region is to create an accurate and genuine inventory map of soil erosion [42]. The revised soil erosion inventory map was developed on the Google Earth Engine (GEE) using a machine learning technique for Sentinel-2 and Sentinel-1 to identify different types of soil erosion in the investigation. Aside from GEE, many types of soil erosion were observed using the Google Earth platform. In the field survey, the inventory produced by the machine learning algorithm employing the GEE platform and Google Earth was confirmed. As shown in Figure 1, we identified and mapped 1752 soil erosion inventory locations as red layers in the study region throughout this research investigation. The inventory data was separated into two datasets: training (70%) and testing (30%).

3.2.2. Production of the Thematic Data Layers

Elevation, slope, aspect, curvature, stream, precipitation, LULC, soil, and road network were all considered in this research to assess the impact of the aforementioned causative factors on soil erosion.

(i) 

Elevation: The elevation of the land influences the type of plants and the pattern of rainfall [43]. SRTM data with a spatial resolution of 30 m was used to obtain the elevation of the research area. Using the natural break technique, as shown in Figure 3a, the elevation was separated into different classes.

(ii) 

Slope: Slope is an important element in soil erosion because the angle of the slope impacts water penetration and velocity; thus, gentle slopes and flat areas have less soil erosion than steep gradients and slope length [44]. As illustrated in Figure 3b, the slop of the research region was estimated using an SRTM DEM with a resolution of 30 m and reclassified into five classes using ArcGIS 10.8.

(iii) 

Aspect: The variable has a major effect on plant growth and health, evapotranspiration, and moisture content, all of which influence soil erosion [12]. As illustrated in Figure 3c, aspects were obtained from SRTM DEM and categorized into nine classes using ArcGIS 10.8.

(iv) 

Curvature: Details of curvature investigation provide useful information on the geomorphology of the investigated area [45]. As shown in Figure 3d, the curvature of the study region was calculated from the SRTM DEM and reclassified into three groups in ARCGIS 10.8.

(v) 

Drainage: Stream channels have a significant impact on the shear strength of rocks and the erosion of sediment [28]. The stream network was extracted from the SRTM DEM and a five-class buffer was applied to the stream network to understand the relationship between soil erosion and drainage network, as illustrated in Figure 4a.

(vi) 

Rainfall: Precipitation is a crucial causal feature of soil erosion because it causes particle dissociation and accelerates soil erosion [46]. The erosivity of precipitation is a primary driving element for sheet and rill erosion. The present research area’s rainfall map was generated from CHIRPS data using a machine learning approach. As illustrated in Figure 4b, the final rainfall map was categorized into five groups using the ArcGIS 10.8 platform.

(vii) 

UIC: LULC change has a significant impact on soil erosion, independent of other predisposing factors such as climate, soil properties, and topography [47]. The research area’s LULC map was generated in GEE using a machine learning approach employing high resolution data. As illustrated in Figure 4c, the LULC map was exported into ArcGIS and reclassified into five groups to examine the various classes of LULC with soil erosion.

(viii) 

Geology: The composition of different types of rocks, such as volcanic, sedimentary, and metamorphic rocks, as well as their geological formations, have varying effects on soil erosion [48]. Figure 4d depicts the lithological map of the research region, which was digitized from the Northern Geological Map of Pakistan.

(ix) 

Fault: A fault is a lithological feature that has deformed the rocks in the studied region. Figure 5a shows a fault map that was digitized from a geological map of North Pakistan. Following digitization, five buffers were used to estimate the association fault with soil erosion.

(x) 

Soil: The texture of the soil surface is an important predisposing factor of soil erosion since the texture and structure of the soil impact soil resistance to erosion because structural stability and plant cover reduce soil erosion [43]. Figure 5b depicts the soil map of the study region obtained from a soil survey in Pakistan.

(xi) 

NDVI: Vegetation cover reduces soil erosion and may affect the sedimentation pattern [49]. In this work, we calculated NDVI in GEE from 2017 to 2022 using a machine learning method. Following the calculation in GEE, we imported the data into ArcGIS and categorized it into low, medium, and high NDVI classes, as shown in Figure 5c.

(xii) 

Road: Road construction has expanded dramatically in recent years to meet the country’s economic needs, but also affects the hydrologic and topographic trend and has a negative impact on soil erosion [50]. The road map was produced in ArcGIS using ground data from the Punjab Highway Authority, as shown in Figure 5d.

3.2.3. Soil Erosion Susceptibility Mapping Techniques

In the present work, we employed WOE and FR approaches to compute LSM for the study region. The following are the specifics of the approaches stated.

(i) 

Weight of Evidence: WOE is a GIS-based quantitative model that uses the Bayes rule to combine data to estimate the likelihood of occurrences [51]. The WOE approach was initially intended to assess mineral potential mapping using geospatial modelling [52]. The WOE model employs statistical techniques to measure the comparative significance established on the log-linear approach of Bayesian probability. In the case of the WOE model, the positive (W⁺) and negative (W⁻) weights are identified as the most important components. The weight of both causative parameters (B) established on the existence or non-existence of soil erosion (C) of research area is evaluated [53] using the equations below.

$${\mathrm{W}}^{+}=\ln \frac{\mathrm{P}\left(\frac{\mathrm{B}}{\mathrm{C}}\right)}{\mathrm{h}\left(\frac{\mathrm{B}}{\mathrm{C}}\right)}$$

(1)

$${\mathrm{W}}^{-}=\ln \frac{\mathrm{h}\left(\frac{ \overline{\mathrm{B}}}{ \stackrel{⃛}{\mathrm{C}}}\right)}{\mathrm{h}\left(\frac{ \overline{\mathrm{B}}}{ \stackrel{⃛}{\mathrm{C}}}\right)}$$

(2)

L = W⁺ − W⁻

(3)

P is the probability, and ln is the natural log, whereas BB and B¯B¯ are the presence and absence of the causal factors, respectively. Similarly, the presence and absence events are represented by CC and C⁻C⁻ [54]. The presence of the predisposing factor at the occurrence sites is represented by W+. Its magnitude demonstrates the positive relationship between the presence of the causative parameter and disaster incidence. A negative weight (W) indicates the absence of the causal component of soil erosion and indicates the degree of detrimental association. The difference between W+ and W is known as the weight contrast, and demonstrates the spatial link between the causative parameter and the event [55]. In this equation, C stands for both positive and negative spatial associations based on positive and negative values. Figure 6a depicts the final LSM of WOE.

(ii) 

Frequency Ration (FR): The FR model has been widely considered as one of the most effective GIS-based models for determining the spatial relationship between two variables [56]. This methodology is a fairly consistent experimental method for producing LSM for the investigated area [57]. The 2.12 algorithm was used to calculate the FR for each predisposing parameter.

$$\mathrm{FR}=\frac{{\mathrm{N}}_{\mathrm{i}}{\mathrm{P}}_{\mathrm{x}}/\mathrm{N}}{{\mathrm{N}}_{\mathrm{i}}{\mathrm{l}}_{\mathrm{Q}}/\mathrm{Nl}}$$

(4)

where FR is the Frequency Ratio, N_iP_x is the number of pixels in each landslide conditioning factor class, N is the total number of pixels in the study area, N_il_P is the number of landslide pixels in each landslide conditioning factor, and Nl is the total number of landslide pixels in the study area. To create LSI for the area of interest, the following mathematical algorithm was used.

$$\mathrm{LSI}={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{FR}}_{\mathrm{ij}}$$

(5)

FR_ij is the frequency ratio value for ‘j’ class of factor ‘i’, and n is the total number of factors. After performing these steps, the LSM map were produced as shown in Figure 6b.

(iii) 

Information Value: This approach was used in the current study to generate LSM of the research region. This model, which is based on GIS and RS, was used to forecast the geographical connection between landslide inventory and several types of predisposing variables [57]. The following computation formula was used to do this evaluation.

$$\mathrm{W}=\log \frac{\frac{\mathrm{MQox}\left({\mathrm{R}}_{\mathrm{o}}\right)}{\mathrm{M}{\mathrm{Q}}_{\mathrm{o}}\mathrm{x}\left({\mathrm{R}}_{\mathrm{o}}\right)}}{{\displaystyle \sum }\frac{{\mathrm{MQ}}_{\mathrm{o}}\mathrm{x}\left({\mathrm{R}}_{\mathrm{o}}\right)}{{\displaystyle \sum }{\mathrm{MQ}}_{\mathrm{o}}\mathrm{x}\left({\mathrm{R}}_{\mathrm{o}}\right)}}$$

(6)

W denotes the weight of causative factor for landslides. $\mathrm{M}\mathrm{o}\mathrm{x}\left({\mathrm{R}}_{\mathrm{o}}\right)$ shows the landslide number of pixels within class “o”$,{ \mathrm{MQ}}_{\mathrm{o}}\mathrm{x}\left({\mathrm{M}}_{\mathrm{o}}\right)$ the umber of all pixels within class “o”, ${\mathrm{MQ}}_{\mathrm{O}}\mathrm{x}\left({\mathrm{R}}_{\mathrm{O}}\right)$ the total number of landslides pixels, and ${\displaystyle \sum }{\mathrm{MQ}}_{\mathrm{O}}\mathrm{x}\left({\mathrm{M}}_{\mathrm{i}}\right)$) the total number of pixels in study area. The LSI can be generated for study area using following formula.

$$\mathrm{LSI}={\mathrm{W}}_{\mathrm{E}}+{\mathrm{W}}_{\mathrm{S}}+{\mathrm{W}}_{\mathrm{A}}+{\mathrm{W}}_{\mathrm{C}}+{\mathrm{W}}_{\mathrm{LULC}}+{\mathrm{W}}_{\mathrm{L}}+{\mathrm{W}}_{\mathrm{P}}+{\mathrm{W}}_{\mathrm{F}}+{\mathrm{W}}_{\mathrm{R}}+{\mathrm{W}}_{\mathrm{D}}$$

(7)

${\mathrm{W}}_{\mathrm{E}}$ Weight of Elevation,

${\mathrm{W}}_{\mathrm{S}}$ = Weight of Slope, ${\mathrm{W}}_{\mathrm{A}}$ = Weight of Aspect, ${\mathrm{W}}_{\mathrm{C}}$ Weight of Curvature, ${\mathrm{W}}_{\mathrm{LULC}}$ = Weight of Landuse Landcover, ${\mathrm{W}}_{\mathrm{L}}$ = Weight of Lithology, ${\mathrm{W}}_{\mathrm{F}}$ = Weight of fault, ${\mathrm{W}}_{\mathrm{R}}$ = Weight of Road, ${\mathrm{W}}_{\mathrm{P}}$ = Weight of Rainfall, ${\mathrm{W}}_{\mathrm{D}}$ = Weight of stream network. The LSM of IV is shown in Figure 6c.

4. Results

Using machine learning techniques and Google Earth imagery, we identified and mapped up to 1252 eroded spots from Sentinel 1, Sentinel 2, and Landsat 8 images. Figure 1 depicts the regional distribution of soil erosion. This soil erosion includes many types of soil erosion, including rill erosion, gully erosion, ravine erosion, and stream-bank erosion. We completed three bivariate models in this research to obtain SEM of the study region.

Table 2. Statistical analysis of soil erosion inventory and soil erosion causative factors.

Parameters	Class	No of Pixels in Class	No of Landslide Pixels in A Class	W+	W−	WoE	% Pixels in Class	% LS Pixels in Class	(FR)	IV = log (A/B)
Elevation	<900	107,632	422	−1.41	0.20	−1.6	22.46	5.52	0.24	−1.40
	900–1100	71,920	900	−0.24	0.03	−0.28	15.01	11.78	0.78	−0.24
	1100–1300	85,410	1216	−0.11	0.02	−0.13	17.82	15.92	0.89	−0.11
	1300–1500	72,886	2800	0.90	−0.29	1.50	15.21	36.65	2.36	1.51
	>1500	141,251	2300	0.02	−0.00	0.03	29.48	30.11	1.021	0.021
Slope	<10	60,344	799	0.60	−0.12	0.73	12.70	10.46	0.82	0.59
	10–20	217,254	3200	0.71	−1.97	2.68	45.73	41.89	0.91	0.70
	20–30	165,432	2940	0.90	−1.46	2.36	34.82	38.49	1.10	0.89
	30–40	30,547	648	1.08	−0.14	1.22	6.43	8.48	1.31	1.06
	>40	1461	51	1.52	−0.01	1.53	0.30	0.66	2.40	1.52
Aspect	F	36,821	572	−0.02	0.00	−0.03	7.75	5.94	0.76	−0.26
	NE	34,546	1201	−0.03	0.00	−0.03	7.75	7.49	0.97	−0.14
	E	52,413	509	0.80	−0.10	0.89	7.27	15.72	2.16	0.41
	SE	74,371	542	−0.51	0.05	−0.55	11.03	6.66	0.60	0.44
	S	63,720	2243	−0.80	0.10	−0.89	15.66	7.10	0.45	−2.04
	SW	42,743	817	0.81	−0.21	1.02	13.41	29.37	2.19	1.198
	W	47,930	640	0.18	−0.02	0.20	9.00	10.70	1.19	0.19
	NW	61,501	388	−0.18	0.02	−0.20	10.09	8.38	0.83	0.42
	N	60,993	688	−0.94	0.09	−1.03	12.95	5.08	0.39	−0.24
Curvature	Concave	118,738	4170	0.81	−0.51	1.32	24.78	54.60	2.20	0.97
	Flat	262,259	962	−1.48	0.67	−2.16	54.74	12.59	0.23	−1.47
	Convex	98,102	2506	0.48	−0.17	0.65	20.48	32.8	1.60	0.47
Distance to Stream	<25	6968	775	1.66	−0.09	1.87	2.05	6.45	3.01	1.60
	25–50	9840	493	1.55	−0.05	1.60	3.77	14.94	2.87	1.49
	50–100	18,071	1141	1.4	−0.13	1.55	1.45	10.15	1.97	1.38
	100–250	47,475	1447	0.66	−0.11	0.77	9.91	18.94	0.93	0.65
	>250	396,744	3782	−0.52	1.11	−1.63	82.81	49.52	0.60	−1.05
Precipitation (mm/year)	1410–1571	47,185	332	−0.84	0.06	−0.90	9.93	4.35	0.44	−0.83
	1571–1681	98,970	972	−0.50	0.10	−0.60	20.83	12.73	0.61	−0.49
	1681–1773	94,900	1041	−0.39	0.08	−0.47	19.98	13.63	0.68	−0.38
	1773–1881	72,489	1042	−0.11	0.02	−0.13	15.26	13.64	0.89	−0.11
	1881–2035	161,494	4251	0.50	−0.40	1.57	34.00	55.66	2.98	1.58
LULC	Forest	210,650	2000	−0.52	0.28	−0.81	43.97	26.18	0.50	−0.52
	Vegetation	141,362	2327	0.03	−0.01	0.05	29.50	30.47	0.73	0.03
	Barren Land	63,576	2010	0.70	1.54	0.93	13.27	26.32	2.46	1.53
	Urban	61,743	200	−1.61	0.11	−1.72	12.89	2.62	0.20	−1.59
	Water	1796	10	−1.06	0.00	−1.06	0.37	0.13	0.35	−0.95
Lithology	Alluvium	1269	17	−0.18	0.00	−0.18	0.26	0.22	0.84	−0.17
	Quaternary	18,192	263	0.24	−0.01	0.25	2.71	3.44	1.27	0.24
	Muree Formation	398,228	6686	0.04	−0.22	0.26	84.58	87.54	1.03	0.04
	Kuldana	20,835	379	0.14	−0.01	0.14	4.35	4.96	1.14	0.13
	Lora	10,620	73	−0.85	0.01	−0.86	2.22	0.96	0.43	−0.84
	Margalla Hill Limestone	6774	177	−0.89	0.03	−0.92	5.59	2.32	0.41	−0.88
	Kuzagali Shale	683	27	0.93	0.00	1.553	0.14	0.35	2.48	1.56
	Mari Limestone	697	10	−0.11	0.00	−0.11	0.15	0.13	0.90	−0.10
Fault Buffer	<25	2468	30	−0.28	0.00	−0.28	0.52	0.39	0.76	−0.27
	25–50	4000	63	−0.24	0.00	−0.24	1.04	0.82	0.79	−0.24
	50–100	7488	94	−0.24	0.00	−0.25	1.56	1.23	0.79	−0.24
	100–250	18,894	377	0.18	−0.01	0.18	4.15	4.94	1.19	0.17
	>250	446,248	7074	0.00	0.02	−0.02	92.73	92.62	1.00	0.00
Soil	Loamy and Clayey non-calcareous soil	292,459	2538	−0.67	0.64	−1.30	64.17	33.23	0.73	−0.66
	Loamy and non-calcareous soil	186,639	5101	0.64	−0.67	1.55	35.83	66.78	2.79	1.56
NDVI	Low	98,102	3306	0.76	−0.34	1.10	20.65	43.28	2.43	1.52
	Medium	262,259	3670	−0.13	0.14	1.46	55.21	48.05	1.53	1.47
	High	118,738	662	−1.06	0.19	−1.25	25.00	8.67	0.05	−0.9
Distance to Road	<20	6437	92	0.68	−0.01	0.69	1.34	0.63	1.53	0.67
	20–40	6401	64	0.32	−0.005	0.32	1.33	0.51	1.37	0.31
	40–100	18,535	157	0.15	−0.006	0.16	3.86	1.55	1.16	0.15
	100–350	67,753	516	0.04	−0.008	0.05	14.14	3.93	1.04	0.04
	>350	379,972	2650	−0.04	0.14	−0.18	79.30	93.36	0.96	−0.04

displays the results of all bivariate approaches in detail. The following is a detailed description of Table 2.

Altitude, as shown in Table 2, has a significant relationship with soil erosion in the study region. The most important elevation class is 1300–1500 m, followed by >1500 m, while the least sensitive class of soil erosion is 900 m. In the present investigation, the slope parameter was regarded as the essential element. The slope parameter was essential because, as shown in Table 2, as the slope increased, so did the rate of soil erosion. According to the data in Table 2, >40° is the most vulnerable class for erosion phenomena, whereas 10° and 10–20° are the most secure and safe zones. The recent survey’s analytical findings revealed that SE was the most important class of aspect, followed by NW. According to the tabulated data, S and SW are not major classes for soil erosion. According to the tabulated outcomes, a concave structure is the important class for soil erosion, as shown in Table 2, with 1.32, 2.20, and 0.79 values for the WOE, FR, and IV models, respectively. While the less vulnerable class of curvature parameters was flat, with values of −2.16, 0.23, and −1.47 for WOE, FR, and IV, respectively. The statistics of soil erosion and drainage showed that both variables had a positive relationship with each other. The bivariate statistical results suggested that a 25 m drainage buffer is the most vulnerable class for soil erosion, followed by the 25–50 and 50–100 m buffers. According to the table, the WOE, IV, and FR model values for the most vulnerable class are 1.60, 2.87, and 1.49, respectively. The stream study showed that the most vulnerable class of stream parameter was >250, with −1.63, 0.60, and −1.05 for the WOE, FR, and IV models, respectively. In the current research, a rainfall map was generated using CHIRPS satellite data and then divided into five classes to assess the relationship of rainfall parameters with landslide incidents. The precipitation data in Table 2 confirmed the notion that precipitation plays a significant role in soil erosion. The analysis indicated that the precipitation class of 1881–2035 mm/year was the most significant for soil erosion, followed by 1773–1881 and 1681–1773 mm/year. The class 1410–1571 mm/year was the least vulnerable class of current causative factors, with −0.90, 0.44, and 0.83 for WOE, FR, and IV models, respectively. The tabulated results clearly show that each class of the LULC parameter had varying effects on soil erosion activities. Table 2 shows that the barren land of the current research region is the most sensitive class to soil erosion when compared to other classes of LULC. WOE, FR, and IV analytical values were 0.87, 1.98, and 0.68, respectively. Soil erosion in the built-up class of LULC parameters can be regarded as minimal. The present study’s tabulated data demonstrated that faults had no direct influence on soil erosion danger. Table 2 describes the corresponding results of all three models. The findings revealed that the fault’s 150–350 m buffer zone was the most vulnerable to soil erosion, while the 25 m buffer was the least vulnerable.

The present study’s results, as shown in Table 2, highlight that lithology is an important parameter in landslides. The analysis also indicates that Kuzagali Shale is the most vulnerable geological formation to soil erosion, whereas Margalla Hill Limestone is the most resistant. It also indicates that loamy and non-calcareous soils are the most prone to soil erosion. Soil erosion and soil variables had an association value of 1.30, 1.86, and 0.62 for WOE, FR, and IV, respectively. According to Table 1, loamy and clayey non-calcareous soils are less prone to soil erosion. The WOE, FR, and IV for the aforementioned soil parameter class and soil erosion were −1.30, 052, and −0.66, respectively. Table 2 shows that the lowest NDVI class was more vulnerable to soil erosion than the highest NDVI class. In comparison to other classes of LULC, the current research region had the highest risk of soil erosion. WOE, FR, and IV analytical values were 0.87, 1.98, and 0.68, respectively. Table 2 also shows that WOE, FR, and IV for road association with soil erosion, which revealed that the road network is not a major contributing element in soil erosion. The 20 m class of road buffer was the most vulnerable to soil erosion, followed by the 20–40 and 40–100 m classes, while >350 was the least sensitive type of road network to soil erosion, with −0.18, 0.96, and −0.04 for WOE, FR, and IV, respectively. A validation method is required in research projects to examine the correctness of models’ analytical outputs since, without it, the research results are not dependable enough to implement [58]. We employed the AUC approach to assess the legitimacy of the models. SRC and PRC were used in this strategy to generate graphs. An AUC >0.90 indicates accurate prediction; 0.80 AUC 0.90 suggests acceptable prediction; 0.70 AUC–0.80 similar efficacy prediction; 0.60 AUC 0.70 denotes bad prediction; and AUC 0.60 reveals unreliability [59,60]. Figure 7 depicts the validation graph of WOE, FR, and IV using 30% of the soil erosion testing data. The greatest AUC value indicated the most dependable model results, and vice versa.

The AUC graph for the WOE model showed that the plotted values for SRC and PRC were 0.88 and 0.92, respectively. Figure 7a,b shows the SRC and PRC for the WOE model, whereas Figure 7c,d shows the SRC and PRC for the FR, which were 0.91 and 0.94, respectively. The SRC and PRC values for IV were 0.87 and 0.90, respectively. Figure 7e,f shows the AUC for the IV model.

5. Discussion

The characteristics of soil have a significant impact on ecology and is linked to the a country’s socioeconomic development. The soil-water relationship is a key factor in numerous socioeconomic trends, with a focus on mitigating ground degradation and renovating territory [58,59]. Researchers have discovered that soil erosion influences LULC changes, soil fertility, and terrain erosion globally. As a result, preventing land degradation is a significant concern for a country’s environment, society, and economy. As a result, proper preparation and procedures for soil erosion susceptibility mapping were used to reduce soil erosion in the present study area. We employed quantitative GIS models to assess the incidence of soil erosion, its relationship with causal variables, and to identify locations prone to soil erosion. GIS-based models are viable ways of producing Soil Erosion Susceptibility Maps to decrease the risk of soil erosion. According to a literature review, 11 characteristics, including elevation, slope, aspect, curvature, stream network, precipitation, LILC, lithology, soil, NDVI, and road network, have been used to quantitatively analyze the relationship of soil erosion with causative variables. In this study, we used and assumed a comparative assessment of three GIS-based statistical models for soil erosion susceptibility mapping: WOE, FR, and IV. Based on AUC validation, we determined that the FR model outperformed the other two strategies, the WOE and IV models, in terms of accuracy.

Elevation has a considerable impact on gully erosion [61]. According to our results, the class of elevations between 916 and 1065 m was the most important for soil erosion, with a 2.95 value. A 0.00 value association of both variables indicates that elevations greater than 1065 m have no effect on soil erosion [62]. The most significant elevation class according to our data was 1300–1500 m, followed by >1500 m, while the least susceptible class for soil erosion was 900 m. Areas of less than 900 m above sea level in our research area are covered in vegetation, whereas areas at 1300–1500 m are barren.

Slope and LULC are important factors in soil erosion [60]. Steep slopes and barren terrain play important roles in soil erosion. Concave curvature, which is associated with soil erosion, is the most vulnerable class [62]. Our investigation also revealed that the concave structure was most susceptible, with values of 1.32, 2.20, and 0.97 for the WOE, FR, and IV models, respectively.

The stream network has a substantial impact on soil erosion [63]. The soil erosion association with the stream network at <313 m is 1.44, whereas research has revealed that an association is <1 for distances more than 1237 m from the stream indicate a lower chance of soil erosion [61]. The stream parameter had a major influence on soil erosion in our study. According to the bivariate statistical data, the 25 m drainage buffer was the most susceptible to soil erosion. According to the table, the WOE, IV, and FR model values for the most vulnerable class were 1.87, 3.01, and 1.60, respectively, while the least vulnerable class of stream parameter was >250, with values for the WOE, FR, and IV models of −1.63, 0.60, and −1.05, respectively. Soil erosion is mostly caused by precipitation, with more precipitation causing more soil erosion [43]. In Pakistan’s Murree, we promote knowledge that highly precipitated areas are more prone to soil erosion. Although precipitation is an important element in soil erosion, it is heavily impacted by the geology of the region [64]. We observed from that rainfall caused soil erosion in areas with loose lithology. Plants have an important role in soil erosion management because their roots compress the soil, reducing the effects of rainfall and weathering on the soil surface. As a result, a region with thick vegetation will be less vulnerable to soil erosion than an area with sparse vegetation and bare ground [65]. Soil erosion decreases when NDVI increases, and vice versa. According to our observations, the class 0.016 had a value of 1.64. This demonstrates the importance of plants in controlling soil erosion [66]. Our current research showed that NDVI is an important parameter in soil erosion analysis, with low NDVI indicating a substantial erosion process and high NDVI indicating reduced vulnerability to soil erosion. The low NDVI class of the current research relationship with soil erosion had values of 1.10, 2.43, and 1.52 for WOE, FR, and IV, respectively. The highest NDVI class, on the other hand, had a lower connection with soil erosion, with −1.25, 0.05, and −0.9 for WOE, FR, and IV, respectively. The road network is an important anthropogenic component that influences soil slope instability [67]. Erosion has a strong connection with the road network, with a greater value of 7.45 at 300 m and a weaker correlation with increasing distance, with a value of 0.39 when the distance exceeds 1200 m [62].

A comparison of the current study with previous study validated results since all predisposing variables associated with soil erosion exhibited the same behavior regarding the significance of soil erosion activity. This research not only identified vulnerable zones but also identified the most important and relevant causal components in the soil erosion process. According to our results in the present research drainage network, precipitation, soil, and geology are the most relevant parameters for soil erosion in our study region. Barren land, NDVI, elevation, and slope are all elements that influence soil erosion in GIS-based models of the research area. The soil erosion susceptibility maps produced could be useful for different national and provincial agricultural and energy-related organizations in reducing the consequences of this hazard in the area, as well as planning for major occurrences such as floods and landslides. To minimize soil loss in the Murree area, the Punjab government has urged the Rawalpindi administration to concentrate on forest and advanced-level conservation strategies.

6. Conclusions

This research aimed to use geospatial modelling to draw a map of soil erosion susceptibility in the study region. Three GIS and RS-based models (WOE, FR, and IV) were examined for soil erosion susceptibility mapping. The causal parameters of soil erosion were investigated based on their relevance in the literature. Elevation, slope, aspect, curvature, stream, precipitation, LULC, lithology, soil, NDVI, and distance to road were all investigated in ArcGIS and analyzed for their relationship with soil erosion utilizing WOE, FR, and IV approaches. Furthermore, the present study’s results were considered using the AUC approach, which revealed that the SRC for WOE, FR, and IV were 88%, 91%, and 87%, respectively. Validation findings showed that the PRC for WOE, FR, and IV were 92%, 94%, and 90%, respectively. These results demonstrate that the FR approach with the specified parameters is the most appropriate for soil erosion. The outcomes show that the three models used in this study, namely WOE, FR, and IV, are extremely suitable techniques. For regional land, data for all current research approaches is easily accessible, cost-effective, and time-consuming. The approach used in this research may be suggested for regions influenced by similar natural (topographic, geologic, hydrologic, and climatic) and man-made activities, such as road networks, and for mapping soil erosion susceptibility. The final susceptibility map of soil erosion developed by three GIS and RS-based models can help policymakers reduce soil erosion in the study region.