Remote Sensing Quantitative Research on Soil Erosion in the Upper Reaches of the Minjiang River

1 Introduction

Soil erosion refers to the impact of water or wind that detaches and removes soil particles. This occurs gradually and causes the soil to be deteriorated, transported, and finally settled down in other places (Greenland and Szabolcs, 1994). Soil erosion is not only the most common geographical phenomenon on the earth’s surface but also the major cause of soil fertility loss and land degradation. Soil deterioration due to erosion and surface runoffs become increasingly severe worldwide. Therefore, mitigating such problems will help restore the local ecosystem. Soil is critical to the ecosystem as it can maintain the system’s stability. In turn, the stability of the ecosystem is also vital to soil conservation. For example, vegetation’s canopy, litter, and roots can reduce the erosive effects from raindrops and surface runoffs. Thereby, the soil will have better erosion resistance, more nutrients, and lower loss (Chen et al., 2012). Universal Soil Loss Equation (USLE) is commonly used to assess soil conservation in different regions (Wei et al., 2021a; Azizian and Koohi, 2021; Fan et al., 2021; Jemai et al., 2021).

China has a vast territory and large population, but its arable land is very small, threatening food security. Soil erosion may lead to soil infertility. Such loss will affect the sustainable development of both economy and society (Qi et al., 2011). Located in the transition zone between the Qinghai–Tibet Plateau and Sichuan Basin, the upper reaches of the Minjiang River are not only a vital ecological barrier in the upper reaches of the Yangtze River and Chengdu Plain but also a critical water resource for the Plain. To protect the ecological security of the headwater systems, quantitative research on soil erosion should be carried out in this area, with economic evaluation of ecological services. Affected by the 2008 Wenchuan earthquake, soil erosion happens more frequently than before in the upper reaches of the Minjiang River (Figure 1). Soil resources are essential for human life, which makes the study of soil erosion significant in this area. One method for the study is USLE, which is a popular practical measurement of soil loss. This empirical model is developed by the U.S. Department of Agriculture in the middle of the 20th century. Its initial mission is to estimate the long-term soil erosion of farmland in the U.S. in a quantitative way. The USLE is developed based on on-site measurement of soil loss and runoffs in more than 10,000 observation sites across the U.S. This model enables a quantified estimation of soil erosion through the multiplication of six factors.

FIGURE 1. Soil erosion phenomenon in the study area (Taken in 2012).

There are two methods to calculate soil erosion, and the traditional method is field measurement and analysis. However, this method requires arduous fieldwork, especially when large areas are involved. As there are many uncertainties in large-scale transformation, it has limited application. Another method is based on a variety of empirical or applied physics professional models, including USLE (Universal Soil Loss Equation), WEPP (Water Erosion Prediction Project), ANSWERS (A Real Nonpoint Source Watershed Environment Response Simulation), and LISEM (Limburg Soil Erosion Model) (Wei et al., 2021b). The WEPP model suits the quantification of sand loss in lands for agriculture and animal husbandry (Flanagan et al., 1995). ANSWERS model works well in analyzing the impact of environmental factors on soil erosion (Beasly et al., 1980). The LESEM model is based on secondary rainfall to simulate sediment yield and abortion, which is a secondary process analysis model with a relatively perfect physical mechanism (De Roo, 1996). The model variables and parameters are numerous, the mechanism is relatively complex, and the calibration is difficult. In comparison, the USLE model has a simple mechanism, easy data acquisition, and fewer limitations in the spatiotemporal scale, which is beneficial to the quantitative study of soil erosion at the regional rainfall scale (Wischmeier and Simith, 1965). The USLE model has been widely applied since the 1980s. It has attracted a great deal of attention in the academic circle (Benavidez et al., 2018; Mondal et al., 2018; Peng et al., 2018; Tung et al., 2018; Alewell et al., 2019; Lin et al., 2019; Zhang et al., 2019). The USLE model is a good choice for soil erosion studies in medium-scale watersheds. However, many studies have ignored the seasonal changes in vegetation coverage and precipitation during the calculation of the precipitation erosivity factor. In such studies, total annual rainfall and single-phase vegetation coverage factors are often adopted, but this may affect the results of the USLE model. Quick estimation shall be accomplished with the RS and GIS techniques. With advanced earth observation technology, satellite remote sensing images about terrestrial vegetation and other evaluation parameters can be obtained in real time. This study adopts the USLE model and has fully considered seasonal changes in vegetation coverage and the physical and chemical properties of the soil. With the scientific method, it begins the quantitative assessment of the soil erosion in the upper reaches of the Minjiang River. By studying soil erosion and related ecological service value in this area, this study provides a reference for soil erosion control and the ecological incentives offered to people in exchange for ecological service.

