Spatiotemporal Variations of Grassland Ecosystem Service Value and Its Influencing Factors in Inner Mongolia, China

Abstract

The services provided by grassland ecosystems are important and irreplaceable in maintaining the balance and stability of ecosystems. The spatiotemporal variations of grassland ecosystem service value (ESV) and its influencing factors in Inner Mongolia from 2000 to 2019 were studied in this paper. Based on the socio-economic data, remote sensing data, geographic data, and meteorological data, a dynamic ESV assessment method based on the equivalent factors was used to calculate the grassland ESV for each year. The spatiotemporal dynamic variation and future trend of grassland ESV were studied by coefficient of variation index (CV), Theil–Sen median trend analysis, Mann–Kendall test, and Hurst index, and the Geodetector was used to determine the main factors affecting the distribution of ESV. The results indicated that (1) the annual average grassland ESV of Inner Mongolia was higher in the northeast than in southwest, the average ESV was 2.0794 million CNY/km², and the pixels were concentrated from 1 to 3 million CNY/km², accounting for 75.46% of the study area; (2) during the study period, the average grassland ESV increased slowly with time at an annual growth rate of 0.2, and the total ESV decreased first and then increased with the change in grassland area; (3) the average volatility was 0.16, and pixels with CV values between 0.1 and 0.2 accounted for 69.2% of the study area, indicating the fluctuation of ESV was relatively stable during the study period; (4) 37.16% of the grassland ESV in Inner Mongolia decreased slightly and 41.77% increased slightly during these years, and the two parts showed opposite trends in the future; and (5) the single factor influencing the spatial distribution of grassland ESV was mainly normalized vegetation index (NDVI) and precipitation, and the multi-factor interactions were NDVI∩slope and NDVI∩precipitation. All influencing factors exhibited a stronger impact through the two-factor interaction. This study can provide reference values for the policymaking of natural resource conservation or restoration.

1. Introduction

Ecosystem services (ES) refer to the life support products and services directly or indirectly obtained through the structure, process, and function of the ecosystem. The valuation of ES is important for environmental protection, ecological function zoning, environmental economic accounting, and ecological compensation decision-making, and also is the material basis of human survival and development [1,2,3,4]. The grassland ecosystem is one of the most widely distributed terrestrial ecosystems. In China, the area of grassland is about four million km² and accounts for 41% of the total land area [5], making grassland the largest terrestrial ecosystem type. The grasslands provide animal husbandry products and plant resources needed for the economic development, and play a key role in maintaining the pattern, function, and process of China’s natural and seminatural ecosystems, especially in arid, alpine, and other harsh habitat areas [6]. The services provided by grassland ecosystems mainly include product supply, climate regulation, gas regulation, soil and water conservation, biodiversity protection, wind prevention and sand fixation, nutrient cycling, and cultural service.

In recent years, a large area of grasslands has degraded due to increased utilization intensity and climate change. Thus, resources and ES that grasslands provide have been seriously affected. It is increasingly recognized that to improve the ecological environment and the sustainable development of human society, we must undertake the valuation of ES to build bridges between ES and market value [4,7]. Such valuation can provide decision-makers with sufficient information to guide the formulation of natural resource conservation or restoration policies, thus reducing or avoiding the misuse and abuse of natural resources. Through the valuation of ES, the important difference and spatial distribution characteristics of regional ecosystems can be clearly defined. This can guide the scientific planning of regional ecological regionalization and ecological protection and realize the rational utilization of resources and regional sustainable development.

At present, the valuation of ES can be roughly divided into two approaches [8,9], the primary-data-based approach and the unit-value-based approach. The primary-data-based approach [6,10,11,12,13,14,15] evaluates the ecosystem service values (ESV) by quantifying ES and evaluating the economic value for each service. This approach requires many input parameters and complex accounting processes [16,17,18,19], thus can only be performed on one or just a few ES at a time. The unit-value-based approach estimates ESV based on economic value per unit area of ecosystem [8,20,21,22,23,24,25,26]. The equivalent factor method is the most widely used unit-value-based approach in China [20,23]. For this method, the economic value of each service in a certain ecosystem is estimated as the product of an equivalent coefficient (dimensionless) and the economic value (CNY/ha) represented by one standard equivalent factor. Then, the total ESV is summed with the value of different ES. Compared with the primary-data-based approach, the equivalent factor method is more intuitive and easier to use, and it requires less data. Furthermore, it is particularly suitable for assessing the spatiotemporal distribution of ESV at regional and global scales [8,21,27]. A dynamic ESV assessment method was proposed by modifying and developing the equivalent factor method and can reflect the spatiotemporal variations of ES that occur in ecosystems [25].

Analyzing the influencing factors of ESV can better understand the synergistic effects of ES [28]. Most of the existing studies used regression analysis to study the changes and influencing factors of ESV in different regions. However, various factors have different effects on ESV, and many factors may act on ESV at the same time, so regression analysis cannot quantify the contribution of a single factor [29,30]. Geodetector can quantitatively analyze the driving force of a single factor and the driving force of multi-factor interactions. It can make up for the deficiency of spatial differentiation of each factor and the shortcomings of traditional analysis methods and has become an effective tool for determining various influencing factors [31].

