Impacts of Land Use Changes on Net Primary Productivity in Urban Agglomerations under Multi-Scenarios Simulation

: Land use is closely related to the sustainability of ecological development. This paper employed a patch-generating land use simulation (PLUS) model for the multi-scenario simulation of urban agglomerations. In addition, mathematical analysis methods such as Theil-Sen Median trend analysis, R/S analysis, Getis-Ord Gi* index and unary linear regression were used to study the temporal and spatial evolution characteristics of net primary productivity (NPP) for the impact of land use changes on NPP in urban agglomerations from 2000 to 2020 and to forecast the future trend of NPP. The results indicate that urban expansion is obvious in the baseline scenario and in the ecological protection scenario. In the scenario of cropland protection, the urban expansion is consistent with the land use plan of the government for 2035. The NPP in Beijing decreased gradually from northwest to southeast. The hot spot areas are concentrated in the densely forested areas in the mountainous areas of northwest. The cold spot areas are mainly concentrated in the periphery of urban areas and water areas. The NPP will continue to increase in forest and other areas under protection and remain stable in impervious surfaces. The NPP of Beijing showed a strong improvement trend and this trend will continue with the right ecological management and urban planning of the government. The study of land use in urban agglomeration and the development trend of vegetation NPP in the future can help policymakers rationally manage future land use dynamics and maintain the sustainable development of urban regional ecosystems.


Introduction
Vegetation changes have recorded the imprint of human activities profoundly [1]. Human activities such as land use transformation directly change the structure, process and function of ecosystems [2][3][4]. In human activity areas, vegetation has been lost or altered for human use (e.g., commercial development and infrastructure construction), and this poses a major threat to biodiversity conservation [5]. Net primary productivity (NPP) is defined as the accumulation rate of organic matter by vegetation, which is equal to the difference between the carbon absorbed by plant photosynthesis and the carbon consumed by plant autotrophic respiration [6,7], and it plays an important role in global carbon balance [8]. Land use changes have had a wide and significant impact on the ecological environment, especially in urban agglomeration where urban expansion is one of the main driving forces of change in NPP [4,9]. So, the temporal and spatial distribution characteristics of dynamic changes of land use and NPP have become an important research content in regional eco-environmental assessment.
The urbanization process in metropolises may offset the positive effects on NPP [3,10]. Vegetation destruction and improper land use caused by human activities have caused serious land degradation and ecological environment problems [11]. However, Li, X. [7] found that the NPP of new urban land in 2019 was better than that of original urban land before 2010. The human activities, mainly ecological restoration projects, made the vegetation coverage on the Loess Plateau show a significant upward trend in the past three decades [12]. As a metropolitan center, Beijing City is the capital of China with great land use transformations which is worthwhile researching. Discussing the land use transformation on the influence of vegetation NPP and whether the future changes are consistent are essential for alleviating the conflicts between sustainable development and environmental protection.
Previous studies of NPP have mainly focused on the two parts. The first part is the establishment, application and optimization of the NPP calculation model [13,14]. The other part is exploring the response mechanism of NPP by taking climate factors and human activities as the main driving factors [15,16]. Meanwhile, the change trend and the future development of NPP are also worth studying. The combination of Theil-Sen Median trend analysis and Mann-Kendall test is used to detect the trend of spatial variation in NPP widely [7,17]. This combination can eliminate the noise of original data and are very robust [17]. Brilli, L. [18] predicted a decrease in net ecosystem exchange and NPP by approaching the most extreme climate conditions. Sun, G. [19] assessed the characteristics of NPP due to future climate change and CO2. These simulation results are realized by predicting the driving factors of NPP instead of the NPP itself. Other scholars simulate NPP in the future through the change trend of NPP in the past [20][21][22]. In addition to analyzing the changing trend of NPP and its driving factors, it is necessary to study the future trend of them during rapid urbanization.
To predict the future land use (one of the driving factors of NPP), many kinds of land use simulation models have been developed, such as cellular automata (CA) [23][24][25][26], conversion of land use and its effects modeling framework (CLUE) model and its improved model [27][28][29], future land use simulation (FLUS) model [30] and patch-generating land use simulation (PLUS) model [31,32]. However, the CA model is weak at revealing the underlying drivers of land use change [33]. CLUE-S model leads to the separation of the macro land use demand projections and the local change allocations [34]. Thus, there is an urgent need for more comprehensive approaches to predict the land use in Beijing. At present, in order to pursue high accuracy, the mainstream of land prediction models has been the coupling model to eliminate the limitations of a single model by carrying out the simulation from two aspects: quantity and spatial distribution. The PLUS model has been proven to be a more effective model, which provides more precise simulation results than the FLUS model [31]. A Markov chain predicts dynamic variation characteristics and is characterized by both high operational precision and high prediction accuracy [30]. This paper used a Markov chain to predict the amount of land use and simulated spatial changes by PLUS model on the basis of physiographic driving factors, transportation driving factors and socioeconomic driving factors. In addition, researchers have evaluated the future spatiotemporal land use under different scenarios combined with simulation models [24,30,31,35]. In this research, different transition matrix restricted different types of land use transformations to obtain different scenarios.
Based on mathematical analysis methods such as unary linear regression, Getis-Ord Gi* index and Hurst index, the objectives of this study were to: (1) explore the spatiotemporal variation trend of vegetation NPP from 2000 to 2020 based on land use; (2) evaluate future land use at 30 m spatial resolution under three different scenarios in 2030; (3) analyze the future trend of NPP. The temporal and spatial variation in land use and NPP were comprehensively analyzed for ecological protection and management in this research. What is more, based on the impact of land use on vegetation NPP, the correctness of NPP change trend in the future can be judged. This research is expected to provide a scientific basis for spatial development plan and represent a valuable reference of the urban agglomeration similar to Beijing.