Study Area and Data

Study Area

The study area (102°35′18.026″E-104°15′14.817″E, 30°45′38.072″N-33°09′12.734″N) is the upper reaches of the Minjiang River in Aba Tibetan and Qiang Autonomous Prefecture of Sichuan Province, China (Figure 2). The Minjiang River, an important source of the Yangtze River, is a tributary with a large discharge in the Yangtze drainage area. It rises from the southern foot of the Minshan Mountain, and the part above the head of the Dujiangyan River is called the upper reaches of the Minjiang River. The upper reaches of the Minjiang River are in the transition zone from the mountainous northeastern edge of the Qinghai–Tibet Plateau to the Sichuan Basin. In terms of the farming method, this area belongs to typical semi-agricultural and semi-pastoral areas (Zhang, 2007). Based on administrative divisions, the study area includes five counties: Wenchuan County, Mao County, Li County, Heishui County, and Songpan County, with an area of 24,783.08 km². The total length of the mainstream in the upper reaches is 337 km (Wang, 1997), with a population of over 400,000.

FIGURE 2. Location of the research area.

Remote Sensing Imagery

This study obtained three phases of Landsat8 OLI images that covered the study area. The three phases had three scenes (130, 37; 130, 38; and 130, 39), with a spatial resolution of 30 m. The time phases were December 2013, June 2014, and February 2015. The cloud cover of all images was less than 3%. Remote sensing images were downloaded from the USGS website.

Meteorological Data

Eight meteorological stations inside and nearby the study area (Malkang, Xiaojin, Songpan, Dujiangyan, Hongyuan, Ya’an, Ruoergai, and Pingwu) were selected for this study. Monthly and annual average rainfall of each meteorological station was collected. The rainfall data ranged from 1980 to 2013, spanning 34 years. Meteorological data were obtained from China Meteorological Administration.

Other Data

Besides remote sensing image and meteorological radiation data, we also obtained other data, including DEM, soil types, vegetation types, land use status, and physical and chemical properties of the soil in the study area. DEM data were downloaded from the USGS website. The resolution of DEM was 30 m. Data about soil types, physical and chemical properties, vegetation types, and land use status were obtained from the Sichuan Academy of Land Science and Technology.

Methods

Soil conservation amount is an indicator of the ecosystem’s ability to conserve the soil. It can be measured by subtracting the actual amount of soil erosion from the potential amount of soil erosion in the same place. The actual amount of soil erosion refers to the amount of eroded soil. The potential amount of soil erosion refers to the amount of eroded soil when there is no vegetation coverage or land management. The actual amount of soil erosion can be directly calculated with the USLE equation. The potential amount can be calculated by canceling the vegetation coverage and management factor C and soil conservation factor P in the USLE equation (Guo et al., 2012). The calculation formula is listed as follows (Eq. 1, Eq. 2):

where A_s denotes the modulus of actual soil erosion, with the unit of t/(hm²·a); and A_q represents the modulus of potential soil erosion, with the unit of t/(hm²*a). When the soil erosion modulus is multiplied by the land area, the amount of soil erosion can be obtained. R indicates the rainfall erosivity factor, with the unit of MJ·mm/(hm²·h·a); K is the soil erodibility factor, with the unit of (t·hm²·h)/(hm²·MJ·mm); LS is the factor of slope length and gradient; C means the factor of vegetation coverage and management; P is the factor of soil conservation; and LS, C, and P are dimensionless factors.

Calculation of Rainfall Erosivity Factor (R)

Rainfall erosivity (R), which is the potential kinetic energy of rainfall to erode soil, reflects the influence of rainfall on soil erosion. As a dominant factor in the soil loss equation, this factor depends on different parameters, such as rainfall amount, precipitation duration, and rainfall kinetic energy (Sun et al., 2011; Liu et al., 2012; Zhang et al., 2014). Traditional methods cannot directly measure rainfall erosivity, so it is often roughly estimated based on rain intensity and rainfall amount. Among the estimation methods, Weschmeier’s empirical formula (McCool et al., 1982) (Eq. 3) is widely acknowledged. This formula derives the erosivity based on the long-term annual and monthly data on the average precipitation of the study area.

where R is the rainfall erosivity factor, with the unit of MJ·mm/(hm²·h·a). P_m represents years of annual average rainfall, with the unit of mm. $P_{m_i}$ denotes years of monthly average rainfall, with the unit of mm.