Few studies have analyzed grassland ESV in China over a long time series. Inner Mongolia is rich in grassland resources, accounting for 20% of China’s grassland area [32]. It is an ecological protective screen in the north of China, and the evaluation and analysis of ESV is of great significance to the local economic development. Therefore, the aims of this study were to (1) adopt the dynamic ESV assessment method to evaluate the grassland ESV on the pixel scale in Inner Mongolia from 2000 to 2019 based on the socio-economic statistical data, remote sensing data, geographic data, and meteorological data; (2) reveal the dynamic spatiotemporal variation in grassland ESV in the study area during these years and explore its future trends; and (3) quantify the influence strength of individual factors on the spatial distribution of grassland ESV and evaluate the interaction effects of two different factors. The findings of this paper are intended to provide theoretical and methodological support for the assessment of grassland ESV and ecological compensation programs in China.

2. Materials and Methods

2.1. The Dynamic ESV Assessment Method

2.1.1. The Equivalent Coefficients Table and the Standard Equivalent Factor

ES are the flows of materials, energy, and information in a certain ecosystem determined by its structure and process. Based on the classification of ES by Costanza et al. [20] (17 categories), Xie et al. [24] reclassified ES into 4 primary categories and 11 secondary categories according to the understanding of Chinese people and decision-makers on ES using the method of MA (Millennium Ecosystem Assessment) [33] (

Table 1. The equivalent coefficients table for ecosystem service value (ESV) per unit area for the grassland ecosystem and four ecosystem services.

Primary Classification	Secondary Classification	Grassland Ecosystem Classification
		Prairie	Meadow
Provisioning services	Food supply	0.10	0.22
	Raw material supply	0.14	0.33
	Water supply	0.08	0.18
Regulating services	Air quality regulation	0.51	1.14
	Climate regulation	1.34	3.02
	Waste treatment	0.44	1.00
	Regulation of water flows	0.98	2.21
	Erosion prevention	0.62	1.39
	Maintenance of soil fertility	0.05	0.11
Habitat services	Habitat services	0.56	1.27
Cultural services	Cultural and amenity services	0.25	0.56

).

The equivalent coefficient is the relative weight of ESV for a certain ecosystem compared to the standard ecosystem (e.g., farmland) [23,24]. The equivalence factors for China’s ES were obtained in three ways [8]: direct comparison with ESV that has been studied in the literature, indirect comparison with ecosystem biomass for ES that lacks the estimation on its value, or based on experts’ knowledge.

In this study, we used the equivalent coefficients table for ESV per unit area constructed by Xie et al. [25] to determine the equivalent coefficients for each service of grassland ecosystem (Table 1).

The standard equivalent factor for ES is the net profit of the annual grain yield per unit area of farmland [24], which is easily traceable through well-functioning markets. The grain yield of farmland is calculated based on three major grain crops: rice, wheat, and corn. Because the prices of the three grain crops and sown area fluctuate in different years, we used the average net profit of the three grains from 2000 to 2019 as the standard equivalent factor. The formula is as follows:

$$D=\frac{1}{20}{\displaystyle \sum }_{n=2000}^{2019}{r}_n\left(S_n^r\times F_n^r+S_n^w\times F_n^w+S_n^c\times F_n^c\right)$$

(1)

where $D$ represents the average standard equivalent factor (CNY/ha), n indicates the year, and $r_n$ is the discount factor. $S_n^r$, $S_n^w,$ and $S_n^c$ denote the percentage (%) of the sown area in the total sown area of all three crops in a certain year for rice, wheat, and corn, respectively. Likewise, $F_n^r$, $F_n^w$, and $F_n^c$ denote the average net profit per unit area (CNY/ha) of all three crops in the study area of the same year for rice, wheat, and corn, respectively. Taking the price level in 2019 as the benchmark, the standard equivalent factor of each year is revised by using the CPI (Consumer Price Index).

2.1.2. Construction of the Dynamic Equivalent Factors

The internal structure and external form of a certain ecosystem in different regions and different periods in the same year are constantly changing, as are the ES from this ecosystem and their economic values. Previous studies [24,34,35] show that ES such as food production, raw material supply, air quality regulation, climate regulation, waste treatment, maintenance of soil fertility, habitat services, and cultural and amenity services are positively correlated with biomass in general. Water supply and water flow regulation are related to precipitation. Similarly, erosion prevention is closely related to precipitation, terrain slope, soil properties, and vegetation coverage. Based on these findings, Xie et al. [25] determined the spatiotemporal dynamic regulation factors as net primary productivity (NPP), precipitation, and erosion prevention (

Table 2. The regulation factors of different types of ecosystem services (ES).

Primary Classification	Secondary Classification	Regulation Factors
Provisioning services	Food production	NPP
	Raw material supply
	Water supply	Precipitation
Regulating services	Air quality regulation	NPP
	Climate regulation
	Waste treatment
	Regulation of water flows	Precipitation
	Erosion prevention	Erosion prevention
	Maintenance of soil fertility	NPP
Habitat services	Habitat services
Cultural services	Cultural and amenity services

).

The spatiotemporal dynamic equivalent factors are constructed as Equations (2)–(5).

