Spatio-Temporal Evolution, Prediction and Optimization of LUCC Based on CA-Markov and InVEST Models: A Case Study of Mentougou District, Beijing

With the rapid advancement of urbanization and industrialization, the contradiction between the social economy and resources and the environment has become increasingly prominent. On the basis of limited land resources, the way to promote multi-objective comprehensive development such as economic, social development and ecological and environmental protection through structure and layout regulation, so as to maximize regional comprehensive benefits, is an important task of current land spatial planning. Our aim is to obtain land-use-change data in the study area using remote-sensing data inversion and multiple-model simulation. Based on land suitability evaluation, we predict and optimize the land use structure of the study area in 2030 and evaluate and compare ecosystem services. Based on remote-sensing images and eco-environmental data from 1985 to 2014 in the study area, land use/land cover change (LUCC) and future simulation data were obtained by using supervised classification, landscape metrics and the CA-Markov model. The ecosystem services were evaluated by the InVEST model. The analytic hierarchy process (AHP) method was used to evaluate the land suitability for LUCC. Finally, the LUCC in 2030 under two different scenarios, Scenario_1 (prediction) and Scenario_2 (optimization), were evaluated, and the ecosystem service functions were compared. In the last 30 years, the landscape in the study area has gradually fragmented, and the built-up land has expanded rapidly, increased by one-third, mainly at the cost of cropland, orchards and wasteland. According to the suitability evaluation, giving priority to the land use types with higher environmental requirements will ensure the study area has a higher ecosystem service value. The rapid development of urbanization has a far-reaching impact on regional LUCC. Intensive land resources need reasonable and scientific land use planning, and land use planning should be based on the suitability evaluation of land resources, which can improve the regional ecosystem service function.


Introduction
Under the influence of natural and human factors, land use/land cover change (LUCC) has had a latent influence on the material cycle and energy conversion in the ecological chain with the continuous advancement of urbanization. These impacts alter habitat InVEST model developed by Stanford University [43]. Among them, the InVEST model has been more widely used because of its accurate quantification, visualization of results and low cost [44]. The InVEST model was mainly divided into marine, land and freshwater modules, and each module was specifically divided into several assessment projects. The assessment mainly included three aspects: The change in ES caused by the change in ecological environment dominated by the change in LUCC and landscape patterns [45]; crop yield evaluation, grassland resource evaluation and forest volume estimation aiming at evaluation and application [46]; and ES driven by ecosystem characteristic parameters, such as the net primary productivity, vegetation coverage and leaf area index [47].
The Mentougou District of Beijing is the western mountainous area of Beijing, the capital of China. As the western green barrier of the capital, it was also the only pure mountainous area in Beijing, providing important ES functions [48,49]. In the mountainous areas of the Chinese capital, land resources have economic, social, natural and ecological attributes [50,51]. This study took the Mentougou District of Beijing as the research object and took seven issues of remote-sensing images and other data (including terrain, landform, climate, soil, etc.) from 1985 to 2014 as the basic data source. Using the methods of supervised classification, landscape metrics, spatial statistical analysis, scenario simulation, grey-linear model, analytic hierarchy process and ES evaluation, the LUCC and ES function in the study area were used, evaluated and predicted. The purpose of this study was to (1) analyze the temporal and spatial transfer characteristics and landscape patterns of land use in the study area from 1985 to 2014; (2) carry out land suitability evaluation and scenario prediction and simulation of future land use; (3) evaluate the differences of ES functions in the study area in 2030 under different scenarios. The intention of the present study was to provide a scientific reference for the formulation of regional ecological environment construction and sustainable development.

Study Area
The Mentougou District is located in the northwest of Beijing, northern China (115 • 25 E-116 • 10 E, 39 • 48 N-40 • 10 N), with a total area of 1447.85 km 2 ( Figure 1). The region has a mid-latitude continental monsoon climate, and 98.5% of the region is mountainous [52]. In 2019, the average annual precipitation in Mentougou District was 405.7 mm and the average temperature was 13.8 • C [53]. Mentougou district has 4 streets and 9 towns, and in 2019, 254,000 people registered for residence and 344,000 permanent residents, of whom 42,000 belonged to the agricultural population and 211,000 to the non-agricultural population [54]. The regional GDP is 391.33 million dollars, and the primary, secondary and tertiary industries accounted for 1.24%, 26

