Relationship of Ecosystem Services in the Beijing–Tianjin–Hebei Region Based on the Production Possibility Frontier

The supply and demand of ecosystem services are affected by land use. Only a few studies have conducted in-depth quantitative analyses. This study adopted the Beijing–Tianjin–Hebei region as the research area. The CLUMondo model was adopted to infer the land-use pattern under protection, development, and natural scenarios in 2035. Moreover, the InVEST model was utilized to evaluate carbon sequestration, water yield, and soil conservation under multiple land-use patterns. The production possibility frontier was drawn to visualize the trade-off relationship further. The trade-off intensity index was calculated to quantify the magnitude of the trade-off. (1) Under the development scenario, the accelerated expansion of urbanized land will occupy a large amount of arable and forest land, which should be planned and controlled. (2) The trade-off and synergistic relationships could be transformed under the different land-use scenarios. (3) The production possibility frontier curve for each ecosystem service trade-off and the optimal value of the trade-off configuration were plotted for the different scenarios. The trade-off intensity of ecosystem services was also calculated. This study combined ecosystem services with land-use regulations and revealed the link between ecosystem services and regional land-use pattern change. The aim is to provide a reference for the synergistic progress of the ecological economy in the Beijing–Tianjin–Hebei region.


Introduction
Ecosystems provide natural environmental conditions and functions that sustain human survival. The benefits that humans obtain directly or indirectly from ecosystems are called ecosystem services [1], which are essential to maintaining the well-being of humans. The concept of ecosystem service was established to improve and protect the ecological environment [2]. Several scholars have conducted extensive research on ecosystem services [3,4]. The Millennium Ecosystem Assessment clarified the contribution of ecosystem services to human well-being. Ecosystem services refer to the full benefits that humans gain from natural ecosystems [5]. Considering the relationship between ecosystem services and ecosystem structure, process and function, ecosystem services usually include four categories: provisioning, regulating, cultural, and supporting services [6,7]. Approximately 60% of ecosystem degradation in the past 50 years was due to the urbanization and population growth worldwide. Therefore, research on ecosystem services has become increasingly urgent. Ecosystems and ecosystem services are widely studied quantitatively, Li et al. utilized Bayesian belief networks to analyze the efficiency of ecosystem services [29]. In line with the distribution of the three ecosystem services in different scenarios, a three-dimensional coordinate system was drawn [30], and the optimized target efficiency curve was obtained.
Land-use change is an important factor that affects trade-offs in ecosystem services [31,32]. Given population growth and urbanization acceleration, the overall land-use pattern is constantly changing, thus affecting the balance of ecosystem services [26,33]. Scenario simulation can be used to select a series of socio-economic and ecological indicators that can be combined with relevant policy planning to find out how human activities and land-use pattern changes influence the balance of ecosystem services [34]. Existing studies have evaluated the relationship among multiple ecosystem services by simulating various land-use scenarios [30,35,36]. Therefore, ecosystem service evaluation under the simulated land-use change scenarios is an effective method to study environmental protection policies. This approach contributes a new idea for identifying the relationship between ecosystem services [37]. It simulates land-use change scenarios on a small scale and fully considers the impact of land-use intensity during the simulation [38]. The final result can reflect the current intensity of human activities and the impact of social and economic development on land-use patterns [39]. The model has been widely used in land-use simulation [40,41].
The introduction of the production possibility boundary provides a new idea for the quantification of the trade-off relationship. On the basis of this method, the strength of the trade-off between different combinations can be quantitatively described, and the trade-off between different services can be displayed directly. A few studies on the trade-off relationship between ecosystem services have used this method.
The Beijing-Tianjin-Hebei region is a highly dense political, economic, cultural, and populated area comprising northern cities in the country. It occupies a place in the country's overall development strategy and economic structure. In 2014, the Chinese government proposed the Beijing-Tianjin-Hebei Coordinated Development Strategy, which organically integrates the advantageous resources of the three places. This integration has increased the diversity of cross-regional cooperation. The region shares the same development goals in ecology, economy, and society. By evaluating ecosystem services in the Beijing-Tianjin-Hebei region under different land-use scenarios, this study analyzes the correlation and trade-off between different services. The PPF curve is adopted to demonstrate the strength of the trade-off between these services visually. Such a depiction provides a reference foundation for the new pattern of coordinated progress in the region.