Study Area
Beijing is located at 115 • 25 -117 • 30 E, 39 • 28 -41 • 05 N and has a total area of 16,410 square kilometers approximately. It is located in the northern part of the North China Plain and surrounded by the Taihang Mountains and Yanshan Mountains (Figure 1). Beijing's terrain is high in the northwest and low in the southeast with an average elevation of 43.5 m above sea level. The mountainous area in this administrative region accounts for 62%, mainly distributed in middle and low mountainous areas. With a typical warm temperate sub-humid semi-arid monsoon climate, Beijing is hot and rainy in summer, cold and dry in winter, and spring and autumn are short. The zonal vegetation type in Beijing is warm temperate deciduous broad-leaved forest with temperate coniferous forest.

Study Area
Beijing is located at 115°25′-117°30′ E, 39°28′-41°05′ N and has a total area of 16410 square kilometers approximately. It is located in the northern part of the North China Plain and surrounded by the Taihang Mountains and Yanshan Mountains (Figure 1). Beijing's terrain is high in the northwest and low in the southeast with an average elevation of 43.5 m above sea level. The mountainous area in this administrative region accounts for 62%, mainly distributed in middle and low mountainous areas. With a typical warm temperate sub-humid semi-arid monsoon climate, Beijing is hot and rainy in summer, cold and dry in winter, and spring and autumn are short. The zonal vegetation type in Beijing is warm temperate deciduous broad-leaved forest with temperate coniferous forest.
As the capital of China, Beijing has strengths in the fields of politics, culture, traffic and science and technology. Impervious surface is one of the land use types which has emerged not only as an indicator of the degree of urbanization, but also a major indicator of environmental quality [36]. The other main land use type in Beijing is forest which covers the largest area and mainly occurs in areas with elevations of 1000-1500 m. The local government is facing a challenge in balancing ecological conservation and local livelihoods under land use change.

