Attribution analysis and multi-scenario prediction of NDVI drivers in the Xilin Gol grassland, China

XU Mengran, ZHANG Jing, LI Zhenghai, MO Yu

College of Enoironment and Bioresources, Dalian Minzu University, Dalian 116600, China

Abstract: Grassland degradation is influenced by climate change and human activities, and has become a major obstacle for the development of arid and semi-arid areas, posing a series of environmental and socio-economic problems. An in-depth understanding of the inner relations among grassland vegetation dynamics, climate change, and human activities is therefore greatly significant for understanding the variation in regional environmental conditions and predicting future developmental trends. Based on MODIS (moderate resolution imaging spectroradiometer) NDVI (normalized difference vegetation index)data from 2000 to 2020, our objective is to investigate the spatiotemporal changes of NDVI in the Xilin Gol grassland, Inner Mongolia Autonomous Region, China. Combined with 12 natural factors and human activity factors in the same period, the dominant driving factors and their interactions were identified by using the geographic detector model, and multiple scenarios were also simulated to forecast the possible paths of future NDVI changes in this area. The results showed that: (1) in the past 21 a, vegetation cover in the Xilin Gol grassland exhibited an overall increasing trend, and the vegetation restoration (84.53%)area surpassed vegetation degradation area (7.43%); (2) precipitation, wind velocity, and livestock number were the dominant factors affecting NDVI (the explanatory power of these factors exceeded 0.4). The interaction between average annual wind velocity and average annual precipitation, and between average annual precipitation and livestock number greatly affected NDVI changes (the explanatory power of these factors exceeded 0.7). Moreover, the impact of climate change on NDVI was more significant than human activities; and (3) scenario analysis indicated that NDVI in the Xinlin Gol grassland increased under the scenarios of reduced wind velocity, increased precipitation, and ecological protection. In contrast,vegetation coverage restoration in this area was significantly reduced under the scenarios of unfavorable climate conditions and excessive human activities. This study provides a scientific basis for future vegetation restoration and management, ecological environmental construction, and sustainable natural resource utilization in this area.

Keywords: normalized difference vegetation index (NDVI); grassland degradation; geographical detector; Cellular Automat (CA)–Markov model; Xilin Gol grassland

1 Introduction

Grassland ecosystems cover 30%–40% of the Earth's total land area, contain about 20% of the global soil carbon pool (Subramanian et al., 2020), and have a large carbon sink potential (Xin et al., 2020). Thus, grassland ecosystems greatly influence the global carbon cycle and regional economy (Wang et al., 2019; Yu et al., 2019). In recent decades, grassland degradation has been influenced by climate change and human activities, and has become a major factor restricting development in arid and semi-arid areas. Grassland degradation poses a series of environmental and socio-economic problems (He et al., 2015; Gu et al., 2018). Dynamic grassland vegetation monitoring and corresponding driver analysis have always been an important topic in globalization research (Pan et al., 2018; Du et al., 2020). Assessing and analyzing the relative roles of climatic factors and human activities in the changes of grassland ecosystem are critical,which has important scientific and practical significance for the restoration and protection of grassland ecosystems, as well as the construction of regional ecological security barriers.

