Ecological Security Pattern Construction in Beijing-Tianjin-Hebei Region Based on Hotspots of Multiple Ecosystem Services

The contradiction between urban expansion and ecological protection in the Beijing-TianjinHebei region (BTH) is increasingly acute, which has become one of the main problems restricting regional development, and sustainable development of ecosystem services is the key to increasing human well-being. Based on GIS platform and multiple models, this paper analyzes the temporal and spatial variation characteristics of four key ecosystem services (water conservation, soil conservation, habitat quality, and plant net primary productivity) in different ecological regions of BTH in recent 20 years, quantifies the impact of different climate factors and land use change on ecosystem services (ESs), and discusses the primary ecosystem hotspots and ecological security pattern. The results showed that the interannual variation of water conservation (WC) and plant net primary productivity (NPP) increased from 2000 to 2020, while the change of soil conservation (SC) was not obvious, which was mainly controlled by climate factors, WC and SC were more affected by precipitation, and temperature was the key factor affecting NPP. Habitat quality (HQ) presented a significant downward trend; it was mainly attributed to the deterioration of ecological environment caused by accelerated urbanization expansion. According to hotspot analysis, it could be found that WC was the fastest-growing ecosystem service function in BTH, and NPP would become the factor with the greatest contribution to ecological importance in the future. The important protected areas and main ecological sources of ecological security pattern were mainly distributed in Yanshan-Taihang mountain area, which was consistent with the key areas of ecosystem services. In this study, the temporal and spatial differences of ecosystem service in BTH were demonstrated in a more intuitive way and provided scientific guidance for decision makers to formulate effective ecological protection policies in different regions.


Introduction
Ecosystem Services (ES) can be defined as the result of positive interactions between humans and the ecosystem environment [1]. Ecosystem services, which can be divided into regulation services, supply services, support services, and cultural and entertainment services, play a fundamental role in helping human getting various benefits with a direct or indirect process [1][2][3]. However, the sustainable supply of ESs has been broken because as people gain the goods and benefits from ecosystem without restraint, ESs face more and more pressures due to the multiple changes of climate and land use, which are exacerbated by socio-economic developments and cause the destruction of ecosystem services such as soil and water loss, habitat conversion, degradation, and fragmentation [4][5][6]. In recent decades, about 60% of ecosystem services around the world have been degraded or unsustainable [2] and global ecosystems are showing declining trends at an unprecedented speed [7]. Scholars have started paying attention to the negative change of ESs model has become the mainstream model to build the ecological security pattern, there is no standardized mode to choose ecological source and to build an ecological security pattern related to the protection of human well-being by combining ecological services [35].
Beijing-Tianjin-Hebei region (BTH) is not only the political, economic, and cultural center in China but also a particularly important ecological area. In the past 50 years, temperature in BTH has increased significantly, while precipitation and potential evapotranspiration have decreased in most sub-regions [36]. Simultaneously, the expansion speed of construction land is accelerated from 6.75% in 1975 to 12.9% in 2020 due to rapid urbanization. In order to alleviate the degradation of ESs due to urbanization, the ecological projects such as Grain to Green Project, Natural Forest Conservation Program, and Beijing-Tianjin Sand Source Control Program have been launched successively. Therefore, reasonable analysis of the impact of climate and land use change on ecosystem services, ecological security pattern construction in BTH, and scientific guidance for decision makers are of great significance to regional development.
The object of this study is to build a comprehensive evaluation to analyze the ecosystem service function of BTH from two aspects: driving force analysis and ecological security pattern construction. The process of model simulation is shown in Section 2. Section 3 mainly analyzes the temporal variation and spatial distribution of four kinds of ecosystem services (WC, SC, HQ, NPP); by different scenarios, the influence degree of driving factors was analyzed, and the ecological security pattern was constructed based on hotspot analysis. Section 4 mainly discusses the improvement of research method compared with the previous research methods, puts forward reference suggestions for the ecological civilization construction, and illustrates future research prospect. The Supplementary

