Spatial Quantification of Cropland Soil Erosion Dynamics in the Yunnan Plateau Based on Sampling Survey and Multi-Source LUCC Data

1. Introduction

Soil, one of the earth’s most precious and threatened resources [1], provides humans with far more than foods but also a large variety of services like biomass production, water filtration, nutrient transformation, carbon storage, habitat and terrestrial biodiversity maintenance [2,3]. However, most of the soil resources worldwide are only in poor health conditions, and accelerated soil erosion induced by inappropriate human activities and related land use change, is the primary driver behind the problem [4,5]. Soil erosion refers to the complex process of soil materials detachment, transportation, and deposition by external erosive forces. It causes on-farm impacts of reduced soil fertility and productivity [6,7], and also leads to greater off-site costs such as muddy flooding, sedimentation, water pollution [8] and at stake are the global biogeochemical cycles [9,10]. Literature and runoff plots observations have demonstrated that cropland is the main source of soil loss [5,11,13], and the soil erosion rates from which can be orders of magnitude greater than natural soil formation speed and those from other land use types [13,14,15], especially in mountainous areas. It is estimated that about 80% of the agricultural land around the world is suffering from serious erosion problems, with an unsustainable mean annual cropland erosion rate of about 30 t/(ha·a). Consequently, more than one-third of cropland has been vanished in the last decades due to soil erosion [16]. To monitor and assess the impacts of soil erosion and make strategies to deal with them, mapping up to date quantitative information on cropland erosion rates at regional scale is essential but also very challenging [17,18], since for most areas worldwide the observed erosion data is woefully inadequate.

On the one hand, it is difficult to gain insight into the spatial pattern of soil erosion without specific soil erosion rates and hotpots information, confusion raises in the allocation of soil erosion mitigation programs and priorities, the formulations of policies, and the effectiveness evaluation of soil conservation measures [19]. Besides, knowledge gaps will be generated in critical fields like climate change, landslides and flood prediction, carbon mitigation scenarios and earth science modelling, and the well-known polices of SDGs, CAP, UNCCD and IPBES will be out of focus [20]. Despite soil erosion modelling and prediction have received considerable attention from governments and scientists for more than seven decades [21,22], with various empirical, conceptual and physically-based models and approaches been developed to measure, estimate and monitor soil erosion from field to landscape scales [23,24,25,26]. Yet, most models are only applicable to micro-scales like field plots, hillslopes and small catchments, and are difficult to applied to large scales due to the spatial heterogeneity of soil erosion affecting factors [28], scale issues [14,27], model limitations and applicability [14,29,30], especially high demand for model input data [20,31]. Since most empirical soil erosion models are established based on the scales of plots and hillslopes with certain applicable conditions and scopes [31]. For example, the most widely applied soil erosion prediction USLE-type models, are originally developed at the plot scale for agricultural lands (gradient less than 18 %) based on the “unit plot concept” of a 22.1 m long, 1.83 m wide plot, with a 9% slope with up-and-down hill tillage [20]. The major limitation of soil erosion modelling for any given area is the microscopic process involved is less considered, and it is difficult to acquire up-to-date soil erosion information like crop rotation, terracing, mulching, contouring and hedgerow planting at large scales, especially in fragmented mountainous landscapes. When upscaling the models to large scales, the input variables or parameters of the models generally simplified, huge uncertainties may lead to extrapolation error, and the reliability of the results is often questioned. Currently, the contradiction between the relatively low resolution of available input data and the high resolution required for runoff-erosion processes is the major obstacle to overcome for large scale dynamic quantification of soil erosion [29]. For mountainous areas, poor data availability, timeliness and data quality has been the biggest obstacle in mapping and visualizing soil erosion rates at large scales [32].

To date, limited by the over parameterization of physical models and poor datasets available, large scale soil erosion assessment methods are generally based on empirical models, and can be divided into two categories of sampling surveys and remote sensing assessments [27,33]. (a) Sampling survey refers to the method of allocating samples within a region according to certain proportion and rules, field investigation on erosion features and parameters is then conducted, soil erosion models will be further applied to quantify soil erosion rates or conditions, and statistical methods will be used to estimate the overall soil erosion patterns of the region finally. The typical examples are the National Resource Inventory (NRI) conducted in the United States [34], the National Soil Erosion Survey in China [35], the EUSEDcollab network [36] and the gully erosion monitoring based on the Land Use/Cover Area frame survey in Europe [37]. (b) Another category is remote sensing based assessments with simplified models. Since large-scale application of complex models is challenging, as the availability of high-resolution remote sensing images increases, large-scale estimation using empirical models such as USLE/RUSLE, becomes feasible due to the relatively simple input data. Compared with field investigation, satellite remote sensing is characterized by timely, affordable data, uniform data over large areas, real-time information acquisition and regular revisit wide view field [38], and has been widely applied in soil erosion modelling and mapping [39]. Particularly, efforts have been put into direct soil erosion detection [19,40], and parameters of rainfall erosivity estimation [41,42], soil related properties derivation [43], topographic factors extraction [44], cover-management (C-factor) and support practices (P-factor) [45,46,47,48], specific soil conservation measures mapping using high-resolution imageries [49,50,51], as well as those large-scale soil erosion assessments with raster layers operations [52,53]. The biggest advantage of sampling survey is it provides reliable soil erosion rates, and large-scale spatial patterns of soil erosion status can be achieved by combining with statistical principles [54,55], but the field measurement and investigation of indices is labor-intensive and high cost. Remote sensing-based methods allow rapid and efficient soil erosion assessment even in areas that intensive field investigation is a challenge, but more in a qualitative or semi-quantitative way. Although high resolution imagery like SPOT 5, IKONOS and Quikbird ensure high quality data in erosion mapping, their utility remain hindered as large area imageries are also unaffordable for most countries [40].

