Spatiotemporal Change Analysis and Multi-Scenario Modeling of Ecosystem Service Values: A Case Study of the Beijing-Tianjin-Hebei Urban Agglomeration, China

Abstract

Ecosystem service value (ESV) assessment is a crucial indicator of regional ecological quality and ecological management effectiveness. Ecosystem services (ES) provide the environmental foundation for human existence and social advancement. However, the future course of land use change (LUC) in urban agglomerations and its implications for human society remains uncertain, which presents a challenge to maintaining a balance between ecological service functions and regional socioeconomic growth. This paper took the Beijing-Tianjin-Hebei (BTH) urban agglomeration as an example and used the future land use simulation (FLUS) model to project the spatial distribution of land use under the natural development scenario (NDS), food security scenario (FSS), and ecological priority scenario (EPS) of BTH in 2030, 2040, and 2050. Next, the changes to ESV under various scenarios were investigated through the equivalent coefficient method. In order to make more targeted recommendations for regional development, the study also used hotspot analyses to explore the impacts of LUCs on ESV. The results showed that: (1) from 2000 to 2020, the LUC in the BTH was dramatic and mainly focused on the interconversions among the three land use categories of cropland, grassland, and built-up land. The total ESV demonstrated the tendency to decrease from CNY 386,859.89 × 10⁶ in 2000 to CNY 371,968.78 × 10⁶ in 2020. (2) Compared with 2020, the ESV in BTH in 2050 under the FSS loses 16,568.78 × 10⁶ CNY, followed by the NDS (CNY 10,960.84 × 10⁶), while the ESV under the EPS increases by CNY 9373.74 × 10⁶. The results of the scenario simulation showed that there was significant variability in ESV under different political orientations. (3) Hotspot analysis indicated that the ESV changes were clustered in the northeastern part and the eastern coastal region of the BTH. On this basis, we identified Chengde, Beijing, Tianjin, and Zhangjiakou as the key cities to focus on and made meaningful suggestions for their future regional environmental protection and sustainable development. This research can serve as a guide in creating sustainable BTH development policies and offer fresh perspectives for investigating how land use patterns affect the ecological environment’s regional quality under various policy trajectories.

1. Introduction

In the 21st century, urban areas have been expanding at an unprecedented rate, significantly impacting natural ecosystems. Ecosystem service (ES) refers to the benefits that humans derive from these ecosystems, which are closely linked to both human well-being and social development [1,2]. Research indicates that over the past fifty years, at least two-thirds of ecosystems have been damaged by human-induced land use changes (LUCs), leading to considerable reductions in ES and significant reduction in biodiversity [2,3]. According to Costanza’s research, the total value of global ecosystems was believed to have declined by USD twenty trillion per year from 1997 to 2011 [4]. This critical situation has compelled society to acknowledge that protecting ecosystems and the services they provide in sustaining human livelihoods is urgent and cannot be postponed [2,5,6]. Concurrently, the demand for reliable information on ES from decision-makers has steadily increased, making the measurement and monitoring of ES status and trends a necessary step to address this need. In addition, assessing ecosystem service functions is essential not only for clarifying regional ecological conditions and maintaining ecological stability [7] but also for optimizing land use patterns [8]. This assessment can raise people’s environmental awareness and underscore the importance of safeguarding natural capital [9,10].

Currently, two mainstream approaches are employed to access ESs: the ecological model and the equivalent coefficient method [11]. Ecological models typically require extensive input data and involve complex computational processes, limiting their application to small-scale ecosystem service calculations [12,13]. Furthermore, these models can only project one or a few types of ES, so they are also unable to fully predict the complete ecosystem service value (ESV) in a region [14,15]. In contrast, the equivalent coefficient method relies solely on land-use/land-cover (LULC) data, facilitating an intuitive assessment of various ESVs and enhancing its applicability for large-scale evaluations [4,16,17]. This method essentially transfers economic benefits from one homogeneous location to another to calculate the regional ESV [18,19]. In 1997, Costanza et al. [1] were the first to apply the equivalent coefficient method to evaluate the total ecosystem service value of global ecosystems, and relevant research has also attracted the attention of scholars in related disciplines [20]. Consequently, this approach has gained popularity in global ESV assessments [4,21,22,23,24]. In response to the specific conditions in China, Xie et al. [24,25,26] developed and updated the equivalence factor table to enhance its relevance for ESV assessments within the country. In 2012, Zhu [27] proposed the inclusion of Gross Ecosystem Product (GEP)—which represents the total value of final goods and services provided by ecosystems—in China’s sustainable development assessment and accounting framework. The introduction of this concept and the launch of the project have once again promoted the development of ecosystem service valuation in China [28]. To date, studies in this area have encompassed comprehensive evaluations across various scales, including global [4], regional [29], watershed [30,31], and specific ecosystems [26]. Each of these scales employs distinct valuation methodologies tailored to the unique characteristics of the ecosystems being studied.