The spatial and temporal distribution of precipitation is uneven in the study area. Therefore, the precipitation data change greatly in different periods and seasons. For example, winter has little precipitation, but summer receives heavy precipitation. Rainfall in this area is high in the wet season. To improve the accuracy of rainfall erosivity, this study has adopted meteorological records from eight meteorological stations. These stations are in the study area and areas adjacent to it (Maerkang, Xiaojin, Songpan, Dujiangyan, Hongyuan, Ya’an, Ruoergai, and Pingwu). This study collects and processes the precipitation data over 34 years from all meteorological stations (1980–2013), including the monthly average rainfall and the annual average rainfall. The rainfall erosivity factor (R) of each meteorological station is calculated based on Eq. 3. According to the results, the rainfall erosivity factor (R) is positively correlated with precipitation. The rainfall erosivity factor (R) is 128.57 in Maerkang, 88.61 in Xiaojin, 85.02 in Songpan 253.23 in the Dujiangyan, 107.71 in the Hongyuan, 588.69 in Ya’an, 98.36 in Zoige, and 156.78 in Pingwu. The interpolation algorithm is applied to calculate the R value based on sample data in meteorological stations (Figure 3).

FIGURE 3. Distribution of rainfall erosivity factor R in the upper reaches of the Minjiang river.

Calculation of Soil Erodibility Factor (K)

Soil erodibility factor (K), which tests soil’s sensitivity to erosion, is a basic factor in the soil loss equation and an important parameter to predict soil erosion (Zhang et al., 2007). In the past, most studies on soil erodibility focused on the relationship between soil properties and soil erosion. It was not until the 1960s, when Wischmeier proposed the soil erodibility index, (Olson and Wischmeier, 1963) that scholars began to estimate soil erodibility. However, the calculation of soil erodibility factor had to wait until the 1990s when Williams et al. proposed a formula in their EPIC (Erosion-Productivity Impact Calculator) model. Therefore, it was not until the 1990s when the calculation of factor K enters became mature. The calculation of soil erodibility factor K in the study area is shown Eq. 4 (Sharpley and Williams, 1990; Guo et al., 2012).

where K is the soil erodibility factor, with the unit of (t·hm²·h)·(hm²·MJ·mm). SAN is the sand content (%), with the particle size of 0.02–2 mm. SIL shows the silt content (%), with particle size of 0.002–0.02 mm. CLA denotes the content of clay particles (%), with a particle size of less than 0.002 mm. C represents the organic carbon content (%).

According to the Soil Type Distribution Map of Sichuan Province, the study area shows strong soil zonality. It includes 11 types of soil, distributed from the North to the South, which are as follows: 1) mountain meadow soil, 2) mountain brown soil, 3) mountain ash soil, 4) mountain red–brown soil, 5) mountain gray–brown soil, 6) mountain gray-cinnamon soil, 7) cold sand yellow mud, 8) plateau black meadow soil, 9) alpine cold desert soil, 10) swamp soil, and 11) mine yellow mud.

This study takes Sichuan Soil Investigation as a reference as it offers detailed documentation of soil attributes in the study area. The study finds parameters about soil types (Table 1), sand content (SAN), silt content (SIL), clay content (CLA), and organic carbon content (C) of the study area. With formula 4, it successfully generates the K value for every soil type (Table 1) and plots the K-value distribution map in the upper reaches of the Minjiang River (Figure 4).

TABLE 1. Main soil physical and chemical parameters and K value in the upper reaches of the Minjiang river.

FIGURE 4. K value of soil erodible factor in the upper reaches of the Minjiang river.

Calculation of Slope Length and Gradient Factor (LS)

The factor of slope length and gradient (LS), which shows the slope’s influence on soil erosion, is important for calculating regional soil erosion. It is expressed as LS during the calculation. Generally, the amount of soil erosion increases with the gradient (Yang et al., 2010). Macool’s research (McCool et al., 1989) indicated that the maximum gradient for applying USLE is 15°. However, most slopes in the upper reaches of the Minjiang River are steeper than 15°, especially those in Wenchuan County and Li County. As the study area belongs to the alpine canyon landform, some slopes in this area are steeper than 75°. Therefore, when calculating the slope factor, Macool’s slope formula is used for less-than-15° slopes, and Liu’s slope formula is used for steeper-than-15°slopes (Liu et al., 1997). The formulas are as follows (Formula 5):