Study Area
The Beijing-Tianjin-Hebei region is the capital economic circle in China, covering Beijing, Tianjin, and Hebei provinces. It has a large area and gradually tilts from the northwest to the southeast. The annual average temperature is approximately 11 • C and the annual average precipitation is approximately 500 mm. The regional land-use types present a significant decrease in arable land and a large increase in urban and rural built-up land. Given the acceleration of urbanization, the land-use pattern in this area has gradually become increasingly complicated.
The study area occupies approximately 218,000 km 2 , with the population more than 100 million. Given the acceleration of urbanization, rapid population growth, and backward ecological infrastructure in recent years, the consumption of regional natural resources has become particularly high, and a series of ecological and environmental issues, such as air pollution, land desertification, and soil erosion, have become increasingly prominent. The region is faced with the status quo of excessive ecological load [42], which restricts further regional development and progress. Land-use and cover changes depend largely on nature, social economy, and human activities. The location of the study area is shown in Figure 1. largely on nature, social economy, and human activities. The location of the study area is shown in Figure 1.

Data Sources
The data are shown in Table 1. Land-use data have been separated into six types: arable land, woodland, grassland, water, urbanized land, and unused land. Infrastructure elements were obtained by importing land-use data into ArcGIS for processing. Data preprocessing and extraction were performed in the ArcGIS 10.2 platform. The spatial resolution of raster data was 1 km × 1 km, and Krasovsky_1940_Albers projection was used.

Data Sources
The data are shown in Table 1. Land-use data have been separated into six types: arable land, woodland, grassland, water, urbanized land, and unused land. Infrastructure elements were obtained by importing land-use data into ArcGIS for processing. Data preprocessing and extraction were performed in the ArcGIS 10.2 platform. The spatial resolution of raster data was 1 km × 1 km, and Krasovsky_1940_Albers projection was used. Ecosystem supply and ecosystem relationships are deeply affected by land-use changes [43]. The CLUMondo model is a new method for design of land-use change patterns; it is driven by the regional commodity demand and considers the factors of the conversion between different regions to simulate changes in the land system [38]. By using the CLUMondo model and with 2015 as the initial year, different driving factors were set, and the conversion parameters between different land-use types were continuously revised according to demand. The land-use conditions of the study area in 2035 under the three scenarios were simulated [44]. This research analyzed and explained land-use problems in the study area through land-use change. It established relevant models to understand the process of land-use change and the driving force for this change. The study also simulated and predicted future scenarios to determine the trend of future land-use changes accurately [45]. Initially, this research analyzed the land-use changes in 2005-2015 and referred to the research results of other scholars [46,47]. Ten natural, socio-economic, and demographic factors were recognized as driving factors. The information on the driving factors is shown in Table 2. The selection of drivers is a key factor in the spatial pattern of land-use types. The area under the receiver operating curve (AUC) is an assessment index to measure the calculation accuracy of a regression model. The results of this study suggested that the AUC values between each category and the driving factors were all over 0.75. Hence, the selected of driving factors can comprehensively explain the spatial patterns of land-use types. The use of forest and built-up land area was adopted as the demand factor. The CLUMondo model was utilized to simulate and predict land-use changes under various scenarios for exploring the trade-off relationship between multiple ecosystem services in each scenario. The first scenario is the natural scenario, which is consistent with the analysis of land-use change in 2005-2015. Land-use change is expected to continue according to the existing natural development pattern until 2035; hence, no special constraints were set for the CLUMondo model. The next scenario is the development scenario, which accelerates the urbanization process based on the overall land-use planning and socio-economic development in the Beijing-Tianjin-Hebei region. It converts a large amount of arable, forest, and unused land into built-up land. In this scenario, the conversion resistance value for all land-use types was set to 1.
In the conservation scenario, the ecological environmental in the Beijing-Tianjin-Hebei region is vigorously protected. The area of arable, forest, and grassland that is converted to built-up land is strictly limited. In addition, the development of urbanization is slowed down. In this scenario, the conversion resistance of arable, forest, and grassland to built-up land was zero.