Previous studies [56,57,58,59] confirmed that land use/cover change (LUCC) is the primary cause of accelerated soil erosion under climate change scenarios, and is the most direct and intuitive reflection of the interaction between human activities and the soils on earth surface [60]. For most cases, the key to prevent soil erosion is to change various unreasonable land uses to a sustainable mode which in line with the principles of sustainable development, such as the projects of returning farmland to forest/grassland and converting slopes into terraces in China. Meanwhile, compared with the inversion of soil erosion indices, remote sensing application in LUCC monitoring is the field with the most complete and mature technology. At large scales, by integrating most commonly usedLandsat series data imagery, relevant studies [61,62] also reveal the long-term impact of land use change on soil erosion, and provides the suitable information necessary for assessing soil erosion intensity, but due to the lack of field-based soil erosion investigation data, the estimated dynamic results are generally potential soil erosion risks without exact dynamic soil erosion rates. Obviously, more detailed field experimental data that accurately quantify soil erosion rates are needed.

In 2010 through 2012, the Ministry of Water Resources of China (MWRC) conducted the first ever and only field-based National Soil Erosion Survey (NSES) in history by using sampling survey and the Chinese Soil Loss Equation (CSLE) [63,64,65]. Those detailed onsite investigated provides abundant information on soil erosion rates at land parcel scale, which reduces the uncertainties in soil erosion modeling and prediction. As soil erosion is a dynamic process demands constant monitoring to obtain up-to-date information on its spatial pattern [40]. The combination of advantages between both sampling survey and remote sensing is definitely the potential solution for the dilemma of large-scale soil erosion rates quantification. The objectives of this paper, is thus to quantitatively assess the cropland soil erosion dynamics induced by the long-term cropland change in mountainous areas, with perspectives from soil erosion field investigation and the LUCC scenarios based on time-series satellite images.

2. Materials and Methods

2.1. Study Area

Yunnan Plateau (20°8′–29°16′N, 97°31′–106°12′E) is a low latitude highland region situated in the southwest border of China, it covers a total area of about 3.83 ×10⁵ km², and borders the Himalayan Range, Myanmar, Laos and Vietnam (Figure 1). Mostly mountainous in character, 94% of the province is dominated by mountains and plateaus, only 6% of which interspersed with small, scattered valley basins. The landscape tilts downward from the northwest to the southeast, and elevation ranges from 76m to 6740 m above sea level, with an average altitude of 2000 m [65]. Affected by the Indian and East Asian monsoons and air masses from the Qinghai-Tibet Plateau, the region has a subtropical plateau monsoon climate with substantial variation, though relatively mild due to the elevation. Mean annual precipitation varies from 600 in dry-hot valleys to 2300 mm in the southern and western mountains, with over half of the rain occurring between June and August, while the dry season (November to April of the following year) accounts for only 20% or less of the 1100 mm annual precipitation. Annual average temperatures in the winter and summer are 6–8 °C and 19–22 °C, respectively [33,35]. Soils in Yunnan are generally rich in clay and formed under high precipitation and temperature conditions such as Acrisols, Cambisols and Luvisols. It is particularly worth mentioning that Yunnan has long been recognized as the hotspot of biodiversity in China. Taking vegetation resources as an example, tropical, seasonal, subtropical evergreen broad-leaved, temperate coniferous forests and meadow steppes can all be found in the region [69]. The landscape is deeply dissected by six major rivers of the Irrawaddy, Nu (Salween), Lancang (upper reach of the Mekong), Jinsha (upper reach of the Yangtze), Honghe and Pearl. The province also owns the largest sloping cropland area in China, as the limited basin areas already been fully utilized, soil erosion pressure on remaining land resources is extremely high.

2.2. Data Sources

2.2.1. Sampling Survey and Primary Sample Units (PSUs)

The sampling survey in the NSES was conducted using a non-equal probability sampling method (Figure 2a). In view of the dominant erosive force of soil loss and the integrity of county boundaries, the whole country was firstly divided into water erosion region, wind erosion region, freeze thaw erosion region and other regions co-occurrence erosion (Figure 2b). Generally, a uniform national gridding was employed for water erosion, which has the greatest implications. Four layers of grids were set according to the grid size, including the county-level (40 km×40 km), township level (10 km×10 km), the control area level (5 km×5 km) and Primary Sample Units (1 km×1 km). The grids were divided according to the Gauss-Krüger Projection zoning method, which divides the country into 3°interval geographical zones (a total of 22 zones). In each zone, the Y-axis direction is divided into grids on both sides based on the central meridian, and the X-axis direction is divided into grids on both sides based on the equator. PSUs for water erosion was allocated with four sampling densities of 4%, 1%, 0.25% and 0.0625% [66,67]. For the plain area, PSU was a single 1 km×1 km grid, while for the mountainous area, PSU was a small watershed of 0.2–3.0 km². PSU is a small geographical area with a fixed location and a certain area that can express the basic characteristics of soil erosion and show spatial heterogeneity in soil erosion factors (especially soil conservation measures). For areas of glaciers, permanent snowfields, deserts, swamps, large lakes and reservoirs where no water erosion occurs, and high altitudes exceeding 4800 m with less human activities, water erosion was generally less considered [66,67]. Based on the erosion characteristics, human disturbance and accessibility of each province, a total of 33, 966 PSUs were determined nationwide in the NSES (Figure 2d). Compared to the NRI, NSES in China is actually an area sampling survey rather than a point sampling survey.