Using current geographic information software to extract the slope length factor is laborious. Many researchers who have contributed to this effort (Sujatha and Sridhar, 2018; Sujatha and Sridhar, 2019) mainly focus on extracting the slope length factor from the DEM. Based on physics knowledge, Moore proposed to calculate LS with the theory of unit flow power (Moore and Wilson, 1992). Wilson, Williams, and Desmet also had in-depth research on calculating slope length in the catchment (Wilson, 1986). Hickey and Van Remortel applied the digital elevation model (DEM) to calculate the slope gradient (Hickey et al., 2001; Van Remortel et al., 2004). Fu et al. applied different formulas according to the ranges of slope gradients. Moreover, they developed a program to calculate the slope length factor in Visual Studio 2010. This is an effective tool to extract the slope gradient factor (Fu et al., 2015). Tang et al. proposed a parallel grid-DEM method to calculate the slope length of the erosion. In their study, the strategy of buffer updated calculation is adopted to overcome data dependence in the parallel calculation (Liu et al., 2015). In this study, the slope length approximate method is used to calculate slope length factor (L). The slope length in the upper reaches of the Minjiang River is calculated based on a 30-m resolution DEM. Eq. 6 and Eq. 7 show the detailed calculation:

where

L indicates the slope length factor, $\lambda$ is the slope length (m), H represents elevation, and $\theta$ is the slope value $(°)$.

Based on the study area’s DEM, the slope length map and the slope map are extracted through the spatial analysis function of the ArcGIS platform. The slope length and gradient factor (LS) diagram of the study area are plotted with Eqs 5, 6, 7 (Figure 5).

FIGURE 5. LS of the upper reaches of Minjiang river.

Calculation of Vegetation Coverage and Management Factors (C)

Vegetation coverage and management factors (C) are vital indicators to evaluate the ability of vegetation cover to resist soil erosion. Under the same condition of soil type, terrain gradient, and precipitation, the vegetation coverage and management factor is the ratio of the amount of soil loss between soil land with high vegetation coverage or field management and bare land where clear plowing is carried out. This value ranges from 0 to 1 and is dimensionless. This value reflects the ability of soil under specific vegetation coverage in resisting resist erosion. The larger the value, the weaker the resistance to soil erosion under the corresponding vegetation coverage will be. That is to say, the larger the value, the worse the soil erosion will be (Zhang et al., 2001). In the USLE model, the creator provides a lookup table to confirm the soil loss ratio. To obtain the C value, the user should select the corresponding vegetation cover type and land use type and then query in the lookup table to obtain the C value (Wischmeier and Smith, 1978); in the Revised Universal Soil Loss Equation (RUSLE), the soil loss ratio is determined by a series of parameters, such as land-use type structure (PLU), vegetation cover (CC), surface roughness (SR), soil moisture content (SM), and ground surface (SC). The soil loss ratio can be obtained by multiplying these five sub-factors. Based on the work schedule and insufficient evidence of the current study area, it is difficult to obtain the abovementioned parameters and related data. Therefore, it is necessary to establish a simple algorithm suitable for a larger scale range of factors. This study takes advantage of remote sensing (RS) technology to extract the C factors effectively. In recent years, many researchers have contributed to applying remote sensing technology to extract large-scale C factors efficiently, and their algorithms are constantly improving (Wu et al., 2012). In terms of applying remote sensing technology to extract C factors, the most commonly acknowledged method is to interpret and classify the types of land cover in the study area by using remote sensing images and then assign factors to different types of land cover based on empirical values or field measurements (Folly et al., 1996; Onyando et al., 2005; Pandey et al., 2009). This method is simple and easy to implement and remains a crucial method for obtaining factors in regional soil erosion assessment. Because each type of land cover corresponds to a fixed C value, it cannot reflect the characteristics of spatial changes in surface vegetation, undermining the accuracy of factor calculation and calculation accuracy of soil erosion (Wang et al., 2002; Li et al., 2008). To accurately reflect the spatial variation characteristics of vegetation coverage, researchers adopt the field measurement data and establish a regression model between the C values and remote sensing image spectral information and then calculate the value corresponding to each pixel according to the remote sensing image of the study area. The most commonly used remote sensing spectral information is the Normalized Difference Vegetation Index (NDVI). For example, De Jong (1994) evaluated the soil erosion risk in southern France by establishing a linear regression model between NDVI and C factors. Karaburun, (2010) used a linear regression model between NDVI and C factors to evaluate the risk of soil erosion in the Buyukcekmece watershed in Istanbul. Cai et al. established the mathematical relationship between vegetation coverage and C factors according to the observation of artificial rainfall and natural rainfall in runoff areas and also established the correlation between slope sediment yield and vegetation coverage (Cai et al., 2000) (Eq. 8). This method comprehensively reflects the relationship between vegetation coverage and the C factor. Based on the review of previous studies, it is evident that vegetation coverage has a significant correlation with the amount of soil erosion. The greater the vegetation coverage, the lower the risk of soil erosion will be (Cai and Tang, 1992). This study uses the vegetation coverage method to calculate the C factor for its easily exercisable mechanism.

where f_c is the vegetation coverage of the study area in the growing season (%)

