Migration of Rural Residents to Urban Areas Drives Grassland Vegetation Increase in China’s Loess Plateau

Human activities are critical factors influencing ecosystem sustainability. However, knowledge on regarding the mechanisms underlying the response of vegetation dynamics to human activities remains limited. To detect the driving factors and their individual contribution to the grassland vegetation dynamics in China’s Loess Plateau, a structural equation model (SEM) and a principal component regression model were built. The SEM showed that population change and urbanization, temperature and humidity, and agriculture and economy accounted for 62.5%, 31.2%, and 7.7%, respectively, of the overall impact directly affecting grassland vegetation dynamics. Furthermore, the principal component regression model demonstrated that the effects of the urbanization rate on the grassland above-ground biomass exceeded those of the other factors. The agriculture population had the maximum negative effect on grassland area. The higher the urbanization rate means the higher the number of residents migrates from rural to urban areas. Following this argument, the disturbances of human activities to grassland vegetation were expected to gradually decrease in rural areas, where the vast majority of the Loess Plateau is located. The migration of rural residents to urban areas promoted the increase in biomass and areas of grassland vegetation. Our findings suggest that the effect of urbanization should be considered when assessing vegetation change.


Introduction
Human activities and climate change are two of the most important factors driving the dynamics of terrestrial ecosystems [1][2][3]. Over the past decades, the spatiotemporal patterns of vegetation growth dynamics have been altered due to the ongoing climate change and intensifying human disturbances [4][5][6][7]. Especially in arid and semi-arid regions, climate change and human activities may have opposite effects on the dynamics of vegetation. On the one hand, overgrazing, deforestation, Thus, in the present study, over 20-year data on 22 counties were analyzed to evaluate the responses of grassland vegetation dynamics to human activities (i.e., population change, urbanization, agricultural production, and economic development) and climate change (i.e., precipitation, temperature, and sunshine) in the Loess Plateau (Figure 1b). We focused on these factors because of the following reasons: (1) temperature and water resources are regarded as the controlling factors driving the seasonal variation of vegetation growth in arid and semi-arid areas [43,44]; (2) urbanization, population density, agricultural production, and economic development are the key human drivers of ecosystem structure and function [45][46][47][48]. Similar to those in other developing countries, rural areas account for the vast majority of China's Loess Plateau. The migration of rural populations into urban areas may decrease human disturbance (such as agriculture cultivation, and grazing) in rural areas and thus increase the grassland vegetation cover. Vegetation growth is more controlled by water resources than by other climate factors in arid and semi-arid areas [49,50]. In the Loess Plateau, the annual precipitation presents a non-significant upward trend [37]. However, vegetation greenness has increased significantly in this region [38]. Thus, climate factors may not be the dominant drivers of grassland vegetation in the region.
The objectives of the current research were to (1) assess the extent to which human activities and climate change factors impact the increase of grassland vegetation and (2) identify the specific human activities that drive the increase of grassland vegetation in the Loess Plateau. Answering these questions will greatly improve our mechanistic understanding of the potential factors affecting long-term grassland vegetation dynamics in the semi-arid regions of the world.

