2.3. Annual soil erosion modeling and mapping (2005 - 2015)
The inherent soil erosion potential depends on rainfall, soil type, LCLU, and terrain characteristics, which are represented by various factors, such as, rainfall erosivity, soil erodibility, slope length and slope steepness, LCLU management, and conservation practices. At the national scale in Pakistan, in this study RUSLE (equation 1) (Renardet al. , 1991) was used to estimate soil erosion at 1 km spatial resolution for 2005 and 2015.
\(A\ =\ R*K*L*S*C*P\) (1)
where, A = soil erosion (ton/ha/year), R = rainfall erosivity factor (MJ mm/ha/hr/year), K = soil erodibility (ton·ha·hr/ha/MJ/mm), L = slope-length factor (dimensionless), S = slope-steepness factor (dimensionless), C = cover management factor (dimensionless), and P = conservation practice factor (dimensionless).
Equation 2 (Renard & Freimund, 1994) was used to calculate rainfall erosivity factor (R) for 2005 and 2015 using annual mean precipitation data of the respective year. Over the diverse geographical and topographical study areas, several researchers utilized equation 2 to calculate rainfall erosivity factor (Howland et al. , 2018; Lamyaaet al. , 2018; Uddin et al. , 2016, 2018).
\(R\ =\ 0.0483P\hat{}1.610\) (2)
where, R = rainfall erosivity factor in Megajoule millimeters per hectare per hour per year (MJ mm/ha/hr/year) and P = annual precipitation in millimeters (mm).
In the FAO soil map, fourteen soil classes were found which behave differently towards soil erosion and hence contain varying soil erodibility (K) factor values that explain the degree of soil erodibility of each soil type. Therefore, for this study, the K factor values were obtained from published research articles, which were assigned to fourteen soil classes to obtain a soil erodibility factor map, used in the RUSLE (Table 1 and Table S1).
In this study, we used equations 3 and 4 to calculate the slope-length (L) and slope-steepness (S) factors respectively which are adopted from variously published studies in the South Asia region (Uddin et al. , 2016, 2018) The combined LS-factor describes the effect of topography on soil erosion. As temporal topography data was unavailable for 2005 and 2015, a one-time computed L & S spatial layers were used for annual soil erosion estimations and mapping for 2005 and 2015 over the entire Pakistan.
\(L\ =([\lambda/22.13]\hat{}m)\) (3)
where, L = slope-length factor (dimensionless), λ =grid size (1 km x 1 km) or field slope length, and m = 0.5 for slopes > 4 %; 0.4 for 4% slope; 0.3 for slopes < 3% (Wischmeier & Smith, 1978)
\(S=\ ((0.43\ +\ 0.30\ s\ +\ 0.043\ s\hat{}2)/6.613)\) (4)
Where, S = slope-steepness factor (dimensionless), s = slope derived from DEM in percentage (%)
In this study, MODIS’s seventeen LCLU classes defined according to the International Geosphere-Biosphere Programme (IGBP) classification scheme were clamped/recoded to six LCLU classes (‘Cropland’, ‘Forest land’, ‘Grassland’, ‘Other land’, ‘Settlements’ and ‘Wetlands/Snow cover’) which were defined by Intergovernmental Panel on Climate Change (IPCC) (Table S2). From several published articles on RUSLE, we synthesized average cover management factor (C) values for each IPCC LCLU class, which were assigned to 2005 and 2015 LCLU maps (Table 2 and Table S3).
In this study to determine the soil cover management (C) and conservation practice factor (P) values, we compiled and used average values from previously published studies in Pakistan, India, China, and on a global scale (Table S4). Each LCLU class was assigned an average C and P factor values given in Table 2, to get cover and use management and conservation practice factor (P) maps for 2005 and 2015.
2.4. Quantification of soil erosion factors andannual rate of soil erosion changes (2005 - 2015)
To generate soil erosion maps for 2005 and 2015, the spatially computed soil erosion factors: ‘rainfall erosivity (R)’, ‘soil erodibility (K)’, ‘slope-length (L)’, ‘slope-steepness (S)’, ‘cover management (C)’ and ‘conservation practice (P)’ were inserted in equation 1.
The final annual soil erosion (ton/ha/year) maps were reclassified and renamed into four classes; Low (<1), Medium (1–5), High (5-20), and Very high (>20). In this study, each factor and final annual soil erosion maps were spatially and statistically compared between 2005 and 2015. For the annual soil erosion change assessment, we used a cross tabular or change matrix method in which, diagonal values show the stability of soil classes, while omission and commission values indicate a shift in area or percentage between the four soil erosion classes (adopted from Gilani et al., (2015) used for LCLUC assessment). In this study, using the change matrix method, soil ‘loss’ (commission) and gain (omission) statistically and spatially were reported between 2005 and 2015.
Latitude wise Pakistan extends from 24° N to 37.5° N. In this study across latitude, mean soil erosion (ton/ha/year) in 2005 and 2015 were plotted to understand eroded intensity linking to spatial distribution of soil erosion areas.