Water Yield
In this work, the annual water yield module of the InVEST model was used to assess the annual water yield of the Beijing-Tianjin-Hebei region in 2035 under multiple scenarios. The water content module combines the land-use type, soil depth, topography, climate, and other reasons based on the Budyko hydrothermal coupling balance hypothesis and the annual average precipitation [48]. The formula precipitation minus the actual evaporation was used to calculate the grid water content.
where Y xj is the annual water yield of grid unit x on land cover type j. It mainly includes surface runoff, soil water content, litter water holding capacity, and canopy interception. AET x is the annual actual evapotranspiration of pixel x, and P x is the annual precipitation of pixel x. The Z value, which is a seasonal constant characterizing rainfall characteristics, must be inputted to the water yield model. This study refers to scholars' research on precipitation characteristics and set the Z value to 3.8 [43]. The biophysical properties of each land-use type involved in the model are illustrated in Table 3.

Soil Conservation
Soil conservation is an important ecosystem service for the Beijing-Tianjin-Hebei region. The sediment delivery ratio module of the InVEST model was used in this research to evaluate the soil conservation under different scenarios in 2035. Given that the SDR module is based on the universal soil loss equation (USLE), which considers the capability of the site itself to intercept upstream sediment and adds sediment retention to the reservoir data, the results of the model are highly realistic and scientifically accurate. Soil conservation, potential soil loss, and actual soil loss are calculated as follows [48]: where USLE is the actual soil loss in the original land-use cover, RKLS is the potential soil loss for bare soil, R is the rainfall erosivity factor, K is the soil erodibility factor, LS is the slope length-gradient factor, C is the crop-management factor, and P is the support practice factor. The values of C and P combined with the characteristics of the study area and information previous studies are shown in Table 4 [49].

Carbon Storage
The carbon storage and sequestration model was used to select distribution of surface land-use and cover types for the calculation and analysis of the four basic carbon pools. Then, the module was utilized to generate a spatial distribution map of carbon storage. In this work, the aboveground biogenic carbon stock and the below-ground biogenic carbon stock were combined into vegetation carbon stock, which was not considered in this study because dead organic carbon stock is difficult to observe, and its stock is unavailable. The carbon storage calculation principle is as follows [48]: where C tot is the total regional carbon storage, C above is the aboveground biological carbon storage, C below is the underground biological carbon storage, C soil is the soil carbon storage, and C dead is the dead organic carbon storage.
With reference to previous studies [50], the carbon density adopted in this research refers to the calculation of the carbon density of each carbon pool in the Beijing-Tianjin-Hebei region, as shown in Table 5.

Plotting the PPF Curve
In this work, the trade-off among the water yield, soil conservation, and carbon storage was expressed as PPF. First, the values of each ecosystem service were extracted, and the trade-offs and synergies between two ecosystem services were analyzed using the correlation analysis tool in R language. Second, the two ecosystem services presenting tradeoffs were normalized to values between 0 and 1 after processing. Lastly, the normalized layers of the two ecosystem services were divided to obtain the ratio layer. Each cell in the ratio layer represented the ratio of the two ecosystem services in the corresponding geographic location. Then, each cell in the ratio layer was arranged in ascending order, and the values of the two ecosystem services in the corresponding geographic location were summed up. The curves were drawn based on the final results. The best trade-off values of the two ecosystem services were visualized [26,51].
The outcome of the changes in the research setting trade-offs reveals the best combination of the two ecosystem services [20]. On this basis and with reference to previous studies, the shortest distance from the mean point to the curve corresponding to the two services is represented as the trade-off intensity. Assuming that the coordinates of the mean point are at the (x 0 , y 0 ) position and the equation of the curve is g = f (x), the coordinates of any point on the curve are (x, g). The distance calculation formula is as follows: Therefore, the strength of the trade-off is H min . A large value of H min indicates a strong trade-off relationship. A small value of H min indicates a weak trade-off relationship.