Study Area
The BTH is located in North China between 36.07 • N-42.65 • N and 113.46 • E-119.79 • E, and the elevation is from 52 m to 2836 m ( Figure 1). The elevation is descending from north to south, which is characterized by the transition from mountain and hill to platform and plain. The climate type is humid and semi-arid continental monsoon climate with four distinct seasons [37]. Affected by elevation and climate factors, vegetation types are diverse, which can be divided into four vegetation zones: Bashang grassland, mountain and hilly deciduous broad-leaved forest, plain deciduous broad-leaved forest, and coastal plain halophyte [38]. With the continuous improvement of urbanization level, the proportion of construction land has risen from 8.17% in 2000 to 12.99% in 2020. About 10.2% of ecological land is converted to construction land, and the occupation of ecological land by urban expansion has become an important issue affecting the construction of ecological civilization.
According to the ecological regionalization by the Research Center for Eco-Environmental Sciences, Chinese Academy of Sciences (http://www.ecosystem.csdb.cn accessed on 20 December 2020), combined with topography, hydrothermal condition, vegetation characteristics, and ecosystem structure of BTH, the study area is divided into seven ecological regions ( Figure 1). According to the ecological regionalization by the Research Center for Eco-Environmental Sciences, Chinese Academy of Sciences (http://www.ecosystem.csdb.cn accessed on 20 December 2020), combined with topography, hydrothermal condition, vegetation characteristics, and ecosystem structure of BTH, the study area is divided into seven ecological regions (Figure 1).

Data Sources
The data required in this study mainly include meteorological, vegetation, land use, soil, DEM, and socio-economic data. The meteorological data was downloaded from China Meteorological Data Service Centre, including daily temperature, precipitation, and potential evapotranspiration data of 123 meteorological stations in BTH and its surrounding areas, and the data was rasterized by Kriging. The remote sensing data of NDVI, Land use, and DEM are raster data, soil and road data are vector data, and Water Resources is textual data from statistical yearbook. All data were consolidated into WGS84 coordinate system and Albert projection. The details are shown in Table 1.

Data Sources
The data required in this study mainly include meteorological, vegetation, land use, soil, DEM, and socio-economic data. The meteorological data was downloaded from China Meteorological Data Service Centre, including daily temperature, precipitation, and potential evapotranspiration data of 123 meteorological stations in BTH and its surrounding areas, and the data was rasterized by Kriging. The remote sensing data of NDVI, Land use, and DEM are raster data, soil and road data are vector data, and Water Resources is textual data from statistical yearbook. All data were consolidated into WGS84 coordinate system and Albert projection. The details are shown in Table 1.

Ecosystem Services Indicators
As an integral part of nature, ecosystem services can be used as an intermediate factor to connect ecosystem with human welfare, which has direct value to human and society, so as to bridge the differences between ecosystem and economic issues. Referring to related research [39], this study comprehensively considers international classification of ecosystem service [40], key detection indicators monitored by local governments, the difficulty of indicators evaluation and the close relationship with human beings [41], and habitat quality (HQ), soil conservation (SC), water conservation (WC), and net plant productivity (NPP) are selected as the key evaluation indicators. The InVEST model and CASA model have been widely used in ecosystem service assessment due to their advantages of simple parameter setting and good applicability [42][43][44]. In the current study, HQ, SC, and WC were quantified by InVEST model, and NPP was simulated by CASA model, the computational methods of parameters are shown in Table 2, and the specific calculation formula of the model is detailed in Supplementary Information.

CA_Markov Model
CA_Markov model is a model which calculates the relationship between initial state and transition probability of different states based on Markov chain [52]. The principle of the model is to predict future state changes by setting rules in the current state of the cell and its surrounding state [53]. This model is divided into CA model and Markov process. The model coupling process is as follows: 1.
The Markov model is used to calculate the land use type under the "possible state", and the area quantity or proportion of mutual transformation of land use types is regarded as a state transition probability. This probability is used to predict the change state of land use structure. The specific formula is as follows: where S(T) is the state of land use structure at time T, and P ij is the state transition matrix.

2.
A 5 × 5 CA filter (is) is constructed to analyze the influence of the rectangular space composed of 5 × 5 cells around the cell on the state change of the cell, and a weight factor with significant spatial significance is established [54], which is expressed as follows: where S is the finite and discrete set of cellular states, N is the domain of cellular, and f is the transformation rule of cellular states in local space.

3.
The starting time of prediction period and the number of CA cycles were determined. Based on the 2000-2015 transfer matrix of BTH region, the number of CA iterations was set as 15 to simulate the land use pattern of BTH region in 2030. Kappa index was used to test the accuracy of the initial transfer pattern changes [53]. A comparison was carried out between the predicted and actual land use in 2015. The Kappa coefficient value is 0.91, which represents a good simulation result. Thus, the verified CA_Markov can be used to predict land use type in 2030.
Three scenarios were set to predict land use change of BTH region in 2030 by CA_Markov model: business as usual (BAU), economic growth (ED), and ecological conservation (EC). For the BAU scenario, we assumed that the transition probability of shifting the Markov chain model from 2015 to 2030 is similar to that from 2000 to 2015: the change patterns of LUCC are sustained. In EC scenario, the area of forestland, grassland, construction land, and wetland transferred into construction land was decreased by 90%, 90%, 50%, and 20% respectively, and the transformation of unused land into construction land was prohibited. Under the ED scenario, the LUCC was based on the development in BTH, the area of construction land, forestland, grassland, and unused land converted to construction land was increased by 3%, and the transformation of water areas into construction land was forbidden. The land requirement of three different scenarios in 2030 is summarized in Table 3.

Hotspot Analysis
Hotspot analysis is a proven technique to identify the region where the value that is of high interest occurs [55,56]. In this study, the grid whose ecosystem service exceeded the respective average values of all grids was defined as the hotspots of this kind of ecosystem service, and the multiple ecosystem services hotspots in BTH were obtained by overlaying the hotspots of four ecosystem services. If the grids did not exceed their corresponding mean values, these raster cells were defined as non-hot area. If the raster cell had only one ecosystem service above the average, then this type of raster was defined as class 1 hotspot. Based on this, Class 2, 3, and 4 hotspots were defined respectively, and the spatial distribution pattern of multiple ecosystem service hot spots was finally obtained [57,58].

Minimum Cumulative Resistance
MCR model considered the source, distance, and landscape interface characteristics, and calculated the total ecological flow when passing through different landscape units, where the ecological flow represents the resistance value of movement on two adjacent ecological sources. By calculating the cumulative cost of ecological flow through the base surface, the minimum cost distance, namely the path of minimum cumulative resistance, was found. The resistance surface reflected the potential possibility and trend of species movement. According to the MCR resistance value, ecological expansion and urban expansion were mutually restricted. To coordinate the process of ecological expansion and urban expansion [59,60], minimum cumulative resistance model is as follows [61]: where F MCR is the minimum cumulative resistance value, f is the positive correlation function between the minimum cumulative resistance and ecological process, and min represents the minimum cumulative resistance value of unit i for different sources; D ij represents the spatial distance from source j to landscape unit i, and R i represents the resistance coefficient of landscape unit i to species movement. Construction of resistance surface: the resistance surface can reflect the resistance of species moving between different landscape units and is the basis of regional ecological security pattern [62]. Based on the current situation of ecological environment in BTH, the resistance surface was constructed from the perspective of nature, economy, and society (Table 4). By analytic hierarchy process (AHP) [27,28], the spatial weighting calculation is carried out in GIS platform to obtain the resistance surface in the expansion process of ecological sources in BTH.       , SC (c), and HQ (d), α represents the change rate of four ecosystem services, which is classified based on the natural partition point method. Red represents the very low change rate, yellow represents the low change rate, green represents the medium change rate, blue represents the fast change rate, and dark blue represents the very fast change rate.

Spatio-Temporal Characteristic of Ecosystem Services in BTH
NPP showed an increasing trend, with the rate of 5.28 gC·m −2 ·year −1 . The lowest NPP was 249.79 gC·m −2 ·year −1 in 2000 and the highest was 249.79 gC·m −2 ·year −1 in 2020 ( Figure  2b). In the last 20 years, NPP in Regions B, D, and F of Yanshan-Taihang Mountain and North China Plain was 308.72~428.66 gC·m −2 ·year −1 ) (Figure 3b), which was much higher than that in Bashang Plateau of Region A (254.01 gC·m −2 ·year −1 ). In addition, the increasing rate of NPP in Region C and D was up to 6.44 gC·m −2 ·year −1 and 8.11 gC·m −2 ·year −1 , respectively, which was much higher than other regions. This was mainly affected by ecological projects, such as the policy of "Three-North Shelter Forest Program" and the program of Grain for Green, which made the productivity of plant communities in mountainous regions raise continuously and the high value area of NPP increase significantly.
From 2000 to 2020, SC decreased from 162.25 t·hm −2 year −1 to 128.8 t·hm −2 year −1 and reduced by 20.62% at a rate of −0.893 t·hm −2 year −1 (Figure 2c). The SC variation was different among all ecological regions (Figure 3c). Except for the slight increase rate of 0.02 t·hm −2 year −1 in Region G, the other six ecological regions showed a decreasing trend. Among them, Northern of Hebei and Yanshan Mountains (Region B) and Bashang Plateau (Region A) demonstrated the most significant decreasing trend, with a rate of 2.09 and 0.91 t·hm −2 year −1 , respectively, which was related to continuous expansion of regional farmland and continuous reduction of forestland with strong root interception ability. In , SC (c), and HQ (d)), α represents the change rate of four ecosystem services, which is classified based on the natural partition point method. Red represents the very low change rate, yellow represents the low change rate, green represents the medium change rate, blue represents the fast change rate, and dark blue represents the very fast change rate.
NPP showed an increasing trend, with the rate of 5.28 gC·m −2 ·year −1 . The lowest NPP was 249.79 gC·m −2 ·year −1 in 2000 and the highest was 249.79 gC·m −2 ·year −1 in 2020 (Figure 2b). In the last 20 years, NPP in Regions B, D, and F of Yanshan-Taihang Mountain and North China Plain was 308.72~428.66 gC·m −2 ·year −1 ) (Figure 3b), which was much higher than that in Bashang Plateau of Region A (254.01 gC·m −2 ·year −1 ). In addition, the increasing rate of NPP in Region C and D was up to 6.44 gC·m −2 ·year −1 and 8.11 gC·m −2 ·year −1 , respectively, which was much higher than other regions. This was mainly affected by ecological projects, such as the policy of "Three-North Shelter Forest Program" and the program of Grain for Green, which made the productivity of plant communities in mountainous regions raise continuously and the high value area of NPP increase significantly.
From 2000 to 2020, SC decreased from 162.25 t·hm −2 year −1 to 128.8 t·hm −2 year −1 and reduced by 20.62% at a rate of −0.893 t·hm −2 year −1 (Figure 2c). The SC variation was different among all ecological regions (Figure 3c). Except for the slight increase rate of 0.02 t·hm −2 year −1 in Region G, the other six ecological regions showed a decreasing trend. Among them, Northern of Hebei and Yanshan Mountains (Region B) and Bashang Plateau (Region A) demonstrated the most significant decreasing trend, with a rate of 2.09 and 0.91 t·hm −2 year −1 , respectively, which was related to continuous expansion of regional farmland and continuous reduction of forestland with strong root interception ability. In addition, construction land in these areas also expanded, resulting in obvious attenuation of regional soil conservation capacity.
The HQ reduced from 0.631 to 0.577 during 2000-2020, with a decrease rate of 8.32% (Figure 2d). Spatially, the decreasing rate of HQ in North China Plain and coastal regions (Regions E, F, and G), greatly affected by human activities, showed the greatest decrease. The expansion of cities such as Beijing and Tianjin and the construction of new areas such as Xiong'an occupied a large amount of surrounding ecological land, which damaged the natural ecosystem and led to a significant deterioration in habitat quality.
3.2. Analysis on the Driving Factor of Ecosystem Service Variation 3.2.1. Impacts of Climate Change on Ecosystem Services Table 5 showed the variation trends of temperature, precipitation, and potential evaporation in different ecological regions of BTH from 2000 to 2020. The regional annual average precipitation and temperature presented an increasing trend, with growth rates of 4.06 mm/year and 0.003 • C/year respectively; the annual average potential evaporation showed a decreasing trend with the linear change rate of −3.56 mm/year. On the whole, the climate in BTH tended to be warm and humid, but there were spatial differences in climate change among each ecological region. The central and southern parts of BTH showed a warming trend, while the northern part showed a weak cooling trend; the increasing trend of precipitation in plain area is significantly higher than that in mountainous area. The temperature rise and decreasing trend of potential evaporation in Region D was the most obvious, while the increase of precipitation was relatively small, which was caused by the strengthened fohn effect by the proximity of Taihang Mountain to the west. The obvious trend of warming and humidification in Region E lied in the accelerated expansion of urbanization in the past 20 years. The proportion of construction land in 2020 had reached 30%, the urban impervious surface area had surged, and urban heat island and wet island effects had been significantly strengthened. To diminish the impact of climatic anomaly in special years, the average climate state during 2010-2020 was set as the baseline. One must quantitatively evaluate the responses of WC, SC, and NPP to variation of 30% and 10% in precipitation and potential evaporation (Figure 4). Under the condition of constant potential evaporation, WC would increase 97.16% or decrease 66.73% with precipitation increase (decrease) 30%, and WC would increase 49.48% or decrease 28.02% with precipitation increase (decrease) 10%. Under the condition of constant precipitation, 30% increase (decrease) in potential evaporation would result in 27.4% decrease or 43.9% increase in WC, and 10% increase (decrease) in potential evaporation would lead to 10.5% decrease or 12.2% increase in WC. It can be seen that precipitation has a greater impact on WC than potential evaporation, and the change amplitude of WC caused by precipitation increase (potential evaporation decrease) was greater than that caused by precipitation decrease (potential evaporation increase), and the drastic change of precipitation presented greater impact on change amplitude of WC. Precipitation was also the main climate factor affecting soil conservation, and there was a very significant positive linear correlation between precipitation and SC (r = 1, p < 0.01). SC increased (decreased) by 36.13 t·hm −2 for every 10% increase (decrease) of precipitation. In addition, under the condition of constant precipitation, a 30% increase (decrease) in temperature could lead to an NPP decrease of 21.68% or increase of 13.04% and a 10% increase (decrease) in temperature would result in a 9.64% decrease or 6.63% increase in NPP. Under the condition of constant temperature, a 30% increase (decrease) in precipitation could cause NPP to increase 4.8% or decrease 8.4% and a 10% increase (decrease) in precipitation would result in a 1.8% decrease or 2.5% increase of NPP. The influence of temperature and precipitation on NPP is nonlinear, which was similar to the relationship between climate factors and WC. Meanwhile, the influence of temperature was greater than precipitation, and the decline of NPP caused by temperature rise (precipitation decrease) was greater than that caused by temperature fall (precipitation increase). increase), and the drastic change of precipitation presented greater impact on change amplitude of WC. Precipitation was also the main climate factor affecting soil conservation, and there was a very significant positive linear correlation between precipitation and SC (r = 1, p < 0.01). SC increased (decreased) by 36.13 t·hm −2 for every 10% increase (decrease) of precipitation. In addition, under the condition of constant precipitation, a 30% increase (decrease) in temperature could lead to an NPP decrease of 21.68% or increase of 13.04% and a 10% increase (decrease) in temperature would result in a 9.64% decrease or 6.63% increase in NPP. Under the condition of constant temperature, a 30% increase (decrease) in precipitation could cause NPP to increase 4.8% or decrease 8.4% and a 10% increase (decrease) in precipitation would result in a 1.8% decrease or 2.5% increase of NPP. The influence of temperature and precipitation on NPP is nonlinear, which was similar to the relationship between climate factors and WC. Meanwhile, the influence of temperature was greater than precipitation, and the decline of NPP caused by temperature rise (precipitation decrease) was greater than that caused by temperature fall (precipitation increase).