There are many methods for estimating vegetation coverage, and the most commonly used method is the pixel dichotomy model. It is a linear pixel decomposition model and is a simple and widely applied method. The NDVI has a strong linear correlation with the density of vegetation. When the vegetation coverage is at a medium and low level, the NDVI increases rapidly as the vegetation coverage increases. The NDVI keeps stable when the vegetation coverage reaches a certain degree. Therefore, the NDVI method is adopted to monitor the dynamic change of vegetation growth. The study uses the NDVI index to calculate the vegetation coverage (Deng et al., 2019) (Eq. 9). To take into consideration the changes in vegetation caused by seasonal changes, the calculation of the C factor is obtained based on the annual average vegetation coverage in the study area, and the annual average NDVI is obtained based on the three-phase NDVI Distribution map and the NDVI pixel binary model (Eq. 9). Therefore, the annual average vegetation coverage f_c in the study area can be obtained, and then the vegetation coverage and vegetation coverage in the study area can be calculated according to Eq. 8. Management factor C (Figure 6).

where NDVI_veg represents the value of pixels with full vegetation coverage. NDVI_soil denotes the value of pixels without vegetation coverage. According to previous results and principles of statistical analysis (Li, 2003), the value that corresponds to 5% of the NDVI statistical histogram is determined as NDVI_soil and the value that corresponds to 95% of the NDVI statistical histogram is determined as NDVI_veg.

FIGURE 6. C value of the upper reaches of Minjiang river.

Calculation of Soil Conservation Factor (P)

The soil conservation factor (P) refers to the ratio of soil loss after soil and water conservation measures to soil loss on the downslope. The value is 1 when there are no soil conservation measures on the soil surface. The value is 0 when there are corresponding measures and no soil loss. The soil conservation measures in this study refer to the measures to mitigate soil erosion caused by rainfall, including the changing terrain or confluence pattern to reduce the surface runoff and runoff rate (Zhang et al., 2015). Soil conservation measures for agricultural land use mainly include contour farming, strip farming, terracing, and other drainage structure construction. The soil conservation measures for grassland and pasture mainly include plowing and ridging along the contour line or its vicinity to moisture the soil surface, reduce the runoff, and resist soil erosion (Yang, 2004). The soil conservation measure factor (P) is calculated by using the comparison method, but the error in different regions is large. The soil conservation measure factor is considered to be a complicated factor in the USLE equation, and it can only be approximately estimated by indirect methods (Lu et al., 2011). Most researchers assign the p values based on in situ measurements in their research area (Yang, 2002; Yang, 2004; Xu and Shao, 2006; Yang et al., 2006). The p value in this study is based on multiple elements, including the investigation of the current situation of soil and water conservation, current research results, empirical p values under different conservation measures (Bo et al., 1997; Cai et al., 2000; Yang, 2004; Zeng and Pan, 2007) (Table 2), USLE manual, the slope percentage calculated through GIS based on DEM data, and overly analysis on the type of land use and slope (Hu et al., 2013; Yaşar Korkanç and Korkanç, 2016).

TABLE 2. P value adopted under different soil and water conservation measures in China.

According to the current land use data upstream of the Minjiang River, there are nine Class-1 land types and 30 Class-2 land types. The Class-2 and Class-1 land types are similar in land attributes, and it can be considered that the soil conservation measures are similar. Therefore, the Class-1 land types can be used as the reference basis for the land type classification of water and soil conservation measures. In the study area, the cultivated land includes paddy fields, dry land, and irrigated land. Among them, dry land accounts for 99.89% of the total cultivated land area that adopts contour belt cultivation. The p value under contour farming can refer to Table 3. The garden land use includes orchards, tea gardens, and other land uses. The garden land use in the study area is mainly planted with fruit trees such as cherries, plum trees, and pear trees. The p value directly adopts contour farming. Forest land includes forest land, shrub forest land, and other land types. Among them, forest land and shrub forest land account for 95.86% of the total forest land, most of which are primitive natural forests. Some areas with lower slopes are planted with artificial secondary forests. For areas where there are no soil conservation measures with higher slopes, the p value is 1. For areas with lower slopes, the p value is 0.8. The grassland in the area includes natural pasture and artificial pasture. Among them, natural pasture and other grasslands account for 99.98% of the total grassland area. In mining land, the water and soil conservation measures are not carried out, and the soil erosion in mining land will be aggravated, so the p value is 1. The Class-2 land use in the study area includes bare land, swampland, and agriculture facility land. There are no soil conservation measures in these areas, and the p value is 1. The Class-2 land-use types in transportation land include railway land, highway land, rural roads, and airport land. The land types are all cement hardened with no water and soil loss, and the p value is 0. The Class-2 land-use types in water surface and water conservancy facility land include the surface of river, lake water, reservoir water, pit pond water, inland beach, hydraulic construction land, glacier, and permanent snow. This type of land has no water and soil conservation measures, and the possibility of soil loss is minor, so the p value is 0. Urban villages and industrial land include Class-2 land types such as towns, villages, scenic spots, and particular land. These land types are generally featured by hardened cement surfaces where the soil loss unlikely to happen, so the value of P is 0. Table 3 shows the division of soil conservation measure factors (P) for each land use type according to different slope gradients. Figure 7 is the distribution map of p value.