Study Region
The Loess Plateau covers an area of approximately 640,000 km 2 , accounting for 6.7% of the total territory of China. The elevation, ranging from 200-3000 m, is high in the west and low in the east [51,52]. The region has a typical semi-arid continental monsoon climate. The mean annual temperature ranges from 3.6 °C in the northwest to 14.3 °C in the southeast. The mean annual precipitation varies from 200 mm in the northwest to 700 mm in the southeast, whereas rainfall mostly occurs as high intensity rainstorms during summer [53]. The region covers 287 counties from 7 provinces and it spans 5 bioclimatic zones from northwest to southeast, namely, arid deserts, arid and semi-arid desert-grasslands, semi-arid typical grasslands, sub-humid to semi-arid forestgrasslands, and sub-humid forest [12,54]. The sub-humid forest zone (FOR Zone) was excluded from the analysis because we focused only on grasslands (Figure 1b). Approximately 8.5% of China's population live in the Loess Plateau, where the population density reaches 168 person/km 2 [41,55]. Thus, in the present study, over 20-year data on 22 counties were analyzed to evaluate the responses of grassland vegetation dynamics to human activities (i.e., population change, urbanization, agricultural production, and economic development) and climate change (i.e., precipitation, temperature, and sunshine) in the Loess Plateau (Figure 1b). We focused on these factors because of the following reasons: (1) temperature and water resources are regarded as the controlling factors driving the seasonal variation of vegetation growth in arid and semi-arid areas [43,44]; (2) urbanization, population density, agricultural production, and economic development are the key human drivers of ecosystem structure and function [45][46][47][48]. Similar to those in other developing countries, rural areas account for the vast majority of China's Loess Plateau. The migration of rural populations into urban areas may decrease human disturbance (such as agriculture cultivation, and grazing) in rural areas and thus increase the grassland vegetation cover. Vegetation growth is more controlled by water resources than by other climate factors in arid and semi-arid areas [49,50]. In the Loess Plateau, the annual precipitation presents a non-significant upward trend [37]. However, vegetation greenness has increased significantly in this region [38]. Thus, climate factors may not be the dominant drivers of grassland vegetation in the region.
The objectives of the current research were to (1) assess the extent to which human activities and climate change factors impact the increase of grassland vegetation and (2) identify the specific human activities that drive the increase of grassland vegetation in the Loess Plateau. Answering these questions will greatly improve our mechanistic understanding of the potential factors affecting long-term grassland vegetation dynamics in the semi-arid regions of the world.

Study Region
The Loess Plateau covers an area of approximately 640,000 km 2 , accounting for 6.7% of the total territory of China. The elevation, ranging from 200-3000 m, is high in the west and low in the east [51,52]. The region has a typical semi-arid continental monsoon climate. The mean annual temperature ranges from 3.6 • C in the northwest to 14.3 • C in the southeast. The mean annual precipitation varies from 200 mm in the northwest to 700 mm in the southeast, whereas rainfall mostly occurs as high intensity rainstorms during summer [53]. The region covers 287 counties from 7 provinces and it spans 5 bioclimatic zones from northwest to southeast, namely, arid deserts, arid and semi-arid desert-grasslands, semi-arid typical grasslands, sub-humid to semi-arid forest-grasslands, and sub-humid forest [12,54]. The sub-humid forest zone (FOR Zone) was excluded from the analysis because we focused only on grasslands (Figure 1b). Approximately 8.5% of China's population live in the Loess Plateau, where the population density reaches 168 person/km 2 [41,55]. Several counties in this region are considered as underdeveloped than those in other regions in China. Nevertheless, the economy of this region started to recover after the implementation of the Reform and Opening up policy in 1978. Agricultural income was no longer the major or sole source of family livelihood. Specifically, an increasing number of rural residents have started to engage in the secondary or tertiary industry in the last two decades [34].

Data Sources and Preprocessing
To assess the impacts of human activities and climate change, we selected 22 study counties with available long-term climate data (Figure 1b). They were distributed exactly across the four bioclimatic zones with decreasing amounts of annual precipitation ( Figure S3), and different economic development levels and population densities. Remotely sensed vegetation indices, grassland cover, climatic and topographic data, and proxies of human activities (e.g., agriculture production, economic development, population, and urbanization) were collected. The data sources and quality control are described, as follows.