Impact of Land Use Change on Ecosystem Services
From 2000 to 2020, cultivated land, grassland, and unused land in BTH decreased, while construction land, forestland, and water area increased ( Table 6). The most remarkable feature of land transfer was the transformation from cultivated land to construction land. In the past 20 years, cultivated land decreased by 9539 km 2 , while construction land increased by 9869 km 2 , and the net transformation area was up to 8614 km 2 . In addition, the increased forestland and water area were mainly artificial economic forests, ponds, reservoirs, and other artificial water bodies, which were mainly converted from farmland and have limited effect on the promotion of ES. The decreased grassland mainly transformed into construction land, which directly led to the shrinkage of ecological land. From the perspective of spatial variation of land transfer (Figure 5), the conversion of ecological land such as forestland and grassland to construction land mainly concentrated in Yanshan-Taihang mountain area, while the transformation of cultivated land to

Impact of Land Use Change on Ecosystem Services
From 2000 to 2020, cultivated land, grassland, and unused land in BTH decreased, while construction land, forestland, and water area increased ( Table 6). The most remarkable feature of land transfer was the transformation from cultivated land to construction land. In the past 20 years, cultivated land decreased by 9539 km 2 , while construction land increased by 9869 km 2 , and the net transformation area was up to 8614 km 2 . In addition, the increased forestland and water area were mainly artificial economic forests, ponds, reservoirs, and other artificial water bodies, which were mainly converted from farmland and have limited effect on the promotion of ES. The decreased grassland mainly transformed into construction land, which directly led to the shrinkage of ecological land. From the perspective of spatial variation of land transfer (Figure 5), the conversion of ecological land such as forestland and grassland to construction land mainly concentrated in Yanshan-Taihang mountain area, while the transformation of cultivated land to construction land mainly occurred in North China Plain. Meanwhile, the intensity of land use transfer in 2010s was much greater than that in 2000s, which was closely related to the large-scale urban expansion in BTH. construction land mainly occurred in North China Plain. Meanwhile, the intensity of land use transfer in 2010s was much greater than that in 2000s, which was closely related to the large-scale urban expansion in BTH.     Under the BAU scenario, WC in Region F and G with great impact of human activities was slightly lower than that in the 2010s. Compared with the 2010s, NPP also showed an increasing trend in 2030 under the three scenarios, with an increase range of 5.16-5.94%, and land use change had little impact on NPP. Spatially, NPP of Region C and A in the northwest of BTH increased most significantly, with an increase of 12% and 17%, respectively, indicating that ecological engineering in Bashang Plateau and the upper reaches of Yongding River had achieved remarkable results. Compared with the base period, SC in 2030 would reduce around 2%, and simulation results under the three scenarios presented little difference. Region F and G distributed in North China Plain, with more than 99% of the areas having the slope of less than 2 • , and the SC was extremely low (2.30-2.54 t·hm −2 year −1 ), which would increase slightly by 3.1-4.6% in 2030. In addition, SC in the other five ecological regions decreased, especially in Yanshan-Taihang Mountain area (Region C and D). Forestland and grassland greatly reduced and the SC decreased significantly, with a rate of 4.6-5.4 t·hm −2 year −1 (1.8-3.2%) under the three scenarios. The HQ in 2030 is 2.6-11.5% lower than that in 2010s under the three scenarios, and the ecological environment is seriously degraded. As a direct reflection of the transformation between ecological land and construction land, HQ is significantly affected by human activities. In Region E, where urbanization degree was the highest, construction land accounted for more than 50%, and HQ decreased significantly, with a reduced rate of 0.11 (30.89%), especially under the ED scenario. However, in Region F, where farmland accounted for more than 70%, even if ecological land (forest land and grassland) was less than 10%, the HQ still increased by 13.49% (0.05) due to the reduction of construction land under the EC scenario. than 99% of the areas having the slope of less than 2°, and the SC was extremely low (2.30-2.54 t·hm −2 year −1 ), which would increase slightly by 3.1-4.6% in 2030. In addition, SC in the other five ecological regions decreased, especially in Yanshan-Taihang Mountain area (Region C and D). Forestland and grassland greatly reduced and the SC decreased significantly, with a rate of 4.6-5.4 t·hm −2 year −1 (1.8-3.2%) under the three scenarios. The HQ in 2030 is 2.6-11.5% lower than that in 2010s under the three scenarios, and the ecological environment is seriously degraded. As a direct reflection of the transformation between ecological land and construction land, HQ is significantly affected by human activities. In Region E, where urbanization degree was the highest, construction land accounted for more than 50%, and HQ decreased significantly, with a reduced rate of 0.11 (30.89%), especially under the ED scenario. However, in Region F, where farmland accounted for more than 70%, even if ecological land (forest land and grassland) was less than 10%, the HQ still increased by 13.49% (0.05) due to the reduction of construction land under the EC scenario.