Normalized difference vegetation index (NDVI) is the best indicator for monitoring terrestrial vegetation status (Alfredo, 2016; Burrell et al., 2017; Ali et al., 2018; Munkhnasan et al., 2018;Rhif et al., 2020), and has been widely used in vegetation productivity estimation (Wu et al.,2020), drought (Gong et al., 2017; He et al., 2018), desertification (Li et al., 2021), and ecological environment monitoring (Xie et al., 2020). Recently, scholars have studied grassland ecosystem changes and their driving mechanisms in different regions based on long-time series NDVI data.Their views can be divided into two categories: Firstly, climate change (including precipitation and temperature) plays the dominant role. For example, Zhao et al. (2012) found that vegetation degradation in the Xilin Gol grassland during 1998–2007 was significantly correlated with decreased precipitation. Temperature and precipitation could account for 50% of the vegetation variation in the Qinghai-Tibet Plateau (Wang et al., 2016). Similar results were found on the Loess Plateau (Naeem et al., 2021) and Ili River valley (Yan et al., 2022). Secondly, grassland degradation is closely related to human activities. For instance, Li et al. (2012) found that animal husbandry, instead of climatic factors, served as the primary vegetation degradation factor. Xilin Gol grassland degradation was affected by urban expansion, road construction, and mining activity, and that the climate change impact was less significant (Sun et al., 2017; Batunacun et al.,2018). Especially, human activity contribution exceeded 90% of NDVI changes in the Sahel region's desert/steppe biome transition zone (Wu et al., 2020). Also, Pan et al. (2017) found that non-climatic drivers played a greater role in the grassland changes of the Qinghai-Tibet Plateau during 1980–2010.

However, the abovementioned studies mostly used methods such as correlation, regression, and trend analysis. In fact, a strict linear relationship may not exist between vegetation growth and its driving factors in the complex response processes, not to mention the possible presence of complex interactions. For example, Shi et al. (2019) found that NDVI increased slowly overall in this region during 2000–2015. In areas where NDVI increased significantly, human activities constituted the dominant driving factor, while in areas where NDVI decreased significantly,precipitation changes served as the primary driving factor. To clarify the complex relationship between drivers and vegetation state, Wang et al. (2017) proposed a geographical detector model(GDM) that detects spatial heterogeneity based on spatial variance analysis, and quantified the relative importance of individual driving factor, and the interaction between them (Wang et al.,2020).

The northern Chinese grassland is ecologically sensitive, which is susceptible to climate change and threatened by land desertification. Natural factors and human activities dominate vegetation changes in this area. During the past 20 a, the average precipitation in the major sand and dust source regions of northern China has increased by about 20% (IPCC, 2013). Meanwhile,a series of ecological projects successively introduced in this area after 2000, such as Grassland Desertification Control, Grain for Green, and Grazing Prohibition (Hu et al., 2013). However,under the dual influence of climate change and human activities, the main factors driving vegetation change may vary and shift, which causes further uncertainty regarding the region's future vegetation changes. Therefore, this study focused on the Inner Mongolian Xilin Gol grassland, and aimed to: (1) analyze temporal and spatial variation in NDVI during 2000–2020; (2)identify the main driving factors of vegetation change, analyze their interactions, and determine the most suitable range or type of each factor in driving vegetation growth; and (3) predict the possible change trajectory of NDVI through multi-scenario analysis.

2 Materials and methods

2.1 Study area

Xilin Gol grassland is located in the central part of the Inner Mongolia Autonomous Region,China (41°35′–46°46′N, 111°09′–119°58′E) with a total area of 2.06×105km2(Fig. 1). It has a mid-temperate arid and semi-arid continental monsoon climate. The area has an annual average temperature of 2.2℃, an average annual precipitation of about 280 mm, and an altitudinal range of 760–1925 m a.s.l. Its terrain is dominated by high plains, and its topography is high in the south and low in the north, with many low mountains and hills in the east and south, and basins scattered in between. The western and northern terrain is flat. Soils from the southeast to the northwest experience a gradual transition from chernozem to light and dark chestnut soils (Hao et al., 2014). The zonal vegetation is composed of grasslands, including typical meadow and desert grassland (Hao et al., 2014). Animal husbandry was the dominant industry, but since 2008, mining has emerged as the dominant economic sector, while animal husbandry has subsequently become the second-largest income source (Yang et al., 2011).

2.2 Data sources

The data used in this study mainly include NDVI, climate, topography, soil type, vegetation type,and human activity (Table 1).

2.3 Data analysis

The research framework is shown in Figure 2. Firstly, Theil–Sen (T–S) trend analysis and Mann–Kendall (M–K) test were performed to analyze the temporal and spatial vegetation changes.Secondly, anthropogenic and natural driving factors affecting NDVI changes were identified using GDM. Finally, multi-scenario prediction of NDVI changes was carried out using the Cellular Automat (CA)–Markov model.