Analysis of Land-Use Simulation Results under Different Scenarios
Under the conservation scenario, in 2035, the area of arable land will decrease by 0.76%, the area of woodland will increase by 12.22%, the area of built-up land will increase by 28.72%, the water area will decrease by 10.91%, and the area of unused land will decrease by 8.13%. Under the development scenario, the rate of urbanization will accelerate sharply in 2035. The area of arable land will decrease by 15.22%, the area of woodland will decrease by 5.94%, the area of built-up land will increase by 95.65%, the water area will decrease by 28.89%, and the area of unused land will decrease by 12.98%. Compared with the situation in 2015, the area of each category shows evident changes. Under the natural scenario, in 2035, the area of arable land will decrease by 2.72%, the area of woodland will decrease by 2.31%, the area of built-up land will increase by 45.96%, the water area will decrease by 8.63%, and the area of unused land will decrease by 7.96%.
Comparison of the three scenarios showed that the increase and decrease in each category in the natural scenario were in the middle. The area of built-up land in the development scenario increased the most, and only the area of woodland in the conservation scenario showed an increasing trend. The kappa values of all land-use types were above 0.8, indicating that the simulation results are highly consistent with the actual land-use and that the accuracy of the model can meet the simulation requirements ( Figure 2).

Spatial Distribution of Ecosystem Services
The distributions of water yield under each scenario change are shown in Figure 3.

Spatial Distribution of Ecosystem Services
The distributions of water yield under each scenario change are shown in Figure 3.  In terms of soil conservation, the spatial pattern changed significantly under the three scenarios. The range of soil conservation in the conservation scenario was 0 t to 8.989×10 7 t. The range decreased significantly, but the density of the high-value distribution area increased. The range of soil conservation in the development scenario, which was 0 t to 318954 t, was partially reduced in the distribution area of its high values compared with the situation in 2015. Soil conservation under the natural scenario ranged from 0 t to 9.147×10 6 t; the high values were mainly concentrated in the western and northern parts of the study area, and the low values were mainly distributed in the southeastern plains of the study area. Soil conservation showed a certain trend of increase under the conservation scenario, indicating that the soil conservation capacity of the study area could be improved by relevant ecological policy constraints (Figure 4). In terms of soil conservation, the spatial pattern changed significantly under the three scenarios. The range of soil conservation in the conservation scenario was 0 t to 8.989 × 10 7 t. The range decreased significantly, but the density of the high-value distribution area increased. The range of soil conservation in the development scenario, which was 0 t to 318,954 t, was partially reduced in the distribution area of its high values compared with the situation in 2015. Soil conservation under the natural scenario ranged from 0 t to 9.147 × 10 6 t; the high values were mainly concentrated in the western and northern parts of the study area, and the low values were mainly distributed in the southeastern plains of the study area. Soil conservation showed a certain trend of increase under the conservation scenario, indicating that the soil conservation capacity of the study area could be improved by relevant ecological policy constraints (Figure 4). The distribution of carbon storage under different scenarios is shown in Figure 5. The carbon storage under the conservation scenario gradually shifted its high-value area from the central area to the western area compared with 2015. In addition, the lowest value was higher than the highest value in 2015. The carbon storage distribution changed considerably under the development scenario, and the high value areas greatly decreased (34,204 Mg to 198,100 Mg). Carbon storage in the natural scenario did not change much compared with 2015. The distribution of carbon storage under different scenarios is shown in Figure 5. The carbon storage under the conservation scenario gradually shifted its high-value area from the central area to the western area compared with 2015. In addition, the lowest value was higher than the highest value in 2015. The carbon storage distribution changed considerably under the development scenario, and the high value areas greatly decreased (34,204 Mg to 198,100 Mg). Carbon storage in the natural scenario did not change much compared with 2015.

Ecosystem Service Trade-Off Relationships
Non-parametric tests were performed on the ecosystem service capacities of water yield, carbon storage, and soil conservation. The tests revealed that all three data series did not follow the normal distribution. Therefore, the Pearson correlation analysis tool in R was used to determine whether a trade-off occurred among the various ecosystem services. The results are shown in Figures 6 and 7.
The results revealed a high correlation among all services. The correlation coefficient (R > 0) for each ecosystem service under the conservation scenario indicated a synergistic relationship among ecosystem services. Under the development scenario, the correlation coefficient between water yield and carbon storage was R = −0.61, and the correlation coefficient between water yield and soil conservation was R = −0.17, indicating a certain trade-off between water yield and these two services. By contrast, the correlation coefficients between water yield and carbon storage and between water yield and soil retention in the natural scenario were negative and had values of −0.19 and −0.086, respectively. These results suggest that the relationships among ecosystem services shifted in the conservation scenario.