PSUs were the main object of field investigation and data collection in the NSES. Each PSU was then divided into pieces of land defined as a land parcel that share the same land type and conservation measure. Data gathers were then trained to investigate each land parcel and collect information concerning soil erosion factors of the CSLE with uniform standards. Specifically, by using high-precision topographic maps as field survey base maps, detailed information such as land parcel number, fraction vegetation coverage, canopy density, land use, vegetation type, engineering measures and crop rotation patterns were acquired. Suffering from serious water erosion problem, Yunnan attracted a lot of attention in the NSES, and both the number of PSUs and land parcels selected were the largest among all provinces in China, with 2,817 and 20,155 respectively (Figure 2c). All these field investigated data was provided by Beijing Normal University (the technical support unit of the NSES) [65,69].

In addition to the data of field investigation, critical data involved in calculating annual average soil erosion rate was obtained nationwide as follows: (a) Daily erosive rainfall (greater than 12 mm) data for three decades; (b) Digital DEMs of 1:10,000 scale for each PSU to extract slope gradient and length; (c) More than 10,000 records of soil profile data and soil types at 1:500,000 scale retrieved from the Second National Soil Survey, and observed unit plot and cropland plot data; (d) time series multi-spectral reflectance data of Sentinel-2, Landsat TM/ETM/OLI images prepared to revise vegetation indices during 2000-2020; (e) High resolution satellite imagery in Yunnan of GF-1 (2m), GF-2 (1 m), GF-7 (0.8 m) images and Beijing-2 (resolution of 0.65 m) at different periods to optimize land use map.

After standardized data processing, each PSUs contain 7 raster layers (spatial resolution of 10 m) of soil erosion factors in the CSLE model, soil erosion modulus was then computed by using raster multiplication operation and statistical methods, soil loss for each land use type was finally evaluated at land parcel, PSU, province and national levels. Figure 3 presents the detailed information on a random PSU of small watershed and each land parcel within it. For this study, the most critical data is the multi-year average soil erosion rate of various land use parcels obtained based on field surveys, which is the basis for comprehensive analysis of soil erosion distribution, area, ratio and intensity.

2.2.2. Land Use/Cover Change (LUCC) Dynamics

Initialized in the latter half of the 1990s, the National Land Use/Cover Database of China (NLUD-C) [70] has been applied to land use/cover related researches as the most well-known LUCC database in China for decades. By drawing boundaries and labelling attributes for each LUCC polygon based on Landsat TM/ETM/OLI images (resolution of 30 m), we have updated the database for several periods using interactive interpretation method. NLUD-C include datasets of land use status and LUCC dynamics of China with a 5-year interval, and land use types are classified into six first-level categories and 25 corresponding second-level classes.

Compared to numerous LUCC products retrieved using automatic classification and change detection methods, the biggest advantage of NLUD-C lies in the professional knowledge, uniform interpretation symbols and unified image acquisition phases before and after, which ensures its high-accuracy information of land attributes. Although the OA of NLUD-C is significantly higher than other similar products, distortion and inaccurate patch boundaries still exists in mountainous areas since manual geometric correction is performed at different periods. At present, remote sensing LUCC data with higher resolution may be a better choice for soil erosion modeling as it reduces model uncertainty, but is also characterized by shortcomings of high image acquisition costs and short time series. To optimize the accuracy of land use data, we employed NLUD-C as the main LUCC basis, supplemented by 7 publicly accessible used non-homogeneous LUCC datasets, which are three long term 30 m resolution products of GLC_FCS30 [71], CLCD [72] and GlobeLand30 [73], as well as four short term 10 m resolution products of ESRI_LandCover [74], ESA_WorldCover [75], CRLC [76] and Dynamic World from Google [77]. The detailed information of those land use datasets is listed in Table 1. In this study, land use maps of Yunnan in 2000, 2005, 2010, 2015 and 2020 were prepared for further analysis.

2.3. Methods

2.3.1. The CSLE Model

By adapting parameters of the Universal Soil Loss Equation (USLE) to China, Liu et al. developed the CSLE model based on measured data from Chinese unit plot and numerous plots modified to unit plot [64]. As the official model of the Ministry of Water Resources of China (MWRC) for soil erosion assessment, the differences of CSLE compared to USLE are modifications that made for crop system, management, practice, soil type, rainfall pattern and topography in China.

CSLE is a model used to estimate annual soil loss by sheet and rill water erosion for a given combination of the soil erosion affecting factors. The major advantage of CSLE is that it is more in line with the topographical conditions and the actual situation of soil conservation measures in China, while the two factors of cover management and support practices (C, P) in USLE were modified into three factors of biological measures (B), engineering measures (E) and tillage measures (T) [64,65]. Five dimensionless factors of slope length, slope gradient, biological measure, engineering measure and tillage measure are used to modify the soil loss determined by the dimensional rainfall erosivity factor and soil erodibility factor in this model.

The principal equation of CSLE can be expressed as follows:

$$A=R\times K\times L\times S\times B\times E\times T$$

(1)

where A is mean annual soil loss with unit of t/(hm²·a); R is the rainfall and runoff factor or rainfall erosivity, MJ·mm/(hm²·h·a); K is the soil erodibility factor, t·hm²·h/(MJ·hm²·mm); L is the slope length factor and S is the slope steepness factor; B is the biological measure factor; E is the engineering measure factor; T is the factor of tillage measures, B, E, T factors have a unitless range of 0–1 and the smaller the value is, the better the soil conservation effect of a certain measure is.