TABLE 3. P values of different land use types in the upper reaches of Minjiang river.

FIGURE 7. P value of the upper reaches of Minjiang river.

Result

Analysis of the Results of Soil Conservation

The potential and the actual amount of soil erosion in the upper reaches of the Minjiang River are quantitatively assessed with the USLE formula. The actual modulus of soil erosion is directly calculated through the USLE formula (Eq. 1), whereas the modulus of potential soil erosion is calculated by canceling the vegetation coverage factor C and the land management factor P in the USLE formula (Eq. 2).

The intensity of soil erosion in the research area is graded by referring to the “Standards for Classification and Grading of Soil Erosion (SL190-2007)”. The standards have stipulated the model for grading the water erosion intensity of soil erosion and specified the grading standards for the intensity of soil erosion accordingly (Table 4). By referring to such grading standards, the actual level of soil erosion is determined based on the intensity of soil erosion, whereas the grading map for the intensity of soil erosion is obtained (Figure 8). As the modulus of potential soil erosion is subtracted from the modulus of actual soil erosion, the modulus of soil erosion is obtained accordingly (Figure 9).

TABLE 4. Grading of the level of soil erosion intensity.

FIGURE 8. Grading map of soil erosion intensity in the upper reaches of the Minjiang river.

FIGURE 9. Modulus of soil erosion in the upper reaches of the Minjiang river.

Judging from the statistical results (Table 5), it shall be noted that the total area of the zone subject to micro-level erosion in the research area amounts to 2,493.07 km², accounting for 10.06% of the total research area. The average modulus of erosion in the research area amounts to 4.57 t/(km²·a), whereas the annual quantity of soil erosion amounts to 0.0114 million tons. The total area of the zone subject to light erosion amounts to 417.31 km², accounting for 1.68% of the research area. The average modulus of erosion in the research area amounts to 1,418.30 t/(km²·a), whereas the annual quantity of soil erosion amounts to 0.5919 million tons. The area of the zone subject to intermediate erosion amounts to 571.85 km², accounting for 2.31% of the research area. The average modulus of erosion in the zone amounts to 3,398.89 t/(km²·a), whereas the annual quantity of soil erosion amounts to 1.9437 million tons. The area of the zone subject to intense erosion amounts to 798.38 km², accounting for 3.22% of the research area. The average modulus of erosion in this zone amounts to 6,067.87 t/(km²·a), whereas the annual quantity of erosion amounts to 4.8445 million tons. The total area of the zone subject to severely intense erosion amounts to 3,200.30 km², accounting for 12.91% of the research area, whereas the average modulus of erosion in this zone amounts to 10,784.13 t/(km²·a), and the annual quantity of soil erosion amounts to 34.5124 million tons. Moreover, the total area of the zone subject to utterly intense soil erosion amounts to 17,302.17 km², accounting for 69.81% of the research area, and the average modulus of erosion in this zone amounts to 34,674.29 t/(km²·a). The annual quantity of soil erosion amounts to 599.9405 million tons.

TABLE 5. Statistics of the grading of soil erosion intensity in the research area.

The primary form of soil erosion in the research area is utterly intense erosion, and the area of soil erosion accounts for 69.81% of the entire research area. The annual quantity of soil erosion in this area amounts to 599.9405 million tons, accounting for 93.47% of the total area of soil erosion in the research area, which is 641.8443 million tons.

The quantity of the potential and actual soil erosion, along with the quantity of soil conservation in the research area, can be further calculated by Eqs 10, 11, 12:

where T_s refers to the quantity of actual soil erosion (unit t), T_q refers to the potential quantity of soil erosion (unit t), T refers to the quantity of soil conservation (unit t), and S refers to the area, unit (hm²).

By multiplying the modulus of soil erosion for each unit pixel in the research area by the area of a single pixel, the actual quantity of soil erosion (T_s) in the upper reaches of the Minjiang River is estimated to reach 641.8443 million tons, whereas the quantity of potential soil erosion (T_q) is estimated to reach 925.2906 million tons, and the quantity of soil conservation (T) is calculated to reach 283.4463 million tons. The rainfall in the research area has imposed a severe impact on soil erosion, and the area of the zone subject to utterly intense erosion amounts to 17,302.17 km², accounting for 69.81% of the research area.