LULC is among the most significant drivers influencing ESs [32,33,34]. The type and extent of LULC can directly affect the type and intensity of ESs [35,36,37], and substantial alterations in LULC have become a primary cause of the degradation and impairment of regional ecosystem functions [38]. Beginning in the 1980s, a series of reforms led to a qualitative leap in China’s economy and a significant increase in the rate of urbanization. The great transformation has exerted enormous pressure on ecosystems, leading to the severe degradation of ecosystems [39,40,41]. This transformation has spurred academic interest in the transitions of land use types and their effects on ecosystems [40]. However, existing research predominantly focuses on the relationship between LUC and ESs, complicating efforts to ascertain the precise amounts and spatial distributions of ESV associated with LUC. In addition, current studies tend to emphasize the implications of historical and existing land use patterns for ecosystems, which, while valuable for assessing past land use policies, often overlook the impacts of regional LUCs on ESV under various future scenarios. This limitation reduces the practical relevance of findings for territorial development planning. In response, many scholars have developed sophisticated spatial models for predicting future land use [42,43,44,45]. These models aid in identifying optimal resource protection strategies by simulating diverse land use scenarios and quantifying the resultant ecosystem changes, thereby providing a foundation for optimizing the spatial configuration of national territory and facilitating rational resource allocation [38].

The Beijing-Tianjin-Hebei urban agglomeration (BTH) is the central growth pole of China and one of the three major urban clusters in China. This area faces significant pressures from the ongoing expansion of built-up land, the reduction of cropland, and the accelerated degradation of ecological spaces. Consequently, addressing the resulting environmental degradation has emerged as a critical issue that threatens the sustainable development of the region. Although some research has been conducted on ecosystem service and ESV changes, few studies have quantified the contributions of various LUCs to ESV in BTH. In addition, earlier studies mainly spotlighted the historical impact of land use on ESV changes, but the way in which ESV modifications adapt to altered future land use patterns has been neglected. Consequently, the equivalent coefficient approach was employed in this study to determine the variations in land use and ESV in BTH across time. Furthermore, we created three scenarios to investigate how the ESV reacts to LUCs in diverse policy orientations. The paper addresses the following issues: (1) What are the spatial and temporal distributions and variations of ESV in urban agglomerations? (2) How do LUCs affect the ESV under different scenarios? (3) Which areas should be prioritized in future regional transformations, and what measures ought to be taken to protect the ecosystem?

2. Materials and Methods

2.1. Study Area

BTH is located in the core area around Bohai Bay, at 30°05′–42°37′ N, 113°11′–119°45′ E. The total land area is approximately 216,485.2 km² and comprises 11 prefecture-level cities in Hebei Province and two municipalities directly under the control of the central government, Beijing and Tianjin (Figure 1). The region consists of the Bashang Plateau (the southern edge of the Mongolian Plateau), the Taihang Mountains and Yanshan Mountains, the northern part of the North China Plain, and the coastal plain of the Bohai Sea. BTH has a temperate continental monsoon climate, with four distinct seasons throughout the year. Spring is cold, dry, and windy; summer is hot and rainy; autumn is clear and cool; and winter is cold with little rain or snow. The average annual temperature in the region ranges from 2.7 to 14.7 °C, while the average annual precipitation is 370–730 mm. The predominant vegetation types include deciduous forests and temperate grasslands.

BTH is the largest and most dynamic economic area in northern China, attracting increasing attention from both China and the entire world. After the adoption of the Outline of Coordinated Development of the Beijing-Tianjin-Hebei Region, BTH’s development was elevated to the pinnacle of national strategy. Since then, a number of policies and measures have been introduced one after the other to provide strong support and a guarantee for the region’s development. The plan’s thorough execution has resulted in rapid growth in BTH’s overall economic output, ongoing improvements to the transportation system, a number of ecological restoration initiatives that have successfully raised the quality of the area’s ecological environment, and the start of a new era in regional development. In 2022, regional development entered a new phase, with a total economic output of over CNY 10 trillion and a population of over 100 million. However, owing to speedy economic development and dramatic population expansion, large-scale city building continues to be carried out, which has certain pessimistic impacts on the regional ecological environment, such as water scarcity [46,47] and air pollution [48], soil erosion, and desertification [49]. The population of the BTH region will show a growing trend in the short term and this will raise demand for urban building sites. There will continue to be a shortage of ecological land, which may lead to more serious environmental problems in the region and not only threaten regional ecological stability but also significantly limit comprehensive and balanced development. By studying the future ecological environment of BTH, regional sustainable development can be supported.

2.2. Framework and Data Sources