Remotely Sensed Vegetation Indices
Considering the severe fragmentation of the land surface in the Loess Plateau [56] data covering all of the 22 study counties were employed to build a seamless mosaic image. The Landsat multispectral data were obtained from the United States Geological Survey (http://glovis.usgs.gov/). Then, we selected cloud-free images during the growing season when the growth of vegetation reaches its peak (July to August), and calculated the NDVI.

Climate and Topography Data
Most areas in the Loess Plateau have an arid and semi-arid climate and thus, precipitation and temperature are the controlling factors of vegetation growth [38]. Precipitation, temperature, and sunshine data at the weather stations located in 22 study counties were downloaded from the China Meteorological Data Network (http://data.cma.cn/). Furthermore, the Standardized Precipitation Evapotranspiration Index (SPEI) with a 12-month time scale was calculated by using the monthly precipitation and temperature data [57].
Elevation was regarded as a key topographical factor influencing the growth of vegetation [58]. The elevation data were acquired from the Shuttle Radar Topography Mission with a spatial resolution of 30 m.

Indicators of Human Activities
The Defense Meteorological Satellite Program (DMSP)/Operational Linescan System (OLS) stable nighttime light (NTL) data have been widely used as proxies of urban extent, socioeconomic development, and population distribution [59][60][61][62][63]. However, due to the lack of on-board calibration, varied atmospheric conditions, satellite shifting, or sensor degradation the NTL datasets need to be calibrated before their application to relevant studies [64,65]. Li and Zhou developed a stepwise calibration for NTL and provided calibrated and temporally consistent NTL time series datasets from 1992 to 2013 [66]. We used their calibrated NTL to represent an integrated indicator of human activities in the 22 study counties.
The net income per capita of rural residents, number of livestock on hand (goats and sheep are the main livestock in the Loess Plateau and thus, the number of goats and sheep on hand is always collected as a livestock indicator), proportion of the added value of the primary industry to the gross domestic product (GDP), total power of agricultural machinery per capita, agricultural population, urban population, and urbanization rate (urban population divided by total population) were used as proxies of human activities. These data were obtained from China's Economic and Social Big Data Research Platform (http://data.cnki.net/) and the statistic departments of local governments in the study counties. As the territory areas are different among the study counties, the number of livestock on hand, agricultural population, urban population, and NTL were divided by the areas of the study counties to derive the density of livestock on hand, density of agricultural population, density of urban population, and density of NTL, respectively.

Extracting Grassland Pixel from Remote Sensing Images
Landsat data are relatively high-resolution among the open access images. The remote-sensing images (cloud-free, clear and easy to identify) were selected to perform the supervised classification for each study county (Table S1). The interval when supervised classification was performed was no more than 5 years. The region of interest, that is, the training samples, was selected from each study county and then compared with images from Google Earth in the same region to improve classification accuracy. Five land use types were classified: (1) forest, (2) grassland, (3) cropland, (4) unused land, and (5) other land use (urban/residential/water). The maximum likelihood classification tool in the Environment for Visualizing Images software was applied to perform supervised classification with a total classification accuracy of ≥ 80% and a Kappa value ≥ 0.8. In this study, accuracy was evaluated on the basis of the collected training samples. Then, the grassland pixels of all images were extracted on the basis of those images with supervised classification at each study county. Lastly, we calculated the above-ground biomass (AGB, see Section 2.3.2 for the equations) and the area of grassland vegetation for each study county.

Above-Ground Grassland Biomass Model
To estimate the AGB of grassland vegetation, we surveyed 21 field sampling sites in the Loess Plateau ( Figure 1b). We selected a 60 m × 60 m plot without human disturbance and measured the AGB in six 1 m × 1 m quadrats from 2013 to 2016. The plants of each quadrat were clipped at the ground level at the peak of the growing season (July or August). The plants that were clipped from the six quadrats of each plot were pooled and oven-dried to a constant weight (65 • C for 72 h). Moreover, the remotely sensed NDVI and elevation at corresponding plots were obtained. The linear relationship between the NDVI, elevation, and AGB of the grassland are shown in Figure S1. Lastly, a satellite-based regression model for evaluating the grassland AGB of the Loess Plateau was built by the following equation: where AGB (x,t) is the grassland AGB (g/m 2 ) of the plot, x is the spatial location, and t is time. ASL (km) is the elevation of the plot. The Statistical Package for Social Sciences software (SPSS 20.0) was used to estimate the parameter values of model (1). The final model is expressed as follows: where R 2 is the coefficient of determination the grassland AGB of the plot, RMSE is the root mean square error, and n is the sample size.
Model (2) was used to calculate the grassland AGB of each study county in ESRI ArcGIS 10. The mean (AGB m , g/m 2 ) and total (AGB T , 10 4 t) grassland AGB were calculated for each county. We also calculated the density of AGB T (AGB T density, t/km 2 i.e., AGB T divided by the area of the corresponding county) to make the AGB T values of the study counties comparable. In this study, AGB T density was different from AGB m . Specifically, the AGB T density was derived from all areas of the counties (including all land use types), whereas AGB m was the average of the grassland AGB derived solely from grassland pixels while also excluding other land use types.

Conceptual Structural Equation Model (SEM)
A conceptual SEM was developed to detect specific driving factors and assess their individual effect on the dynamics of grassland vegetation in the Loess Plateau ( Figure S2). Five blocks (each block represents a measurement model) are marked by dashed boxes in Figure S2. Herein, the block in the SEM was M (manifest)-to-L (latent) block type, which is referred to as formative, and causation was presumed to flow from the manifest variables to the latent variable [67]. Thus, in each measurement model in Figure S2, the arrows flow from the manifest variables (rectangle) to the latent variable (ellipse).
Several distinct approaches are available to perform SEM. The partial least squares (PLS) approach focuses on the analysis of variance. It is a soft modeling approach to SEM. Moreover, the PLS approach is recommended as an operative analytical tool to reduce errors [68]. We used Smart PLS 3.2.6 software to perform the measurement model and the structural model determination ( Figure S2). The software is a PLS-based SEM technique for studying causal models that involve multiple constructs with multiple indicators [69].
From the spatial aspect, we used the data of all study counties to perform the SEM firstly. And then, we performed the SEM while using the data of study counties at each bioclimatic zone. From the temporal aspect, the data of all study counties were divided into different time periods, namely, 1992-1997, 1998-2002, 2003-2007, and 2008-2013. The SEMs of these time periods were also developed.

Principal Component Regression Model
To complement the findings of the SEM, we also investigated the driving factors for the dynamics of grassland vegetation in the Loess Plateau by using principal component analysis (PCA). A total of 14 original variables were selected to reflect human activities and climate change. SPSS 20.0 was employed to perform the PCA and extract the principal components (PCs) from the original variables. As a result, the AGB T density, AGB m , and proportions of grassland area to total area (GA/TA) were taken as response variables. Moreover, all PCs were treated as explanatory variables. The PCs with the most significant influence on response variables can be identified by implementing the principal component stepwise regression model.

Trends in Annual Climate
The average precipitation from 1986 to 2013 was approximately 177 mm in the desert zone, 300 mm in the desert-rassland zone, 406 mm in the typical grassland zone, and 485 mm in the forest-grassland zone. The average temperature in the same period was around 9 • C in the desert and the desert-grassland zones, 7 • C in the typical grassland zone, and 10 • C in the forest-grassland zone. Temperature exhibited an increasing trend in the four bioclimatic zones from 1986 to 2013. However, the increasing trend was not significant in the desert-grasslands ( Figure S3). Moreover, no significant trend in precipitation was found ( Figure S4).
The Mann-Kendall trend test for each study county also showed similar patterns for temperature and precipitation (Figure 2a,b). Specifically, the mean annual precipitation did not show a growing trend in all 22 study counties (Figure 2a), whereas the mean annual temperature showed an upward trend in over 60% of the study counties (Figure 2b).

Trends in Grassland AGB and Area
The AGBT demonstrated an increasing trend in over 50% study counties (Figure 2c). AGBm exhibited a growing trend in the desert, typical grassland, and forest-grassland zones (Figure 2d). No decreasing trend in AGBT and AGBm was observed in any county (Figure 2c,d). Moreover, the grassland area exhibited an increasing trend in all of the study counties in the desert-grassland and typical grassland zones and in over 70% of the study counties in the desert and forest-grassland zones ( Figure 2e). Therefore, a significantly increasing trend was also noted in grassland AGB and the area in the four bioclimatic zones of the Loess Plateau (Figure 2c-e).

Correlation Between the Dynamics of Grassland Vegetation and Human Activities
The seven proxy variables of human activities had a clear correlation with AGBT density, AGBm, and the proportions of grassland area to the total area (Table 1). Among these proxy variables, the urbanization rate and the density of urban population, both of which reflect the urbanization development level, had the highest positive correlation with AGBT density and AGBm, with coefficient values of 0.49 and 0.57, respectively.

Trends in Grassland AGB and Area
The AGB T demonstrated an increasing trend in over 50% study counties (Figure 2c). AGB m exhibited a growing trend in the desert, typical grassland, and forest-grassland zones (Figure 2d). No decreasing trend in AGB T and AGB m was observed in any county (Figure 2c,d). Moreover, the grassland area exhibited an increasing trend in all of the study counties in the desert-grassland and typical grassland zones and in over 70% of the study counties in the desert and forest-grassland zones (Figure 2e). Therefore, a significantly increasing trend was also noted in grassland AGB and the area in the four bioclimatic zones of the Loess Plateau (Figure 2c-e).

Correlation Between the Dynamics of Grassland Vegetation and Human Activities
The seven proxy variables of human activities had a clear correlation with AGB T density, AGB m , and the proportions of grassland area to the total area (Table 1). Among these proxy variables, the urbanization rate and the density of urban population, both of which reflect the urbanization development level, had the highest positive correlation with AGB T density and AGB m , with coefficient values of 0.49 and 0.57, respectively. Note: r is the Pearson correlation coefficient; P is the p-value and the sample size n is 401; significance levels are as follows: ** P < 0.01, * P < 0.05. Human activities variables include: net income per capita of rural residents (NI per capita); the density of livestock on hand (LOH density); the proportion of the added value of the primary industry to GDP (AGDP/GDP); total power of agricultural machinery per capita (TPOAM per capita); the density of agricultural population (AP density); the density of urban population (UP density). Grassland vegetation dynamics variables are: AGB T density, AGB m and the proportions of the grassland area to the total area (GA/TA).

Spatial Effects of Human and Climate Factors on Grassland Vegetation Dynamics
Over the entire study region (the union of the four bioclimatic zones), the increase in population and urbanization directly promoted the increase of grassland vegetation (Figure 3a). Population and urbanization showed the largest direct effect (0.625), followed by temperature and humidity (0.312) (Figure 3a). Agriculture production and economic development has no significant direct effect, but indirectly impacts the dynamics of grassland vegetation by interacting with population and urbanization in a significantly negative way. Similarly, the population increases and urbanization development could reduce the AGB and area of grassland vegetation indirectly when the interaction with agriculture and economic was considered (Figure 3a). Overall, the indirect effect of population and urbanization was far less marked than the direct effects.
The effects of human activities and climate change on the dynamics of grassland vegetation were unequal across bioclimatic zones. First, increases in population and urbanization directly increased grassland vegetation AGB and area in the desert and typical grassland zones, and the direct effect was close to that in the entire study region (Figure 3b,d). On the contrary, increases in population and urbanization had significant and negative direct effects in the forest-grassland zone ( Figure 3e) and no direct effect in the desert-grassland zone (Figure 3c). Second, temperature and humidity exerted a nonsignificant effect on the grassland vegetation dynamics in all zones (Figure 3b-e). Finally, agriculture production and economic development had a significant and positive direct effect on grassland vegetation dynamics only in the desert-grasslands zone (Figure 3c), and no significant and direct effects in other zones (Figure 3b,d,e). Moreover, the indirect effects (interacting with population and urbanization) of agriculture and economy on grassland vegetation dynamics was only significant in typical grassland zone (Figure 3d). In summary, population and urbanization remained as the dominant factor directly driving the dynamics of grassland vegetation in all the bioclimatic zones, except in the desert-grassland zone (Figure 3b-e).

Spatial Effects of Human and Climate Factors on Grassland Vegetation Dynamics
Over the entire study region (the union of the four bioclimatic zones), the increase in population and urbanization directly promoted the increase of grassland vegetation (Figure 3a). Population and urbanization showed the largest direct effect (0.625), followed by temperature and humidity (0.312) (Figure 3a). Agriculture production and economic development has no significant direct effect, but indirectly impacts the dynamics of grassland vegetation by interacting with population and urbanization in a significantly negative way. Similarly, the population increases and urbanization development could reduce the AGB and area of grassland vegetation indirectly when the interaction with agriculture and economic was considered (Figure 3a). Overall, the indirect effect of population and urbanization was far less marked than the direct effects.
The effects of human activities and climate change on the dynamics of grassland vegetation were unequal across bioclimatic zones. First, increases in population and urbanization directly increased grassland vegetation AGB and area in the desert and typical grassland zones, and the direct effect was close to that in the entire study region (Figures 3b and 6d). On the contrary, increases in population and urbanization had significant and negative direct effects in the forestgrassland zone (Figure 3e) and no direct effect in the desert-grassland zone (Figure 3c). Second, temperature and humidity exerted a nonsignificant effect on the grassland vegetation dynamics in all zones (Figures 3b-e). Finally, agriculture production and economic development had a significant and positive direct effect on grassland vegetation dynamics only in the desert-grasslands zone (Figure 3c), and no significant and direct effects in other zones (Figure 3b,d,e). Moreover, the indirect effects (interacting with population and urbanization) of agriculture and economy on grassland vegetation dynamics was only significant in typical grassland zone (Figure 3d). In summary, population and urbanization remained as the dominant factor directly driving the dynamics of grassland vegetation in all the bioclimatic zones, except in the desert-grassland zone (Figure 3b-e).  The path coefficients and standardized total effects for deserts zone, desert-grasslands zone, typical grasslands zone, and forest-grasslands zone, respectively.

Temporal Changes in the Effects of Human and Climate Factors on Grassland Vegetation Dynamics
The direct effect size of population and urbanization was larger than those of the other driving factors (temperature and humidity, agriculture, and economy) in all time periods (Figure 4). The effect of population and urbanization decreased annually. On the contrary, the effect of temperature and humidity increased annually (Figure 4).  The path coefficients and standardized total effects for deserts zone, desert-grasslands zone, typical grasslands zone, and forest-grasslands zone, respectively.

Temporal Changes in the Effects of Human and Climate Factors on Grassland Vegetation Dynamics
The direct effect size of population and urbanization was larger than those of the other driving factors (temperature and humidity, agriculture, and economy) in all time periods (Figure 4). The effect of population and urbanization decreased annually. On the contrary, the effect of temperature and humidity increased annually (Figure 4).  The path coefficients and standardized total effects for deserts zone, desert-grasslands zone, typical grasslands zone, and forest-grasslands zone, respectively.

Temporal Changes in the Effects of Human and Climate Factors on Grassland Vegetation Dynamics
The direct effect size of population and urbanization was larger than those of the other driving factors (temperature and humidity, agriculture, and economy) in all time periods (Figure 4). The effect of population and urbanization decreased annually. On the contrary, the effect of temperature and humidity increased annually (Figure 4).

Effects Based on PCA
The first six PCs explained 87.2% of the total variance, and thus, they were extracted to reduce the dimensionality of the 14 original variables ( Table 2). The first PC was heavily loaded on sunshine percentage, sunshine hour, days with daily precipitation >0.1 mm, mean annual precipitation, and mean relative humidity. For ease of interpretability, we termed these aspects of climate change as "precipitation and sunshine" (PC1). The second PC was loaded heavily on urbanization rate, the density of urban population, and the density of livestock on hand. Thus, we termed these aspects of human activities as "urban population and urbanization" (PC2). The third PC was termed as "agricultural population" (PC3). The fourth PC brought together the total power of agricultural machinery per capita, the net income per capita of rural residents, and the proportion of the added value of the primary industry to GDP. Thus, PC4 was termed as "agriculture and economy". The fifth PC (PC5) and sixth PC (PC6) were loaded on SPEI and mean annual temperature, respectively. Note: TPOAM per capita is the total power of agricultural machinery per capita; AGDP/GDP is the proportion of the added value of the primary industry to GDP; SPEI is the standardized precipitation evapotranspiration index. Table 3 shows that PC2 had the largest and positive effect on grassland vegetation biomass (AGB T density and AGB m ). The standardized effect sizes of PC2 on AGB T density and AGB m were 0.47 and 0.60, respectively. PC3 had a negative and the largest effect (−0.40) on the grassland vegetation area. PC1 had a negative effect on grassland vegetation biomass and area. As a humidity factor, PC5 had a positive effect on grassland vegetation biomass. PC6 had a positive effect on AGB m . PC4, which represents agriculture and economy, had the weakest influence on grassland vegetation biomass and area compared with the population and urbanization components (PC2 and PC3) and climate components (PC1 and PC5). Table 3. Principal component stepwise regression results between six principal components (PCs) and grassland above-ground biomass (AGB) and area in the Loess Plateau (n = 401).

Principal Component Regression Model
Note: y 1 is AGB T density. y 2 is AGB m . y 3 is the proportions of the grassland area to the total area (GA/TA). R 2 is coefficient of determination. RMSE is the root mean squared error. r is correlation coefficient.

The Night-Time Lights in Different Bioclimatic Zones
The mean value of NTL in the other three bioclimatic zones (deserts, desert-grasslands, and typical grasslands zone) was lower than that in the forest-grassland zone ( Figure 5). Meanwhile, the increasing rate of NTL accelerated since the full implementation of the GFG in 2000. Table 3. Principal component stepwise regression results between six principal components (PCs) and grassland above-ground biomass (AGB) and area in the Loess Plateau (n = 401). is the proportions of the grassland area to the total area (GA/TA). R 2 is coefficient of determination. RMSE is the root mean squared error. r is correlation coefficient.

The Night-Time Lights in Different Bioclimatic Zones
The mean value of NTL in the other three bioclimatic zones (deserts, desert-grasslands, and typical grasslands zone) was lower than that in the forest-grassland zone ( Figure 5). Meanwhile, the increasing rate of NTL accelerated since the full implementation of the GFG in 2000.

The Role of Urban Population and Urbanization
As shown in Table 2; Table 3, PC2 was positively correlated with the urbanization rate and density of the urban population. Meanwhile, PC2 had the most important and positive effect on grassland vegetation biomass. It suggested that the increase in urbanization rate drove the increase in grassland vegetation biomass in the Loess Plateau. Urbanization is an important driver of ecosystem function change [18,47]. Furthermore, precipitation, temperature, and economy development may also be the drivers [24,38]. Comparison of our findings with previous studies indicated that urbanization played a dominant role in the increase in grassland vegetation biomass. Most previous studies did not consider comparing the relative important among different drivers of ecosystem function change [70,71]. In this study, we found that urbanization was the most important driver of increase in grassland vegetation biomass in the Loess Plateau compared with other drivers. Therefore, urbanization had the positive influence on vegetation restoration in the Loess Plateau.

The Role of Urban Population and Urbanization
As shown in Table 2; Table 3, PC2 was positively correlated with the urbanization rate and density of the urban population. Meanwhile, PC2 had the most important and positive effect on grassland vegetation biomass. It suggested that the increase in urbanization rate drove the increase in grassland vegetation biomass in the Loess Plateau. Urbanization is an important driver of ecosystem function change [18,47]. Furthermore, precipitation, temperature, and economy development may also be the drivers [24,38]. Comparison of our findings with previous studies indicated that urbanization played a dominant role in the increase in grassland vegetation biomass. Most previous studies did not consider comparing the relative important among different drivers of ecosystem function change [70,71]. In this study, we found that urbanization was the most important driver of increase in grassland vegetation biomass in the Loess Plateau compared with other drivers. Therefore, urbanization had the positive influence on vegetation restoration in the Loess Plateau.
PC3 was positively correlated agricultural population (Table 2) and it had a negative effect on grassland vegetation area. It indicated that agricultural population decreases led to increases in the grassland vegetation area. Urbanization means that more and morerural residents moved to cities and towns [72,73]. Meanwhile, farmland is more inclined to be abandoned in those areas with higher levels of urbanization [74]. Hence, the decrease in agricultural population caused by increasing urbanization rate promoted the increase of grassland vegetation area in the Loess Plateau. Collectively, these results indicated that urbanization was the critical influence factors in the sustainable management of ecosystems in the Loess Plateau.

Dominant Driving Factors at Spatiotemporal Scales
The findings from the SEM and PCA revealed that the spatiotemporal variations of grassland vegetation were mainly regulated by population (urban and agricultural) and urbanization in the Loess Plateau. Moreover, the driving factors on the dynamics of grassland vegetation showed a spatial difference among the four bioclimatic zones (Figure 3b-e).
In the forest-grassland zone, the increase in urbanization rate led to the decrease in grassland vegetation AGB and area ( Figure 3e); this finding is consistent with the results of studies on forest biomass carbon or terrestrial NPP [46,75]. However, the increase in urbanization rate would promote the increase of grassland vegetation in the other three bioclimatic zones (Figure 3b-d). DMSP/OLS NTL data could be used to explain this phenomenon: the more brightness of night-time lights indicated that more intensive population and higher development of urbanization [76,77] (Figure 5 and Figure  S5). Thus, we argue that urbanization plays an important role in gathering population in the other three bioclimatic zones with sparse population, where grassland vegetation increased due to the eased ecological pressure. On the contrary, the increase in population and urbanization would squeeze the space of grassland vegetation in the forest-grassland zone with greater brightness of NTL and higher NTL density ( Figure 5 and Figure S5). Thus, the effects of the increase in urbanization rate on the dynamics of grassland vegetation were negative in the forest-grassland zone and positive in the other three bioclimatic zones.
We identified the increase in urbanization rate as the dominant driving factor, although the influences of temperature and humidity on the dynamics of grassland vegetation gradually increased with time ( Figure 4). The brightness of NTL has rapidly increased in the Loess Plateau in recent years ( Figure 5 and Figure S5). Thus, most areas of the Loess Plateau have entered the stage of rapid urbanization. At the same time, an increasing number of rural residents have moved to cities and towns [73]. Thus, we speculated that the disturbances of human activities to grassland vegetation gradually reduced in the rural region, where the vast majority of the Loess Plateau is located. Correspondingly, the grassland vegetation increased naturally, and its growth was largely driven by climate factors in recent years [13].

The Response Mechanisms of Grassland Vegetation Dynamics to Urbanization
Growing evidence suggested that the vegetation in the Loess Plateau has increased after the implementation of the GFG Program [3,78]. Our study revealed that the increase in urbanization rate is the underlying and dominant factor driving the increase of grassland vegetation in the Loess Plateau. Many rural households have had a surplus of labor forces who might become migrant workers or be involved in other economic sectors since the implementation of the GFG Program in the region [72,73,79]. On the one hand, GFG can be considered as a "trigger" that accelerated the pace of rural residents' migration to urban areas. On the other hand, most rural migrant workers have adapted to urban life in search for high income and good education for their offsprings, and few people are intent on destroying vegetation that has been restored [80,81].
As shown in Figure S6, the driving mechanism of the urban population and urbanization that is related to the increase of grassland vegetation is as follows: the urbanization process promotes the population shift from rural to urban areas, thereby deeply reducing the ecological pressure of the vast countryside areas in the Loess Plateau and facilitating the increase of grassland vegetation. Moreover, urbanization ensures the success of the GFG program. Specifically, urbanization leads to more farmlands being abandoned. Thus, grassland vegetation increases accordingly [82].

Conclusions
In this study, we investigated the specific factors driving the increase of grassland vegetation and assessed their impacts on grassland AGB and area in the Loess Plateau of China. Both SEM and PCA analyses both showed that the increase in population and urbanization is the most important driving factor for the increase of grassland vegetation in the region. Specifically, the increase in population and urbanization exerts positive effects on the increase of grassland vegetation. The urbanization process promotes the population shift from rural areas to urban areas, thus reducing the ecological pressure of the vast countryside areas in the Loess Plateau and facilitating the increase of grassland vegetation. The influences of temperature and humidity have gradually increased in recent years, but human activities, especially population change and urbanization, were still critical and indispensible influence factors in the Loess Plateau. Our study suggested that the effect of urbanization should be considered when assessing vegetation change.
Supplementary Materials: The following are available online at http://www.mdpi.com/2071-1050/11/23/6764/s1, Figure S1: Relationships between NDVI, elevation and AGB of sampling sites plot., Figure S2: Conceptual model depicting the proposed relations among agriculture production and economic development (AP&ED), population changes and urbanization (PC&U), interaction of AP&ED and PC&U (High-order interaction), temperature and humidity degree and dynamics of grassland vegetation. The measurement model relates the indicator to the latent variables. The structural model relates all latent variables. Manifest variables: net income per capita of rural residents, the density of livestock on hand, the proportion of the added value of the primary industry to GDP (AGDP/GDP), total power of agricultural machinery per capita; the density of agricultural population, the density of urban population, urbanization rate; density of stable nighttime light (NTL density); mean annual precipitation, mean annual temperature, sunshine percentage; AGB T density, AGB m , the proportions of the grassland area to the total area (GA/TA)., Figure S3: Inter-annual average temperature( • C) of all study county changes in: (a) deserts zone; (b) desert-grasslands zone; (c) typical grasslands zone; (d) forest-grasslands zone., Figure S4: Inter-annual average precipitation (mm) of all study county changes in: (a) deserts zone; (b) desert-grasslands zone; (c) typical grasslands zone; (d) forest-grasslands zone., Figure S5: The DMSP/OLS night-time light intensity images of the Loess Plateau in different year., Figure S6: The inherent relation between grassland vegetation increase and urbanization in the Loess Plateau. Table S1: The specific year when the remote sensing images were used to perform the supervised classification for each study counties.
Author Contributions: J.-Z.W., K.Z., F.Z. and J.-S.Y. conceived and designed this study, J.-Z.W. wrote the manuscript. All authors revised and approved the manuscript.