Identification of Ecosystem Services Hotspots from 2000 to 2030
The spatial distribution of the primary ecosystem hotspots in BTH is shown in Figure  7. WC and NPP have gradually become the primary ecosystem service of BTH from 2000 to 2020. The primary hotspots of WC gradually spread from coastal areas (Region G) to mountainous areas (Regions B, C, and D). NPP hotspots widely distributed in North China Plain (Regions E and F) in 2030, and their area increased by about 1.06 times in 30 years. At the same time, the proportion of HQ hotspots decreased significantly and was expected to shrink by about 65%. The SC hotspots were very small: less than 0.001% of the total area. In the future, NPP will become the biggest contributor for the ecological importance of BTH, and WC is the fastest-growing ecosystem service function at present. Therefore, more attention should be paid to the protection and development of WC and NPP in future ecological protection policies of BTH.
According to the statistical results of multiple service hotspots ( Figure 8, Table 7), the area of Class 0 service hotspots accounted for the largest proportion (38~41%), and the area with Class 2 service hotspots accounted for a relatively small proportion (8~11%) during 2000-2030. There were significant differences in distribution of ecosystem service classes among different ecological regions. The ecosystem service of North China Plain (Regions E, F, and G) was poor, and Class 0 and Class 1 service hotspots were mainly distributed in this region (Supplementary Information, Table S2). The Class 0 hotspots area always exceeded 50%, and the proportion even reached 95% in Region G. The ecosystem service became better in Bashang Plateau and mountainous areas, and Class 2, 3, and 4 service hotspots were widely distributed. The area of Class 4 service hotspots in Region B always accounted for more than 30%, while the area of Class 4 service hotspots in Region D kept increasing, which was expected to reach 33% by 2030 under the EC scenario. Therefore, the core areas of future ecological protection in BTH are mainly concentrated in the northern Bashang Plateau and Yanshan-Taihang Mountain and improving the existing ecological protection policies in this region is the key point of ecological protection. The ecosystem service of urban and agro-ecological regions is poor, the contradiction between urban expansion and ecosystem service protection is acute, and their coordinated development is the main task of ecological protection in this region.