The research projected future land use patterns using historical data for 2000, 2010, and 2020 and created a methodology to look into how LUC affected ESV under various scenarios (Figure 2). This paper involved four types of data, namely land use data, socio-economic factors, and climate and environmental data (

Table 1. Data and sources.

Categories	Name	Source
Land use data	Land use/land cover data	RESDC (https://www.resdc.cn/ (accessed on 1 October 2023))
Climatic and environmental data	Precipitation
	Temperature
	Soil texture
	DEM	USGS (https://www.usgs.gov (accessed on 1 October 2023))
	Soil organic carbon content (SOC)	ISRIC Data Hub (https://www.isric.org/explore/isric-soil-data-hub (accessed on 1 October 2023))
	Slope	Computed using the ArcGIS 10.4 spatial analysis and hydrological analysis modules.
	Watershed
Socioeconomic data	GDP	RESDC (https://www.resdc.cn/ (accessedon 1 October 2023))
	Population
	Administrative center	NGCC (http://www.ngcc.cn/ (accessed on 1 October 2023))
	Road
	Railroads
	Crop production	Derived from the statistical yearbooks of BTH.
	Unit price of Crop production

). In particular, land use data were classified into six main categories based on the actual situation of BTH: cropland, forestland, grassland, water, built-up land, and unused land. Otherwise, we adjusted all data to 500 m × 500 m spatial resolution.

2.3. Methods

2.3.1. Future Land Use Trajectory Prediction

Model Principles and Scenarios Setting

The future land use simulation (FLUS) model was improved and developed by Li [50] and is based on the traditional cellular automatic (CA) principle. The idea is to calculate the development possibilities of every class in a region by combining the development probability with the domain influence factor, the adaptive inertia coefficient, and the conversion costs. This is undertaken by using the artificial neural network (ANN) algorithm to calculate the baseline era’s land use data and the information for each driving component. By using the roulette wheel competition mechanism, the simulation result is ultimately obtained after determining the cell’s overall conversion probability.

In particular, this model comprises three blocks: (1) an ANN, with its ability to map intricate non-linear correlations between historical land use and diverse causes, which has been used to determine the likelihood of suitability of various types of land [36,51]; (2) a discrete and abstract computational system with self-adaptive inertia and a competitive principle based on CA, which makes it possible to accurately predict the long-term spatial pathways of LULC. It may also be utilized to evaluate the particular relationships between various land use types, which are typically represented as interactions and competition, and (3) a Markov model, a stochastic model used to simulate pseudo-random systems and forecast demand in the future using past data [52].

Referring to related studies and considering the existing development situation, the study designed three scenarios: the natural development scenario (NDS), the food security scenario (FSS), and the ecological priority scenario (EPS). The scenarios were developed in accordance with BTH’s social development plan, following the three strategic objectives of economic development, food security, and ecological protection. Under the NDS, the land use changes were exactly the same as in the past, and the encroachment of built-up land into cropland was strictly controlled under the FSS, while under the EPS, nature reserves and ecological reserves in the study area were set as prohibited areas for the flow of land, and the tendency of the flow of ecological land to non-ecological land was also limited. Combined with the goal of China’s ecological civilization construction, the study projected how land will be used in BTH in 2030, 2040, and 2050 over a period of ten years.

Parameter Setting

(1)

Conversion constraint matrices

The model makes use of a conversion constraint matrix to indicate whether a certain kind of land usage can be converted to another [53]. Values of 1 and 0 indicate whether or not this kind of land can be converted. Referring to the previous study, we set the transformation constraint matrix of BTH in three scenarios (Figure 3).

(2)

Effect of neighborhood

The ease with which one type of land usage can be replaced by another within the neighborhood is known as the neighborhood effect. It represents the capacity to expand, which is mainly due to external factors, and is usually expressed in terms of a range of critical values (ranging between 0 and 1) [54]. The larger the value, the more expandable the terrain type is. Referring to the relevant literature [54,55,56] and expert knowledge, after many experiments, the neighborhood effects were obtained in

Table 2. The neighborhood effect.

	Cropland	Forestland	Grassland	Water	Built-Up Land	Unused Land
NDS	0.3	0.6	0.5	0.4	0.7	0.1
FSS	0.5	0.5	0.5	0.2	0.6	0.1
EPS	0.2	0.8	0.7	0.5	0.6	0.1

.

(3)

Model validation

Based on the changing laws of land use patterns between 2000 and 2010, this research simulated the land use of BTH in 2020 and compared it to the land use as it actually occurred in 2020. To verify the correctness, the paper adopted overall accuracy, the Kappa coefficient, and the Figure of Merit (FOM) coefficient for verification. In this paper, the overall accuracy, Kappa coefficient, and FOM of the analog result are 0.881, 0.841, and 0.3958, respectively. The precision of the FLUS model in land use trajectory prediction experiments has been verified and it is equally well suited in this experiment.

2.3.2. Assessment of Ecosystem Services

This paper is based on the equivalent factor table of service value per unit area of the Chinese terrestrial ecosystem prepared by Xie [26]. In considering the realities of BTH, five land types other than built-up land were estimated, taking into account the convenience of calculation (

Table 3. Ecosystem service value per unit area in BTH (unit: CNY/hectare/year).

Service Type		Coefficient Value (CNY ha−1 yr−1)
		Cropland	Forestland	Grassland	Water	Built-Up Land	Unused Land
Provisioning services	Grain production	1.11	0.23	0.23	0.66	0.00	0.01
	Production of material	0.25	0.54	0.34	0.37	0.00	0.02
	Water resource supply	−1.31	0.28	0.19	5.44	0.00	0.01
Regulating services	Gas regulation	0.89	1.76	1.21	1.34	0.00	0.07
	Climate regulation	0.47	5.27	3.19	2.95	0.00	0.05
	Cleaning the environment	0.14	1.57	1.05	4.58	0.00	0.21
	Hydrological regulation	1.50	3.81	2.34	63.24	0.00	0.12
Supporting services	Soil conservation	0.52	2.14	1.47	1.62	0.00	0.08
	Nutrient cycling	0.16	0.16	0.11	0.13	0.00	0.01
	Biodiversity	0.17	1.95	1.34	5.21	0.00	0.07
Cultural services	Aesthetic landscape	0.08	0.86	0.59	3.31	0.00	0.03
Sum		3.95	18.57	12.06	88.82	0.00	0.65

). We calculated the ESV of the research area in historical periods (2000, 2010, 2020) and simulation periods (2030, 2040, 2050) under three scenarios. Considering the size of the BTH, the fishing net tool was used to create a fishnet with a side length of 1 km in ArcGIS 10.4 (218,253 grids in total), and the ESV corresponding to each category in each grid was calculated separately. The process of calculation is as follows:

$${\mathrm{V}\mathrm{C}}_{\mathrm{i}}={\mathsf{\Sigma }}_{\mathrm{j}=1}^{\mathrm{m}}{\mathrm{E}\mathrm{C}}_{\mathrm{i}\mathrm{j}}{\times \mathrm{E}}_{\mathrm{a}}$$

(1)

$$\mathrm{E}\mathrm{S}\mathrm{V}={\mathsf{\Sigma }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{A}}_{\mathrm{i}}{\times \mathrm{V}\mathrm{C}}_{\mathrm{i}}$$

(2)

where ESV indicates the ecosystem service value in CNY/year, VC_i is the ESV per unit size of land use type i in CNY/hectare/year, EC_ij is the equivalent value of ESs type j for land use type i, m is the quantity of ES categories, E_a is the economic value of one unit of ESs, expressed in CNY/hectare/year, and Ai is the size of land use type i in hm².

A standard equivalent factor for ESV is described as the economic value of the annual natural grain yield of one hectare of the national average wheat production in cropland [4,57]. It is generally assumed that, calculated at market prices, it corresponds to one-seventh of the economic value of the grain yielded per hectare of cropland. In order to better match the research results with the real situation in BTH, this study selected the annual yields of the main grain crops wheat, rice, and corn in the research area over the past 20 years to calculate the grain yield situation. To counteract the effects of food price volatility, national food price data for 2020 were used. The result of this calculation is that the economic value for every piece of ESs of BTH is CNY 2037.06/hectare/year. The concrete calculation process is as follows:

$${\mathrm{E}}_{\mathrm{a}}={\displaystyle \frac{1}{7}}{\mathsf{\Sigma }}_{\mathrm{i}=1}^3{\displaystyle \frac{{\mathrm{m}}_{\mathrm{i}}{\mathrm{p}}_{\mathrm{i}}{\mathrm{q}}_{\mathrm{i}}}{\mathrm{M}}}$$

(3)

where E_a represents the economic value for every piece of ESs, expressed in CNY/hectare/year, i indicates the three cereals, m_i is the average price of cereals, expressed in CNY/kg, p_i is the farming area of cereals i, q_i is the yield of cereals i, expressed in kg/hm², expressed in hm², and M is the entire cultivated area of cereals i, expressed in hm².

2.3.3. Determine Hotspots and Cold Spots of Change in ESV

Hotspot analysis is a commonly used tool for identifying the space orientation of cold spots and hotspots. In this paper, using this tool, the spatial cluster analysis of highly valuable metropolitan areas (hotspots) and poor-value metropolitan areas (cold spots) was carried out based on the changes in ESV of the BTH [58,59]. The accumulation of changes in ESV will identify areas where research needs to be focused and where more appropriate solutions need to be applied. The study used area statistical tools to add ESV changes per hectare to the attribute table of fishing nets. The concrete calculation process is as follows:

$${\mathrm{E}}_0={\mathrm{E}}_{\mathrm{m}1}-{\mathrm{E}}_{\mathrm{m}2}$$

(4)

$${\mathsf{\Delta }\mathrm{E}\mathrm{S}\mathrm{V}=\mathsf{\Sigma }}_{\mathrm{k}=1}^{23}{\mathrm{E}}_0{\mathrm{A}}_0$$

(5)

where E₀ indicates the influence coefficient of land use type conversion on the ESV per unit area, E_m1 represents the ESV per unit of the converted land use, and E_m2 represents the ESV per unit area before transformation. ΔESV is the extent of change in ESV caused by each modification under different scenarios and A₀ refers to the area of transformation types, with a total of 23 transformation types.

$${\mathrm{G}\mathrm{i}}^{\mathrm{*}}={\displaystyle \frac{{\mathsf{\Sigma }}_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i},\mathrm{j}}{\mathrm{x}}_{\mathrm{j}}-\mathrm{X}{\mathsf{\Sigma }}_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i},\mathrm{j}}}{\sqrt[\mathrm{s}]{\left[\frac{{\mathrm{n}{\mathsf{\Sigma }}_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i},\mathrm{j}}^2-\left({\mathsf{\Sigma }}_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i},\mathrm{j}}\right)}^2}{\left(\mathrm{n}-1\right)}\right]}}}$$