Ecosystem Service Trade-Off Relationships
Non-parametric tests were performed on the ecosystem service capacities of water yield, carbon storage, and soil conservation. The tests revealed that all three data series did not follow the normal distribution. Therefore, the Pearson correlation analysis tool in R was used to determine whether a trade-off occurred among the various ecosystem services. The results are shown in Figures 6 and 7.
The results revealed a high correlation among all services. The correlation coefficient (R > 0) for each ecosystem service under the conservation scenario indicated a synergistic relationship among ecosystem services. Under the development scenario, the correlation coefficient between water yield and carbon storage was R = −0.61, and the correlation coefficient between water yield and soil conservation was R = −0.17, indicating a certain trade-off between water yield and these two services. By contrast, the correlation coefficients between water yield and carbon storage and between water yield and soil retention in the natural scenario were negative and had values of −0.19 and −0.086, respectively. These results suggest that the relationships among ecosystem services shifted in the conservation scenario.

Ecosystem Service Trade-Off Strength and Optimal Allocation
The results of the correlation between ecosystem services showed trade-offs between water yield and carbon storage and between water yield and soil storage in the natural and development scenarios. The optimal configuration curve between the two ecosystem services was plotted to illustrate the optimal configuration of ecosystem services. In the PPF curve, the maximum productivity of service supply was represented by the points on the curve. The situation where the service could be provided but not at the optimal allocation at this point was indicated by the points inside the curve. The situation where the service could not be provided with limited land resource allocation was indicated by the points outside the curve.
The results showed trade-offs between carbon storage and water yield services and between soil conservation and water yield services only in the natural and development scenarios. Therefore, trade-off curves were drawn according to the normalization of ecosystem services. The trade-off curves among the different ecosystem services in the two scenarios were different. The PPF curves between carbon storage and water yield for the natural and development scenarios are shown in Figure 8. The slope of the curve for the natural scenario decreased, indicating that the trade-off was gradually decreasing. Meanwhile, the slope of the curve increased gradually under the development scenario, indicating that the trade-off was gradually increasing. The PPF curves between soil conservation and water yield for the two scenarios are shown in Figure 9. The lower part of the curve in the development scenario is concave. By contrast, the upper part of the curve in the natural scenario is concave, showing that the trade-off relationship between the two services is gradually becoming significant. The mean points in the development and natural scenarios represented by points A and B, respectively, are located below the corresponding curves. In both cases, the relations between carbon storage and water yield and between soil retention and water yield do not reach the PPF curve. This scenario shows that the two services are not in an optimal configuration, indicating the potential for optimization. In Figure 8, under the development scenario, the cumulative water yield at point A1 is similar to that at point A, and the cumulative carbon storage at point A2 is similar to that at point A. The trade-off relationship between water yield and carbon storage under the development scenario is the best at points A1 and A2. Similarly, the trade-offs between water yield and carbon storage under the natural scenario are optimally configured at points B1 and B2. In Figure 9, under the development scenario, the cumulative water yield of point A1 is similar to that of point A. In addition, the cumulative soil conservation of point A2 is similar to that of point A. Points A1 and A2 are the best trade-off relationships between water yield and soil conservation under the development scenario. Similarly, points B1 and B2 are the optimal allocation points for the trade-off relationship between water yield and soil conservation under natural scenarios.