$$F_{nij}=\{\begin{array}{l}{P}_{ij}\times F_{n1} \mathrm{or}\\ R_{ij}\times F_{n2} \mathrm{or}\\ S_{ij}\times F_{n3}\end{array}$$

(2)

where $F_{nij}$ refers to the dynamic equivalent factor per unit area for ecosystem service n of a grassland ecosystem in region i in year j; $P_{ij}$, $R_{ij}$, and $S_{ij}$ are the spatiotemporal regulation factors of NPP, precipitation, and erosion prevention in region i in year j, respectively. $F_{n1}$, $F_{n2}$, and $F_{n3}$ represent the equivalent coefficient for the services regulated by NPP, precipitation, and erosion prevention, respectively. The values of $F_{n1}$, $F_{n2}$, and $F_{n3}$ are listed in Table 1.

The equation of the spatiotemporal regulation factor of NPP ($P_{ij}$) for grassland ecosystem is as follows:

$$P_{ij}=\left(B_{ij}/\overline{B}\right)$$

(3)

where $B_{ij}$ refers to the NPP (t/hm²) of a grassland ecosystem in region i in year j, and $\overline{B}$ represents the annual average NPP (t/hm²) of a grassland ecosystem in the study area.

The equation of the spatiotemporal regulation factor of precipitation ($R_{ij}$) for a grassland ecosystem is as follows:

$$R_{ij}=\left(W_{ij}/\overline{W}\right)$$

(4)

where $W_{ij}$ refers to the average precipitation (mm/hm²) of a grassland ecosystem in region i in year j, and $\overline{W}$ refers to the annual average precipitation (mm/hm²) of grassland ecosystem in the study area.

The Universal Soil Loss Equation (USLE) is the most widely used equation to estimate soil erosion in the world [36]. It simulates the difference between potential soil erosion and actual soil erosion to derive the quantity of erosion prevention using precipitation, terrain slope, soil properties, and vegetation coverage [37]. Based on the USLE, the spatiotemporal regulation factor of erosion prevention ($S_{ij}$) for grassland ecosystems is estimated as,

$$S_{ij}=\left(E_{ij}/\overline{E}\right)$$

(5)

where $E_{ij}$ refers to the quantity of soil conservation (t/hm²) of a grassland ecosystem in region i in year j, and $\overline{E}$ represents the annual average quantity of soil conservation (t/hm²) of a grassland ecosystem in the study area.

The grassland ESV for each pixel i in the study area in year j, ${\mathrm{ESV}}_{ij}$, is calculated as,

$${\mathrm{ESV}}_{ij}={\displaystyle \sum }_{n=1}^{11}D\times F_{nij}$$

(6)

The total grassland ESV in the study area in year j, ${\mathrm{ESV}}_j$, is calculated as,

$${\mathrm{ESV}}_j={\displaystyle \sum }_{i=1}^m{\mathrm{ESV}}_{ij}$$

(7)

where m is the number of grassland pixels in the corresponding year.

2.2. Spatiotemporal Dynamic Analysis of ESV

2.2.1. Coefficient of Variation Method

Coefficient of variation (CV) is a common index to measure the dispersion degree of time series [38]. In this study, the pixel-based CV was used to analyze the fluctuation characteristics of grassland ESV. Larger CV values indicate greater volatility of grassland ESV and thus a more unstable grassland ecosystem. On the contrary, smaller CV values indicate more concentrated distribution of time series data and thus a more stable grassland ecosystem. The formula for the CV is as follows:

$$CV=\frac{\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left({\mathrm{ESV}}_i-\overline{\mathrm{ESV}}\right)}^2}{n}}}{\overline{\mathrm{ESV}}}$$

(8)

where ${\mathrm{ESV}}_i$ is the grassland ESV in the i-th year, $\overline{\mathrm{ESV}}$ is the average grassland ESV on the pixel scale in the study area from 2000 to 2019, and n is the number of statistical years.

2.2.2. Theil–Sen Median Trend Analysis and Mann–Kendall Test

Trend analysis is widely used in time dynamic analysis to explore the characteristics of interannual variation [39]. The combination of Theil–Sen median trend analysis and Mann–Kendall test could quantitatively analyze the variation trends of ESV on the pixel scale for long time series.

The Theil–Sen median trend analysis is a robust non-parametric statistical method [40,41,42], which does not require the data to follow a certain distribution and has strong resistance to individual outliers. This method estimates trend magnitudes by calculating the median slope of n(n − 1)/2 pairs of combinations. We assume the time series of ESV is {${\mathrm{ESV}}_i$} (i = 1, 2, 3, …, n), and the calculation formula is as follows:

$$\beta =Median\left(\frac{{\mathrm{ESV}}_j-{\mathrm{ESV}}_i}{j-i}\right)$$

(9)

where ${\mathrm{ESV}}_i$ and ${\mathrm{ESV}}_j$ represent the ESV in the i-th and j-th year on the pixel scale, respectively; 1 ≤ i < j ≤ 20 in this study; β is the change trend of the ESV. If β > 0, the ESV shows an upward trend and vice versa. If β = 0, the ESV shows a stable trend.

Mann–Kendall test is widely used to judge the significance of long time series [43], and this method is also applicable to non-normally distributed data series. For the same time series of ESV, the test statistic S is calculated as follows:

$$S={{\displaystyle \sum }}_{i=1}^{n-1}{{\displaystyle \sum }}_{j=i+1}^{n}sgn\left({\mathrm{ESV}}_j-{\mathrm{ESV}}_i\right)$$

(10)

where n represents the length of the time series, and $sgn$ is a symbolic function, which is defined as:

$$sgn\left({\mathrm{ESV}}_j-{\mathrm{ESV}}_i\right)=\{\begin{array}{cc}1& ({\mathrm{ESV}}_j>{\mathrm{ESV}}_i)\\ 0& \left({\mathrm{ESV}}_j={\mathrm{ESV}}_i\right)\\ -1& ({\mathrm{ESV}}_j<{\mathrm{ESV}}_i)\end{array}$$

(11)

Mann [44] and Kendall [45] proved that when n ≥ 8, S basically obeys normal distribution, its mean value is 0, and its variance is:

$$Var\left(S\right)=\frac{n\left(n-1\right)\left(2n+5\right)}{18}$$

(12)

The formula for standardizing S is:

$$Z_c=\{\begin{array}{cc}\frac{S-1}{\sqrt{Var\left(S\right)}}& (S>0)\\ 0& \left(S=0\right)\\ \frac{S+1}{\sqrt{Var\left(S\right)}}& (S<0)\end{array}$$

(13)

In the formula, $Z_c$ is the standardized test statistic, and its range is (−$\infty$, +$\infty$), which follows the standard normal distribution. We combined Theil–Sen median trend analysis with Mann–Kendall test, and when β = 0, the null hypothesis was true, which indicates that the time series data have no monotonic trend. When β ≠ 0, the null hypothesis was rejected. Under a given significance level $\alpha$, if |$Z_c$| > $Z_{1-\alpha /2}$, it indicates that a significant change occurs in the time series data. Generally, the value of $\alpha$ is 0.05, and the Z value for the 95% confidence interval is 1.96.

2.2.3. Hurst Index

The Hurst index based on the Rescaled Range analysis (R/S) is one of the main methods used to quantitatively analyze the sustainability of long time series data [46], and is widely utilized in hydrology, meteorology, and geography [47,48,49]. In this study, the pixel-based Hurst index is used to demonstrate the future trend of ESV. The calculation principle is as follows.

For the same time series of ESV {${\mathrm{ESV}}_{H\left(t\right)}$} (t = 1, 2, 3, …, n), given any positive integer $\tau$ ≥ 1, the mean value series is defined as:

$$\overline{{\mathrm{ESV}}_{H\left(\tau \right)}}=\frac{1}{\tau }{{\displaystyle \sum }}_{t=1}^{\tau }{\mathrm{ESV}}_{H\left(t\right)}\hspace{1em}(\tau =1, 2, 3, \dots ,n)$$

(14)

The cumulative deviation is formulated by

$$X\left(t,\tau \right)={{\displaystyle \sum }}_{t=1}^{\tau }\left({\mathrm{ESV}}_{H\left(t\right)}-\overline{{\mathrm{ESV}}_{H\left(\tau \right)}}\right)\hspace{1em}(1\le t\le \tau )$$

(15)

The extreme deviation sequence is calculated by

$$R\left(\tau \right)=\underset{1\le t\le \tau }{\max }\left(X\left(t,\tau \right)\right)-\underset{1\le t\le \tau }{\min }\left(X\left(t,\tau \right)\right)\hspace{1em}(\tau =1, 2, 3, \dots ,n)$$

(16)

Further, the standard deviation series is formulated by

$$S\left(\tau \right)={\left[\frac{1}{\tau }{{\displaystyle \sum }}_{t=1}^{\tau }{\left({\mathrm{ESV}}_{H\left(t\right)}-\overline{{\mathrm{ESV}}_{H\left(\tau \right)}}\right)}^2\right]}^{\frac{1}{2}}\hspace{1em}(\tau =1, 2, 3, \dots ,n)$$

(17)

Define R($\tau$)/S($\tau$) ≌ R/S, if R/S∝${\tau }^H$, which indicates that the Hurst phenomenon exists in the analyzed time series data. H is called the Hurst index, the values of which range from 0 to 1. It can be obtained by least squares in the double logarithmic coordinate system (lnt, lnR/S). If 0.5 < H < 1, the time series of ESV has sustainability, that is, the future trend of ESV is consistent with the past. The closer H is to 1, the stronger the positive sustainability. If 0 < H < 0.5, the time series of ESV has anti-sustainability, that is, the future trend of ESV is opposite to the past trend. The closer H is to 0, the stronger the reverse sustainability [50,51]. If H = 0.5, the time series of ESV is random and the future trend is unknown. According to the range of H, its sustainability is usually divided into four grades: strong sustainability (0.75 ≤ H < 1), weak sustainability (0.50 ≤ H < 0.75), weak anti-sustainability (0.25 ≤ H < 0.50), and strong anti-sustainability (0 ≤ H < 0.25).

2.2.4. Geodetector Model

Geodetector (http://www.geodetector.cn/, accessed on 1 January 2020) is a spatial analysis model for quantitatively detecting spatial differentiation and identifying the related factors [52,53,54]. The key principle of its theory is to detect the consistency between spatial distribution patterns of dependent and independent variables, which can quantitatively determine the explanatory power of individual factors and two-factor interactions [52,55]. It has no presuppositions and constraints for data, effectively overcoming the limitations of traditional statistical analysis methods when handling categorical variables, and has been applied in natural sciences, social sciences, environmental sciences, and human health [52,53,56].

The results of Geodetector include four parts: a factor detector, an interaction detector, a risk detector, and an ecological detector. This paper focuses on the driving force of ESV, so the factor detector and interaction detector were selected for quantitative elaboration and analysis. The factor detector can detect the explanatory power of individual factor X on the spatial differentiation of variable Y with the q value. The equation of the q statistic is as follows:

$$q=1-\frac{SSW}{SST}$$

(18)

$$SSW={\displaystyle \sum }_{h=1}^L{N}_h{\sigma }_h^2$$

(19)

$$SST=N{\sigma }^2$$

(20)

where SSW and SST represent the within sum of squares and the total of sum of squares, respectively; h = 1, …, L is the strata of variable Y or factor X, that is the classification or partition of Y or X; $N_h$ and $N$ are the number of units in layer h and the whole region, respectively; and ${\sigma }_h^2$ and ${\sigma }^2$ are the variances of variable Y in layer h and the whole region, respectively. The range of q is [0, 1], and the larger the value of q, the more explanatory power of factor X has for variable Y. The q value means that X explains 100 × q% of Y. The interaction detector can identify the interaction effect between two different factors (X1 and X2) and evaluate whether the explanatory power of the two-factor interaction is enhanced, weakened, or independent (

Table 3. Types of interaction mode between the two factors.

Interaction Mode	Criterion
Weaken, nonlinear	q(X1∩X2) < Min(q(X1), q(X2)) 1
Weaken, univariate	Min(q(X1), q(X2)) < q(X1∩X2) < Max(q(X1), q(X2))
Enhance, bivariate	q(X1∩X2) > Max(q(X1), q(X2))
Independent	q(X1∩X2) = q(X1) + q(X2)
Enhance, nonlinear	q(X1∩X2) > q(X1) + q(X2)

¹ X1 and X2 represent the driving factors of ESV distribution. The symbol ∩ denotes the interaction between X1 and X2.

). We evaluated whether the explanatory power of the of two-factor interaction is enhanced or weakened compared with their independent effects.

2.3. Study Area

Inner Mongolia (97°12′~126°04′ E, 37°24′~53°23′ N) is in the northern border region of China, covering an area of 1.183 million km². The study area is dominated by the Inner Mongolia Plateau, and the average altitude is 1000–1200 m (Figure 1). The terrain extends obliquely from northeast to southwest in a long and narrow shape, with the Great Khingan Mountains in the east and the Yinshan and Helan Mountains in the south. Affected by geography and topography, most of the region is dominated by temperate continental climate. Winter is cold and long, summer is short, the precipitation is low and uneven, and the rain and heat are in the same period. From northeast to southwest, the precipitation gradually decreased, while the temperature gradually increased, and gradually transited from humid and semi-humid areas to semi-arid and arid areas. Influenced by temperature and precipitation, the vegetation also showed a near-meridional spatial differentiation, with forest, grassland, and desert, in turn, from northeast to southwest. Grassland in the study area accounts for about 67% of the total area. The ecological environment of grassland is very fragile and is one of the most sensitive areas to climate change [57].

2.4. Data Source and Preprocessing

The land cover of Inner Mongolia was derived from the MODIS Land Cover Type Product (MCD12Q1, NASA (https://lpdaac.usgs.gov/, accessed on 1 January 2020), which maps global land cover at annual time steps and 500 m spatial resolution since 2001. The MCD12Q1 product provides 5 legacy classification schemes, among which the International Geosphere Biosphere Program (IGBP) [58] scheme is the most widely used. Tiles for MCD12Q1 over the study area were mosaicked and remapped using the Modis Reprojection Tool (MRT). ArcGIS (ESRI Inc., Redlands, CA, USA) was used to clip to the study area, and then the grassland region was extracted. The land cover map of Inner Mongolia in 2019 is shown in Figure 2.

NPP data were obtained from the MOD17A3HGF dataset product provided by EOS/MODIS data of NASA (https://lpdaac.usgs.gov/, accessed on 1 January 2020), with a temporal resolution of 1 year and a spatial resolution of 500 m. The scale factor of this dataset is 0.0001. MRT was used for data preprocessing at first, then the unit conversion and clipping were performed by Python (https://www.python.org/, accessed on 1 January 2020).

Meteorological data were obtained from the monthly precipitation dataset with a resolution of 1 km for China provided by the National Tibetan Plateau Data Center (https://data.tpdc.ac.cn/, accessed on 1 January 2020) [59,60].

Soil, elevation, and NDVI data were required for calculating the quantity of soil conservation in the USLE model. The soil data were derived from the Harmonized World Soil Data Base (HWSD). By connecting soil spatial distribution with the soil attribute table according to corresponding fields, the soil texture and soil organic matter content of different soils were obtained. Digital elevation model (DEM) data were derived from the Space Shuttle Radar Topographic Map Mission (SRTM) (http://srtm.csi.cgiar.org/, accessed on 1 January 2020), and the spatial resolution was resampled from 90 m to 1 km. The NDVI data were obtained from MOD13A2 (https://lpdaac.usgs.gov/, accessed on 1 January 2020), the 16-day maximum value composite (MVC) vegetation index (VI) product at the resolution of 1 km. Pixel dichotomy model was used to calculate the vegetation coverage based on NDVI.

The socio-economic statistical data were obtained from the China Statistical Yearbook and the National Agricultural Products Cost Return Assembly Yearbook from 2000 to 2019. The main variables considered were the net profit and sown area of three major crops (e.g., rice, wheat, and corn), and CPI that was used in Equation (1). The revised average standard equivalent factor of the study area from 2000 to 2019 is 2542 CNY/ha.

Five natural factors, DEM, NDVI, slope, temperature, and precipitation, and five social factors, gross domestic product (GDP), population density, distance to road, railway, and waterway were selected to detect the driving factors of grassland ESV distribution [61,62] (Figure 1 and Figure 3). The spatial distribution grid data of GDP and population density were supplied by the Research Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC, https://www.resdc.cn, accessed on 1 January 2020). The data of road, railway, and waterway were obtained from OpenStreetMap (https://www.openstreetmap.org, accessed on 1 January 2020), and the Euclidean distance was calculated for each factor using ArcGIS.

The Geodetector can only handle discrete variables, so the spatial grids with the same resolution of all variable X (influencing factors) are needed to be discrete to match the variable Y (ESV) to acquire the q value. With the help of the “GD” package in R [63], equal breaks, natural breaks, quantile breaks, geometric breaks, and standard deviation breaks were used to classify the influencing factors, and the number of classifications was set to 3 to 10. The parameter combination with the largest q value was screened out for spatial discretization. After that, the discretized values of influencing factors and ESV were imported into the “GD” package to conduct the analysis [52].

3. Results

3.1. Spatial Distribution and Temporal Change in Grassland ESV

The distribution of average grassland ESV in Inner Mongolia has strong spatial heterogeneity in terms of geographic distribution, which is higher in the northeast along the Great Khingan Mountains and lower in the southwest at the Inner Mongolia Plateau desert region (Figure 4). Overall, the average grassland ESV from 2000 to 2019 was 2.0794 million CNY/km², most of them are concentrated from 1 to 3 million CNY/km², and pixels with values in this range accounted for 75.46% of the study area (

Table 4. Area ratio of average grassland ESV in Inner Mongolia from 2000 to 2019.

ESV (Million CNY/km2)	Area Ratio (%)
<1	10.73
1–2	42.74
2–3	32.72
3–4	9.92
4–5	2.43
>5	1.46

). Most of the areas ranging from 1 to 3 million CNY/km² are distributed around Hulun Lake and the southern part of the Inner Mongolia Plateau.

The regional mean ESV showed a slowly increasing trend over time, from 2,100,758 CNY/km² in 2000 to 2,108,773 CNY/km² in 2019, with an annual growth rate of 0.2‰ (Figure 5). The total ESV is the product of the regional mean ESV and the total area of the grassland in that year. The grassland distribution was extracted from MCD12Q1 year by year. From 2001 to 2019, the grassland area in Inner Mongolia first decreased and then increased, and the lowest value occurred in 2016. Therefore, the total ESV also showed the same trend, and the turning point was CNY 1359.41 billion in 2016.

3.2. Spatiotemporal Dynamic Variation and Future Trend of Grassland ESV

3.2.1. Fluctuation Characteristics of Grassland ESV

The CV value was divided into five grades to describe the spatial fluctuation characteristics of the grassland ESV in Inner Mongolia: low volatility (CV < 0.1), middle low volatility (0.1 ≤ CV < 0.15), medium volatility (0.15 ≤ CV < 0.2), middle high volatility (0.2 ≤ CV < 0.25), and high volatility (CV ≥ 0.25) (Figure 6). In terms of spatial distribution, the interannual fluctuation characteristics of grassland ESV in the study area were spatially heterogeneous. Compared with Figure 4, the pixels of grassland ESV with high volatility are mainly distributed in the area where grassland ESV is low, and the grassland ESV with low volatility is mainly distributed in the area where grassland ESV is high. That is, pixels with low grassland ESV fluctuated greatly, while pixels with high grassland ESV fluctuated little. The results showed that the volatility of grassland ESV was negatively correlated with its value during the study period.

The volatility of grassland ESV in Inner Mongolia from 2000 to 2019 was mainly middle low and medium volatility, accounting for 34.30% and 44.90%, respectively (

Table 5. Area ratio with different volatility level of grassland ESV in Inner Mongolia from 2000 to 2019.

The Volatility of ESV	Area Ratio (%)
Low volatility (CV < 0.1)	6.72
Middle low volatility (0.1 ≤ CV < 0.15)	34.30
Medium volatility (0.15 ≤ CV < 0.2)	44.90
Middle high volatility (0.2 ≤ CV < 0.25)	12.52
High volatility (CV ≥ 0.25)	1.56

). The average volatility was 0.16. This means that the ESV of most areas in the study area fluctuated little from 2000 to 2019, and the grassland ecosystem is relatively stable.

3.2.2. Variation Trends of Grassland ESV

The grassland ESVs in Inner Mongolia are mainly increased slightly or decreased slightly: 40.59% of the grassland ESVs showed a slight decreased trend, and 45.51% showed a slight increased trend (Figure 7,

Table 6. Area ratio with different variation trends of grassland ESV in Inner Mongolia from 2000 to 2019.

Variation Trends of ESV	Area Ratio (%)
Decreased significantly (β < 0, |$Z_c$| > 1.96)	2.97
Decreased slightly (β < 0, |$Z_c$| ≤ 1.96)	40.59
Increased slightly (β > 0, |$Z_c$| ≤ 1.96)	45.51
Increased significantly (β > 0, |$Z_c$| > 1.96)	10.93

). In space, the areas with slight and severe decrease are in Xilin Gol League and Chifeng City. It is worth noting that in the southwest junction of these two cities, the annual average ESV is 2–3 million CNY/km², but the ESV is seriously decreasing. Measures should be taken to strengthen the ecological protection of these areas. The areas with slight and significant increase are distributed relatively scattered: most of them are distributed around Hulun Lake, Hinggan League, Tongliao City, etc.

3.2.3. Future Trends of Grassland ESV

The Hurst index of ESV ranges from 0.11 to 0.85, with an average value of 0.40. Pixels with Hurst index ranges (0.25, 0.50) account for 90.5% of the study area. This shows that the future trend of grassland ESV in the study area is mainly weak anti-sustainability, while the other sustainability can be basically ignored. Therefore, the sustainability situation is re-divided into three classes based on the actual result of the Hurst index: weak anti-sustainability (H < 0.35), medium anti-sustainability (0.35 ≤ H < 0.5), and weak sustainability (H > 0.5). In addition, the spatial overlay analysis of the Hurst index and change trend of grassland ESV (Figure 7) was carried out in ArcGIS to obtain the spatial distribution of future trends of grassland ESV in Inner Mongolia (Figure 8).

In the future, 24.89% of grassland ESVs will show a medium anti-sustainability trend, which are dispersedly distributed in Xilin Gol League, Ulanqab City, and Baotou City (Figure 8,

Table 7. Area ratio with different future trends of grassland ESV in Inner Mongolia.

ID	Future Trends of ESV		Area Ratio (%)		Change in the Future
1	Medium anti-sustainability	decreased significantly	0.42	24.89	From decrease to increase
2		decreased slightly	11.86
3		increased slightly	11.04		From increase to decrease
4		increased significantly	1.57
5	Weak anti-sustainability	decreased significantly	1.94	66.46	From decrease to increase
6		decreased slightly	25.30
7		increased slightly	30.73		From increase to decrease
8		increased significantly	8.49
9	Weak sustainability	reduced significantly	0.62	8.65	Continuously decrease
10		reduced slightly	3.43
11		increased slightly	3.73		Continuously increase
12		increased significantly	0.87

). The area with weak sustainability accounted for 8.65%, and the distribution was relatively sporadic. Except for medium anti-sustainability and weak sustainability areas, most of the study area showed weak anti-sustainability, accounting for 66.46%. Among them, 30.73% (ID 7) showed weak anti-sustainability and slight increase, that is, this region will show a slight decreased trend in the future; 25.30% (ID 6) showed weak anti-sustainability and slight decrease, which means that this region will show a slight increased trend in the future. More attention should be paid to the ID 7 region so that its ecological environment will not change in the opposite direction at all.

3.3. Analysis of Driving Forces of Grassland ESV

3.3.1. Individual Effect of Influencing Factors

To better understand the influences of individual factors, spatial analysis is required. Geodetector was applied to perform the spatial analysis (Figure 1, Figure 3 and Figure 4) to determine the individual effect of influencing factors. As shown in Figure 9, the q values of each influencing factor were ranked as follows: NDVI (0.39) > precipitation (0.26) > slope (0.14) > population density (0.13) > distance from waterways (0.11) > DEM (0.10) > GDP (0.06) > temperature (0.05) > distance from railways (0.02) > distance from roads (0.01). The spatial distribution of ESV may not be the result of a single factor, as the explanatory power of individual factors is usually limited. It is necessary to further explore the interaction effects of different factors.

3.3.2. Interaction Effects of Influencing Factors

The interaction detector determines whether two influencing factors acted independently or interacted, and if they interacted, whether they are enhancing or weakening. Based on the 10 influencing factors, the associated strengths of two different influencing factors on the distribution of grassland ESV in Inner Mongolia were performed using the interaction detector in Geodetector. The result (Figure 10) shows that the interaction effects of two different influencing factors has a greater impact on ESV in the study area compared with those of an individual factor. The study demonstrated that multiple driving factors comprehensively influenced ESV in a region.

Among all the interactions, the strength of NDVI∩slope reached the highest value (q = 0.434), followed by NDVI∩precipitation (q = 0.431) (Figure 10). This indicates that the combined influence of natural factors (such as NDVI, slope, and precipitation) exhibits a higher impact than social factors on grassland ESV distribution in Inner Mongolia. In particular, the interaction between NDVI and other factors has a great impact, with all the q values greater than 0.4 except the combination of NDVI∩distance from roads. In addition, the interactions between precipitation and eight other factors also dramatically enhanced the influencing strength. The q value of the combination of meteorological factors precipitation∩temperature reached 0.391, which is next to the NDVI combination.

4. Discussion

4.1. Spatiotemporal Dynamic of Grassland ESV

At present, the dynamic ESV assessment method based on equivalent factors is widely used in ESV studies to estimate the ESV of various terrestrial ecosystems in China [62]. This method has the advantages of less data requirement, simple application, easy operation, comprehensive evaluation, unified method, and easy comparison of results, and can be used as a rapid accounting tool for ecosystem service value assessment. It is more suitable for the actual situation of ESV evaluation in China [25]. However, due to data availability and technical limitations, the secondary detailed classification of the same land use type was not carried out in the commonly used land use and cover datasets, therefore it is difficult to assess ESV at the secondary ecosystem level [62]. Despite its shortcomings, the dynamic ESV assessment method is often the best or only option for resource managers and policy analysts to evaluate multiple ES under time constraints when raw data at large geographic scales, such as regional and national scales, are not available [64]. In most cases, it is valuable to dynamically evaluate ESV in different regions at pixel scale [25,65].

4.2. Influencing Factors of Spatial Distribution of ESV

Vegetation is one of the main factors influencing the spatial distribution of ESV [66]. The results of individual effects of influencing factors showed that NDVI contributed the most to the spatial differentiation of ESV and was in a dominant position. Previous studies also found the same conclusion in different regions. For example, NDVI plays a key role in the spatial distribution of ESV in the Yellow River Basin [67], Yinchuan City [61], and the “Silk Road Economic Belt” [62]. The results of the interaction effect of influencing factors showed that the combination of NDVI∩slope and NDVI∩precipitation exhibited a higher impact than other combinations on ESV distribution. Similar results have been found in karst areas [68] and the “Silk Road Economic Belt” [62]. Most of the studies have found that interaction between each two factors has a greater impact on ESV compared to a single factor [61,62,67,68]. Thus, the ESV spatial distribution is the result of multiple factor interactions.

Precipitation in the study area decreased from east to west (Figure 3c). The slope of the study area is gentle, and most areas have slopes of less than 10° and are suitable for grass growth (Figure 3a). The spatial distribution of NDVI is consistent with that of grassland ESV (Figure 3b and Figure 4). In other words, the area with larger NDVI also has higher ESV. On the contrary, the area with smaller NDVI also has lower ESV. Therefore, better growing grassland can provide greater ESV for human beings. In order to maintain the ESV of a grassland ecosystem, we should continue to strengthen the protection and management of grasslands, moderate grazing, and actively carry out comprehensive control of soil erosion and grassland desertification in the central part of Inner Mongolia. At the same time, we should improve the quality of returning farmland to grassland. Thus, the grassland ecosystem can provide abundant ES for human beings continuously.

4.3. Uncertainties and Future Directions

The results of the total ESV of grassland in Inner Mongolia is approximate with the results of Xie et al. [8], which total the ESV of Inner Mongolia at USD 271.7 billion in 2010, and our result is CNY 1383.72 billion only for grassland. There are still some problems. Firstly, the spatial resolution and accuracy of MCD12Q1 is not very high compared with other land use and cover datasets [69]. However, the temporal resolution and span of this product met the requirements of this study on land use and cover datasets. In addition, the land cover obtained by MCD12Q1 did not subdivide grassland types, while Xie et al. [8,25] has classified grassland types into prairie, shrub grass, and meadow. There are several grassland types in the study area [70], but the shrub grass rarely appears. Therefore, the equivalent coefficients of prairie and meadow were averaged as the equivalent coefficients of the whole grassland in the study area [61]. Simply averaging the equivalent coefficients of different grassland types cannot reflect the role of different grassland types in ES. Therefore, in subsequent studies, grassland types should be subdivided, and then use the corresponding equivalent coefficients to calculate the ESV of the corresponding region.

5. Conclusions

This study calculated the grassland ESV in Inner Mongolia from 2000 to 2019 and analyzed the spatiotemporal variations and its main influencing factors. The main findings of this work are as follows: (1) The annual average grassland ESV of Inner Mongolia has strong spatial heterogeneity which was higher in the northeast and lower in the southwest; 75.46% of the grassland ESV in the study area was concentrated from 1 to 3 million CNY/km². (2) During the study period, the average grassland ESV increased slowly with time at an annual growth rate of 0.2‰, and the total value decreased first and then increased with the change in grassland area. (3) The fluctuation of grassland ESV was relatively stable during the study period and the volatility was negatively correlated with its value; 37.16% of the grassland ESV in Inner Mongolia decreased slightly and 41.77% increased slightly during these years, and the two parts showed opposite trends in the future. (4) The influencing factor detection shows that the single influencing factor is mainly NDVI and precipitation, and the multiple influencing factor interaction shows NDVI∩slope and NDVI∩precipitation contributing most to the spatial distribution of grassland ESV in Inner Mongolia. The spatial distribution of ESV is influenced by the combination of natural and social factors. These results can provide a scientific basis for policymaking in grassland ecosystems based on the values of ES and site-specific decisions for ecosystem management in different regions. In the future, applying the detailed grassland types, adding more regulation factors to evaluate the spatial distribution of ESV, and analyzing more influencing factors of grassland ESV are some directions for research.