(6)

$$\mathrm{X}={\displaystyle \frac{1}{\mathrm{n}}}{\mathsf{\Sigma }}_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{x}}_{\mathrm{i}}$$

(7)

$$\mathrm{S}=\sqrt{{\displaystyle \frac{1}{\mathrm{n}}}}{\mathsf{\Sigma }}_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{x}}_{\mathrm{j}}^2-{\left(\mathrm{X}\right)}^2$$

(8)

where n is the quantities of spatial grid cells, x_i and x_j are the observed values of cells i and j, respectively, and w_ij is the spatial weight matrix constructed based on the spatial k-adjacency relationship.

3. Results

3.1. Spatiotemporal Evolution of LULC

Through statistics and analyses of regional LULC between 2000 to 2020, we found that the major land use type in BTH was cropland, accounting for more than half of the entire region. However, it showed a continuous downward trend (Figure 4). Over the past 20 years, it has decreased by 11,932.25 km², and the dynamic land use index was 0.47%. The grassland area and water area continued to decline, decreasing by 2436.75 km² and 372.50 km², respectively, between 2000 and 2020. However, the percentage of forestland and built-up land kept rising. The individual dynamic settings were 0.38% and 10.69%, respectively. The total growth in built-up land area was 12,751 km², from 2.80% in 2000 to 8.79% in 2020, but the growth rate showed a first fast and then slow trend. The proportion of unused land was lowest within 80 km² at less than 1%.