2.3.1 T–S trend analysis

T–S trend analysis and M–K test can be beneficially combined to determine long-time series data trends, and have been successfully applied to the long-time series vegetation analysis (Meng et al.,

2020; Sun et al., 2021). T–S trend analysis, a non-parametric statistical method effectively avoids outliers and measurement errors without being interfered by outliers. It is therefore suitable for time series trend analysis with outliers (Wu et al., 2015). The calculation equation is as follows:

Fig. 2 Research framework in this study

wherexiandxjare the time series values in yeariand yearj, respectively;nis the time series data length (21 a in this case); median is the time series median; andSNDVIis the NDVI time series variation. WhenSNDVI>0, NDVI increases; and whenSNDVI<0, NDVI decreases. The M–K test (a nonparametric test) is used to test the variation of series and the occurrence time of sudden changes (Mann, 1945; Kendall, 1990).

Referring to the classification methods used by Yang et al. (2019) and Jiang et al. (2015), this study introduced three area types, i.e., improved area (SNDVI>0.0005), unchanged area(–0.0005≤SNDVI≤0.0005), and degraded area (SNDVI< –0.0005). The significance level (α=0.05)was tested by M–K test, and the results were divided into significant (|Z|>1.64) and non-significant changes (|Z|≤1.64).

2.3.2 GDM

GDM is a statistical approach that reveals the drivers behind factors by detecting their hierarchical spatial heterogeneity (Chen et al., 2020; Liu et al., 2020), and was developed by Wang et al. (2020). This model includes factor, interaction, risk, and ecological detectors. The first three detectors were applied in this study, including factor detector, interaction detector, and risk detector. Factor detector was used to detect dependent variable spatial heterogeneity (NDVI in this study), and to analyze the influence of each driving factor on NDVI spatial distribution through aqvalue comparison:

whereLis the number of classes of dependent variableYor factorX;NhandNare the numbers of units in classhand the entire area, respectively; σ2hand σ2are the variances ofYin classhand the entire area, respectively; SSW and SST are the intra-group sum of squares and the total sum of squares, respectively. The valueqis the degree to which the spatial heterogeneity of vegetation changeYcan be explained by factorX, and has a range of 0–1. A greaterqvalue means that the driving factor has a greater influence on vegetation change (Zheng et al., 2021). In this study, theqvalue indicated the spatial distribution consistency between NDVI change and their drivers.

Interaction detector was used to identify the interaction between every two factors (independent variables), and to evaluate the explanatory power of human activities and natural factors, jointly(enhanced or weakened) or separately, in explaining the spatial NDVI distribution. The determination method is shown in Table 2.

Table 2 Interaction relationship of each factor

Risk detector was used to identify suitable and unsuitable driving factor ranges to NDVI(dependent variable). A test was utilized to determine whether the mean NDVI value would vary across different zones for the same factor (Wang et al., 2019).

Analysis scale differences may affect GDM results to some extent. In this study, 12 driving factors were selected, namely, temperature, precipitation, wind velocity, elevation, aspect, slope,soil type, vegetation type, land use type, population density, GDP per capita, and livestock number (Table 1; Fig. 3). The highest driving factorqvalue served as the criterion for establishing the optimal spatial scale (Fu et al., 2018; Table S1). After scale comparison, a value of 7000 m was the best spatial scale for the hierarchical spatial heterogeneity analysis (Fig. S1).

2.3.3 Simulation and prediction of NDVI change

(1) CA–Markov model

CA–Markov model effectively combines the advantages of CA and Markov models. CA usually includes four basic elements, i.e., unit, state, adjacent range, and transformation rule(Wang et al., 2018). It can be expressed as follows:

whereSis the finite set of discrete cell states;fis the transfer function that defines the cell transition from timetto timet+1; andNis the neighborhood of the cell.