Due to the impact imposed by the geographical environment (i.e., alpine canyon landform), soil erosion commonly takes place in the research area. The 2008 Wenchuan earthquake had a significant impact on the geological environment in the region, and the surface soil has loosened ever since. With the arrival of the rainy season, the soil loss in the area becomes particularly severe. Especially in high altitudinal areas, ecological vegetation has yet to return to the state prior to the earthquake. On the other hand, low altitudinal areas are subject to the influence of human factors, including but not limited to building of hydropower stations, deforestation, construction of roads and railways, and mining of mineral sands.

Estimation of Economic Value of Soil Conservation

For the purpose of calculation, the economic value of soil conservation is measured from three aspects, namely, fertility preservation, soil fixation, and silt reduction. The specific value can be estimated by the shadow price method, opportunity cost method, and alternative engineering method (Zhao, 2014).

(1) Value of Fertility Conservation

The soil is capable of absorbing and preserving nutrients. Such capacity of fertility conservation is mainly achieved by adsorbing molecular nutrients on the surface of the soil cavity but without altering the material structure of the nutrients. This approach to conserving fertilizer through physical adsorption helps prevent nutrients from being leached, while enabling the nutrients to reveal a certain gradient of concentration in the soil. Besides physical adsorption, the soil is able to maintain fertility through chemical fixation, which refers to the chemical reaction of chemical ions in the soil with nutrient segregates to form compounds with utterly low solubility in the soil. Such sort of chemical reaction not only facilitates the conservation of nutrients to a certain extent but also helps reduce the poisoning level of certain harmful substances to growing plants. In the meantime, the erosion of rainfall into the soil surface will take away a large quantity of soil parent materials, in addition to a large amount of organic matter and nutrients. In general, the total nitrogen (N), total phosphorus (P), total potassium (K), and organic matter content of the soil are adopted as parameters of soil fertility when estimating the value of soil fertility preservation. According to the results of total N, total P, total K, and organic matter content of varying sorts of soils (Table 6), the soil nutrients lost in the research area can be calculated on an annual basis.

TABLE 6. Physico-chemical properties of primary types of soil in the upper reaches of Minjiang river (%).

Based on the map on the distribution of the modulus of soil erosion (Figure 9) and the map on types of soil and the physicochemical properties of each type of soil (Table 10), the total N, total P, total K, and soil conservation in the upper reaches of Minjiang River are measured accordingly. The contents of total N, total P, total K, and organic matter amount to 1.5306 million tons, 0.0021 million tons, 3.2596 million tons, and 37.9898 million tons, respectively. Moreover, using the shadow price method, the total N, total P, total K, and the content of organic matter with soil conservation are converted into chemical fertilizers and organic fertilizers of the same quality. Through the survey, the price of chemical fertilizer is estimated to be 6.6 yuan per kilogram on average, and the price of organic fertilizer is estimated to be 2 yuan per kilogram on average. On this basis, it is estimated that the economic value of nutrients retained by soil in the upper reaches of the Minjiang River amounts to 107.5934 million yuan per year.

(2) Value of Soil Fixation

To measure the value of soil fixation of the ecosystem, the total amount of soil conservation in the research area is generally converted into the land area based on bulk density and soil thickness. Such land area represents the area of land fixed by the function of soil and water conservation in the ecosystem, and the opportunity cost method is adopted (Eq. 13) to calculate the annual income to be yielded by the land area.

where V_g refers to the value of soil fixation of the ecosystem in the research area, T refers to the total amount of soil conservation,$\rho b$ refers to the soil bulk density, h refers to the thickness of the soil layer, and v refers to the annual income per unit area of land.

The amount of soil conservation is estimated to reach 283.4463 million tons. Referring to the “Standards for the Classification and Grading of Soil Erosion (SL190-2007),” the soil bulk density $\left(\rho b\right)$ is estimated to reach 1.35 g/cm³. In addition, the thickness of the soil layer is estimated based on the research findings of Ouyang et al., and the average soil thickness in China is estimated to be 50 cm. The upper reaches of Minjiang River are characterized by land cover types of woodland and grassland. Moreover, the annual income per unit area of land (v) is estimated based on the research findings on China’s forestry income in 1999, and the inflation has also been taken into account, whereas the income per unit area of land (v) is estimated to reach 50,000 yuan*km²/a.

Based on Eq. 19, the value of soil fixation of the ecosystem in the upper reaches of the Minjiang River is estimated to reach 20.996 million yuan.