Through two land use transfer maps, we have more clearly delineated the features and regularity variations of the spatial distribution of LUCs. Between 2000 and 2020, the entire area of LUCs was 27,365.75 km², with notable changes occurring among cropland, built-up land, and grassland (Figure 5). The three LUCs with the largest area were the conversion from cropland to built-up land (10,275.25 km²), cropland to grassland (4587.5 km²), and grassland to cropland (3419 km²). The analysis showed that cropland and grassland have the highest proportion of erosion. The total share of transfer out of both of them is more than 90% of the total removal amount (cropland accounted for 60.91% of the total removal amount and the share of grassland was 29.29%). The largest amount transferred into is built-up land, which comprised 47.61% of the overall amount of transferred land; moreover, the vast majority of the transferred area of water and unused land was also transferred into built-up land. Overall, incoming and outgoing flows of cropland were the major forms of land transfers. Cropland was primarily sacrificed in favor of the expansion of built-up land, but the total area of cropland remained relatively stable. This finding also corroborated the policy of requiring balanced occupation of and compensation for cropland.

Regarding the arrangement of spaces, the primary location with significant LUCs is the eastern portion of the BTH urban region (Figure 6). From 2000 to 2010, the most dramatic changes happened in the metropolitan areas of Beijing, Tianjin, and Tangshan. The main changes in these cities center on the transformation of cropland. There is a clear downward trend in the intensity of the shift in the last decade compared to the first decade.

3.2. Change of ESV in BTH

Between 2000 and 2020, BTH’s ESV declined from CNY 386,859.89 × 10⁶ to CNY 371,968.78 × 10⁶ (Figure 7a). Particularly, from 2000 to 2010, the ESV change was considerable, with a total decline of CNY 10,971.72 × 10⁶. After 2010, the ESV decreased slightly, which is inseparable from the various environmental protection and ecological restoration measures continuously carried out by the government for the past few years.

Regarding the portion of a single land type that contributes to the total ESV, grassland accounts for over 30%, the highest contribution of any one land type, followed by cropland and forestland (Figure 7a). Otherwise, from the viewpoint of one particular service, the hydrological regulations and the climate regulations are the largest contributors, which is in line with the trend of changes in the total ESV, taking into account all secondary categories except water supply, which shows a decreasing trend (Figure 7b). Specifically, water supply was negative, which is due to the status of land use in BTH with more cropland and fewer water bodies. However, over the last two decades, the water area has increased and the regional water supply capacity has grown, which has increased the ESV offered. These findings showed that the prime reasons behind ESV changes in BTH were support and regulatory services, and their share is even more than 70% of the total.

3.3. Scenario Simulation

3.3.1. Model Simulation and Analysis

Figure 8 illustrates land use in BTH across three scenarios during the period 2030–2050, revealing similar spatial patterns. Cropland remained the paramount land use type, widely distributed across eastern and southern BTH, while forestland and grassland were separated into the north. Under the NDS, a significant expansion of built-up land was evident, primarily at the expense of cropland in the southern plains. In contrast with 2020, the fragmentation of cropland under the NDS continued to increase. In contrast, land use patterns under the FSS have not changed significantly since 2020. The trend of urban sprawl taking up cropland has been more curbed in the last 20 years, mitigating the downside of urban sprawl resulting from urbanization on the issue of farmland protection. LUCs under the EPS concentrate on forestland and grassland in the north. Overall, there was an increase in forestland and a greater concentration of grassland in the north of BTH.