LUCC Classification and Quantification of Landscape Patterns
In this study, we used remote-sensing image data of the Landsat Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM) and Operational Land Imager (OLI) for nearly 30 years (one image per year from 1985,1990,1995,2000,2005, 2010 and 2014) [56,57]. Based on the Chinese land use classification system and the characteristics of the

LUCC Classification and Quantification of Landscape Patterns
In this study, we used remote-sensing image data of the Landsat Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM) and Operational Land Imager (OLI) for nearly 30 years (one image per year from 1985, 1990, 1995, 2000, 2005, 2010 and 2014) [56,57]. Based on the Chinese land use classification system and the characteristics of the study area, land use is divided into nine types through supervised classification [58,59]. LUCC includes cropland (CL), orchard (OC), forestland (FL), shrubland (SL), grassland (GL), bareland (BL), waterbodies (WB), wasteland (WL) and built-up land (BUL), forming a grid image with a spatial resolution of 30 m (Table A1). Using the method of stratified random sampling, we randomly selected f 270 points (30 points for each land use type) in the study area, using higher-resolution remote-sensing images and related planning data (before 2005) and field survey (after 2005), and tested the accuracy of interpretation, and the overall accuracy of data interpretation for each year was over 85%.

LUCC Transfer Matrix
Through the analysis of the land use transfer matrix, two different periods of transformation could be obtained. In order to clearly express the data of land use transfer from 1990 to 2015, we made a graph called circos and listed the data as an Appendix ( Figure 2 and Tables A2-A8). The mapping method using the template is from Canada's Michael Smith Genome Sciences Centre (CMSGSC) (http://mkweb.bcgsc.ca/tableviewer/ accessed on 14 August 2021) [60,61].