Ecosystem Service Trade-Off Strength and Optimal Allocation
The results of the correlation between ecosystem services showed trade-offs between water yield and carbon storage and between water yield and soil storage in the natural and development scenarios. The optimal configuration curve between the two ecosystem services was plotted to illustrate the optimal configuration of ecosystem services. In the PPF curve, the maximum productivity of service supply was represented by the points on the curve. The situation where the service could be provided but not at the optimal allocation at this point was indicated by the points inside the curve. The situation where the service could not be provided with limited land resource allocation was indicated by the points outside the curve.
The results showed trade-offs between carbon storage and water yield services and between soil conservation and water yield services only in the natural and development scenarios. Therefore, trade-off curves were drawn according to the normalization of ecosystem services. The trade-off curves among the different ecosystem services in the two scenarios were different. The PPF curves between carbon storage and water yield for the natural and development scenarios are shown in Figure 8. The slope of the curve for the natural scenario decreased, indicating that the trade-off was gradually decreasing. Meanwhile, the slope of the curve increased gradually under the development scenario, indicating that the trade-off was gradually increasing. The PPF curves between soil conservation and water yield for the two scenarios are shown in Figure 9. The lower part of the curve in the development scenario is concave. By contrast, the upper part of the curve in the natural scenario is concave, showing that the trade-off relationship between the two services is gradually becoming significant. The mean points in the development and natural scenarios represented by points A and B, respectively, are located below the corresponding curves. In both cases, the relations between carbon storage and water yield and between soil retention and water yield do not reach the PPF curve. This scenario shows that the two services are not in an optimal configuration, indicating the potential for optimization. In Figure 8, under the development scenario, the cumulative water yield at point A1 is similar to that at point A, and the cumulative carbon storage at point A2 is similar to that at point A. The trade-off relationship between water yield and carbon storage under the development scenario is the best at points A1 and A2. Similarly, the trade-offs between water yield and carbon storage under the natural scenario are optimally configured at points B1 and B2. In Figure 9, under the development scenario, the cumulative water yield of point A1 is similar to that of point A. In addition, the cumulative soil conservation of point A2 is similar to that of point A. Points A1 and A2 are the best trade-off relationships between water yield and soil conservation under the development scenario. Similarly, points B1 and B2 are the optimal allocation points for the trade-off relationship between water yield and soil conservation under natural scenarios.

Discussion
Nowadays, the method of combining land-use scenario simulation and ecosystem service assessment is becoming widely used [52][53][54]. The future relationship between land-use and ecosystem services can be explored by this method [55]. In this study, landuse data in 2005 and 2015 were selected to analyze the 10-year land-use changes in the Beijing-Tianjin-Hebei region, and the land-use transfer matrix was derived. The CLUMondo model was adopted to simulate the land-use map under different scenarios  Table 6 presents the intensity of the trade-off between carbon storage and water yield, that is, natural scenario (0.051) < development scenario (0.199). The strength of the trade-off between soil conservation and water yield follows the order natural scenario

Discussion
Nowadays, the method of combining land-use scenario simulation and ecosystem service assessment is becoming widely used [52][53][54]. The future relationship between landuse and ecosystem services can be explored by this method [55]. In this study, land-use data in 2005 and 2015 were selected to analyze the 10-year land-use changes in the Beijing-Tianjin-Hebei region, and the land-use transfer matrix was derived. The CLUMondo model was adopted to simulate the land-use map under different scenarios in 2035 in accordance with the transfer matrix. The CLUMondo model's setting for land demand avoided the limitation of the land-use area. At the same time, it combined land-use change with multiple land-use demands so that spatial and temporal land-use changes in multiple scenarios could be simulated accurately [56]. Under the conservation scenario, the expansion of built-up land was restrained to a certain extent, and the proportion of forest land increased. The highest amount of carbon storage services was found in this scenario, and the high-value areas were mainly distributed in forest land, indicating that the supply benefits of carbon storage services could be enhanced by woodland. This outcome is consistent with the results of Deng et al. [57]. The highest values of water yield were found in the development scenario, and most of the high-value areas were located on woodland, which is consistent with the research results of Gao et al. [58]. Water yield can be affected by evapotranspiration from the surface, soil root depth, and soil porosity [59]. In addition, the soil porosity of built-up land is small [60]. Since built-up areas are usually regarded as areas with zero water retention capacity, where artificial building growth leads to an increase in the impervious surface, which reduces the amount of precipitation evaporation, and infiltration, leading to an increase in water yield.
The dominant functions in each region in the development process of the Beijing-Tianjin-Hebei region can be determined through the division of ecological function zones. The blind expansion of urban built-up land must be managed and controlled, and the protection of arable land and natural habitat should be strengthened to reduce the loss of ecosystem services. Then, the overall arrangement and deployment of the land-use structure layout can be implemented to avoid planning beyond the maximum range that the ecosystem services can tolerate. The integrity of the ecosystem structure and the continuity of the process can be improved and maintained by planning urban land-use, adhering to ecological priorities, promoting the development of economical and intensive land-use, and effectively enhancing ecosystem service functions.
The synergy and trade-off among various ecosystem services can be studied quantitatively by introducing PPF. Apart from macro analyses, such as correlation, the degree of influence among various ecosystem services was analyzed [61]. PPF can be combined with other model frameworks or decision-making systems to determine the optimal combination of ecosystem services and ultimately provide a scientific basis for improving the environment and ecological restoration. In most previous studies, production possibility boundaries were used primarily to evaluate the trade-offs or synergies between ecosystem services qualitatively [62]. In essence, they only compare the benefits of different programs or goals, weigh the advantages and disadvantages of various service combinations, and identify the best solution. However, the strength of the balance has not been determined yet. In the future, the strength of the trade-off at each point can be determined through PPF, and the nature of the trade-off can be expressed visually [20]. The coefficients of the objective function for the new optimization simulation can be adjusted according to preference information, and the operation can be repeated to obtain the most ideal solution [63]. A wide range of land-use policies can be evaluated through PPF, which can accurately extract key factors in ecological, agronomic, and economic coupling systems. Its application closely links multiple fields, thus allowing scholars to conduct interdisciplinary research and analysis. In addition, the use of PPF can prevent the complicated conceptualization of problems.
Most models have certain limitations due to many factors. First, the research and development of the InVEST model entail a certain geographical orientation, and a certain degree of inapplicability emerges when evaluating the characteristics of ecosystem services in different regions. Second, to reduce the difficulty of using the model, developers simplified the calculation process of the various ecological services of the simulated ecosystem, leading to the lack of scientific calculation results. Lastly, the carbon sequestration module does not fully consider the dynamic changes in the carbon sequestration rate of different vegetation types and the impact of the carbon cycle process between different carbon pools, resulting in uncertainty in the final assessment results [64]. The shortcoming of this work is that the number of research scenarios was insufficient. Moreover, the conversion resistance parameter was set in accordance with the results of the land-use transfer matrix of the study area in 2005-2015. To make the model run successfully and achieve the expected results, the parameters were debugged several times until the final results were obtained. In this process, several human factors that could interfere with the simulation results were added, which also exerted some influence. In future research, a larger number of scenarios will be set, additional ecosystem service trade-offs will be considered for expansion, and scientific and authoritative methods will be explored further to convert resistance parameter values and make the final results more convincing.