For the purpose of examining the quantitative characteristics of various land use types in three situations, their area changes were computed (Figure 9). In contrast to the 115,269.5 km² in 2020, the cropland under the NDS decreased slightly to 100,550 km² during 2030–2050, representing an obvious declining trend under the EPS, with a decrease of 14,996.25 km² (

Table 4. The area of each land use type under three different scenarios (unit: km²).

		Cropland	Forestland	Grassland	Water	Built-Up Land	Unused Land
	2020	115,269.5	27,790.75	48,035.75	3100.25	18,717.25	78
NDS	2030	110,036.75	28,705.75	47,564	3015.75	23,602.25	67
	2040	105,135	29,535.5	47,059	2923.75	28,271.75	66.5
	2050	100,550	30,274.25	46,527	2834.75	32,739.5	66
FSS	2030	116,587.5	27,492.75	47,022	2963.75	18,849	76.5
	2040	118,252.5	27,254.25	45,561.75	2815.5	19,031.75	75.75
	2050	119,945.5	27,114	43,978	2668	19,212.25	73.75
EPS	2030	110,036.75	28,705.75	49,500.75	3128	21,553	67.25
	2040	105,041.5	29,610.5	50,868	3144.5	24,260	67
	2050	100,273.25	30,505.25	52,142	3160.25	26,844	66.75

). By comparison, cropland under the FSS expanded greatly from 116,587.5 km² to 119,945.5 km², an increase of 3358 km². As shown in Figure 8, the trend of occupation of cropland by built-up land has been effectively curbed under the FSS, and the amount of cropland has grown steadily and maintained its original spatial shape. The development of forestland in the period 2030–2050 is opposite to that of cropland. The forestland area increased by 2714.50 km² and 2483.50 km² under the EPS and NDS, respectively, but decreased by 378.75 km² under the FSS. Under the EPS, there was a trend of forestland expanding northward along the lines of the Yanshan Mountains and Taihang Mountains.

Compared to 18,717.25 km² in 2020, built-up land increased to different extents in different scenarios from 2030 to 2050. Under the NDS and ESS, built-up land expanded at rates of 467.41 km²/a and 270.89 km²/a, respectively. By contrast, under the FSS, because the expanse of built-up land encroaching on cropland was restricted, built-up land increased by only 495 km² between 2020 and 2050. The sprawl of built-up land under the NDS and FSS focused on the city construction areas in Beijing, Shijiazhuang, Xingtai, and Handan. Grassland and water showed the same change trend, decreasing to varying degrees under both the FSS and NDS and increasing under the EPS. Specifically, the grassland area increased from 48,035.75 km² in 2020 to 52,142 km² in 2050 under the EPS, a rise of 4106.25 km², but under the FSS it decreased by 4057.75 km² and by 1508.75 km² under the NDS. The number of water bodies decreased the most under the FSS, which is consistent with the strategy of cropland protection. The variations of unused land were insignificant owing to the small percentage (only 0.02%).

3.3.2. Change of ESV Under Multi-Scenarios in BTH

As expected, the total ESV under the EPS would increase significantly from 2030 to 2050, increasing by CNY 371,968.78 × 10⁶ to CNY 381,342.52 × 10⁶ (

Table 5. Ecosystem service values for the different service types under three scenarios, negative values represent decreases (unit: CNY 10⁹).

Services	FSS			EPS			NDS
	2030	2040	2050	2030	2040	2050	2030	2040	2050
Grain production	30.07	30.21	30.35	28.85	27.83	26.86	28.74	27.64	26.61
Production of material	12.33	12.30	12.27	12.31	12.26	12.21	12.17	11.97	11.78
Water resource supply	−24.04	−24.30	−24.55	−22.23	−20.78	−19.40	−22.43	−21.20	−20.06
Gas regulation	43.44	43.34	43.25	43.29	43.05	42.83	42.78	42.04	41.32
Climate regulation	73.16	72.55	71.96	75.28	76.68	78.02	73.95	74.00	73.96
Cleaning the environment	25.02	24.82	24.63	25.71	26.17	26.61	25.19	25.13	25.04
Hydrological regulation	119.22	118.87	118.52	119.67	119.72	119.77	117.30	115.03	112.80
Soil conservation	39.50	39.30	39.09	23.82	40.30	40.58	39.41	39.07	38.72
Nutrient cycling	5.70	5.71	5.71	5.60	5.50	5.41	5.55	5.41	5.28
Biodiversity	31.17	30.93	30.69	32.04	32.62	33.18	31.40	31.32	31.22
Aesthetic landscape	14.38	14.27	14.16	14.77	15.03	15.27	14.46	14.41	14.34
Grain production	30.07	30.21	30.35	28.85	27.83	26.86	28.74	27.64	26.61
Production of material	12.33	12.30	12.27	12.31	12.26	12.21	12.17	11.97	11.78

). In contrast, the total ESV under the FSS would decrease the highest, with a decline of CNY 16,568.78 × 10⁶ compared to 2020. The total ESV under the NDS would be between the EPS and NDS, so the ESV would be effectively regulated in different simulation scenarios.

Of the several kinds of ESs, hydrological regulation and climate regulation services supported a significantly higher ESV, representing more than half of the overall ESV in each scenario (Figure 10). Most ESVs declined from 2020 to 2050 under the NDS and FSS, except for grain production and nutrient cycling under the FSS and water resource supply and climate regulation under the NDS. Under the EPS, grain production, production of material, gas regulation, and nutrient cycling all declined to varying degrees, but soil conservation, water resource supply, and climate regulation all increased to a greater extent, thereby raising the overall ESV.

To provide greater clarity on the spatial heterogeneity of the ESV changes under various scenarios, this paper conducted a hotspot analysis of changes in the ESV. From Figure 11, the significant region of ESV changes was mainly focused on the northeast and east coast of BTH. Under the NDS, the proportion of both hotspots and cold spots showed an increasing trend (the proportion of hotspots increased from 2.70% to 4.14% and the proportion of cold spots increased from 2.76% to 3.11%), while the FSS scenario showed that the opposite is the case (the share of hotspots decreased from 8.88% to 3.72% and the share of cold spots decreased from 16.74% to 3.33%), suggesting that the ESV has declined most significantly. The EPS scenario showed a trend where the proportion of hotspots and cold spots decreases and then increases, with the overall change being smaller.