The specific calculation methods of R, K, L and S factors are described in detail in our previous work [60,78] and related literatures [79,80] and listed as follows:

$$R={\displaystyle {\sum }_{k=1}^{24}{R}_k}$$

(2)

$$R_k=\frac{1}{N}{\displaystyle {\sum }_{i=1}^N{\displaystyle {\sum }_{j=0}^m(\alpha \cdot P_{i,j,k}^{1.7265})}}$$

(3)

$$W{R}_k=\frac{R_k}{R}$$

(4)

where k represents the 12 months in a year, R_k is the average rainfall erosivity in the k-th month (MJ∙mm∙ha⁻¹∙h⁻¹∙a⁻¹), N refers to the time series, the term α is a value of 0.3937 for the warm season and 0.3101 for the cold season, P_i,j,k is the actual erosive rainfall (≥12 mm) of the j-th day in the k-th month in the i-th year, m is the number of days with erosive rainfall in the corresponding month. WR_k is the ratio of average rainfall erosivity in the k-th month to the average annual rainfall erosivity, which reflects the seasonal distribution of rainfall erosivity.

$$K=\left[2.1\times {10}^{-4}{M}^{1.14}(12-OM)+3.25(S-2)+2.5(P-3)\right]/100$$

(5)

$$M=N_1(100-N_2)$$

(6)

$$M=N_1(N_3+N_4)$$

(7)

where N₁ (particle size: 0.002–0.1 mm) is the percent of silt (0.002–0.05 mm) plus very fine sand (0.05–0.1 mm), N₂ (< 0.002 mm) is the clay fraction, (100–N₂) (0.002–2 mm) represents all soil fractions other than clay, OM is the soil organic matter content (%), S is the soil structure code and P is the soil permeability code.