Soil erosion can impose a siltation effect on rivers, lakes, reservoirs, etc., which will weaken the capacity of water storage in reservoirs and undermine their function of flood control and drought resistance. With respect to pit-pond reservoirs, curbing the soil erosion is equivalent to increasing the storage capacity of reservoirs. Due to the special topography of alpine and gorge regions in the upper reaches of the Minjiang River, the majority of the soil and silt lost from soil erosion will ultimately enter the river system. Based on previous analytical and research findings on the orientation and movement of soil loss and sediment accumulation caused by soil erosion in major river basins in China, it is shown that 24% of the soil loss derived from soil erosion will be deposited at the bottom of reservoirs and rivers, and 33% will be stagnated in river channels, whereas 37% will end up in the sea (Bo et al., 1997). The cost of dredging silt deposited in reservoirs and rivers in the research area with engineering and technical approaches is measured to assess the value loss of desilting in the ecosystem. To calculate the dredging cost, researchers refer to the engineering cost of reservoir dredging specified in the “Specification for Assessing the Service Function of Forest Ecosystem,” which amounts to 6.11 yuan/m³. Finally, the value of sediment reduction in the upper reaches of the Minjiang River is estimated to reach 307.8857 million yuan.

In a nutshell, the amount of soil conservation in the upper reaches of the Minjiang River is estimated to reach 283.4463 million tons, and the total economic value of such soil conservation is estimated to reach 436.4751 million yuan.

(3) Recommendations

Since the ecological service value of soil conservation in the research area is huge, it is necessary to enhance the quality of the ecological environment in the region. In lower altitudinal regions of the research area, human activities such as building of hydropower stations, deforestation, construction of roads and railways, and mining of mineral sands may lead to deterioration of vegetation, thus posing threats to soil erosion. Recommendations are provided as follows: we should constantly launch a series of projects of ecological restoration and reconstruction so as to contribute to the restoration of the ecological environment. Moreover, we need to impose restrictions on human activities, such as building of hydropower stations, deforestation caused by man-made factors, construction of roads and railways, and mining of mineral sands.

Conclusion

(1) Soil conservation is deemed as one of the fundamental functions of ecosystem services and one of the critical functions of supporting services. By quantitatively measuring the amount of soil erosion and conservation in the upper reaches of the Minjiang River, it is possible to quantify the value of ecological service functions of soil conservation in the region. In this study, the remote sensing technology and USLE model are adopted to quantitatively measure the amount of soil conservation in the upper reaches of the Minjiang River so as to further assess the economic value of soil conservation in the ecosystem of the research area based on the shadow price method, opportunity cost method, and alternative engineering method. The research on soil erosion in the upper reaches of the Minjiang River aims to provide evidence for the study of ecological service value and is also expected to provide a sound reference for projects of soil conservation and early warnings against soil erosion in the region.

(2) Judging from the final calculation, the quantity of actual soil erosion in the upper reaches of Minjiang River amounts to 641.8443 million tons, whereas that of potential soil erosion amounts to 925.2906 million tons, and thus the quantity of soil conservation is measured to be 283.4463 million tons. The rainfall in the research area can impose a prevalent impact on soil erosion. The area of the zone subject to utterly intense erosion is 17,302.17 km², accounting for 69.81% of the research area. Such land is primarily concentrated in the southern part of the research area, namely, Li County and Wenchuan County. The valleys in these regions are rough and feature huge altitudinal differences as well as steep slope, and thus it is necessary to adopt measures to enhance water and soil conservation.

(3) The value of ecological services is monetized. The quantity of soil conservation in the upper reaches of the Minjiang River amounts to 283.4463 million tons, and the total economic value of soil conservation amounts to 436.4751 million yuan. In particular, the economic value of annual soil conservation nutrients amounts to 107.5934 million yuan, and the value of soil fixation amounts to 209.996 million yuan, whereas the value of ecological silt reduction amounts to 307.8857 million yuan.

(4) Quick estimation can be accomplished through the RS and GIS techniques. However, considerable uncertainties are brought forward by such assessment techniques. For starters, the basic facts are not completely compatible with the actual land surface. In addition, it is still hard for certain empirical models to reveal the physical process of soil erosion.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Author Contributions

DB and SL jointly conceived and designed the experiments, developed the processing program and wrote the study, jointly performed the experiments and analyzed data, and reviewed and edited the study as supervisors. DB and SL contribute equally to the article. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Young Science and Technology Talents Development Project (Qian Jiao He KY NO. (2022)198; NO. (2020)149), and supported by the open research funded projects of Guizhou Province Key Laboratory of Ecological Protection and Restoration of Typical Plateau Wetlands (NO. (2021)06; NO. (2020)149).

Conflict of Interest

Author SL is employed by Space Star Technology Co., Ltd.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References