4. Discussion

4.1. ESV Response to LUCs

Different land use types produced tangibly different ESVs. In BTH, grassland and forestland made up over 50% of the total ESV. Between 2000 and 2020, the ESV from forestland increased from CNY 9,769,511 × 10⁶ to CNY 10,512,457 × 10⁶ (from 25.25% to 28.26%), while the ESV from cropland decreased from CNY 10,234,928 × 10⁶ to CNY 9,274,852 × 10⁶ (from 26.46% to 24.93%). The changes in the total ESV from 2020 to 2050 were: a decrease of CNY 10,959.27 × 10⁶ under the NDS, an increase of CNY 9373.74 × 10⁶ under the EPS, and a decrease of CNY 16,568.48 × 10⁶ under the FSS. This shows that LUCs greatly decreased the ESV under the NDS and FSS, but increased it sharply under the EPS. It also shows that differences in LUC can lead to variations in the rise or fall in ESV. Land use transformation from forestland to cropland, from forestland or grassland to built-up land, and from water to another land use type led to larger damage to the ESV, whereas transformation from cropland to forestland or water and from cropland, forestland, or built-up land to water greatly increased ESV (

Table 6. The VC change of each LUC type (unit: CNY/ha/yr).

LUC Types	ΔVC	LUC Types	ΔVC
Cropland→Forestland	29,781.82	Water→Cropland	−172,875.10
Cropland→Grassland	16,520.56	Water→Forestland	−143,093.28
Cropland→Water	172,875.10	Water→Grassland	−156,354.54
Cropland→Built-up land	−8046.39	Water→Built-up land	−180,921.48
Cropland→Unused land	−6722.30	Water→Unused land	−179,597.39
Forestland→Cropland	−29,781.82	Built-up land→Cropland	8046.39
Forestland→Grassland	−13,261.26	Built-up land→Forestland	37,828.20
Forestland→Water	143,093.28	Built-up land→Grassland	24,566.94
Forestland→Built-up land	−37,828.20	Built-up land→Water	180,921.48
Grassland→Cropland	−16,520.56	Built-up land→Unused land	1324.09
Grassland→Forestland	13,261.26	Unused land→Cropland	24,566.94
Grassland→Water	156,354.54	Unused land→Grassland	180,921.48
Grassland→Built-up land	−24,566.94	Unused land→Water	1324.09
Grassland→Unused land	−23,242.85	Unused land→Built-up land	1324.09

).

Under the NDS, the ESV reduced by CNY 19,566.83 × 10⁶ from 2030 to 2050, primarily owing to the transition of water and cropland to built-up land. Although the LUCs from cropland to forestland and grassland enhanced the ESV by CNY 19,280.49 × 10⁶, it was not enough to offset the decline caused by other conversions. Under the EPS, from 2020 to 2050, the ESV increased by CNY 18,029.57 × 10⁶, mainly due to the LUCs from cropland to forestland and grassland, while the LUC converting cropland to built-up land decreased by CNY 9379.07 × 10⁶. In contrast to the EPS and NDS, the shape of LUC in BTH from 2020 to 2050 under the FSS is more homogeneous, with a smaller area subject to transformation (only 3.31%). From Figure 12, according to FSS, the primary causes for the decline in ecosystem services are the transformation from forestland and grassland into cropland and the transformation from water into built-up land, which caused a loss of CNY 15,363.60 × 10⁶, accounting for 92.62% of the total loss of ESV. It is worth mentioning that although the water area converted into built-up land was 34,350 km² (only 4.86 of the total area changed), this resulted in a reduction in ESV of CNY 6214.65 × 10⁶. These results suggested that the transition from cropland to forestland and grassland caused a rise in ESV, such as in the EPS, whereas changes from alternative land types to built-up land or from forestland and grassland to cropland brought about a loss of ESV, such as in the NDS and FSS.

4.2. Territory Development Plan Based on Impacts of LUC on ESV

Over the years, research on ESV has been widely used in regional planning and construction [60,61], and scholars have selected appropriate land use patterns for regional development depending on how various land use scenarios affect ESV [62,63]. In this paper, we investigated the ESV in BTH based on simulated land use patterns from 2030 to 2050 and analyzed the impact of the projected changes on the ESV. To assist decision-makers in making more rational logical land use decisions, our research explores the following aspects, which provide important information for regional territorial spatial planning:

Differences in the spatial distribution of hotspots and cold spots can provide a scientific basis for establishing environmental management in accordance with local conditions, so the study further analyzed the hotspot and cold spot areas in each city. The results revealed the regions where land use planners need to pay particular attention, such as Chengde, Zhangjiakou, Beijing, and Tianjin. Among them, Chengde occupied the largest area of hotspots and cold spots, which need to be monitored and controlled (Figure 13).

The types of LUCs in the areas that need to be focused on indicate the specific causes of changes in ESV. We examined the various forms of LUC in these regions and found that the paramount reason for the ESV decline was the spread of built-up land onto ecological land including cropland, forestland, and grassland (Figure 14a). The major form of LUC that resulted in a rise in ESV of BTH was the transformation from cropland into forestland and grassland (Figure 14c). To promote coordinated development, the study proposed corresponding ecological and economic management policies and environmental protection countermeasures for different regions.