Conclusions
Three ecosystem services, namely, water yield, soil conservation, and carbon storage, in 2015 were estimated for the Beijing-Tianjin-Hebei region. The estimation was used to simulate the three ecosystem services under three different scenarios (conservation, development, and natural scenarios) in 2035 and determine the relationship between two ecosystem services. The trade-offs were expressed visually through PPF, which also revealed changes in the ecosystem services under different land-use patterns. This research simulated the analysis of the trade-off intensity of ecosystem service value under three land-use scenarios and can be extended to explore the trade-off relationship between ecosystem services under additional scenarios in the future. The main conclusions are as follows.
(1) Under the development scenario in 2035, the sprawl of urbanized land in the Beijing-Tianjin-Hebei region will accelerate. This expansion will lead to the occupation of a large amount of arable land and woodland, which is not conducive to the construction of ecological civilization and the sustainable and healthy development of cities in the future. Hence, experts need to plan and control the conversion of a certain amount of arable land and woodland into built-up land. (2) A trade-off relationship exists between water yield and soil conservation and between water yield and carbon storage services under development and natural scenarios. Synergistic relationships exist among ecosystem services under the conservation scenario. In addition, the trade-off and synergy relationships can be transformed under different land-use scenarios. (3) In this study, the conservation scenario was found to have the highest value of carbon storage and soil conservation, and the water yield was also at a high level. Wood-land is an important supply area for carbon storage. It has high soil and water conservation capacity. Hence, appropriate measures should be implemented to improve the soil and water conservation function in the region while enhancing the value of ecosystem services by returning the land to the forest to a certain extent. (4) The PPF curves between ecosystem services under different scenarios in 2035 were plotted, and the trade-off intensity of ecosystem services was calculated. In particular, visual representation of the results of agricultural planning decisions through PPF can help guide the direction of subsequent interdisciplinary research and policy. The findings showed that by controlling the conversion of land-use types, the planning goals of ecological service functions can be achieved. This work provides an important theoretical basis for the sustainable development of land resources in the Beijing-Tianjin-Hebei region.