Identification of Ecosystem Services Hotspots from 2000 to 2030
The spatial distribution of the primary ecosystem hotspots in BTH is shown in Figure 7. WC and NPP have gradually become the primary ecosystem service of BTH from 2000 to 2020. The primary hotspots of WC gradually spread from coastal areas (Region G) to mountainous areas (Regions B, C, and D). NPP hotspots widely distributed in North China Plain (Regions E and F) in 2030, and their area increased by about 1.06 times in 30 years. At the same time, the proportion of HQ hotspots decreased significantly and was expected to shrink by about 65%. The SC hotspots were very small: less than 0.001% of the total area. In the future, NPP will become the biggest contributor for the ecological importance of BTH, and WC is the fastest-growing ecosystem service function at present. Therefore, more attention should be paid to the protection and development of WC and NPP in future ecological protection policies of BTH.
According to the statistical results of multiple service hotspots ( Figure 8, Table 7), the area of Class 0 service hotspots accounted for the largest proportion (38~41%), and the area with Class 2 service hotspots accounted for a relatively small proportion (8~11%) during 2000-2030. There were significant differences in distribution of ecosystem service classes among different ecological regions. The ecosystem service of North China Plain (Regions E, F, and G) was poor, and Class 0 and Class 1 service hotspots were mainly distributed in this region ( Supplementary Information, Table S2). The Class 0 hotspots area always exceeded 50%, and the proportion even reached 95% in Region G. The ecosystem service became better in Bashang Plateau and mountainous areas, and Class 2, 3, and 4 service hotspots were widely distributed. The area of Class 4 service hotspots in Region B always accounted for more than 30%, while the area of Class 4 service hotspots in Region D kept increasing, which was expected to reach 33% by 2030 under the EC scenario. Therefore, the core areas of future ecological protection in BTH are mainly concentrated in the northern Bashang Plateau and Yanshan-Taihang Mountain and improving the existing ecological protection policies in this region is the key point of ecological protection. The ecosystem service of urban and agro-ecological regions is poor, the contradiction between urban expansion and ecosystem service protection is acute, and their coordinated development is the main task of ecological protection in this region. Sustainability 2021, 13, x FOR PEER REVIEW 15 of 22