Chengde serves as a crucial ecological barrier within the study area, playing a vital role in water conservation, air purification, and ecological protection [64]. To enhance the functionality of the ecosystem and ensure it can provide more natural products and services to surrounding regions in the future, it is essential to strengthen environmental protection and maintenance efforts. As illustrated in Figure 11, the variations in cold spots and hotspots in Chengde across different scenarios exhibited significant differences. Under the NDS and EPS, the restoration of cropland to forestland was consistently implemented, resulting in a continuous increase in ecological land. In contrast, cropland encroached on ecological areas under the FSS scenario, which led to the emergence of extensive cold spots. We recommend the following: (1) Safeguard regional biodiversity by implementing stringent regulatory measures for ecological reserves within nature reserves. (2) Continue promoting afforestation and greening initiatives, expand greening projects on degraded mining sites, sustain the cultivation of existing plantation forests, and establish an ecological network of waterways and green corridors to create a regional ecological security barrier.

Beijing and Tianjin are densely populated, economically developed megacities characterized by significant human–land conflicts, particularly concerning water scarcity and air pollution, which serve as critical constraints on urban development [65,66]. Given the unique characteristics of the region, the population is projected to continue growing at a high rate in the coming years, leading to urban expansion that will inevitably encroach upon natural land and disrupt the stability of ecosystem service provision. Therefore, we recommend the following: (1) Strictly regulate the boundaries of urban expansion, taking into account regional ecological supply and demand when planning urban development, to ensure the continuity and integrity of remaining natural areas. (2) Enhance the capacity for regional ecosystem service provision by strengthening green infrastructure, establishing farmland protection forest networks, promoting modern agricultural practices, and reducing the demand for water resources.

From a regional perspective, we recommend the establishment of a compensation mechanism to clarify the roles of compensation providers and recipients based on the supply and demand of ecosystem services within BTH. The northern part of the region, characterized by a strong capacity to supply ecosystem services, should be designated as the recipient of ecological compensation. In contrast, Beijing, Tianjin, and the southern regions, which experience significant human–land conflicts, should contribute compensation funds as beneficiaries of ecological protection. Relevant government departments in the region must collaborate to implement environmental construction and restoration initiatives, thereby promoting coordinated regional development.

4.3. Limitations

This work analyzed the ESV and generated the trajectory of potential land uses. However, there were also some uncertainties and limitations. When projecting future land use patterns, the analysis only took into account the individual land mass to examine the ESV while ignoring other ESV drivers. For instance, the water flow in the water bodies of the study area may decrease with increasing forest cover due to increased evapotranspiration, which contradicts the equivalence coefficient table [67,68]. In the process of forecasting, due to a series of objective factors such as data limitations, the study neglected some factors that may affect the accuracy of long-term forecasts, such as climate change [69,70,71,72]. In future research, it will be important to fully consider the various factors affecting the accuracy of long-term forecasts, so as to improve the accuracy and credibility of the forecast results.

Moreover, this study used the equivalent coefficient method, recorded the yields of major grain crops and the acreage area of BTH from 2000 to 2020 to calculate the grain production situation, and utilized national grain price data in 2020 to remove the impact of variations in grain prices. However, due to the limitations in the data collection process, the study calculated the worth of every ES category with a uniform standard, and the calculated results have the problem of neglecting regional spatial differences, so the results were analyzed and discussed mainly from a static point of view.

5. Conclusions

A long history of rapid industrialization and urbanization has had serious impacts on BTH’s ecological environment. Given the objective requirements of an expanding population and the advancement of the economy, policymakers urgently need to integrate regional economic development needs and ecological protection strategies to formulate a reasonable sustainable development plan. In this study, BTH was used as an example to successfully simulate and evaluate the ESV using the FLUS model and the equivalent coefficient method. The following conclusions were drawn from the study:

(1)

Cropland was the primary land use category in BTH between 2000 and 2020, and the prominent features of regional LUCs were the decline in cropland (decrease by 11,932.25 km²) and the increase in built-up land (addition of 12,750.75 km²). The biggest contributors to ESV were grassland and forestland, and there has been a noticeable decrease in regional ESV (total CNY 14,634.22 × 10⁶).

(2)

The study successfully simulated the land use patterns according to the NDS, FSS, and EPS in 2030, 2040, and 2050, and measured their ecological impacts. Among them, the ESV lost under the EPS from 2020 to 2050 is CNY 16568.78 × 10⁶, the NDS was the second largest (loss of CNY 10960.84 × 10⁶), and the ESV under the EPS increased by CNY 9373.73 × 10⁶. This shows that the EPS was the optimal choice for the future development of BTH.

(3)

The form of BTH land use transformation varies in different scenarios, but the main changes were concentrated in Zhangjiakou, Chengde, Beijing, and Tianjin, which were the areas to focus on for future urban development. Meanwhile, it has been discovered that the regional ESV will decrease as ecological land is transferred to built-up territory, while returning cropland to forestland and grassland will greatly enhance the natural environment’s quality, which will aid in achieving sustainable development in the area.