$$S=\{\begin{array}{l}10.8 sin \theta +0.03 \theta \le 5°\\ 16.8 sin \theta -0.50 5°<\theta \le 10°\\ 21.9 sin \theta -0.96 \theta >10°\end{array}$$

(8)

$$L_i=\frac{({\lambda }_{out}^{m+1}-{\lambda }_{in}^{m-1})}{[({\lambda }_{out}-{\lambda }_{in})\times {22.13}^m]}\{\begin{array}{l}m=0.2 \theta \le 1°\\ m=0.3 1°<\theta \le 3°\\ m=0.4 3°<\theta \le 5°\\ m=0.5 \theta >5°\end{array}$$

(9)

where L_i is the slope length factor of the i-th pixel, λ_out and λ_in are the pixel exit and entrance slope lengths and m is the slope length exponent depending on the slope.

$$FVC=\frac{NDV{I}_{max}-NDV{I}_{soil}}{NDV{I}_{veg}-NDV{I}_{soil}}$$

(10)

$$B=\frac{{\displaystyle {\sum }_{i=1}^{12}{B}_i{R}_i}}{{\displaystyle {\sum }_{i=1}^{12}{R}_i}}$$

(11)

$$B_i=\{\begin{array}{l}1 FVC=0\\ 0.6508-0.3436\lg FVC\times 100 0<FVC\le 0.783\\ 0 FVC>0.783\end{array}$$

(12)

where NDVI_max refers to the regional maximum NDVI; NDVI_veg is the NDVI value of the pure vegetation pixels; NDVI_soil is the NDVI value of the pure bare soil pixels; B_i is the B-factor of the i-th month. The relationship between the FVC and B value was compiled using Equation (12). The ET-factors in the CSLE are mainly collected based on the field survey, with reference to corresponding values from runoff plot experiments by local experts and an extensive literature review analysis on a national effort.

2.3.2. Non-Homogeneous Voting and LUCC Optimization

The rapid development of remote sensing technology provides a key technical approach to obtain comprehensive information on large-scale land use/cover distribution and changes. In recent years, scientists worldwide have incorporated image processing methods to interpret and analyze remote sensing images, and produced numerous LUCC products with different spatial resolutions.

Our accuracy evaluation of datasets and extensive literature review demonstrate that: (a) automatic image classification and change detection approaches can provide satisfactory results only when applied to certain land use types with homogeneous color and texture, such as water bodies, built-up land, and bare rock [70]. Due to the difference in satellite sensors, processing methods and classification systems, the description ability of earth surface conditions of different LUCC products are also different, especially in fragmented mountainous areas like Yunnan, affected by cloudy and rainy weather, high quality image availability is poor and the difference can be even more significant, and the reliability of data is often questioned. For example, despite the 8 datasets we mentioned above share similar definition of cropland, the cropland area of Yunnan in 2020 provided by them are totally different. The cropland areas of the 4 long time series 30m-datasets are 6.74×10⁴ km² NLUD-C, 3.24×10⁴ km² for GLC_FCS30, 8.39×10⁴ km² for CLCD and 10.99×10⁴ km² GlobeLand30, and the cropland areas of the 4 recent 10m-datasets are 3.58×10⁴ km² for ESRI_LC, 5.09×10⁴ km² for ESA_WC, 9.77×10⁴ km² for CRLC and 2.50×10⁴ km² for Dynamic World, respectively. For Yunnan, based on our field surveys, CLCD significantly underestimates the impervious area, GLC_FCS30 generally underestimates the cropland area, GlobeLand30 significantly overestimates the cropland area, ESRI_LC and Dynamic World significantly overestimates the impervious area, CRLC misclassifies woodland into cropland and ESA_WC misclassifies divides impervious area into bare land. Misclassification, omission and high confusion degree between grassland and shrub can be found in almost all the datasets. High resolution dataset does not always provide more reliable information, significant differences can be found at the same resolution. Obviously, the overall accuracy of global or national-scale LUCC datasets in local areas needs to be verified, since the verification methods and reference data for datasets are different, the independently completed accuracy assessments cannot be directly compared [81]. Therefore, critical evaluation of the suitability and optimization of LULC products based on application purposes should be made before use.

For plateau mountain areas like Yunnan, long term high-resolution (higher than 10 m) LUCC data is inaccessible and impractical so far. Here, we proposed a non-homogeneous voting method to modify and optimize the NLUD-C data. Specifically, we selected time series geometrically fine-corrected Landsat images from October to February in each year to ensure the uniformity of multiple phases. Based on NLUD-C, the interpretation symbols and the professional knowledge, the non-homogeneous voting was then performed to analyze the classification consistency. For high consistency degree area, we kept the first-level type to reduce uncertainty, while for high confusion degree area, especially for basin valleys and urban periphery areas with frequent human activities and dramatic changes, high-resolution remote sensing images were used to assist in precise visual interpretation. Finally, the obtained time series LUCC data combines the advantages of multiple datasets and avoids the respective shortcomings. Figure 4 presents the workflow of the non-homogeneous voting method for modifying the LUCC data, and Figure 5 shows our revised NLUD-C results in a typical inter-mountain basin, by using high-resolution remote sensing images, we also compare them to three other time series land use products of different underlying surfaces.

3. Results

3.1. Soil Erosion Pattern of Yunnan Based on Sampling Survey and Field Investigation

3.1.1. Investigated Land Parcel Basics of Yunnan in the NSES

According to the sampling survey, a total of a total of 20,155 land parcels in 2,781 PSUs were investigated on site in Yunnan. The average area of the PSUs was 34.83 ha, and the average patch area of the land parcels was 4.79 ha, which is in consistency with the model requirements. The measured mean slope gradient for all the land patches was 20.23° and with a range of 0−82.8°, while the mean slope length is 47.73 m and a range of 0−233.9, the average soil erosion rate of all the land use parcels is 17.02 t/(ha·a). Table 2 lists the basic information in land parcels of the first-level NLUD-C types. Obviously, in terms of spatial distribution, woodland has the largest average patch area and number (a total of 10015 land parcels) and is the dominant land use in the landscape. The average patch area of cropland is 2.77 ha, and the average slope reaches 17.88°, which is very prone to soil erosion.

3.1.2. Soil Erosion Rate Variations under Different Land Use Types and Topography

Table 3 lists the various soil erosion factors in the CSLE and the multi-year average soil erosion rates of the NLUD-C second-level land types. Woodland has an overall higher average rainfall erosivity than other types, with a value of 3570.37 t·hm²·h/(MJ·hm²·mm), followed by grassland, cropland, built-up land, water bodies and unused land. The highest rainfall erosivity was found in paddy fields and garden plantations (classified as woodland in NLUD-C) as they are mainly distributed in southern Yunnan. Additionally, the lowest lower R values were found in sparse grass and bare land, which are mostly distributed in the northern part of Yunnan and the dry-hot valleys of the six major rivers. Since the main soil types are highly sticky, the K values are relatively close and small in all regions. In terms of terrain factors, except for built-up land, water bodies, and irrigated land (basically cropland in flat areas with irrigation conditions), all other land types are characterized by great slope steepness and short slopes, LS values are generally high for sloping cropland, gardens and unused land. In the CSLE model, the vegetation factor B is mainly for woodland and grassland. The vegetation coverage of woodland in Yunnan is significantly higher than that of grassland (lower B factor value), while the impact of vegetation on soil erosion is incorporated into the tillage measure factor T, and intercropping and rotation are the dominant tillage measures. Out of 6714 cropland parcels investigated, 52.25% of them were adopted with engineering measures. Overall, cropland without ECMs contributes 83.51% of the total soil loss from cropland with only a land area of 47.75%.

As Yunnan is the only province with a fully plateau mountainous landscape in China, the impact of terrain on soil erosion rate is crucial, and it is often the most important factor in water erosion models. With the change of terrain conditions, the soil erosion rates of the five major land use types except water bodies show varying degrees of difference in different slope gradient and length zones, and soil erosion is much more sensitive to changes in slope steepness than to changes in slope length (Figure 6). For cropland, the average annual soil erosion rate is 40.47t/(ha·a) and the erosion ratio of 70.11% (4707 out of 6714 land parcels with soil erosion rate higher than soil loss tolerance of 5 t/(ha·a)), far exceeding other land types (28.57% for woodland, 35.42% for grassland, 10.27% for built-up land and 60.53% for unused land). Rain-fed dry land is the main type of cropland (5644 out of 6714 land parcels), and the soil erosion rate is 45.34 t/(ha·a) and erosion ratio of 81.48% (4599 out of 5644 land parcels with soil erosion rate higher than soil loss tolerance), which are also significantly higher than those of paddy fields and irrigated cropland. For woodland and grassland parcels, higher values ​​are mostly found in sparse vegetations and garden plantations (divided as woodland but retains the attributes of cropland). Due to low coverage and intensive disturbance, high soil erosion rates of built-up land were mainly found in mining areas. Bare soil suffers from some of the highest soil erosion rates among all NLUD-C second-level land types, but it occupies a small area and does not contribute much to the total soil.

To better present the spatial pattern of cropland erosion rates, we mapped the spatial distribution pattern of cropland erosion rates at the PSUs scale (Figure 7) according to the classification standards of the MWRC. As can be seen from the figure, soil erosion rates in downstream areas of the six major river basins are generally higher than upstream areas, and aggregation occurred. Cropland land erosion rates in the central flat area was much lower. Considering the influence of terrain, and the similarity in soil properties, planting system and rotation patterns, the difference can be attributed to the fact that rainfall erosivity in the downstream areas are significantly higher than those in the upstream areas.

3.1.3. Impact of Engineering Conservation Measures on Cropland Soil Erosion

For a long time, scientists and soil conservancy departments have attached great importance to soil conservation in cropland, but limited by data availability, little is known about the effectiveness of engineering conservation measures (ECMs) at large scales. ECMs refer to the measures applied by changing micro terrain conditions to intercept runoff, increase soil infiltration or crop production, such as the level terrace, sloping terrace, fruit tree pit, check dam, intercepting drains, diversion canal, etc. In the NSES, one of the major tasks is to conduct detailed field surveys on the type, distribution, quantity and area of ECMs. We also quantify the impact of ECMs on soil erosion based on literature, standard runoff plot observations and data collected under natural or artificial simulated rainfall conditions. To further understand the impact of ECMs on cropland erosion, we also mapped the spatial distribution of ECMs on PSUs that contain cropland parcels, and analyzed their impact on soil erosion rates at different slope steepness intervals (Figure 8a). Obviously, except for the cultivated areas with flat terrain and low soil erosion rate in the central part, ECMs are spatially distributed throughout the province. The proportion of cropland parcels with ECWs shows a decreasing trend as the slope steepness increases. This is mainly because large scale adoption of remediation measures on steep slopes in plateau mountainous region is difficult, unaffordable, and may cause disaster like landslides. Besides, farmers are more encouraged to return steep slope cropland to forest/grass rather than terraced fields based on the Grain for Green Policy.

Figure 8b reveals the effectiveness of soil conservation measures in mitigating cropland erosion beyond just a reduction in the magnitude of erosion rate. For the sloping cropland in the province, soil erosion rate of both cropland with and without ECMs increases as the slope gradient increase. Based the field investigation, cropland parcels can be found in 1863 out of 2871 PSUs, the average cropland soil erosion rate with ECMs is 12.14 t/(ha·a), while the average cropland soil erosion rate without ECMs is 67.25 t/(ha·a), and the rate difference is more than five times, which is based on the premise that low erosion rate cropland in the central part does not equipped with ECMs. For steep slope zones, this difference can be even bigger and reaches a gap of 100 t/(ha·a) for steepness class of 30−35°, indicating that the effect of ECMS also decreases as the slope gradient increases. Apparently, the allocation of ECMs largely affects the spatial distribution pattern of cropland erosion rates. In terms of soil loss prevention, almost all the croplands in Yunnan demand ECMs to control the soil erosion within tolerant rate. But more attention to focused on assessing the difficulty, cost and effectiveness of the treatment (including soil productivity) to determine priorities where projects should really be adopted.

3.2. Land Use Change Dynamics in Yunnan from 2000 to 2020

By following the update mode of NLUD-C data, we optimized land use change dynamics data of Yunnan from 2000 to 2020. Figure 9 shows the corresponding land use maps of year 2000, 2005, 2010, 2015, and 2020. Woodland (67.79%−68.44%), cropland (16.56%−18.51%), and grassland (10.43%−11.64%) are the dominant land types. Woodland can be found in the whole province, cropland is mainly distributed in the central, northeast and southeast parts, and grassland is concentrated in the Jinsha River Basin and the dry-hot valleys of other river basins. The proportions of water bodies (1.19%−1.22%), built-up land (0.82%−3.33%) and unused land (0.05%−0.06%) are relatively low.

Based on the land use transfer matrix, the land use transformation process and change dynamics from 2000 to 2020 were presented in Figure 10. For the 2000-2005 and 2005-2010 periods, the total land change areas were 1.62×10⁴ km² and 1.73×10⁴ km² respectively, of which the conversion between cropland and woodland were the dominant change. For the 2010-2015 and 2015-2020 periods, the total land change areas were 2.06×10⁴ km² and 1.97×10⁴ km² respectively, and the cropland-grassland transformation dominated the periods, the one-way conversion from cropland to built-up land was also noticeable.

The main LUCC characteristics of Yunnan Province from 2000 to 2020 can be summarized as follows. The built-up land area continues to rise, with a net increase of 9451.57 km² and a rate of 300.39%. The cropland area continues to decrease, with a net decrease of 7461.83 km² and a rate of −10.55%. The grassland area continues to decrease, and the net decrease was 4603.39 km², and a rate of −10.36%. The woodland area increased by 2482.63 km², but because of the large area base, the change rate was the lowest at 0.95%. The water bodies and unused land were relatively stable, with increasing areas of 92.75 km² and 36.75 km², respectively. Although the change trend of LUCC in each period is relatively consistent, significant differences were also found in region, quantity and main change scenarios. It should be highlighted that the land conversion area related to cropland accounts for 74.02% of all the transformation scenarios, and is the most significant type of land use conversion.

To finally achieve dynamic quantitative monitoring of cropland erosion rate and soil loss, as a focus, we analyzed the conversions between cropland and other land types spatially. For each time period, most of the cropland change area was converted to woodland, grassland or built-up land, with less conversion to water bodies and unused land, while the cropland reclamation area was also mostly from grassland and woodland (Figure 11). The cropland−woodland conversions were mainly found in the river basins of Lancang, Nu, Irrawaddy and southern parts, while the cropland−grassland conversions mainly occurred in the Jinsha River Basin and the central grassland−cropland−built-up transition zones. In the past 20 years, the cropland loss area of the six major river basins was ordered as Lancang(21.35 km²)＞Honghe(13.86 km²)＞Nu(13.45 km²)＞Pearl(10.76 km²)＞Jinsha(9.70 km²)＞Irrawaddy(5.50 km²), while the cropland loss ratio was ordered as Nu(28.76%)＞Irrawaddy (22.95%)＞Lancang (20.34%)＞Honghe(10.20%)＞Pearl(5.74%)＞Jinsha(4.66%). Apparently, northwestern Yunnan are suffering serious cropland degradation and loss problem.

As can be seen from Figure 11e that the areas with the largest net decrease in cropland are concentrated in the western Yunnan, and mainly distributed in the Irrawaddy River Basin, Nu River Basin and Lancang River Basin. The spatial pattern of cropland decrease trend in the three river basins are basically consistent with each other. Specifically, the net cropland decrease from 2000 to 2005 was mainly found in the upper reaches of the basins, as time changes, the center of the net cropland decrease gradually moves from the upstream to the downstream areas. Nearly 90% of the basin area is experiencing cropland loss. Especially in the Irrawaddy River Basin, the net decrease was mainly concentrated in the middle and upper reaches before 2015. From 2015 to 2020, the entire basin was experiencing a net decrease in cropland.

3.3. Cropland Soil Erosion Dynamics in Yunnan from 2000 to 2020

In previous work, we optimized the input parameters of the CSLE with the annual average values of vegetation coverage and erosive rainfall data from 2000 to 2020, and recalculated annual average soil erosion rate for each land type of the PSUs. Since NLUD-C update and NSES originate from different programs and with different application purposes, the land use/cover classification system of NLUD-C is much simpler than field investigation. Thus, we unified the classification system and retained the six first-level land types of NLUD-C, and determined the average soil erosion rate using an area percentage weighted average method in PSUs. The high-precision soil erosion rate for different land types in PSUs was further interpolated to the corresponding control area of PSUs using a nearest neighbor interpolation method. Finally, net soil erosion rate and soil loss changes for different land use conversion scenarios was calculated (Table 4), and then up-to-data quantitative soil erosion dynamic information on cropland can be acquired by incorporating with the revised LUCC dynamics area during the past two decades (Figure 12). Considering the average soil erosion rate of each land type, for different LUCC scenarios in Yunnan, cropland change induced increases and decreases in soil erosion rate and soil loss exists during the four periods. Conversions from cropland generally reduced the soil erosion intensity, while slope cropland reclamation was the main LUCC type that intensifies of soil erosion.

Figure 12 presents a representative area located in Jinsha River basin, each grid is the control area (grid size of 5 km × 5 km ) of a certain Primary Sample Unit. Through sampling surveys, we calculated the multi-year average soil erosion rate of all the land use types in the corresponding grids, as well as the erosion rate change of various conversion scenarios. As can be seen that both cropland abandonment and expansion were found in the region, there is a significant soil erosion rate variation for the same LUCC conversion among grids despite they are very close in spatial distance. During this period, 1.66 km² and 0.91 km² of cropland abandonment and expansion occurred in the representative area, with the largest area of cropland converted to built-up land, followed by area converted to water bodies. In the dynamic zones, soil erosion rate caused by different cropland transformation types are also different. Returning cropland to woodland and grassland, converting cropland to built-up land reduced the soil erosion rate dramatically.

When land use type changes, the soil erosion rates and intensity will change accordingly. From 2000 to 2020, as cropland related LUCC transformations dominant the landscape, changes of soil erosion rates and soil loss are also significant in corresponding areas. In the past two decades, the amount of soil loss due to the cropland transfer out was −1.28×10⁸ t, while the amount of soil loss due to the cropland reclamation was 0.96×10⁸ t, and the net change in soil erosion caused by cropland transformations was −0.32×10⁸ t, with a decrease ratio of 12.12% of the total cropland soil loss (Figure 13).

The cropland area continues to decline over the past two decades, and most of it was converted into woodland, grassland and built-up land. The conversions from cropland to woodland and grassland account for the largest proportions of cropland erosion reduction, with 43.75% and 34.64% respectively. The conversion from cropland to built-up land accounts for 21.36% of the reduction. Conversion from cropland to woodland has been a major contributor to the decline in soil erosion in cropland areas. As for the soil erosion increase, about 99% of which was caused by cropland reclamation, and conversion from grassland to cropland contributed 56.18% of the total soil erosion increase.

The conversions between cropland and other land uses also have different effects on soil erosion amount in different periods. Among them, compared with the mutual conversions with woodland and grassland, the conversion between cropland and built-up land is more often a unidirectional change. Due to the larger net change area, the net reduction in erosion generated from cropland to built-up land is also the largest, and the reduction in soil erosion during the 2010-2015 period accounts for nearly 40% of the total cropland-built up land change scenario, as it is also the period with the fastest urbanization speed. For the four periods of 2000-2005, 2005-2010, 2010-2015 and 2015-2020, the mutual conversions between cropland and woodland are relatively stable in each period. Most of the soil erosion reduction caused by the conversion from cropland to grassland is concentrated in 2010-2020, while the soil erosion increase caused by the conversion from grassland to cropland shows a decreasing trend year by year.

4. Discussion

Constant soil erosion rates measurement and observation at large scales have proven to be extremely challenging and unrealistic. Based on field sampling surveys, the CSLE model and LUCC data, we proposed a rapid monitoring method to extrapolate cropland soil erosion rates and soil loss from point to surface in mountainous areas. The field investigated 20,155 land parcels share same standards in data quality and all of them meet the USLE-type empirical model requirements in size and scale (less than 150 ha). The LUCC data was further improved using a non-homogeneous voting method, with steps of accuracy assessment, consistency analysis and standardization of the classification system. To facilitate decision-making, we provided continuous distribution information on cropland erosion rates, hotspots and soil loss amounts. The soil erosion rates of each land type are in good consistency with the reported values in literatures [6,30]. Apparently, when choosing a soil erosion model, one should pay more attention to model strengths, limitations and application scope. If the input data does not meet the requirements, the results produced by over-parameterization and scaling extrapolation are often less reliable than those given by a simple model.

Under climate change and land use change scenarios, cropland erosion and degradation are a mutually promoting process. In areas with extremely high biodiversity like Yunnan, the implementation of policies such as returning farmland to forest/grass is of great value in controlling soil erosion and protecting habitats and biodiversity. However, our research shows that a cropland area of 7461.83 km² (−10.55%) has vanished during the past 20 years. which is an extremely shocking number, and the originally small cropland area per head of population is continuing to shrink, the newly reclaimed slope lands are often accompanied with severe soil erosion rates, directly threatening local food security. The real threat should be noted here is that more and more land is becoming unfarmable due to high soil erosion rates. It is estimated that if the current soil erosion rate in China continues, the food production will decrease by 40% in the next 50 years [14]. Moreover, the rising global population demands intensification of agricultural production to meet food demand, which is expected to increase by 50% in 2030 and possibly a doubling in 2050 [1]. If current population growth speed and soil erosion rates continue unchecked, humankind may eventually lose the ability to feed itself in the future barring unforeseen scientific advances [5]. The regulation of sloping croplands is extremely difficult in mountainous areas, as the croplands are fragmentedly distributed on steep slopes. According to local statistics [82], the annual average cropland land area with newly treated erosion control measures in Yunnan is 31.69 km². It will take more than 1,000 years and 180 billion yuan to complete the regulation of unmeasured sloping cropland and, and this is assuming that each cropland can be managed without considering the difficulty of governance. Considering the cropland loss speed, urgent action is needed to face the threat of cropland soil erosion with shared understanding by considering collaboration and interrelationships among stakeholders, different roles (e.g. scientists, governments, farmers, environmentalists).

Field observation of soil erosion is always closer to the truth than the modelling results, and it is the most vital part of scientific investigation. However, most regions around the world have the problem of under-representation of observational data. Currently, remote sensing is instrumental for investigating, evaluating, monitoring and understanding the spatial extent and rate of soil erosion due to the advantages of large coverage area, short revisit period monitoring. High resolution imageries provide high quality data and less uncertainties in soil erosion mapping, but their utility remain hindered due to the acquisition cost. As the spatial, hyperspectral and temporal resolution continuously increase, it sheds more and more light on small scale heterogeneity, and most of the limitations of large-scale soil erosion modelling may eventually dissipate in the future. With a robust framework of sample density and samples, remote sensing applications in large scale dynamic soil erosion mapping and monitoring will be very promising.

We proposed a combination method of point (PSUs) and surface (LUCC data) for quantitative soil erosion assessment in a large region, the work depended greatly on the detailed data collection in the field. The NSES was the first ever national soil erosion investigation using based field investigation, which ensures the accuracy of the input data. However, the quality and representativeness of the data for areas with low sampling density and missing sample information requires more in-depth evaluation.

5. Conclusions

Long-term, quantitative large-scale cropland erosion rates information is vital for agricultural planning and management, but long been hindered by data availability and model limitations. Taking the CSLE as monitoring tool, by integrating a large number of field sampling surveys and LUCC remote sensing data in the national surveys, we proposed a long-term time series dynamic method of monitoring cropland soil erosion rates and soil losses, and conducted an application research in the Yunnan Plateau with complex terrain conditions. Different from previous studies, this study was conducted based on a large number of field surveys and remote sensing for improving model input data and reduces the uncertainties. The results showed that:

(1) The average soil erosion rate and erosion ratio of cropland are significantly higher than other land use types, and huge spatial erosion differences were within each land use type. In addition, soil erosion rates are generally more sensitive to slope than slope length for all land uses. Soil conservation measures adopted in croplands are highly effective in controlling soil erosion and changed the spatial pattern of soil erosion significantly.

(2) In the past 20 years, due to the Grain for Green Policy, population growth and rapid urbanization expansion, the area of cropland and grassland in Yunnan continue to decrease, with the reduction ratios both exceeding 10%, while the built-up impervious land has increased by 300% in land area. The conversions between cropland and grassland is mainly concentrated in the Jinsha River Basin and northern parts, while the conversion between cropland and woodland is widely distributed throughout the province, especially in the southern region. Cropland related conversions account for 74.02% of all LUCC scenarios and show significantly different transformation intensities for each period.

(3) Significant land use changes in landscape scale pose huge impacts on cropland erosion in Yunnan. During 2000−2020, the amount of cropland soil loss has decreased by 0.32×10⁸ t, with a decrease rate of 12.12%. Net soil loss change varies significantly in the six major river basins for different periods and LUCC scenarios. Except for the reclamation of cropland in the lower reaches of river basins and southern Yunnan, which bring a large amount increase in net soil loss, soil erosion in other areas significantly reduced due to the sharp reduction in cropland area. It is the first long-term quantitative study of cropland soil erosion in the area with multiple national investigation efforts, and is of great significance in under-standing the soil erosion patterns of cropland, clarifying the direction and focus of prevention, as well as protecting precious cropland resources to ensure food security in mountainous areas.