Identification of Ecological Security Pattern
Based on ESs, nature, economy, and society factors, the ecological security pattern of BTH was built ( Figure 9). The ecological sources are the concentrated distribution area of ecological land, of which forestland accounted for 80.13%, followed by grassland (15.69%). The ecological corridor mainly covered the forestland (72.53%) and grassland (7.51%) of Yanshan-Taihang Mountain. There were eighteen important ecological nodes in the region, of which fifteen were distributed in the Yanshan Mountain. The ecological buffer zones were distributed in a belt along the ecological corridor, and the medium level buffer zone was the largest (60.36%), which was widely concentrated in Yanshan-Taihang mountain area. The high-level buffer zone was small (17.63%), which was banded around the middle level buffer. The low-level buffer zone (the last guarantee zone to ensure the safety of ecological source area) was mainly distributed in the north of Yanshan Mountain and the east of Bashang Plateau. Based on the distribution of various components of the ecological security pattern in BTH, the Northern Yanshan Mountain and the Intermountain basin in the upper reaches of Yongding River were important parts of the ecological security pattern. The key ecological corridors connected the ecological sources in a counterclockwise direction from northeast to southwest to form an ecological security network with a semi encirclement situation, which integrated the typical grassland area in the east of Bashang Plateau, the deciduous broad-leaved forest area in Yanshan Mountain, the forestland, farmland, and grassland of the Intermountain basin in the upper reaches of Yongding River, and the deciduous broad-leaved forest in the Taihang Mountains. These regions formed important ecological barriers in the agricultural areas of BTH and North China Plain due to high coverage of natural vegetation, low impact of human activities, and high ecosystem services.