Data
The border location data regarding the geographical limits of the research region were downloaded from the Chinese Academy of Sciences Resource and Environmental Science Data Center (http://www.resdc.cn, accessed on 30 October 2021). As the capital of China, Beijing has strengths in the fields of politics, culture, traffic and science and technology. Impervious surface is one of the land use types which has emerged not only as an indicator of the degree of urbanization, but also a major indicator of environmental quality [36]. The other main land use type in Beijing is forest which covers the largest area and mainly occurs in areas with elevations of 1000-1500 m. The local government is facing a challenge in balancing ecological conservation and local livelihoods under land use change.

Data
The border location data regarding the geographical limits of the research region were downloaded from the Chinese Academy of Sciences Resource and Environmental Science Data Center (http://www.resdc.cn, accessed on 30 October 2021). Physiographic factors, transportation factors and socioeconomic factors were determined to be development potential as spatial driving factors of the future land use changes in Beijing. The driving factors between 2000 to 2010 were used to simulate the land use in 2020 for accuracy verification. Correspondingly, the driving factors between 2010 to 2020 were used to simulate the land use in 2030. Referring to the previous studies [31,41], we selected several forces (Table 1) as driving factors to predict the impacts of natural and socioeconomic conditions on land use changes. Among them, DEM was used as the main data source for more driving forces (slope, aspect). The data set was provided by Geospatial Data Cloud site, Computer Network Information Center, Chinese Academy of Sciences (http://www.gscloud.cn, accessed on 30 October 2021), processed by using the data of the first version (V1) of Aster GDEM in 2009. It is a digital elevation data product with a global spatial resolution of 30 m.
(2) Vegetation net primary productivity data. The MOD17A3HGF NPP products [42] came from the Terra Moderate Resolution Imaging Spectroradiometer (MODIS) which is a satellite-based sensor used for earth and climate measurements. MODIS data and products acquired through the LP DAAC. The temporal resolution is 1 year, and the spatial resolution is 500 m. The MOD17A3HGF is derived from the sum of all 8-day net photosynthesis (PSN) products (MOD17A2H). The PSN value is the difference of the Gross Primary Productivity (GPP) and the Maintenance Respiration (MR). We cut and mosaicked the band of NPP from 2000 to 2020 and the pixel value is multiplied by the scale factor of 0.0001 to obtain the annual real value of NPP. The pixel sizes of all data are unified to 30 m, and the applied geographic coordinate system and map projection were the WGS-84 datum and the Universal Transverse Mercator, respectively. The PLUS model contains a land expansion analysis strategy (LEAS) and a CA model that includes a patch-generation mechanism based on multi-type random seeds (CARS). The LEAS transforms the mining of transition rules of each land use type into a binary Remote Sens. 2022, 14, 1755 6 of 21 classification problem (obtain the change and inertia probabilities of each land use type), which can be solved with many data mining methods [31]. The LEAS transform needs land use data of two periods to overlay and extract changing grids from later data, then selected random sampling points. Then the random forest classification (RFC) algorithm analyzes the relationship between the growth of each kind of land use types and various driving factors to obtain the transition rules by calculating the growth probability. In this research, 16 driving factors (Table 1) were chosen to calculate the probability of occurrence for each landscape type among these factors. The CARS module is a Cellular Automata (CA) model base on multi-type random patch seeds. The important property of CA is that they demonstrate the spatial and dynamic process [24]. The CA model in PLUS model has been improve with multi-type random patch seeds and threshold descent rules.
Neighborhood weight represents the expansion intensity of the land type, that is, the expansion ability of each land type driven by external factors. Arnold Jr, C.L. [43] conducted dimensionless processing on the variation of Total Area (one of the landscape metrics, TA) of each land type, making its threshold between 0 and 1. Combined with the requirements of neighborhood weight on data structure, it is concluded that the dimensionless value of TA variation conforms to the requirements of neighborhood weight both in terms of parameter meaning and data structure. Formula is as follows: where X = (X 1 , · · · , X i ); X * is the i th normalized value. The Markov chain model uses the Markov process method to project future land which has been widely used in the prediction of land demands [30,31,43]. This method is characterized by both high operational precision and high prediction accuracy. It can be presented using the following equation as: where S t+1 and S t are the landscape composition at times t and t + 1, respectively; P ij is the transition probability that land use type i is converted to land use type j at time t.
Firstly, the PLUS model was used to simulate the land use in 2020 with the land use in 2000 as the starting year and 2010 as the end. The simulation results were compared with the real 2020 land use map via a Kappa test, as follows: where p o represents overall accuracy which is calculated by dividing the sum of the correctly classified samples in each category by the total number of samples; p e represents the correct proportion of the model in the random case; p p represents the proportion of the correct simulation in the case of ideal classification which is 1 normally. The verification results revealed that the overall accuracy and kappa coefficient reached 0.8268 and 0.747233, respectively. This demonstrates that the accuracy of the model is substantial, and the PLUS model can be used for future simulation reliably in Beijing. Secondly, data from 2010 to 2020 were used to simulate future land use under different scenarios. In this study, different scenarios included the baseline scenario, cropland protection scenario and ecological protection scenario [35]. For the Chinese government's policy, these scenarios were all simulated over the planning constraints zones and development zones. The baseline scenario is the basis of considering other constraint conditions for simulation of land use change in urban agglomerations. Based on the rule of land use change from 2010 to 2020, no restriction on conversion between different types of land is set. The transition matrix of the baseline scenario is unit matrix. The cropland protection scenario protects the cropland in urban agglomeration and strictly controls the conversion of cropland to other types of land which is set as 0 in the transfer matrix. The ecological Remote Sens. 2022, 14, 1755 7 of 21 protection scenario is similar to the cropland protection scenario. The factors of ecological protection are added into the baseline scenario by setting forest, wetland and water body not to be converted into other types of land use. Neighborhood weights and land demands are reported in Table 2. For other specific input parameters and simulation processes, see reference [31]. The overall framework of scenario-based future land use and NPP is illustrated in Figure 2.   trend analysis does not require data to obey specific distribution characteristics and is less affected by outliers. This method quantifies the change trend by calculating the median slope of n(n − 1)/2 data combinations: where β is the changing trend of NPP, β > 0 indicates an upward trend, while β < 0 indicates a downward trend; NPP i and NPP j represent the NPP values of pixels in year i and year j, respectively. The Mann-Kendall test method is a nonparametric statistical test used to judge the significance of trends, which does not require samples to follow a certain distribution and is not disturbed by a few outliers [44]. This method is used to evaluate the significance of the NPP trend, and the formula is as follows: where S is the test statistics; sgn is a sign function; n is the length of the time series; Z is the standardized statistics. The trend test method is the null hypothesis. Under a given significance level α, ±Z 1−α/2 is normalized normal deviation. When |Z| > Z 1−α/2 , refuse the null hypothesis and it indicates that there are changes at level α. On the contrary, the change is not significant. In this paper, significant change was defined when the trend was significant at α = 0.05. When α = 0.01, it was extremely significant.

Rescaled Rang Analysis (R/S Analysis) and Hurst Index Method
The Hurst index method can predict the future development trend according to the past trend of long time series data. The Hurst index method based on the R/S analysis method was first developed by a British hydrologist [45] in the study of Nile reservoir flow. It is widely used in hydrology, economics and climatology [46,47], and has recently been applied to the study of vegetation cover change [21,44,47]. The calculation steps are as follows: where NPP i is the NPP value in year i; k is constant; H is the Hurst index which can be calculated by linear fitting with the least squares in double logarithmic coordinate system (lnm, ln(r/s)). The norm of the Hurst index is 0~1. When 0 < H < 0.5, it indicates that the NPP will reveal a trend that is opposite of that in the past; when H = 0.5, it indicates that the change of NPP is mutually independent which is a random event; when 0.5 < H < 1, it indicates that there are long-term correlations and the NPP will be consistent with the changing trend before.

Hot Spot Analysis
Prior to the implementation of hot spot analysis, the NPP changes of vegetation in Beijing from 2000 to 2020 has used Moran's Index to analyze spatial autocorrelation ( Figure 3). The result shows that it has the characteristics of spatial clustering, which can be carried out for spatial cold hot spot analysis.
Remote Sens. 2022, 14, x FOR PEER REVIEW 9 of 21 the change of NPP is mutually independent which is a random event; when 0.5 < H < 1, it indicates that there are long-term correlations and the NPP will be consistent with the changing trend before.

Hot Spot Analysis
Prior to the implementation of hot spot analysis, the NPP changes of vegetation in Beijing from 2000 to 2020 has used Moran's Index to analyze spatial autocorrelation (Figure 3). The result shows that it has the characteristics of spatial clustering, which can be carried out for spatial cold hot spot analysis. Hot spot and cold spot analysis is performed to identify the states having the high or low values of cluster spatially by Getis-Ord Gi* statistic. The local sum for a feature in the neighborhood of a grid cell is compared proportionally to the sum of all features using Getis-Ord Gi* statistic which can obtain the Z-score [48]. For statistically significant, the positive and larger Z-score indicated the more intense the clustering of high values (hot spot) and negative and the smaller the Z-score signified the more intense the clustering of low values (cold spot) [49,50]. The formula is as follows: where is the attribute value for feature ; ( ) is the spatial weight between feature and ; ( * ) and ( * ) is the mathematical expectation and coefficient of variation of * , respectively. If ( * ) > 0 and significant, it indicates that the value around the re- Hot spot and cold spot analysis is performed to identify the states having the high or low values of cluster spatially by Getis-Ord Gi* statistic. The local sum for a feature in the neighborhood of a grid cell is compared proportionally to the sum of all features using Getis-Ord Gi* statistic which can obtain the Z-score [48]. For statistically significant, the positive and larger Z-score indicated the more intense the clustering of high values (hot spot) and negative and the smaller the Z-score signified the more intense the clustering of low values (cold spot) [49,50]. The formula is as follows: where X j is the attribute value for feature j; W ij (d) is the spatial weight between feature i and j; E G * i and Var G * i is the mathematical expectation and coefficient of variation of G * i , respectively. If Z G * i > 0 and significant, it indicates that the value around the research area is higher than the mean value and the area belongs to the hot zone. On the other hand, the value around the research area is low and the area belongs to the cold spot area.

Future Spatial-Temporal Land Use Changes in Beijing
During 2000-2020, the land use type of urban agglomeration was mainly forest, which accounted for more than 40% of the entire area, followed by cropland and impervious surface. Water body, wetland and bare area accounted for only a small portion of the study area ( Figure 4). The total area of grassland and cropland decreased, among which the area of cultivated land decreased by 1764.68 km 2 in 20 years, the proportion of grassland decreased by about 2%. The land use of water body first decreased and then rose, because the closer it is to the water, the more suitable it is for production and living, and the easier it is to be exploited and utilized by human beings. The water area of Beijing is mainly distributed in Miyun, and the area around the reservoir was transformed into cropland for planting chestnut, plum and other cash crops in early years. People gradually realized the importance of water body, therefore, the area of water body increased by 87.14 km 2 from 2010 to 2020. In 2010, the construction land area of urban agglomeration was 3475.76 km 2 , increasing by 1885.80 km 2 compared with 2000 which tended to show clear changes. The expansion area of impervious surface was mainly concentrated in the urbanization area of urban agglomeration, including Haidian district, Chaoyang District, Daxing District and Tongzhou District.
The mountain terrain in the northwest region is not suitable for large-scale development, except Yanqing district, which is relatively flat and has more cropland and impervious surface with frequent human activities. Additionally, people's requirements for urban environment have been enhanced with the continuous development of urbanization and the continuous improvement of residents' living standards, so large amounts of mountainous land and farmland have been listed as northwest ecological conservation area. Similarly, cropland is vital. Therefore, land use simulation is divided into the following three scenarios: the baseline scenario, cropland protection scenario and ecological protection scenario.
Under the baseline scenario in 2030 (Figure 4d), the area of cropland will decrease by 955.04 km 2 compared with 2000 and the area of other land use types will increase. The area of impervious surface, water body, grassland and forest will increase by 862.43 km 2 , 48.67 km 2 , 23.91 km 2 and 19.07 km 2 , respectively, and the area of bare land and wetland will hardly rise. The baseline scenario is the basis for the simulation of urban agglomeration land use change considering other constraints, while the cropland protection scenario is to add the concept of farmland protection to the baseline scenario to prevent the basic farmland from being occupied by other land in the process of urbanization. As a result, the land use under cropland protection scenario (Figure 4e) has a little change than 2020. Specifically, the area of forest will change the most (from 7388.60 km 2 in 2020 to 7382.57 km 2 in 2030) and other land use areas will change by less than 3 km 2 under cropland protection scenario. The simulation results under the ecological protection scenario (Figure 4f) are similar to that under the baseline scenario. The areas of forest, grassland, wetland, impervious surface and bare land will be 7407.67 km 2 , 1395.65 km 2 , 12.04 km 2 , 4338.19 km 2 and 4.89 km 2 , respectively, under these two scenarios. Under the ecological protection scenario, more of the area will be converted to water body (288.61 km 2 , accounting for approximately 1.76% of the total inundation area), as compared to the baseline scenario in 2030, the area of cropland will reduce by another 1.61 km 2 .
cropland for planting chestnut, plum and other cash crops in early years. People gradually realized the importance of water body, therefore, the area of water body increased by 87.14 km 2 from 2010 to 2020. In 2010, the construction land area of urban agglomeration was 3475.76 km 2 , increasing by 1885.80 km 2 compared with 2000 which tended to show clear changes. The expansion area of impervious surface was mainly concentrated in the urbanization area of urban agglomeration, including Haidian district, Chaoyang District, Daxing District and Tongzhou District.

Spatial Variation Characteristic
As shown in Figure 6a and Table 3, changes of vegetation NPP showed an obvious upward trend from 2000 to 2020 in Beijing. The NPP rising area accounted for 86.77% of the total area, showing a ring shape. The extremely significant increase (p < 0.01) accounted for about 70.62% of the total area, mainly distributed in the west of Fangshan District, Mentougou District, north of Changping District, Yanqing District, Huairou District, Miyun District, Pinggu District, Shunyi District, Tongzhou District, Daxing District. Among the above areas, Shunyi District, Tongzhou District, Daxing District and Pinggu District are flat and far from the urban agglomeration, and the land use type is given priority to cropland and grassland. In other areas, the land use type is mainly forest and grassland with high terrain and less negative human disturbance, however, the changes of vegetation NPP of the mountain ridge area belongs to slightly rise and significantly rise. The main types of land use in the stable area are mainly impervious surface and water body which has no vegetation. The decrease of NPP accounted for 0.4% of the total area, among which, the extremely significant decline (p < 0.01) and significant decline (0.01 ≤ p < 0.05) accounted for about 10.04% of the total area, mainly distributed in the central area of urban agglomeration with flat terrain. Due to the negative influence of social and economic development, human activities and population density, vegetation activity weakened and NPP decreased.

Spatial Variation Characteristic
As shown in Figure 6a and Table 3, changes of vegetation NPP showed an obvious upward trend from 2000 to 2020 in Beijing. The NPP rising area accounted for 86.77% of the total area, showing a ring shape. The extremely significant increase (p < 0.01) accounted for about 70.62% of the total area, mainly distributed in the west of Fangshan District, Mentougou District, north of Changping District, Yanqing District, Huairou District, Miyun District, Pinggu District, Shunyi District, Tongzhou District, Daxing District. Among the above areas, Shunyi District, Tongzhou District, Daxing District and Pinggu District are flat and far from the urban agglomeration, and the land use type is given priority to cropland and grassland. In other areas, the land use type is mainly forest and grassland with high terrain and less negative human disturbance, however, the changes of vegetation NPP of the mountain ridge area belongs to slightly rise and significantly rise. The main types of land use in the stable area are mainly impervious surface and water body which has no vegetation. The decrease of NPP accounted for 0.4% of the total area, among which, the extremely significant decline (p < 0.01) and significant decline (0.01 ≤ p < 0.05) accounted for about 10.04% of the total area, mainly distributed in the central area of urban agglomeration with flat terrain. Due to the negative influence of social and economic development, human activities and population density, vegetation activity weakened and NPP decreased.

Future Trend
The average Hurst index of NPP of vegetation in Beijing is 0.64, which shows a continuous change on the whole. About 85.12% of the total area show continuous change, which is much higher than the area with anti-continuous change (2.05%).
The results of the R/S analysis and Hurst index method were reclassified and performed spatial overlay analysis in order to reveal the future NPP in Beijing, so the persistence characteristics of future vegetation change trends can be obtained (Figure 6b and Table 4). According to statistical analysis, the vegetation NPP of most of the total area (84.8%) will rise or continuous rise in the future, indicating that the ecological protection of Beijing will achieve certain results on the basis of consolidating the existing protection achievements. The decline and continuous decline areas account for 2.37%, and the areas where the NPP will continue to decline in the future are mainly the areas where the Daxing Airport, Beijing Capital Airport, flyovers and other traffic facilities are located. The areas of NPP decline in the future are not only around the areas of continuous decline, but also at the ridges in the west and north of Beijing where has shown a slowly improvement in the past two decades.  Figure 7 shows the spatial distribution information of high value (hot spot) and low value (cold spot) of NPP change trend. When p < 0.01, the hot spots of NPP change are

Future Trend
The average Hurst index of NPP of vegetation in Beijing is 0.64, which shows a continuous change on the whole. About 85.12% of the total area show continuous change, which is much higher than the area with anti-continuous change (2.05%).
The results of the R/S analysis and Hurst index method were reclassified and performed spatial overlay analysis in order to reveal the future NPP in Beijing, so the persistence characteristics of future vegetation change trends can be obtained (Figure 6b and Table 4). According to statistical analysis, the vegetation NPP of most of the total area (84.8%) will rise or continuous rise in the future, indicating that the ecological protection of Beijing will achieve certain results on the basis of consolidating the existing protection achievements. The decline and continuous decline areas account for 2.37%, and the areas where the NPP will continue to decline in the future are mainly the areas where the Daxing Airport, Beijing Capital Airport, flyovers and other traffic facilities are located. The areas of NPP decline in the future are not only around the areas of continuous decline, but also at the ridges in the west and north of Beijing where has shown a slowly improvement in the past two decades.   Land use change is one of the important driving factors of NPP. According to the cold hot spot area of vegetation NPP in Beijing as the mask extraction range, the transfer matrix of land use change in cold spot and hot spot area was extracted, so as to understand the impact of land use change on vegetation NPP change. It can be seen from Table 5 that the main land types in the hot spot of vegetation NPP in Beijing are forest and grassland, and the area of wetland and bare area is extremely few (0.019 km 2 and 0.022 km 2 in 2020, respectively). In terms of land use change, there is little transfer. The area of forest and grassland without land transformation are 1148.492 km 2 and 43.427 km 2 . The transformation of land use types in hot spots mainly occurs in the transformation from grassland to forest about 21.816 km 2 . Although other lands are still transferred to water body and impervious surface, the change area is small which is only 6.364 km 2 after the transformation. The main land use types in this area are still forest land and grassland.  Land use change is one of the important driving factors of NPP. According to the cold hot spot area of vegetation NPP in Beijing as the mask extraction range, the transfer matrix of land use change in cold spot and hot spot area was extracted, so as to understand the impact of land use change on vegetation NPP change. It can be seen from Table 5 that the main land types in the hot spot of vegetation NPP in Beijing are forest and grassland, and the area of wetland and bare area is extremely few (0.019 km 2 and 0.022 km 2 in 2020, respectively). In terms of land use change, there is little transfer. The area of forest and grassland without land transformation are 1148.492 km 2 and 43.427 km 2 . The transformation of land use types in hot spots mainly occurs in the transformation from grassland to forest about 21.816 km 2 . Although other lands are still transferred to water body and impervious surface, the change area is small which is only 6.364 km 2 after the transformation. The main land use types in this area are still forest land and grassland. As can be seen from Table 6, the area of impervious surface without land use transformation in cold spot (p < 0.01) area is 708.79 km 2 , and the area of other land use types transformed into impervious surface is 1532.14 km 2 in total, leading to a low level of NPP in this region. Land use transformation in cold spot area mainly occurs between cropland to impervious surface, forest to grassland and grassland to other land use types. In 2000, 1305.30 km 2 cropland was converted to impervious surface in 2020, while 37.92 km 2 forest land, 183.51 km 2 grassland along with 5.26 km 2 water body was also turned into impervious surface. With the continuous strengthening of urban development in Beijing, the rapid expansion of construction land areas is the main reason for the cold spot areas of NPP.

Spatial Development Plan
In the twenty-first century, preparing cities for urban land use planning is a critical challenge, and the relevant authorities require guidance from the simulation results to define adaptation strategies and policies. In accordance with "The Master Plan of Beijing (2016-2035)" by The People's Government of Beijing City, natural resources such as green space, water system and wetland are important ecological barriers and great attentions is attached to the protection of this famous historical and cultural city, therefore, these zones are the planning constraints zones in this study. The development area is also obtained from the centralized construction area in the released "planning map of two lines and three districts in the city". Compared with the current land use situation, areas of the north of Daxing District, areas around Beijing Daxing International Airport and so on will be intensively constructed.
According to "The Master Plan of Beijing (2016-2035)", the total scale of construction land (including urban and rural construction land, special land, external transportation land and some water conservancy facilities land) should be controlled within 3720 km 2 by 2020; the amount of cropland would be no less than 1.66 million mu by 2020; the city's forest coverage rate would increase from the current 41.6% to 44% by 2020. In the actual land use in 2020, the area of impervious surface was 3475.76 km 2 ; the area of cultivated land was 3902.51 km 2 (5.85 million mu); and the area of forest land was 7388.60 km 2 (45.07%). The current regional development met the planning requirements, that is to say, the future development guidance of this study can refer to the land development planning.
From 2000 to 2010, the NPP value of vegetation in Beijing fluctuated greatly ( Figure 5). Then China has implemented a series of ecological forestry projects, such as Grain for Green Program, the Three-North Shelterbelt Project, Wildlife Protection and Nature Reserve Development Program around Beijing which have made remarkable achievements in vegetation restoration and ecological system construction [21,52].The future work plan is to establish an ecological restoration support mechanism by strengthening the ecological restoration of cropland forest network and river lake wetland in plain areas, and go on carrying the ecological forestry projects to continuously improve the quality of ecological space.

Review of Land-Use Simulation Results
The change of impervious surface area is the most obvious under three different scenarios, which is the main feature of land use change in urban agglomeration in the future and shows the mode of expanding from the existing construction land to the periphery.
"The Master Plan of Beijing (2016-2035)" requires the total scale of the city's construction land (impervious surface, including urban and rural construction land, special land, external transportation land and some water conservancy facilities) will be controlled about 3670 square kilometers by 2035. The area of impervious surface under cropland protection scenario (Figure 4d) meets the requirements, nevertheless, there is still a certain gap with the centralized construction area planned in the future, such as the northernmost areas around Beijing Daxing International Airport which should be developed and constructed. Comparing the three scenarios, the baseline scenario has the most severe urban expansion and has more negative impact on NPP; cropland protection scenario is the only scenario in which cropland is not transformed into impervious surface in a large area, but all kinds of land basically maintain the current situation, lacking the space for urbanization development; the ecological protection scenario has good protection for forest, wetland and water body, but there is a situation that a number of cropland is turned into impervious surface. Therefore, the land use planning of urban agglomeration should comprehensively consider the scenarios of ecological protection and cropland protection, which can not only ensure regional ecological security, but also improve the quality of cultivated land and protect food security. The change of land use in the process of urbanization in Beijing is representative. The land use transfer in Beijing is similar to that in Europe from 1950 to 1990: the proportion of forest and semi-natural areas increased slightly; artificial surfaces and agricultural areas changed greatly, in which the proportion of artificial surfaces increased and the proportion of agricultural land decreased [53]. The results suggested different processes are taking place in different kinds of land use: the mountains northwest of Beijing are dominated by forest management; intensification and urbanization are mainly encountered on the plain. However, the international community has increasingly recognized the importance of maintaining "sustainable" ecosystems to protect the environment on which we depend. The researchers explored the methods and tools of integrated landscape management for nature conservation, physical/territorial planning, watershed management, land arrangement projecting and forestry planning [54,55].

Effect on Vegetation Net Primary Productivity (NPP)
On the one hand, the change of land use from one type to another is often accompanied by a large amount of carbon exchange. For example, when forests with high biomass are converted into farmland, grassland or cities with low biomass, large amounts of CO 2 will be released into the atmosphere which may cause a decrease in NPP of vegetation. The expansion of global croplands, pastures, plantations, and urban areas in recent decades generated considerable losses of biodiversity and also potentially undermined the capacity of ecosystems to sustain food production [56]. in the process of industrialization in Beijing, a large amount of land is being changed into construction land, resulting in the increase of impermeable water surface. According to the result of hot spot analysis (Figure 7), when p < 0.01, the cold spots are mainly concentrated in the periphery of urban area and water area which may taking negative effects to NPP. On the other hand, in addition to natural factors, policy regulation and economic driving are also important causes of land use change and its spatial and temporal disparities. Typical examples of this are the abandonment of the former republics following the upheaval in Eastern Europe and the collapse of the Post-Soviet Union [57]. With the continuing development of the city, urban planning priorities have gradually changed, the urban greening intensity has increased, and the vegetation condition has steadily improved [58]. It can be seen from Table 5 that the main land types in the hot spot of vegetation NPP are forest and grassland which are mainly distributed in the northwest of Beijing. What is more, the transformation of land use types in hot spot mainly occurs in the transformation from grassland to forest about 21.816 km 2 . This is also a partial confirmation of the gradual improvement of vegetation cover in China, which may play a role in relieving the two main problems in Beijing: soil erosion and air pollution.
The results of this study (Figures 4 and 5) are similar to those of Peng, J.'s study [59]: from 2001 to 2009, at first, with the expansion of urban areas, ecological land was gradually occupied which led to a decrease in NPP, and finally the NPP increases occurred in highly developed urban zones, which depended to some extent on the growth of green land, local eco-environmental management and urban planning. Therefore, NPP in Beijing showed an obvious upward trend with a series of ecological forestry projects from 2010 to 2020. As Beijing has been highly developed urban zones, the future land use will not change greatly in the next 10 years and the promotion of NPP will exist which is consistent with the predicted results ( Figure 6b). In addition, continuous decline areas of the future development of NPP belong to the development area in this study and continuous rise areas account for 84.76%.
However, unlike Beijing, land cover and land use changed in some areas have reduced vegetation coverage and net ecosystem production [60,61]. The reason may be the local inaction to the forest which resulting in the decline of its area under the action of human activities. Therefore, the simulation of future development can provide a new basis for local planning, especially in similar countries and regions in the process of urbanization.
In 2019, Nature [62] released research reports pointing out that China alone accounts for 25% of the global net increase in leaf area with only 6.6% of global vegetated area and China is engineering programs to conserve and expand forests with the goal of mitigating land degradation, air pollution and climate change. Although the area of cold spot of NPP changes in Beijing is large (Figure 7), it has more forest cover and the government attaches importance to the protection of ecological environment. The overall NPP has maintained a stable rising level and will continue to rise ( Figure 6).In the process of urbanization, the protection of forests and other natural resources in Beijing made the proportion of vegetation coverage increase, and the vegetation NPP also showed an upward trend, which has a certain reference significance for the future urban planning of developing countries. For example, in the development process of some cities, the proportion of vegetation and net ecosystem production decreased [60,61]. Beijing, which has a relatively rapid urbanization process, can provide important information for land use planning and sustainable urban development in similar urban expansion areas, especially in less developed countries. Enoguanbhor, E.C. and others observed mismatches between past/current land cover and the existing land use plan [61]. The mismatch was also between Beijing's regional land use plan and unregulated urban expansion ("The Master Plan of Beijing (2004-2020)" requires that the scale of construction land in 2020 be controlled at 1650 square kilometers, but the actual area is 3475.76 square kilometers). Therefore, it is more necessary to carry out multi scenario future land use simulation to provide a scientific basis for urban planning. This paper can provide a way for cities in China and around the world to carry out spatial planning for sustainable development.

Limitations and Future Work
At present, some achievements have been made in the research of NPP estimation model, but there are a few studies on NPP prediction. Lixia, W. [22] coupled CASA and CA-Markov model, finding that comparing the estimate and predict values of NPP in 2015, the kappa coefficient reaches 0.8776 which indicates that the coupled model has high accuracy and good applicability for NPP prediction. This study also attempted to use the coupled model (CASA and PLUS model) to predict the values of NPP. NPP was divided into 5 types, 7 types and 10 types as land use types, and the kappa coefficients are 0.584571, 0.524487 and 0.508822, respectively, which only achieved the level of moderate consistency. In the next step, we will try to improve the accuracy of the coupling model or study a better model to predict the value of NPP. Moreover, the influencing factors of NPP are relatively complex, such as the combination of water and heat, the physiological and ecological characteristics of vegetation, and there are also interactions among various climate factors. Therefore, it is required in the following study to take multiple driving factors into account predicting spatiotemporal variations of NPP.
Our future work plan will also appropriately expand the spatial scale of the research. The research will not only be limited to the better developed Beijing area, but also include the areas with large differences in development, such as Hebei Province and Tianjin city, China.

Conclusions
Taking Beijing as an example, this research first used the PLUS model and Markov chain model to project the land use for 2030 under three different scenarios based on physiographic driving factors, transportation driving factors and socioeconomic driving factors. These three different scenarios are the baseline scenario, cultivated land protection scenario and ecological protection scenario. Then, the Theil-Sen Median trend analysis and Mann-Kendall test were used to quantify and analyze the spatiotemporal changes in NPP (2000-2020). The R/S analysis and Hurst index method predict the future development trend of NPP according to the past trend of long time series data. Finally study the relationship between the cold and hot spots of land use and the change of NPP and verify whether it is consistent with the simulation results. The conclusions are divided into four specific sections: (1) From 2000 to 2030, the main land use type in Beijing is forest whose proportion has been increasing. The future land use change will focus on cropland and impervious surface. Among the three scenarios, the cropland protection scenario is the only scenario in which cropland will not be transformed into impervious surface in a large area; the ecological protection scenario has good protection for forest land, wetland and water body, and the area basically does not change; in the baseline scenario, the land use change is the largest, and the area of impervious surface increases the most. (2) NPP of vegetation fluctuated and increased from 2000 to 2020 and the rising area accounted for 86.77% of the total area. The future change of the vegetation NPP showed a rise trend on the whole (rise areas and continuous rise areas account for 0.04% and 84.76%, respectively). The NPP of Beijing showed a strong improvement trend and will maintain. (3) The main land types in the hot spot areas of vegetation NPP are forest and grassland, and the main transformation of land use types is from grassland to forest. The cold spot areas are mainly concentrated in the periphery of urban areas and water areas. (4) The NPP increases occurred in highly developed urban zones with the growth of green land and the right ecological management and urban planning.