LUCC Classification and Quantification of Landscape Patterns
In this study, we used remote-sensing image data of the Landsat Thematic M (TM), Enhanced Thematic Mapper Plus (ETM) and Operational Land Imager (O  nearly 30 years (one image per year from 1985,1990,1995,2000,2005,2010 and [56,57]. Based on the Chinese land use classification system and the characteristics study area, land use is divided into nine types through supervised classification [ LUCC includes cropland (CL), orchard (OC), forestland (FL), shrubland (SL), gra (GL), bareland (BL), waterbodies (WB), wasteland (WL) and built-up land (BUL), fo a grid image with a spatial resolution of 30 m (Table A1). Using the method of str random sampling, we randomly selected f 270 points (30 points for each land use ty the study area, using higher-resolution remote-sensing images and related plannin (before 2005) and field survey (after 2005), and tested the accuracy of interpretatio the overall accuracy of data interpretation for each year was over 85%.

LUCC Transfer Matrix
Through the analysis of the land use transfer matrix, two different periods of formation could be obtained. In order to clearly express the data of land use transfe 1990 to 2015, we made a graph called circos and listed the data as an Appendix (Fi  and Tables A2-A8). The mapping method using the template is from Canada's M Smith Genome Sciences Centre (CMSGSC) (http://mkweb.bcgsc.ca/tableviewer/ ac on 14 August 2021) [60,61].  An explanation of the circos chart. Each stripe represents the occurrence of a transformation process at two different times, and the color represents the land use types. (a) The total amount of A 1 transferred in and transferred out; (b) the proportion of each land-use type transferred into A 1 ; (c) the proportion of each land-use type in transferred from A 1 ; (d) the total area of A 1 transformed out and transformed into other types; (e) the types of land use transformed from A 1 ; (f) the area of A 1 transformed into other land-use types; (g) the types of land use transformed into A 1 ; (h) the area of A 1 transformed from other land-use types.

Statistical Analysis
To assess changes in the structural characteristics of landscape patterns at the subwatershed level, we selected the number of patches (NP), mean patch size (MPS), patch density (PD) and patch cohesion index (COHESION), to characterize landscape patterns (Table 1) [62,63]. All calculations were extracted from the FRAGSTATS 4.2 software (Developed by the Clark Labs, Clark University). It can be used to characterize landscape fragmentation.
Aggregation COHESION It reflects the aggregation degree of patches in the landscape.
NOTE: a ij is the area of patch ij, p ij is the common border length of patch ij, and n i is the number of patches in the landscape of patch type (lass) i.

AHP and CA-Markov
According to ES demand and natural climate factors, a land suitability evaluation map was constructed ( Figure 3). The analytic hierarchy process of determining the importance of each indicator was conducted by more than 20 experts in ecology-related fields. All results were tested for consistency using Yaahp software (developed by the Shanxi Yuan Decision Software Technology Co., LTD in China). Spatial simulation is based on the CA-Markov module in IDRIS 17.0 software (developed by the Clark Laboratory of Clark University in the United States) to simulate LUCC changes in the study area [64][65][66].  (Table A9). According to the results of the suitability analysis, priority should be given to the land use types requiring higher environmental conditions, and these types with higher land-use requirements should be allocated to the appropriate space, so as to obtain the optimization of the LUCC spatial distribution in 2030 (Scenario_2, regarding which landuse types are allocated to more suitable areas). (7) Lastly, the climate factors in 2014 were brought into the InVEST model to calculate and compare the ES of the predicted and optimized LUCC in 2030.
The kappa coefficient is calculated as follows [62]: where is the observed consistency rate between the reference data and the simulation data, is the expected correct simulation proportion in the case of randomness and is the correct simulation proportion in the case of ideal classification, which is generally 1.  The process was as follows ( Figure 4): (1) Determine the transition area matrix and transition probability matrix of LUCC in the study area from 1985 to 2014, which were used to test the accuracy and predict the results with the CA-Markov model. (2) Then, we established the atlas of LUCC change suitability rules; the spatial distribution probability atlas of each LUCC was obtained by analyzing the limiting factors and driving factors as independent variables. (3) Next, the CA filter was constructed; in this study, a 5 × 5 mole neighborhood [21] (four adjacent cells above and below a cell are domains) was used as the filtering parameter of the CA-Markov model. (4) We then determined the starting time and the number of iterations. Firstly, the period from 1990 to 2010 was taken as the starting time of prediction, and the number of CA iterations was set as 10 to simulate the spatial distribution of land use types in the study area in 2010. Compared with the real results, kappa coefficients of all types of land use were > 0.70%, which was feasible. (5) Taking 1990 to 2014 as the starting point of prediction, the number of iterations was set as 16 to predict the LUCC in the study area in 2030. The prediction result was obtained (Scenario_1). (6) According to the land quantity structure and land suitability evaluation results, the land use transformation was restrained and controlled. The gray-linear model was used to solve the quantity of land use with equal emphasis on ecology and economy (Table A9). According to the results of the suitability analysis, priority should be given to the land use types requiring higher environmental conditions, and these types with higher land-use requirements should be allocated to the appropriate space, so as to obtain the optimization of the LUCC spatial distribution in 2030 (Scenario_2, regarding which land-use types are allocated to more suitable areas). (7) Lastly, the climate factors in 2014 were brought into the InVEST model to calculate and compare the ES of the predicted and optimized LUCC in 2030.    The kappa coefficient is calculated as follows [62]: where P o is the observed consistency rate between the reference data and the simulation data, P c is the expected correct simulation proportion in the case of randomness and P p is the correct simulation proportion in the case of ideal classification, which is generally 1.

Water Yield
The water yield (WY) module of the InVEST model can estimate the WY function of different LUCC. The core algorithm is used to calculate the WY of each grid by using the water balance method combined with climate, terrain and LUCC types. The WY is the precipitation of each grid in the region minus the actual evapotranspiration. The calculation formula of WY is as follows [67].
Y x is the annual water yield of grid x (mm); T AEx is the actual annual evapotranspiration of grid x (mm); P x is the precipitation (mm) of grid x (mm); T AEx is approximated from the Budyko curve [67]; R x is the Budyko dryness index of grid x, which is the ratio of potential evapotranspiration (T PEx ) to precipitation (P x ); C Kx is the evapotranspiration coefficient of vegetation, which is different in different vegetation types, and represents the ratio of potential evapotranspiration to reference evapotranspiration (T E0x ) of plants at different growth stages; T AEx is evaluated directly by T E0x , and its value is determined jointly by T E0x and P x ; C Awx is the ratio of soil water availability to precipitation; C Awx is the soil available water content of grid x (mm); Z is the seasonal constant.

Soil Conservation and Soil Loss
The Sediment Delivery Ratio model of InVEST is used to calculate soil loss and soil conservation based on the soil loss equation. The original soil loss equation is added to the model to intercept the sediment and grid upstream sediment, and on this basis, the soil conservation amount is calculated [68].
where RUSLE n is the actual amount of soil erosion, USLE n is the potential soil loss under no vegetation coverage, SEDRET n is the soil conservation, R n is the rainfall erosivity, K n is the soil erodibility factor, LS n is the slope length-gradient factor, C n is the crop management factor, P n is the support practice factor, and SED n is the amount of sediment and sediment intercepted upstream.

Carbon Stocks
The Carbon stocks model of InVEST quantifies carbon storage and sequestration based on four carbon pools: Aboveground biomass, underground biomass, dead organic matter and soil organic matter [69,70]. Carbon storage in vegetation is estimated by multiplying the evegetation carbon density by the area of each LUCC based on local research results [71].
where C total is the total carbon stocks (t/hm 2 ), C above is the aboveground carbon stocks t/hm 2 ), C below is the underground carbon stocks (t/hm 2 ), C soil is the soil carbon stocks (t/hm 2 ), and C dead is the dead organic carbon stocks (t/hm 2 ).

Temporal and Spatial Transfer Characteristics of LUCC
During the study period, the total BUL increased by 21.11 km 2 (32.40%). This increased BUL was mainly the cost of CL (15.49 km 2 ), OC (2.58 km 2 ) and WL (2.42 km 2 ). During the study period, the CL decreased rapidly, with a total decrease of 66.97% (20.13 km 2 ). Before 2000, the reduced area of CL was mainly transformed into BUL and WL. After 2000, CL was mainly transformed into OC and BUL. Before 2000, the transferred CL area was 18.80 km 2 , of which 53.14% (9.99 km 2 ) was converted to BUL and 43.24% (8.15 km 2 ) was converted to WL. From 2000 to 2010, the total area of CL transferred out was 14.94 km 2 . Of the reduced CL, 32.20% (4.81 km 2 ) was transformed into BUL and 42.24% (6.31 km 2 ) into OC. In general, in the past 30 years, the urbanization development of the study area was mainly manifested in the transformation of agricultural land (CL and OC) into BUL. CL accounts for 73.38% of the total land transformed into BUL and OC accounts for 12.22%. In addition, the WL also decreased significantly by 80.87%, mainly transformed into BUL, FL, SL and OC, of which the area converted to FL accounted for 2.58%, BUL accounted for 2.36%, OC accounted for 10.85% and SL accounted for 82.03%. The transfer into SL comprised 224.94 km 2 and the transfer out SL was 74.82 km 2 . The area transferred into FL comprised 77.43 km 2 and the area of FL transferred out was 120.46 km 2 ( Figure 5).
creased BUL was mainly the cost of CL (15.49 km ), OC (2.58 km ) and WL (2.42 km ). During the study period, the CL decreased rapidly, with a total decrease of 66.97% (20.13 km 2 ). Before 2000, the reduced area of CL was mainly transformed into BUL and WL. After 2000, CL was mainly transformed into OC and BUL. Before 2000, the transferred CL area was 18.80 km 2 , of which 53.14% (9.99 km 2 ) was converted to BUL and 43.24% (8.15 km 2 ) was converted to WL. From 2000 to 2010, the total area of CL transferred out was 14.94 km 2 . Of the reduced CL, 32.20% (4.81 km 2 ) was transformed into BUL and 42.24% (6.31 km 2 ) into OC. In general, in the past 30 years, the urbanization development of the study area was mainly manifested in the transformation of agricultural land (CL and OC) into BUL. CL accounts for 73.38% of the total land transformed into BUL and OC accounts for 12.22%. In addition, the WL also decreased significantly by 80.87%, mainly transformed into BUL, FL, SL and OC, of which the area converted to FL accounted for 2.58%, BUL accounted for 2.36%, OC accounted for 10.85% and SL accounted for 82.03%. The transfer into SL comprised 224.94 km 2 and the transfer out SL was 74.82 km 2 . The area transferred into FL comprised 77.43 km 2 and the area of FL transferred out was 120.46 km 2 ( Figure 5).

Evolution Characteristics of Landscape Patterns of LUCC
During the study period, BUL became more fragmented, and landscape heterogeneity increased. The landscape heterogeneity of OC first increased and then decreased. Before 2000, the landscape patterns of OC showed a trend of fragmentation, and after 2000, landscape connectivity increased. CL and FL became fragmented, reducing the connectivity. There were no obvious changes in WB, GL and BL during the study period. The PD of BUL increased from 1.35 patches/km 2 to 1.59 patches/km 2 . The NP and MPS of CL decreased from 1985 to 2014, reaching 340 patches (8.84 × 10 −2 km 2 ) and 224 patches (4.43 × 10 −2 km 2 ), respectively ( Figure 3). The NP of WL increased first and then decreased, from

Evolution Characteristics of Landscape Patterns of LUCC
During the study period, BUL became more fragmented, and landscape heterogeneity increased. The landscape heterogeneity of OC first increased and then decreased. Before 2000, the landscape patterns of OC showed a trend of fragmentation, and after 2000, landscape connectivity increased. CL and FL became fragmented, reducing the connectivity. There were no obvious changes in WB, GL and BL during the study period. The PD of BUL increased from 1.

ES of Mentougou District in 2014
The total amounts of SC, SLO, WY and CS in 2014 were 0.50 × 10 8 t, 0.66 × 10 6 t, 2.02 × 10 8 m 3 and 1713.54 × 10 4 t, respectively. In terms of spatial characteristics, CS, SC and SLO were mainly distributed in some hilly areas in Mentougou District, followed by hilly areas with relatively fragmented LUCC along the river. This part of LUCC was mainly composed of CL, OC and WL. The smallest area was mainly concentrated in the southeast Mentougou District plain area and part of the BL. In general, the CS and WY of different LUCC vary greatly. FL and SL had the strongest ES capacity, with a total CS of 1.73 × 10 8 m 3 (85%), SC of 0.43 × 10 8 t (80%), SLO of 0.46 × 10 6 t (70%), and WY of 0.16 × 10 8 t (96.48%) (Figure 7).
Mentougou district has a total area of 107.59 km 2 with a slope between 0° and 5°, accounting for 7.40% of the whole region. An area of 386.62 km 2 has a slope between 5° and 15°. An area of 475.07 km 2 has a slope between 15° and 25°, accounting for 32.65%. The area with a slope of >= 25° was 485.60 km 2 , accounting for 33.38% of the county area. The sunny slope accounts for 25.01% of the total area of the district, and the shady slope for 38.62%. Soil pH was nearly neutral with values between 6.5 and 7.5, covering a total of 1027.39 km 2 , accounting for 70.62% of the whole region. The acidic soil (pH 5.5~6.5) covers an area of 395.5 km 2 , covering 27.18%. The area of alkaline soil (pH 7.5~8.5) was 31.99 km 2 , covering 2.20%.
The area with SLO less than or equal to 10 t/hm 2 had a total of 1405.77 km 2 (96.63%). The area with SC less than or equal to 200 t/hm 2 was 140.06 km 2 , accounting for 9.63% of the total area of the whole region. The area of SC between 200 t/hm 2 and 600 t/hm 2 was 1281.19 km 2 , covering 88.06%. The area with WY less than or equal to 500 m 3 /hm 2 , 500~1000 m 3 /hm 2 , 1000~1500 m 3 /hm 2 and greater than or equal to 1500 m 3 /hm 2 were 84.63 km 2 (5.82%), 361.03 km 2 (24.82%), 375.98 km 2 (25.84%) and 633.24 km 2 (43.53%), respectively. The evaluation was rated on a scale of 0 to 2, 2 to 3 and 3 to 4, which were suitable, relatively suitable and unsuitable, respectively (Table A7 and Figure A1).

Spatial Distribution Characteristics of Suitable LUCC
The results of the land suitability grade evaluation showed that the area suitable for CL was 514.36 km 2 (35.35%), mainly distributed in the gully and hilly regions ( Figure A1 and Table A10). The area unsuitable for CL was 940.52 km 2 , accounting for 64.65% of the total area of the whole region ( Figure 9 and Table 2). The area suitable for OC was 1130.64 km 2 (77.71%), mainly distributed in the eastern plain and part of the central hilly region, Mentougou district has a total area of 107.59 km 2 with a slope between 0 • and 5 • , accounting for 7.40% of the whole region. An area of 386.62 km 2 has a slope between 5 • and 15 • . An area of 475.07 km 2 has a slope between 15 • and 25 • , accounting for 32.65%. The area with a slope of >= 25 • was 485.60 km 2 , accounting for 33.38% of the county area. The sunny slope accounts for 25.01% of the total area of the district, and the shady slope for 38.62%. Soil pH was nearly neutral with values between 6.5 and 7.5, covering a total of 1027.39 km 2 , accounting for 70.62% of the whole region. The acidic soil (pH 5.5~6.5) covers an area of 395.5 km 2 , covering 27.18%. The area of alkaline soil (pH 7.5~8.5) was 31.99 km 2 , covering 2.20%.
The area with SLO less than or equal to 10 t/hm 2 had a total of 1405.77 km 2 (96.63%). The area with SC less than or equal to 200 t/hm 2 was 140.06 km 2 , accounting for 9.63% of the total area of the whole region. The area of SC between 200 t/hm 2 and 600 t/hm 2 was 1281.19 km 2 , covering 88.06%. The area with WY less than or equal to 500 m 3 /hm 2 , 500~1000 m 3 /hm 2 , 1000~1500 m 3 /hm 2 and greater than or equal to 1500 m 3 /hm 2 were 84.63 km 2 (5.82%), 361.03 km 2 (24.82%), 375.98 km 2 (25.84%) and 633.24 km 2 (43.53%), respectively. The evaluation was rated on a scale of 0 to 2, 2 to 3 and 3 to 4, which were suitable, relatively suitable and unsuitable, respectively (Table A7 and Figure A1).

Spatial Distribution Characteristics of Suitable LUCC
The results of the land suitability grade evaluation showed that the area suitable for CL was 514.36 km 2 (35.35%), mainly distributed in the gully and hilly regions ( Figure A1 and Table A10). The area unsuitable for CL was 940.52 km 2 , accounting for 64.65% of the total area of the whole region ( Figure 9 and Table 2). The area suitable for OC was 1130.64 km 2 (77.71%), mainly distributed in the eastern plain and part of the central hilly region, mainly in most plain and hilly areas, and some mountainous areas. The area of suitable FL and relatively suitable FL was 1350.94 km 2 (92.86%). The area of suitable GL and relatively suitable GL was 1350.94 km 2 (92.85%). The area unsuitable for SL was 103.94 km 2 (7.14%).

. LUCC Prediction and Optimization Results in 2030
The comparison between Scenario_1 in 2030 and the situation in 2014 showed that BUL will mainly expand into adjacent areas in the future. The difference was that the patch size and expansion rate of each BUL were different. The spread degree of BUL in plain areas was higher than that in the hilly areas, and the spread degree of BUL was the slowest in mountainous areas ( Figure 10). Compared with the LUCC situation in 2014, all the WL in the Scenario_2 (Optimization of 2030) was changed into other LUCC with higher comprehensive benefits, and the WL in the prediction results increased by 1.22 km 2 (5.07%). The area of FL and WB increased by 22.91 km 2 (4.12%) and 5.88 km 2 (35.74%), respectively, while the area of CL and SL had little change. The FL, OC and WB decreased by 53.47 km 2 (9.61%), 11.78 km 2 (22.14%) and 0.28 km 2 (1.70%), respectively, and the BUL increased in both scenarios. In terms of spatial distribution, Cl and OC migrate from an unsuitable area to a suitable area ( Figure 10 and Table 3).

Comparison of the ES of Present, Prediction and Optimization
The predicted and optimized land use map of 2030 was used to predict the model, together with the meteorological data and soil data of 2014 (precipitation, potential evapotranspiration, water available to vegetation, soil depth and rainfall erosivity, etc.). The predicted and optimized land use ES was calculated and compared with the ES in 2014.

Comparison of the ES of Present, Prediction and Optimization
The predicted and optimized land use map of 2030 was used to predict the model, together with the meteorological data and soil data of 2014 (precipitation, potential evap-otranspiration, water available to vegetation, soil depth and rainfall erosivity, etc.). The predicted and optimized land use ES was calculated and compared with the ES in 2014.

Impacts of Human Activities and Policies on LUCC
The landscape patterns of Mentougou District have undergone great changes. Since the promulgation of the Land Management Law in 1986, the market economy system (an economy in which social resources were allocated through the market) have gradually occupied a dominant position [72]. Mentougou District lost more than half of its CL in the nearly 30 years from 1985 to 2014. People tend to choose more fertile land for farming and thus the productivity of the land decreased distinctly during this period [73][74][75]. Meanwhile, BUL increased steadily during the study period. Such changes are common in rapidly urbanizing suburbs, and the Los Angeles and California suburbs of Orange County

Impacts of Human Activities and Policies on LUCC
The landscape patterns of Mentougou District have undergone great changes. Since the promulgation of the Land Management Law in 1986, the market economy system (an economy in which social resources were allocated through the market) have gradually occupied a dominant position [72]. Mentougou District lost more than half of its CL in the nearly 30 years from 1985 to 2014. People tend to choose more fertile land for farming and thus the productivity of the land decreased distinctly during this period [73][74][75]. Meanwhile, BUL increased steadily during the study period. Such changes are common in rapidly urbanizing suburbs, and the Los Angeles and California suburbs of Orange County and San Bernardino County have witnessed a dramatic loss of prime farmland and orange groves to suburban development [76]. The study area implemented the policy "Grain to Green" in 1999, prompting further losses of cropland [49]. Future economic development should take into account the unique natural conditions of mountainous areas and their landscape patterns. The government should adopt development and land management policies so that new urban development can be sustainable.

The Impact of Landscape Patterns Changes
The WL in Mentougou District was gradually fragmented, especially after 2005, and landscape heterogeneity increased rapidly. The BUL, CL and BL showed a trend of fragmentation in this period, and the degree of landscape aggregation decreased. The plain sub-region of Mentougou District accounted for less than 10% of the region's area but contained more than 60% of the population [50]. People preferred to build houses in the plains, near the existing residential development [17]. Population pressure in the plains decreased the area of cropland as it was converted to built-up land. Over time, small patches of FL increased, and SL formed large patches, which became the main LUCC in the study area. The landscape heterogeneity of OC and WB increased first and then decreased. Before 2000, the landscape patterns of OC were fragmented, and after 2000, landscape connectivity increased. The habitat was divided into fragmentary patches, which reduced ecosystem stability and accelerated the invasion of alien species due to the edge effect. Since 2005, invasive alien species have been found in Mentougou District, including Solanum rostratum Dunal. and Pueraria phaseoloides (Roxb.) Benth. [77]. These alien species, with their strong ability to survive, will take over the resources of other plants and animals and cause the extinction of local species [78,79].

Spatial Suitability of Various Types of LUCC in the Region
Through the classification and evaluation of Mentougou's basic eco-environmental factors (including slope, aspect, elevation, soil type, soil pH and ES), the suitable spatial distribution positions and quantities of various LUCC in the study area were selected. According to the Chinese policy of Grain to Green, areas with a slope greater than 25 • were strictly not allowed to be CL, so the suitable area for CL accounted for only one-third of the total area [80,81]. The suitable area of CL was mainly distributed in the eastern plain area and some parts of the central hilly area. These places were usually rich in water and experienced less soil loss [82]. The CL-suitable area was mostly brown soil and tidal soil, with a small slope, complete water conservancy facilities, convenient irrigation and high management level. It was a basic farmland area with high and stable output in Mentougou District. The unsuitable area of OC accounted for 22.29%. Cinnamon accounts for more than two-thirds of the study area, and most areas were suitable for FL, GL and SL. After returning CL to FL, the FL in the study area would further increase. The FL could improve the regional microclimate, reduce soil erosion and water evaporation, prevent wind and consolidate soil and improve soil physical and chemical properties. However, our study did not compare different suitability evaluation methods. Future studies should add other evaluation methods, such as machine learning and sensitivity analysis, to increase the objectivity of suitability evaluation [83,84].

Impacts of Future LUCC on ES
Based on the existing spatial transfer matrix of LUCC, the CA-Markov model was used to predict LUCC in 2030, the control constraints were set and the Scenario_2 spatial distribution map of LUCC in 2030 was obtained. Compared with the LUCC in 2014, Sce-nario_2 was more in line with the protection needs of the local government for CL. From the goal of common development of ecological benefits (ecological space) and economic benefits (built-up land space), Science_2 is more suitable than Science_1. Linking the InVEST and CA-Markov models has promising applications for guiding ecological engineering.
First, the CA-Markov model was applied to simulate LUCC by combing eco-environmetat, which provides a method for simulating LUCC in other areas [17,26]. More importantly, the predicted land cover map aided in the identification of areas where ecological space is more likely to occur, which has great significance for site selection for investment in ecological engineering. Second, the InVEST model can quantify a number of ecosystem services [44,46].
Against the background of the south-to-north water transfer and the government's encouragement of ecological restoration of mountains, rivers, forests, farmland and lakes, Scenario 2 would provide important guidance for planning. In addition, the OC of Sce-nario_1 was reduced by 11.78 km 2 (22.14%) and the OC of Scenario_2 remained unchanged. In Scenario_1, SL increased, while in Scenario_2, SL remained unchanged. Because the economic and ecological benefits of the GL in the study area were not high, the area was reduced in both Scenario_1 and Scenario_2 (Table 3). In terms of spatial distribution, Sce-nario_2 was more conducive to the future development of the research area. ES (WY, CS, SC and SL) was evaluated using environmental and climatic factors in 2014, and each index proved that Scenario_2 had higher ES capacity.

Conclusions
From the perspective of landscape patterns, during the past 30 years (1985-2014), LUCC in Mentougou District had undergone drastic changes, with rapid urbanization. BUL increased by one-third, CL decreased by more than half, and WL decreased by more than two-thirds. At the same time, the landscape patterns of various LUCC also changed, generally showing a trend of gradual fragmentation, decreasing landscape connectivity and increasing heterogeneity. In terms of ES and spatial suitability, the total CS, SC, SLO and WY in the study area in 2014 were 1713.54 × 10 4 t, 0.50 × 10 8 t, 0.66 × 10 6 t and 2.02 × 10 8 m 3 . The area suitable for CL and OC in the study area accounted for one-third and two-thirds of the whole study area, and the area suitable for FL, GL and SL accounted for more than 90%. In 2030, the areas of FL, WB, OC and CL in the study area would decrease in Scenario_1, and CL would especially decrease by about half. In Scenario_2, FL and WB would increase by about 40% in total, and BUL also would develop well. In scenario_2, ES would be also better than scenario_1, WY and SC would increase by more than one-third and SLO would decrease. According to the land-use configuration of Scenario_2, the vegetation would grow under more suitable environmental conditions, and the spatial distribution of land use would be more reasonable. This study is helpful for policymakers, planners and landscape designers to determine urban land use schemes and promote the rational allocation of urban land.     The data in each row add up to the total area of land use in the same category in 1995, and the data in each column equal the total area of land use in the same category in 2000. The data in each cell represent the area of the row land-use type transferred to column land-use type from 1995 to 2000. The total area of the study area is 64,500 km 2 .   The data in each row add up to the total area of land use in the same category in 2010, and the data in each column equal the total area of land use in the same category in 2015. The data in each cell represent the area of the row land-use type transferred to column land-use type from 2010 to 2015. The total area of the study area is 64,500 km 2 .    Figure A1. Map of land use type suitability evaluation in Mentougou District; (a-e) the spatial distribution of suitability for cropland, orchard, forestland, grassland and shrubland; (f) areas that are not included in model prediction. Figure A1. Map of land use type suitability evaluation in Mentougou District; (a-e) the spatial distribution of suitability for cropland, orchard, forestland, grassland and shrubland; (f) areas that are not included in model prediction.