Identification of Ecological Security Pattern
Based on ESs, nature, economy, and society factors, the ecological security pattern of BTH was built ( Figure 9). The ecological sources are the concentrated distribution area of ecological land, of which forestland accounted for 80.13%, followed by grassland (15.69%). The ecological corridor mainly covered the forestland (72.53%) and grassland (7.51%) of Yanshan-Taihang Mountain. There were eighteen important ecological nodes in the region, of which fifteen were distributed in the Yanshan Mountain. The ecological buffer zones were distributed in a belt along the ecological corridor, and the medium level buffer zone was the largest (60.36%), which was widely concentrated in Yanshan-Taihang mountain area. The high-level buffer zone was small (17.63%), which was banded around the middle level buffer. The low-level buffer zone (the last guarantee zone to ensure the safety of ecological source area) was mainly distributed in the north of Yanshan Mountain and the east of Bashang Plateau. Based on the distribution of various components of the ecological security pattern in BTH, the Northern Yanshan Mountain and the Intermountain basin in the upper reaches of Yongding River were important parts of the ecological security pattern. The key ecological corridors connected the ecological sources in a counterclockwise direction from northeast to southwest to form an ecological security network with a semi encirclement situation, which integrated the typical grassland area in the east of Bashang Plateau, the deciduous broad-leaved forest area in Yanshan Mountain, the forestland, farmland, and grassland of the Intermountain basin in the upper reaches of Yongding River, and the deciduous broad-leaved forest in the Taihang Mountains. These regions formed important ecological barriers in the agricultural areas of BTH and North China Plain due to high coverage of natural vegetation, low impact of human activities, and high ecosystem services. Sustainability 2021, 13, x FOR PEER REVIEW 17 of 22 Figure 9. Ecological security pattern of BTH.

The Discussions on Driving Factor Analysis
At present, scenario simulation is widely used in analyzing the impact of land use change on ecosystem services, but it has hardly been used in analyzing the impact of climate change on ecosystem services. Previous studies mainly used multiple regression analysis or correlation analysis to explore the impact of climate on ecosystem services. They can only judge the dominant driving factors of ecosystem services through correlation or simply analyze the increase or decrease of ecosystem services caused by climate change. For example, precipitation is the main positive factor affecting water yield, but it is difficult to assess the variation magnitude of water yield caused by precipitation [63,64]. Indeed, the results obtained through climate change scenario simulation in this study not only fully demonstrate the main driving force of different ecosystem services but also further quantitatively analyze the difference of change rate between driving force and ESs. For example, temperature is the main factor leading to NPP change. Compared with temperature decrease, the change rate of NPP is more obvious when temperature increases.

The Discussions on Hotspots Analysis
In recent years, more and more studies have focused on classifying ecosystem services through hotspots analysis to explore the relationship between hotspot areas and urbanization. They believe that the research on this issue can determine the key areas of ecosystem service supply and explore the impact of land use change on ecosystem services [56,65,66]. However, previous studies mostly identified areas with high or low concentration of ecosystem service indicators through hotspots analysis (Getis-Ord Gi*) tool, which could not further identify the median-value area of ecosystem service level, resulting in the neglect of ecological protection [66][67][68]. This study comprehensively analyzes the ecosystem service hotspots in BTH from the perspective of multiple ecosystem service hotspots and primary ecosystem service hotspots and obtains five types of ecosystem service hotspots and two primary ecosystem services. Compared with other studies, the

The Discussions on Driving Factor Analysis
At present, scenario simulation is widely used in analyzing the impact of land use change on ecosystem services, but it has hardly been used in analyzing the impact of climate change on ecosystem services. Previous studies mainly used multiple regression analysis or correlation analysis to explore the impact of climate on ecosystem services. They can only judge the dominant driving factors of ecosystem services through correlation or simply analyze the increase or decrease of ecosystem services caused by climate change. For example, precipitation is the main positive factor affecting water yield, but it is difficult to assess the variation magnitude of water yield caused by precipitation [63,64]. Indeed, the results obtained through climate change scenario simulation in this study not only fully demonstrate the main driving force of different ecosystem services but also further quantitatively analyze the difference of change rate between driving force and ESs. For example, temperature is the main factor leading to NPP change. Compared with temperature decrease, the change rate of NPP is more obvious when temperature increases.