A Markov model is a special stochastic motion process with non-aftereffect property and stability (Wang et al., 2018). According to Wang et al. (2018), we used a state transition matrix to simulate the dynamic future NDVI change (Culik et al., 1990; Wang et al., 2018; Aguejdad, 2021).Its expression is as follows:

wherePis the state transition matrix; andPijis the transition probability fromithlevel tojthlevel in one cycle. According to the study of Wang et al. (2018),Sois the NDVI distribution of a given pixel at the initial time, and the cycle transfer probability matrix is obtained by multiplying withPn. The NDVI distribution afterncycles then becomes:

CA–Markov model is implemented by GeoSOS-FLUS software (http://www.geosimulation.cn/flus.html), which is based on the suitability probability calculation module of neural networks,and can quickly obtain the suitability probability of various land distribution types.

(2) NDVI classification

Since raster data with a discrete space and state must be used as input for CA–Markov model,we divided NDVI data into five levels according to the classification standard proposed by Wang et al. (2018): low vegetation coverage area (NDVI≤0.3), relatively low vegetation coverage area(0.30.6), with assignments of 1–5, respectively.

(3) Modeling and multi-scenario settings

Taking NDVI spatial distribution maps in 2010 and 2015 as the basic data, this study screened the driving factors based on GDM analysis results, and calculated NDVI transition matrix and suitability map using GeoSOS-FLUS software. According to transformation rule, we simulate NDVI distribution in 2020 by CA–Markov model, and used Kappa precision and Figure of Merit(FOM) (the smaller the value, the higher the simulation accuracy) analysis to compare simulated values with actual values, thus testing model reliability in terms of dynamic vegetation change prediction. Taking 2030 as the target year, the created CA–Markov model was used to design the following seven development scenarios by adjusting the range size of dominant driving factors(Table 3).

Table 3 Parameter settings in different scenarios

Fig. 3 Map of driving factors of NDVI change in the Xilin Gol grassland. (a), annual mean temperature; (b),mean annual precipitation; (c), mean wind velocity; (d), elevation; (e), aspect; (f), slope; (g), soil type; (h),vegetation type; (i), land use type; (j), population density; (k), GDP per capita; (l), livestock number.

3 Results

3.1 Spatiotemporal dynamic changes in NDVI

Inter-annual NDVI varaitaion of the Xilin Gol grassland fluctuated between 0.25 and 0.35,showing an overall increasing trend (Fig. 4). The minimum value occurred in 2001, while the maximum values occurred in 2016. As seen from inter-annual changes, NDVI gradually increased after 2000, and fluctuated greatly during 2007–2013 and 2015–2018.

Fig. 4 Inter-annual NDVI variation of the Xilin Gol grassland during 2000–2020

NDVI variation presented a spatially heterogeneous distribution. During 2000–2020, the areas with improved vegetation coverage accounted for 84.53%, while the degraded areas accounted for 7.43% (Table 4). The areas with significantly improved vegetation were distributed in the southeastern and northeastern regions of the Xilin Gol grassland, as well as some areas of central region. The areas with severely degraded vegetation were scattered in the eastern, northern, and central regions of the Xilin Gol grassland (Fig. 5).

Table 4 Statistics of NDVI trend

Fig. 5 Trend of NDVI of the Xilin Gol grassland during 2000–2020

3.2 Drivers of vegetation change

3.2.1 Effect of single driving factor on vegetation change Factor detector revealed the driving effects of natural and human activity factors on NDVI change.These single factors were ranked as follows in a descending order: wind velocity>precipitation>soil type>livestock number>temperature>population density>GDP per capita>vegetation type>slope>land use type>DEM>aspect (Fig. 6). Among them, the effects of wind velocity (q=0.6236)and precipitation (q=0.6190) on NDVI exceeded 0.6, and were therefore the dominant factors affecting vegetation status of the Xilin Gol grassland. The explanatory power of livestock number in explaining NDVI change in the Xilin Gol grassland exceeded 0.4 (q=0.4490), suggesting that grazing is also an important factor for NDVI change.

Fig. 6 Single factor affecting vegetation change in the Xilin Gol grassland during 2000–2020 and its q value.TEM, temperature; PRE, precipitation; WIND, wind velocity; ELEV, elevation; ASP, aspect; SLP, slope; SOL,soil type; VEG, vegetation type; LAND, land use type; POP, population density; GDP, GDP per captital; LIV,livestock number. The abbreviations are the same in Figure 7.

3.2.2 Effects of interactions among driving factors on vegetation change

Interaction detector evaluated the interaction effects among driving factors on NDVI change(Table S2; Fig. 7). First, the interactionqvalues among driving factors all exceeded that of single factor (Table S2), showing mutual or nonlinear enhancement. Among natural factors, average annual precipitation and wind velocity had the strongest interaction (q=0.7456). Secondly, the interaction explanatory power between average annual precipitation and temperature in explaining the spatial NDVI distribution (q=0.7192) also exceeded 0.7. Significant interactions also existed among factors such as wind velocity, temperature, precipitation, and soil use type(α=0.05, withqvalues above 0.5000 in each case). The interactionqvalues between human activity and other factors all showed an increasing trend, and the strongest interaction was observed between average annual precipitation and livestock number (q=0.7182). In this study,after two-factor interactions, nonlinear enhancement was uniformly found in the interactions between slope and other factors (Fig. 7). This result indicated that slope, as an indirect driver of vegetation growth, affects NDVI by affecting other factors.

3.2.3 Suitability analysis on the influencing factors of NDVI

Risk detector can be used to detect the suitable range or type of each driving factor in driving NDVI increase or decrease. The greater the mean NDVI value of different zones, the more favorable the range for vegetation growth. Detection results showed that precipitation (suitable range of 303.89–393.36 mm), elevation (suitable range of 1256–1695 m), slope (suitable range of 8.22°–16.18°), and population density (suitable range of 3.6–13.9 person/km2) were all positively correlated with NDVI (Fig. 8; Table 5), while temperature (suitable range of 2.74℃–3.33℃) and wind velocity (suitable range of 2.48–2.84 m/s) exhibited significant negative correlations with NDVI. The vegetation, soil, and land use types most suitable for vegetation growth were deciduous broad-leaved forest, gray forest soil, and forest land, respectively.

Fig. 7 Interactions among driving factors of NDVI

Fig. 8 Risk detection results of driving factors of mean NDVI. (a), temperature; (b), precipitation; (c), wind velocity; (d), elevation; (e), aspect; (f), slope; (g), soil type; (h), vegetation type; (i), land use type; (j), population density; (k), GDP per capita; (l), livestock number.

Table 5 Suitable range/type of each driving factor (95% confidence level)

3.3 Simulation of spatiotemporal change of NDVI

3.3.1 Simulation process and precision

Based on mean NDVI distribution map in 2010 and 2015 (Fig. 9a and b), we used the GeoSOS-FLUS v.2.3 software to generate NDVI transition matrix and transition probability(Tables S3 and S4). According to Figure 6, we selected driving factors with an explanatory power exceeding 0.3 for NDVI changes (temperature, precipitation, wind velocity, soil type, population density, GDP per capita, and livestock number) for predicting NDVI changes. NDVI transformation suitability probability map was created, and NDVI spatial distribution map in 2020 was acquired through simulation (Fig. 9d). A comparison between actual and simulated mean NDVI spatial distribution maps in 2020 (Fig. 9c and d) resulted in a Kappa coefficient of 0.7023 and an FoM coefficient of 0.1526, which suggested a high simulation precision.

Fig. 9 Real spatial distribution of NDVI in 2010 (a), 2015 (b), 2020 (c) and its simulated result in 2020 (d)

3.3.2 Simulated spatial distribution of NDVI under different scenarios Figure 10 showed simulated spatial change of NDVI under different scenarios in 2030. NDVI transition path was low vegetation coverage area–relatively low vegetation coverage area–medium vegetation coverage area–relatively high vegetation coverage area–high vegetation coverage area under the scenarios of reduced wind velocity (WIN-, Fig. 10c), increased precipitation (PRE+, Fig. 10d), and ecological protection (PD-, Fig. 10g). Compared with the business as usual (BAU, Fig. 10a), the above scenarios showed the increases in high vegetation coverage area by 18,195, 11,206, and 16,174 km2, respectively (Fig. 11). Overall, NDVI showed an increasing trend (vegetation restoration), and the increased areas located in the northeast (West Ujimqin Banner and East Ujimqin Banner) and southeast (Duolun County and Taipusi Banner) of the Xilin Gol grassland.

Under the scenarios of increased wind velocity (WIN+, Fig. 10b), decreased precipitation(PRE-, Fig. 10e), and economic priority (PD+, Fig. 10f), NDVI transition path showed an opposite trend. Compared with the business as usual (BAU, Fig. 10a), these scenarios showed decreases in high vegetation coverage area by 1000, 854, and 783 km2, respectively, and decreases in relatively high vegetation coverage area by 1895, 2630, and 1938 km2, respectively(Fig. 11). Overall, NDVI showed a decreasing trend, suggesting that vegetation coverage was significantly reduced under climate conditions that were not conducive to vegetation growth.

Fig. 10 Simulated spatial distribution of NDVI in the Xilin Gol grassland under different scenarios in 2030. (a),BAU; (b), WIN+; (c), WIN-; (d), PRE+; (e), PRE-; (f), PD+; (g), PD-. The detailed scenarios are presented in Table 3.

Fig. 11 Area change of NDVI at each level under different scenarios (WIN+, WIN-, PRE+, PRE-, PD+, and PD-)in 2030. The detailed scenarios are presented in Table 3.

4 Discussion

4.1 Dynamic changes in NDVI

T–S trend analysis found that, since 2000, vegetation in the Xilin Gol grassland has been in a state of overall recovery and partial degradation, and that vegetation changes have shown obvious spatial heterogeneity, which is consistent with the results of Batunacun et al. (2018) and Shi et al.(2019). Variation in NDVI demonstrated a distribution pattern of high in the northeast and low in the southwest, and the areas where the vegetation coverage has increased significantly locate mainly in the northeastern (East Ujumqin Banner), central (Xilinhot City and Abag Banner), and southeastern regions of the Xilin Gol grassland (Duolun County, Taipusi Banner, and Zhenglan Banner). Although the Xilin Gol grassland has been restored to a certain extent during the past 21 a, serious degradation has still occurred in some areas. These areas mainly located in West Ujimqin Banner and East Ujimqin Banner, and scattered in the Xilinhot City and Sonid Left Banner. Grassland degradation in these areas is associated with local land use changes, especially the development of mining industry (Hu et al., 2013). For example, during open-pit mine (coal mines and quarries) exploitation, wastage is stripped off in large quantities during the mining process to form numerous dumps, which occupies grasslands and destroys vegetation, resulting in a degraded ecological environment (Dai et al., 2014). Moreover, rapid development of industrial infrastructure related to mining has further aggravated grassland fragmentation, eventually leading to grassland degradation (Qian et al., 2014).

4.2 Driving mechanisms of NDVI changes

Factor detection revealed that precipitation and wind velocity were the main natural factors affecting vegetation changes, while temperature did not play a significant role. This finding is in accordance with Zhang et al. (2006), Zhao et al. (2012), Wang et al. (2016), and Naeem et al.(2021). Precipitation is positively correlated with NDVI (Zhao et al., 2012). Precipitation is the leading factor of vegetation growth and productivity improvement (Liu et al., 2019), and determines the distribution space and growth state of vegetation. The warming and humidification change of climate in northern China in recent decades (IPCC, 2013) has provided a more favorable growth environment for vegetation in the Xilin Gol grassland. Among them, Duolun County, Taibus Banner and the northeast of East Ujimqin Banner have the most significant improvement. A wind velocity increase leads to increased soil evapotranspiration. Consequently,decreased soil water content induces water stress, which is detrimental to net primary grassland productivity, and ultimately affects grassland NDVI (Gardiner et al., 2016). Excessive wind velocity blows soils away and damages plant seedlings, resulting in decreased land productivity and serious ecological problems (Wang et al., 2016; Zou et al., 2018). In contrast, topographic factors and vegetation types affect vegetation changes to lesser extents (withqvalues less than 0.2000).

Human activities also affected the spatial distribution of NDVI. A very strong positive correlation existed between grassland degradation and livestock number (Xie et al., 2012; Hao et al., 2014), and residual trend analysis showed that animal husbandry was the main driver of vegetation change in the Xilin Gol grassland during 1981–2006 (Li et al., 2012). In this regard,Sun et al. (2017) and Batunacun et al. (2018) reported similar results. This study also concluded this, that is, livestock number (q=0.4490) and population density (q=0.3220) are the main human factors affecting vegetation change. Nonetheless, risk detection indicated that NDVI was positively correlated with population density, and that local large-scale vegetation restoration(increased NDVI) might be related to ecological restoration projects and policies, intending to balance grasslands and livestock (Dong et al., 2007; Shi et al., 2019). Particularly, a series of ecological restoration projects, such as Grassland Desertification Control, Grain for Green,Ecological Immigration Policy, and Beijing-Tianjin Sandstorm Source Control (Hu et al., 2013)had promoted local vegetation restoration (Shi et al., 2019) since 2000. Additionally, grazing has been prohibited in severely degraded grasslands by executing a policy of grazing rest for moderately degraded grasslands, and by implementing rotational grazing for mildly degraded grasslands since 2002 (Dong et al., 2007; Shi et al., 2019). The implementation of the above-mentioned projects and policies effectively reduced the pressure imposed by human activities on NDVI changes. Consequently, the grassland degradation drivers gradually shifted from anthropogenic (such as grazing) to natural (climate change) factors. Correspondingly, the effects of human activities were weakened, while those of climatic factors were strengthened.

NDVI changes in this area were caused by the superposition of multiple factors. The two-factor superposition of driving factors increased their explanatory power in the spatial distribution of NDVI (Fig. 7); that is, the interactions between climatic factors, as well as between climatic factors and other factors, enhanced their explanatory power. For example, theqvalues of average annual wind velocity and precipitation were both ≤0.6500 in explaining NDVI separately (Fig. 6).However, the strongest interaction occurred between them (q=0.7456), followed in succession by the interaction between average annual precipitation and annual average temperature (q=0.7192).The combined action of multiple factors ultimately affects vegetation growth (Hao et al., 2017;Qiao et al., 2018). In this study, soil type combined with other factors also exerted important effects on vegetation growth and rainwater re-use (Piao et al., 2006; Otgonbayar et al., 2017).

4.3 Policy implications

The factors driving the changes in grassland ecosystem services are distinctively different at different time scales. Policy makers should therefore formulate suitable and time-specific measures for grassland ecosystem management (Hao et al., 2017). Although the Xilin Gol grassland is progressing towards recovery, its vegetation status is still threatened by climate change and human activities, mainly due to its fragile ecological environment. The multi-scenario simulation results provide reference and data support for designing grassland restoration policies to coordinate regional socio-economic development and ecological governance. For example,under the scenarios of reduced wind velocity, increased precipitation, and reduced livestock population, overall vegetation recovery exists. In contrast, under the scenarios of increased wind velocity, reduced precipitation, and increased livestock population, significant vegetation degradation occurs.

Grasslands are a renewable resource. Under reasonable management, grasslands can be continuously regenerated and reproduced, and be put into sustainable use by human beings. If this is not done correctly, such grasslands will degrade or even disappear. Thus, for future ecological management, governments should prioritize the determination of grassland areas with high water availability, and correlate grass availability with livestock population density. Together with this,grassland management and animal husbandry should be deployed scientifically. Depending on the specific water resource availability in the study area, grassland carrying capacity should be determined, and the avoidance of livestock overload should be taken as the primary measure to curb grassland degradation. Simultaneously, animal husbandry development should be adopted as the breakthrough point for agriculture and animal husbandry industrialization, and artificial grass planting should be vigorously implemented to promote such activities and to balance grassland and livestock biomass.

Detailed recommendations are as follows: according to the variation range of annual precipitation and wind speed, the range of grassland carrying capacity should be comprehensively evaluated. At the same time, according to the natural and socio-economic conditions of different regions, the restoration and management of severely degraded grassland should be focused on,and the ecological protection of moderately and mildly degraded grassland should be appropriately strengthened to realize the sustainable development of agriculture and animal husbandry.

4.4 Limitations of the study

This study utilized the MODIS NDVI dataset to analyze the distribution of, and changes in vegetation of the Xilin Gol grassland during 2000–2020. The spatial resolution of this dataset is relatively low, which means that it may not be an ideal choice for some areas, and may end up exaggerating the research results of these areas. For example, in areas with dense vegetation coverage, the data quality was acceptable, while in areas with sparse vegetation coverage, NDVI was greatly affected by soil. Therefore, NDVI data of Landsat or Sentinel-2 can be considered in the future. Additionally, due to the rapid fluctuation of annual plants with precipitation, using NDVI as an indicator of vegetation status may produce erroneous estimates of vegetation restoration status in areas with low vegetation coverage (Tao et al., 2020). Therefore, net primary productivity and enhanced vegetation index can be introduced to improve future research.

The selection of driving factors did not consider the influence of groundwater, soil water content, market, or policy factors. In the future, more remote sensing climatic data (such as groundwater dynamic changes, soil moisture, evapotranspiration, and drought index) and policy factors can be taken into account for an in-depth analysis of NDVI spatial distribution drivers.

The CA–Markov model could predict the changes in vegetation in the Xilin Gol grassland, but due to the interactions among data, models, driving factors, and other aspects, there were errors in prediction results. Limited by the model transformation rule, it was necessary to convert continuous NDVI data into five levels (discrete type), which led to precision loss in the data input into the model, affecting simulation results. Subsequent research can incorporate long short-term memory (LSTM) into the Markov transformation rule, and repeatedly consider the phase of driving factors to improve model precision.

5 Conclusions

This study revealed the temporal and spatial distribution characteristics of NDVI in the Xilin Gol grassland during 2000–2020, its natural and human activity drivers, and the effects of interactions between drivers on NDVI changes. The effects of dominant driving factors on the spatial distribution of NDVI under multiple scenarios were also simulated using CA–Markov model.Vegetation of the Xilin Gol grassland had an overall improvement, however, it was still being degraded locally. Temporally, NDVI in the Xilin Gol grassland generally showed a gradual increase during 2000–2020. Spatially, NDVI exhibited a distribution pattern of high in the northeast and low in the southwest. After 2000, natural factors were the dominant factors affecting vegetation coverage in the Xilin Gol grassland, but the influence of human activities could not be ignored. Coordinating the relationship between human activities and natural environment is greatly significant to vegetation restoration and ecological protection.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (31500384, 31971464), the Young Science and Technology Talents Support Program in Inner Mongolia Autonomous Region (NJYT-19-B31), and the Liaoning Province Joint Fund Project (2020-MZLH-11).

Appendix

Table S1 Zoning effect of geographical detector

Fig. S1 Scale effect of geographic detector result (q value and 90% quantile of q value). x1–x14 are the detailed driving factors.

Table S2 Interaction of driver factor

Continued

Table S3 Conversion matrix of NDVI

Table S4 Conversion probability of NDVI