The Discussions on Hotspots Analysis
In recent years, more and more studies have focused on classifying ecosystem services through hotspots analysis to explore the relationship between hotspot areas and urbanization. They believe that the research on this issue can determine the key areas of ecosystem service supply and explore the impact of land use change on ecosystem services [56,65,66]. However, previous studies mostly identified areas with high or low concentration of ecosystem service indicators through hotspots analysis (Getis-Ord Gi*) tool, which could not further identify the median-value area of ecosystem service level, resulting in the neglect of ecological protection [66][67][68]. This study comprehensively analyzes the ecosystem service hotspots in BTH from the perspective of multiple ecosystem service hotspots and primary ecosystem service hotspots and obtains five types of ecosystem service hotspots and two primary ecosystem services. Compared with other studies, the distribution trend of Class 4 hotspot is consistent with the other study [56], which is mainly distributed in the northern mountainous areas with high vegetation coverage, and the distribution of Class 0 and 1 hotspots basically coincides with other studies, mainly distributed in North China Plain. However, in our study the non-hotspots are further divided into Class 2 and 3 hotspots areas. Meanwhile, the primary ecosystem service hotspots analysis showed that NPP gradually took the place of WC and became the regional dominant ecosystem services. Therefore, in the process of formulating ecological security planning measures, we can refer to different levels of ecosystem service hotspots for ecological protection with different intensities and strengthen the key protection of WC and NPP, which will greatly improve the efficiency of ecological protection.

Linking Ecological Security Pattern and Planning Strategy in BTH
BTH is the largest and the most developed economic area in North China. The unreasonably rapid economic development at the cost of ecological loss in the early stage makes the region face problems such as water shortage, biodiversity loss, soil and water loss, ecological land degradation, and so on [69]. In recent years, more and more studies have realized the necessity and urgency of ensuring the ecological security pattern for the future development planning of BTH. However, previous studies on the ecological security pattern mostly focused on calculating the ecological environment carrying capacity from the administrative division scale, and reflecting the evolution of ecological security estimation through ecological security index, the evaluation results cannot reflect the impact of regional geographical environment on ecological security pattern [70,71]. Some other studies focus on the impact of the natural environment on the regional ecological security pattern, but the research mainly focuses on the construction of ecological security pattern in small regions, such as Xiong'an and the mountainous area in the north of Beijing [72][73][74]. Few studies have conducted a comprehensive analysis of the regional ecological security pattern of BTH. Therefore, based on the resistance surface in previous studies, this study constructs a set of ecological security patterns that can reflect regional natural and social factors by integrating natural factors, socio-economic factors, and ecosystem services [24].
A reasonable ecological security pattern can provide scientific guidance for optimizing the layout of ecological land. The northern low-level buffer zone of ecological security should reduce human intervention on ecology and is designated as an extremely important zone of ecological security. Developing economic activities on the basis of protection and restoration is the key issue in middle-level buffer zones, for example, forest parks can be properly developed. Other examples are taking the Beijing Winter Olympic Games as an opportunity to promote the construction of white snow and ice and green industries [75], preventing the construction of stadiums for the Winter Olympic Games from destroying the landform and landscape, reducing biodiversity, irrational use of water resources, environmental pollution, and other ecological problems [76][77][78]. High level buffer zone is an ideal area to protect regional ecological security. Moderate economic development activities will not block the ecological corridor and affect the ecological security pattern, so economic construction should be actively played to promote and protect ecological construction.

Conclusions
In summary, this study established a comprehensive ecosystem service analysis and regional security pattern evaluation system. The main conclusions are as follows: the inter-annual changes of WC and NPP showed an increasing trend from 2000 to 2020, and the increased rate in Yanshan-Taihang mountain area was significantly higher than that in the other area. Precipitation and temperature were the key climate factors controlling WC and NPP, respectively, and the impact of future land use change caused by policy orientation on WC cannot be ignored. At present, WC is the fastest-growing ecosystem service function in BTH, and NPP would become the greatest contributor to ecological importance in the future. Besides, SC was mainly affected by precipitation and showed an insignificant decreasing trend, while HQ was significantly reduced by the accelerated urbanization expansion, and the ecological environment tended to deteriorate in the past 20 years. Yanshan-Taihang mountain area has always been the key area of ecosystem service and the core region of ecological security in BTH. Due to high natural vegetation coverage and little effect by human activities, the ecological source was mainly distributed here, so the protection value here is pretty important. The result provides an effective direction and a scientific reference for the future territorial space planning and ecological restoration in BTH; the policy maker should reasonably control urban expansion in accordance with the ecological security pattern and pay attention to the development of areas with high NPP and WC.
However, there are still some problems to be improved in future research. The maximum light energy value required for CASA model simulation refers to the research results of light energy utilization calculation in China without better considering the special situation of the region. How to obtain the real data in the study area is the focus of future research. Meanwhile, our study only studies four supply services; more ecosystem services, such as supervision, culture, and supporting services, are not comprehensively selected to conduct a more comprehensive assessment of ecosystem services in BTH, which might be addressed in further studies.