Multi-Scenario Simulation of Ecosystem Service Value in Wuhan Metropolitan Area Based on PLUS-GMOP Model

Abstract

Rapid construction and development activities in large urban areas have significantly impacted land use and land cover (LULC). They have brought great pressure to urban sustainable development. Current studies have shown that changes in LULC structure significantly affect regional ecosystem service functions but lack the sufficient scientific basis to provide reasonable strategies for the future development of urban areas. Based on land use and related data for the Wuhan metropolitan area (WMA) in 2000, 2010, and 2020, in this study, we construct a land use and ecosystem service value (ESV) simulation method based on a coupled PLUS-GMOP (patch generation land use simulation and grey multi-objective optimization) model and find that the changes in LULC structure from 2000 to 2020 are mainly reflected in the decrease in farmland and water area and the increase in built-up land, which are spatially reflected in Wuhan city center and other surrounding urban centers. The ESV also exhibits a slight increase and then a significant decrease, and a consistent overall pattern of high in the west and low in the east. By presupposing three scenarios for 2030 (ND, natural development; EFD, ecological first development; EECD, ecological and economic coordinated development), the analysis shows that although the ecological service value is still decreasing, the EECD scenario achieves a relatively high economic value (+90.134 billion yuan) by losing less ecological service value (0.27 million yuan) than EFD, which is the development model advocated in this study. The PLUS-GMOP coupling model proposed in this study provides a scientific reference for coordinating regional economic development and ecological protection in large cities, and provides a new technical path for metropolitan area sustainable development and planning.

1. Introduction

The relationship between ecosystem services and human wellbeing is very close [1,2]. Frequent human activity continues to cause significant changes in LULC (land use and land cover) patterns that affect the supply and distribution of ecosystem services, and imbalances in the LULC structure severely limit ecosystem service value (ESV) development [3,4]. Therefore, the LULC structure needs to be optimized to achieve optimal allocation of land resources in different sectors, divisions, and regions [5]. Metropolitan areas are the most economically active, open, and densely populated regions of a country [6,7,8]. Many countries have made metropolitan areas the centers of national development [9]; for example, China proposed building about 30 modernized metropolitan areas in 2019 to participate in global competition [10,11]. ESV is a more mature index for judging regional ecological quality at present, and is of great significance for optimizing land use in metropolitan areas [12].

There is a strong link between land use and ESV [13,14,15]. However, to understand the specific impact relationship between LULC and ESV, it is necessary to simulate and predict this relationship. For example, some studies use multi-objective programming or system dynamics (SD) to optimize the LULC structure [16,17], but these methods do not adequately represent the spatial variation of LULC and ESV [18]. The optimal allocation methods based on model construction mainly include the multi-agent system (MAS), the minimum cumulative resistance (MCR) model, cellular automata (CA), or the Markov model [19,20]. These models lack accuracy and scientific validity for the simulation of future sites [21]. In this study we combine the advantages of grey multi-objective optimization (GMOP) and patch generation land use simulation (PLUS) models in quantitative and spatial optimization, and establish a coupled model based on both [22]. On the one hand, PLUS is a new model used to improve the accuracy of regional land use simulation, and on the other hand, it reveals the relationships between land use change and various driving factors [23]. The use of a PLUS-GMOP coupling model to study land use change and ecological service value can exploit the model’s accuracy advantages in terms of area quantity and spatial placement, and improve the scientific value of this type of future scenario simulation research. The optimization method proposed in this study can help to explore the changes in LULC and ESV in urban areas, and can provide scientific guidance and advice on ecological protection and urban planning in metropolitan areas.

2. Study Area

WMA (Wuhan metropolitan area) is located in the eastern part of Hubei Province (112°29′–116°08′ E, 31°49′–33°20′ N) in China. The metropolitan area is centered on the provincial capital city of Wuhan, with a radius of about 100 km, roughly exhibiting a radial extension, and including Wuhan, Huangshi, Qianjiang Xianning, Ezhou, Xiantao, Huanggang, Xiaogan, and Tianmen. The area is about 58,000 km², and the population was 31.86 million in 2020. Figure 1 shows the study area location. WMA belongs to the regional economic association, which is the most dense and dynamic area of industry and production in Hubei Province [24,25]. As the largest economic development association in central China, WMA is also a pilot area for the construction of “two types of society”, and the contradiction between economic and social development progress and ecological resource constraints is increasingly prominent.

3. Methodology

3.1. Data and Processing

The datasets studied in this paper comprise the following:

(1) WMA 2000, 2010, and 2020 land use data. According to the national land use classification standard and the research purpose of this paper, the land use types were reclassified into six categories: farmland, forest land, grassland, water, built-up land, and unused land;

(2) Climate and environmental data. The main natural factors driving land use change are topography, climate, precipitation, river distance, etc.;

(3) Socioeconomic data. This is divided into raster data that drive land use change and statistical data used to calculate the quantitative structure of multi-objective land use. Raster data include population distribution, GDP distribution, etc. The statistical information mainly includes grain output and crop area.

The specific data sources can be found in

Table 1. Data type and source.

Category	Data	Year	Original Resolution	Data Source
Land Use Data	LULC	2000, 2010, 2020	30 m	National Geomatics Center of China (http://www.globeland30.org/ accessed on 1 August 2021)
Climate and Environmental Data	Elevation	-	90 m	Geospatial Data Cloud (https://www.gscloud.cn/ accessed on 1 August 2021)
	Slope	-	90 m
	Average annual precipitation	2020	1000 m	CAS Resource and Environmental Science and Data Center (https://www.resdc.cn/ accessed on 1 August 2021)
	Average annual temperature	2020	1000 m
	River system	2020	-	National Earth System Science Data Center (http://www.geodata.cn/ accessed on 1 January 2021)
Socioeconomic Data	GDP	2020	1000 m	CAS Resource and Environmental Science and Data Center (https://www.resdc.cn/ accessed on 1 January 2021)
	Road network	2020	-	Open Street Map (https://www.openstreetmap.org/ accessed on 1 January 2021)
	Railroad stations	2020	-
	Government locations	2020	-
	Grain sown area and production	2010–2020	-	Government Statistical Yearbook (http://tjj.hubei.gov.cn/ accessed on 10 October 2021)
	Crop area	2010–2020	-

.

3.2. The Coupled PLUS-GMOP Model

Using WMA land use data and natural socioeconomic data for three years, in this paper we explore the LULC change pattern from 2010 to 2020, fit and optimize the land use demand using Markov and GMOP models, and simulate and predict the land use distribution and ESV distribution pattern of WMA in 2030 under ND (natural development), EFD (ecological first development), and EECD (ecological and economic coordinated development) scenarios using a PLUS model. Figure 2 shows the research framework.

3.2.1. Plaque Number Prediction Using GMOP Model

GMOP was combined and developed from the gray prediction theory of accurate linear delineation. Its main idea is to define the uncertain objective function and constraints in a specific determined range of satisfactory regions. The model can solve various uncertainties of the objective process and conditions in land use and solve multi-objective problems in land use structure optimization. The construction of the model involves four main aspects: decision variable selection, setting the objective function, determining the gray constraints, and determining the solution method. In this paper, we use three scenarios to predict the changes in the number of land use patches in WMA in 2030.

(1) Set decision variables

The decision variables are set based on the consideration of land use in WMA and the availability and operability of data. Six decision variables are selected: farmland, forest, grassland, water, built-up land, and unused land.

Table 2. Area of each land use type in WMA from 2000 to 2020 (km²).

Year	Farmland	Forest	Grassland	Water	Built-Up Land	Unused Land
2000	32,525.60	17,270.79	1482.60	5120.65	1673.91	32.25
2010	32,899.56	17,165.46	1417.07	4503.89	2068.98	47.43
2020	32,031.27	16,659.22	1323.99	3641.50	4358.19	91.74

shows the area of each type of site.

(2) Multiple-scenario setting and objective function

The land system is a complex system composed of natural, social, and economic subsystems, and its objective system includes economic development objectives, social benefit objectives, ecological protection objectives, and so on [26]. The objective function is established in two respects: economic benefits and ecological benefits.

ND: The ND scenario uses a Markov model directly to simulate and forecast based on historical land use data.

EFD: In the EFD scenario, the optimization objective is set to maximize the ecological benefits, where the ecological benefits are represented by the ESV indicator.

$$F_1\left(x\right)=max{\displaystyle \sum }_{i=1}^6{c}_i{x}_i$$

(1)

In formula (1), $F_1\left(x\right)$ is ESV (unit: ten thousand yuan), $c_i$ is the ESV per unit area of each decision variable (unit: ten thousand yuan/km²), and $x_i$ is the area of the decision variable (unit: km²).

We adopted the “Table of Ecological Service Value per Unit Area of China’s Terrestrial Ecosystem” constructed by Xie et al. [27], analyzed the relevant statistical yearbook data to obtain the grain economic output value per unit area from 2010 to 2020, and predicted the value of grain per unit area from 2010 to 2020. The ecosystem service value coefficient corresponding to the land use type in the circle in 2030 (

Table 3. Table of ecological benefits per unit area of different land use types in WMA (ten thousand yuan/km²).

First Class Type	Secondary Type	Farmland	Forest	Grassland	Water	Built-Up Land	Unused Land
Supply services	Food production	5.17	1.71	2.22	2.74	0.00	0.10
	Raw material	2.02	15.40	1.86	1.81	0.00	0.21
Conditioning services	Gas regulation	3.72	22.32	7.75	2.64	0.00	0.31
	Climate regulation	5.01	21.03	8.06	10.65	0.00	0.67
	Water flow regulation	3.98	21.14	7.86	97.00	0.00	0.36
	Waste treatment	7.18	8.89	6.82	76.74	0.00	1.34
Support services	Soil maintenance	7.60	20.77	11.58	2.12	0.00	0.88
	Biodiversity protection	5.27	23.31	9.66	17.73	0.00	2.07
cultural services	Aesthetic view supply	0.88	10.75	4.50	22.95	0.00	1.24
Total		40.83	145.32	60.31	234.36	0.00	7.18

), and finally the ecological benefit objective function of the WMA, were obtained as follows: $F_1\left(x\right)=$40.83$x_1$ + 145.32$x_2$ + 60.31$x_3$ + 234.36$x_4$ + 0$x_5$ + 7.18$x_6$.

EECD: The EECD scenario involves neither ecological maximization nor economic maximization, but aims to promote a coordinated and balanced development model of ecological protection and economic development. We constructed the objective function through the multi-objective decision-making of the GMOP model:

$$Max\left\{F_1\left(x\right),F_2\left(x\right)\right\}$$

(2)

$$F_2\left(x\right)=max{{\displaystyle \sum }}_{i=1}^6{d}_i{x}_i$$

(3)

In formula (3), $F_2\left(x\right)$ is the total economic benefit (unit: ten thousand yuan), $d_i$ is the economic benefit coefficient per unit area of each decision variable (unit: ten thousand yuan/km²), and $x_i$ is the area of the decision variable (km²).

In this paper we use the gross product output per unit of land area to represent its economic benefit coefficient. Combined with the statistical yearbook data for cities in WMA, the output values of agriculture, forestry, animal husbandry, and fishery are regarded as the output benefits of farmland, forest, grassland, and water area, respectively. The economic benefits of built-up land are represented by the GDP of the secondary and tertiary industries. We calculated and obtained the economic benefit coefficient of each category in the period 2010–2020, used the grey model for prediction, and determined that the financial benefit coefficients for farmland, forest, grassland, water, and built-up land in 2030 were 134.56, 18.80, 10.26, 15.69, and 18,531.26, respectively. Based on related research [28], we set the economic benefit of unused land as one thousand yuan per km². Finally, the objective function of economic benefit of WMA was obtained as: $F_2\left(x\right)=$ 134.56$x_1$ + 18.80$x_2$ + 10.26$x_3$ + 15.69$x_4$ + 18,531.26$x_5$ + 0.01$x_6$.

(3) Determine the gray constraint

Whether the selection of constraint parameters is accurate or not determines whether the simulation results are feasible, and also determines whether a static model close to the discrete-time dynamic system can truly reflect the laws of the system itself.

(a) Total land area constraint: The sum of the land use area of each decision variable is the total area of the study area (58,105 km²):

$${\sum }_{i=1}^6{x}_i=58,105$$

(4)

(b) Constraints on the total population: The total population carried by agricultural land and urban land should be controlled within the total population in the forecast year of 2030.

$$D_1({\sum }_{i=1}^3{x}_i)+D_2\left(0.71x_5\right) \le P$$

(5)

In formula (5), $D_1$ and $D_2$ are the average population densities of agricultural land and urban land in WMA in the target year (person/km²); the coefficient of 0.65 is the average of the proportion of urban area to the total area of built-up land from 2010 to 2020; and $P$ is the predicted value of the total population of WMA in 2030. According to the population and LULC structure data from 2010 to 2020 and using the grey model for prediction: 71 ($x_1+x_2+x_3$) + 6300 ($0.65x_5$) $\le 37,000,000.$

(c) Food security constraints: Jianghan Plain is an important food base for safeguarding WMA and even Central China, and ensuring regional food security is an important measure by which to achieve the sustainable development of WMA.

$$x_1\times f_{3 }\times f_4\times f_5\ge f_1\times f_2\times P$$

(6)

In Formula (6), $x_1$ is the area of arable land; $f_3$ is the planting proportion of grain crops in the target year; $f_4$ is the multiple planting index; $f_5$ is the grain yield per unit of arable land (unit/km²); $f_1$ is the grain self-sufficiency rate; $f_2$ is the standard grain consumption per capita (according to China’s 2018 basic target of food consumption and nutrition, per capita grain consumption is 517.30 kg/a); $P$ is the predicted total population in 2030 in the target year; obtained from the statistics of the economic data of cities in the WMA: $f_3$ = 0.55, $f_4$ = 2, $f_5$ = 550,000, $f_1$ = 1, $f_2$ = 517.30, $P$ = 37,000,000, $x_1$ ≥ 31,636.53.

(d) Cultivated land area constraints: In addition to the food security guarantee mentioned above, due to the expansion of urban built-up land, the farmland area in the WMA has slightly decreased from 2000 to 2020, thus it is unlikely that the farmland area will increase in 2030. The highest limit can be obtained at the constraint level:

$$31,636.53\le x_1\le 32,031.27$$

(7)

(e) Forest area constraints: WMA is implementing a more stringent strategy of returning farmland to forests and nature reserves, which can effectively slow down the reduction of forest area. Therefore, we set the upper limit of forest area in 2030 as the current area and the minimum as the current area. According to the forest area simulated by Markov:

$$16,168.23\le x_2\le 16,659.22$$

(8)

(f) Grassland area constraints: According to the trend of substantial reduction of grassland in WMA over the years and the need to protect grassland and strengthen the construction of artificial grassland in the “WMA Territorial Spatial Plan (2020–2035)”, we set a multi-objective scenario for grassland. The area should be smaller than the status quo, but larger than the grassland area under the natural development scenario. However, under the influence of the policy, the reduction rate of grassland area is decreasing. Markov prediction suggested that the various land areas in 2030 are based on the change trend from 2010 to 2020. Therefore, the reduction rate for grassland area under the multi-objective scenario simulation should be lower than that of the Markov prediction. The simulated grass reduction rate is as follows:

$$1261.06\le x_3\le 1323.99$$

(9)

(g) Water area constraints: In the past 20 years, the water area of WMA has been drastically reduced, from 5120.65 km² in 2000 to 3641.5 km² in 2020. This is also a major problem faced by many large cities in China. At present, in order to curb the encroachment of water resources, Hubei Province and Wuhan City have introduced relevant management measures to strictly protect the rational use of water resources. Therefore, although the water area will continue to shrink under the multi-scenario model, the rate will be significantly slower than before, given as follows:

$$3099.83\le x_4\le 3641.5$$

(10)

(h) Constraints on the area of built-up land: From 2000 to 2020, the rapid expansion of built-up land in the WMA promoted economic growth on the one hand, and aggravated environmental damage on the other. Therefore, in response to the Chinese government’s incremental development turning to the proposal of stock development, the future growth of built-up land in the WMA must be more intensive and economical. Therefore, we set the upper limit of the built-up land area in 2030 as the built-up land area based on the Markov simulation, and the lower limit is represented by the current built-up land area:

$$4358.19\le x_5\le 6299.09$$

(11)

(i) Unused land area constraints: Since the WMA is still in the stage of rapid development, some unused land may still be developed as built-up land. Therefore, we determined that the unused land area in 2030 should be smaller than the current situation, but will not disappear:

$$0\le x_6\le 91.74$$

(12)

(j) Decision variable non-negative constraint: the area of each land use type is not negative:

$$x_i\ge 0, i=1,2,3,4,5,6$$

(13)

(4) GMOP model solution

On the basis of using the predicted value to whiten the gray number in the GMOP, the function is solved using Lingo software, and finally the area of 6 land use types is obtained, representing the land demand quantity of the multi-objective optimization scenario in 2030.

3.2.2. Using PLUS Model to Simulate Spatial Distribution of Patches

The PLUS model is a new simulation model of land use change generated by patches. Its principle comes from the CA model and has been dramatically improved based on the CA model. The PLUS model has been greatly improved on the basis of the traditional CA model. It combines transformation and pattern analysis strategies to improve the accuracy of regional land use simulation. It has been widely used in the related research practice of land use structure simulation.

(1) Selection of driving factors for land change

LULC change comprises the all-round performance of internal and external factors such as social economy and nature. With the rapid development of WMA, land use change is not only affected by natural elements but also by a combination of driving factors such as the social economy and spatial location. Based on the research results and data availability for the relevant driving factors of land use change, we selected elevation, slope, temperature, precipitation, GDP, distance from water, and the government from the terrain, accessibility, and socioeconomics. Twelve driving factors, such as distance from the railway station and spread to different grades of a road, are shown in Figure 3. These 12 drivers were rasterized and unified to a resolution of 30 m × 30 m to ensure smooth model operation.

(2) Cost matrix and limit expansion area setting

The cost matrix represents the conversion rules between various types of land and reflects whether each type of land can be converted into other types. When a specific type of land cannot be transformed into different types of land, the corresponding value of the matrix is 0; when the transformation is allowed, the corresponding value of the matrix is 1. Combined with the actual situation of the change of land use types in the WMA, with the development of economic and technological levels, it is entirely possible to convert any land into built-up land, but the difficulty and cost of converting built-up land to other types of land will be great, and the actual situation is perfect. Therefore, in this study we set the constraint that built-up land could not be converted into other types of land. Secondly, it is impossible to directly judge whether the mutual transformation between water, forest land, farmland, grassland, and unused land is allowed, thus the model cost matrix parameters in different scenarios need to be set according to their constraints.

Due to the requirements of relevant policies, the basic needs of some land use types remain unchanged. To truly reflect the simulation process of different scenarios, we set restricted development areas for ecological protection scenarios and ecological and economic coordination scenarios. Among these, the environmental protection scenario masks the ecological protection red line in the area and sets it as a restricted expansion area. During the simulation process, the conversion of land within the red line to other land types is prohibited.

(3) Accuracy verification

The development of WMA has expanded the city scale; built-up land is rapidly growing towards the border counties (cities) of various cities. A critical zone appears between areas of built-up land. The GMOP and PLUS model starts from “top-down” and “bottom-up” perspectives, comprehensively considering natural, economic, and social factors, as well as limiting factors such as local ecological control areas. The data are simulation benchmark data. The PLUS model was used to simulate the land use situation of the WMA in 2020. The kappa coefficient and the Fom coefficient were used to verify the accuracy of the simulation results. The kappa value was set to 0.80, and the Fom value was set to 0.29, indicating that the model simulates urban agglomeration. The accuracy of land use change was high, indicating that the model is reliable and stable and can be used to simulate land use change in the WMA in 2030.

3.2.3. ESV Accounting Method

ESV accounting was carried out on the basis of research scholars in China and abroad [16]. In order to reduce the impact of regional differences on the ESV assessment results, the national ESV coefficient table was revised to construct an ESV accounting model for the WMA. The ESV equivalent factor coefficient was determined based on the relative contribution rate of ecosystem services that different ecosystems can produce, and its value is about 1/7 of the average grain yield market value in the study area. According to the “China Agricultural Product Price Survey Yearbook”, the average rice yield in the study area from 2000 to 2020, as the main rice producing area, was 5701.44 kg/ha. Calculated based on the average market price of rice in the same period of 1.89 yuan/kg, for the revised ESV in the study area, the unit factor value amount is 1539.39 yuan/ha. On this basis, the ESV coefficient per unit area of WMA can be calculated (Table 2), and the ESV of WMA can be obtained by combining formula (14). Due to the low ESV of built-up land, it was not considered in this study. The trade-off synergy model (Equations (15) and (16)) was introduced to explore the relationships between different ecosystem service functions from 2000 to 2020. In order to ensure the reliability of the ESV coefficient correction results for ESV estimation, a sensitivity test was carried out on the correlation of ESV and value coefficients of various land use types in WMA:

$$V_{ES}={\sum }_{i=1}^n{A}_i\times C_{Vi}$$

(14)

In Formula (14), $V_{ES}$ is the ecosystem service value; $A_i$ is the area of the $i$ land use type; and $C_{Vi}$ is the ESV coefficient per unit area of the $i$ land use type:

$$I_i=\frac{S_{ib}-S_{ia}}{S_{ia}}$$

(15)

$$D_{ij}=\frac{1}{2}\left(\frac{I_i}{I_j}+\frac{I_j}{I_i}\right)$$

(16)

In Formulas (15) and (16), $S_{ia}$ and $S_{ib}$ are the values of the $i$ ecological service function at the beginning and end of the research period, respectively; $I_i$ is the change rate of the $i$ ecological service function value; $D_{ij}$ denotes the $i$ and $j$ ecological services; the trade-off of functionality is the degree of synergy. $D_{ij}$ < 0 means that the $i$ and $j$ ecosystem services are negatively correlated, showing a trade-off relationship; conversely, $D_{ij}$ > 0 means that the $i$ and $j$ ecosystem services are negatively correlated, showing the same increase or the same reduced synergy.

3.2.4. ESV Spatial Analysis

To explore the spatial distribution characteristics and evolution rules of ESV in WMA, the ArcGIS 10.2 software Create Fishnet function was used to create grids of different scales such as 30 m × 30 m, 1 km × 1 km, 5 km × 5 km, 10 km × 10 km, etc., to characterize the ESV spatial distribution. According to the experiment, using a 5 km × 5 km grid scale in the study area can not only ensure that the research demand is met but also ensure that the calculation amount is relatively low. ArcGIS 10.2 was used to analyze the spatial correlation of ESV in WMA, and the global Moran index was used for the spatial autocorrelation analysis. The value range of the Moran index is [−1, 1]. The space is positively correlated when the Moran index is greater than 0. The higher the index, the stronger the correlation of the ESV spatial distribution and the more pronounced the spatial agglomeration. When the space is negatively correlated, and the index is lower, the weaker the correlation and the more significant the spatial difference. On this basis, the hotspot detection tool of the ArcGIS 10.2 software was used to analyze the dynamic evolution law of ESV hotspots in WMA from 2000 to 2020, and the spatial distribution map of regional ESV hotspots and cold-spots was drawn.

4. Results

4.1. Characteristics of Temporal and Spatial Evolution of Land Use

In the past 20 years, with the rapid socioeconomic development of the WMA, significant changes have occurred in terms of land use [29]. According to

Table 4. List of changes in land use area in WMA from 2000 to 2020 (km²).

Land Use Type	2000	2010	2020	2000–2010	2010–2020	2000–2020
Farmland	32,525.60	32,899.56	32,031.27	373.96	−868.29	−494.33
Forest	17,270.79	17,165.46	16,659.22	−105.33	−506.24	−611.57
Grassland	1482.60	1417.07	1323.99	−65.53	−93.08	−158.61
Water	5120.65	4503.89	3641.50	−616.76	−862.39	−1479.15
Built-up land	1673.91	2068.98	4358.19	395.07	2289.21	2684.28
Unused land	32.25	47.43	91.74	15.19	44.30	59.49

, built-up land and water area have changed the most in the past 20 years. Built-up land has increased by 2684.28 km², accounting for 61.59% of the total built-up land in 2020, and the water area has decreased by 1479.15 km², accounting for 40.62% of the water area in 2020. From 2000 to 2010, the most significant area changes were still represented by built-up land and water areas, which increased by 395.07 km² and decreased by 616.76 km², respectively. The increase of 2289.21 km² and the decrease of 868.29 km², respectively, show that the past ten years have been the most intensive time for urban and rural construction in WMA, and a large amount of farmland and water area have also been occupied in the process.

According to Figure 4, the spatial area of land use has changed in the WMA over the past 20 years; the area with the most significant change in built-up land is located in the central part of the WMA; that is, the peripheral area of Wuhan’s main site, and the urban regions of surrounding cities. In the rapid expansion of built-up land, the surrounding farmland and waters have also been continuously occupied, which is especially obvious from 2010 to 2020. Therefore, the urban areas of various cities, mainly Wuhan, have experienced the most significant changes in built-up land, farmland, and water. Analysis of the Sankey diagram shows that in the period 2000–2010, the slight shift in farmland was mainly due to the conversion of some water areas and forest land into farmland, while in the period 2010–2020, a large amount of forest land, farmland, and water areas were converted into built-up land. Other than that, grassland and unused land have not changed much in general over the past 20 years, which is partly related to their small size and partly because unused land is generally more remote in location and therefore has changed less.

4.2. ESV Spatiotemporal Evolution Characteristics

Based on the 5 km × 5 km grid scale, the unit grid ESV of WMA was divided into the low-value area (0–800), low-value area (800–1600), and medium-value area (1600–2500) using the natural discontinuous point method. This method was also used to obtain two high-value areas: (2500–4000) and (>4000). From

Table 5. Variation in the value of individual ecological services in WMA from 2000 to 2020 (billion yuan).

First Class Type	Secondary Type	ESV			ESV Variation
		2000	2010	2020	2000–2010	2010–2020	2000–2020
Supply services	Food production	21.10	21.66	20.70	0.57	−0.96	−0.39
	Raw material	34.10	34.27	33.03	0.17	−1.24	−1.07
Conditioning services	Gas regulation	52.76	53.00	51.09	0.24	−1.91	−1.67
	Climate regulation	57.69	59.18	56.03	1.49	−3.15	−1.66
	Water flow regulation	85.94	100.17	84.33	14.22	−15.84	−1.61
	Waste treatment	67.67	79.15	66.63	11.48	−12.52	−1.04
Support services	Soil maintenance	63.08	63.39	61.26	0.30	−2.13	−1.82
	Biodiversity protection	65.29	67.81	63.47	2.51	−4.34	−1.83
Cultural services	Aesthetic view supply	30.46	33.74	29.69	3.29	−4.05	−0.76
Total		478.10	512.38	466.24	34.28	−46.14	−11.86

and Figure 5, it can be seen that from 2000 to 2020, the ESV in WMA showed a trend of rising first and then falling, from 478.10 in 2000. One hundred million yuan rose to 51.238 billion yuan in 2010 and dropped to 46.624 billion yuan in 2020, an overall decrease of 1.186 billion yuan. The values of nine secondary types of ecological services decreased, as shown in Figure 6, among which soil conservation and biodiversity maintenance decreased the most. From the perspective of spatial distribution (Figure 5), the areas where food production, climate regulation, soil conservation, biodiversity maintenance, and providing aesthetic landscape increased and decreased overlap to a certain extent, and show a global crossover phenomenon; the areas where raw material production and gas regulation functions increased and decreased are mainly in the eastern part of WMA, and there is little change in Tianxianqian (Tianmen Xiantao and Qianjiang). The increase in water flow regulation and waste treatment mainly occurs in the Yangtze River Basin and near the East Lake in Wuhan, indicating that the implementation of the protection of the Yangtze River and the environmental protection measures of the East Lake are relatively good. From Figure 5, it can be seen that the spatial distribution pattern of ESV in WMA is relatively stable. The impression concept of the urbanization process is the most significant. The large-scale transformation of natural landscapes by human activities has resulted in a large amount of farmland being converted into built-up land, which has eventually led to a continuous reduction in the value of ecosystem services; high and high-value areas are mainly distributed in the northeast of WMA. The mountainous and hilly areas in the southeast are dominated by woodlands, grasslands, and water. In general, the ESV distribution pattern of WMA from 2000 to 2020 is high in the east and low in the west. The global spatial autocorrelation analysis showed that the Moran index was above 0.6 each year, indicating that the spatial distribution of ESV in WMA showed a significant positive spatial correlation and maintained a high degree of spatial aggregation. According to the hotspot detection results in Figure 7, the cold-spot agglomeration areas are distributed in the western area of WMA, including Qianjiang, Xiantao, Tianmen, and other locations. The hotspot agglomeration areas are mainly concentrated in the Dabie Mountain area and Mufu Mountain area in the east, including Huanggang, Xianning, and other locations. From 2000 to 2020, the changes in the agglomeration areas of cold- and hotspots were not significant. The number of cold-spots showed a slow growth trend. In terms of the number of hotspots, this showed a decreasing trend from 2000 to 2010, and the number continued to increase after 2010.

4.3. Analysis of Multi-Scenario Simulation Results in 2030

4.3.1. Multi-Scenario Land Use Simulation

(1) ND scenarios

In this study, we assumed that the development law of the land use types employed remained unchanged in the ND scenario; as shown in Figure 8 and

Table 6. LULC area and change in WMA under each scenario in 2030 (km²).

Year	Index	Farmland	Forest	Grassland	Water	Built-Up Land	Unused Land
2020	Area	32,031.27	16,659.22	1323.99	3641.50	4358.19	91.74
2030 (ND)	Area	31,158.80	16,168.23	1261.06	3099.83	6299.09	115.12
	Variation	−872.47	−490.99	−62.93	−541.67	1940.9	23.38
2030 (EFD)	Area	31,626.16	16,552.68	1297.57	3299.34	5298.08	28.29
	Variation	−405.11	−106.54	−26.42	−342.16	939.89	−63.45
2030 (EECD)	Area	31,356.162	16,282.6839	1279.9998	3344.3352	5786.6778	52.272
	Variation	−675.11	−376.54	−43.99	−297.16	1428.49	−39.47

, the original 2030 forecasts were generated by the CARS module. Under the natural development scenario, WMA’s farmland, forest land, grassland, and water area have all declined. The farmland area in 2030 will be 31,158.80 km², a decrease of 872.47 km²; the forest land area will be 16,168.23 km², a reduction of 490.99 km²; the grassland area will be 1261.06 km², a decrease of 62.93 km²; and built-up land and unused land will increase by 6299.09 km² and 115.12 km², respectively. It can be seen that in the natural development scenario, the expansion of built-up land is the most obvious, and in the process, farmland and water are mainly occupied.

(2) EFD scenarios

Under the EFD scenario, although the farmland, forest, grassland, water area, and unused land in WMA will still decrease, according to Figure 8 and Table 6, the decline is moderated compared with the realistic scenario. The farmland area in 2030 will be 31,626.16 km², a decrease of 405.11 km²; the forest area will be 16,552.68 km², a reduction of 106.54 km²; the grassland area will be 1297.57 km², a decline of 26.42 km²; the water area will be 16,552.68 km², a decrease of 342.16 km²; and the built-up land area will be 5298.08 km², an increase of 939.89 km². It can be seen that under the EFD scenario, the expansion of built-up land is limited, and ecological land, such as farmland, water area, and forest land, is primarily protected to avoid the phenomenon of environmental land encroachment caused by construction.

(3) EECD scenario

Under this scenario, compared with 2020, according to Figure 8 and Table 6, other types of land in WMA are still decreasing except for built-up land. Among them, farmland has reduced the most, with an area of 675.11 km², followed by forest land, with an area of 376.54 km², and built-up land has increased, with an area of 1428.49 km². Compared with the ND and EFD scenarios, this scenario presents an obvious compromise mode. Although the ecological land is still being eroded, the degree can be controlled within a reasonable range.

4.3.2. Simulation of Temporal and Spatial Changes in ESV

Table 7. Comparison of indicators under various scenarios in WMA from 2020 to 2030.

Scene	Index	Farmland	Forest	Grassland	Water	Built-Up Land	Unused Land
2020 (current situation)	Area (km2)	32,031.27	16,659.22	1323.99	3641.50	4358.19	91.74
	ESV (104 yuan)	1,307,836.75	2,420,917.85	79,849.84	853,421.94	0.00	658.69
	Economic benefit	4,310,127.69	313,193.34	13,584.14	57,135.14	80,762,752.02	0.92
2030 (ND)	Area (km2)	31,158.80	16,168.23	1261.06	3099.83	6299.09	115.12
	ESV (104 yuan)	1,272,213.88	2,349,567.44	76,054.59	726,477.09	0.00	826.53
	Economic benefit (104 yuan)	4,192,728.38	303,962.75	12,938.48	48,636.39	116,729,993.00	1.15
2030 (EFD)	Area (km2)	31,626.16	16,552.68	1297.57	3299.34	5298.08	28.29
	ESV (104 yuan)	1,291,296.19	2,405,436.02	78,256.64	773,232.20	0.00	203.15
	Economic benefit (104 yuan)	4,255,616.36	311,190.46	13,313.10	51,766.57	98,180,155.43	0.28
2030 (EECD)	Area (km2)	31,356.16	16,282.68	1279.99	3344.33	5786.67	52.27
	ESV (104 yuan)	1,280,272.09	2,366,199.62	77,196.78	783,778.39	0.00	375.31
	Economic benefit (104 yuan)	4,219,285.15	306,114.45	13,132.79	52,472.61	107,234,430.8	0.52

and Figure 9 show the values of ecosystem services in WMA in 2000, 2010, 2020, and 2030 under the three scenarios. The ecological service values of WMA in 2030 are 442.51 billion yuan (ND scenario), 454.84 billion yuan (EFD scenario), and 450.78 billion yuan (EECD scenario). It can be found that the ecological service value shows a downward trend. In the three scenarios, farmland is the land type with the most significant transfer area. In the environmental priority scenario, the transfer of forest, grassland, and water area to other land types is restricted. Compared with the ND scenario, the ecosystem service values of farmland and forest land under the EFD scenario have increased by 4.67 million yuan and 3.84 million yuan, respectively. In the EECD scenario, farmland is protected, the area of farmland is reduced, and other land types are allowed to transfer to built-up land. There is a 1.97 million yuan increase in the ecological service value of farmland in this scenario compared with the ND scenario. Compared with the EFD scenario, the ecological service value is reduced by 2.70 million yuan. However, from an ecological economics perspective, the total economic value under the scenario of EECD is 111,825.43 billion yuan, an increase of 901.33 billion yuan compared with the EFD model. The EECD scenario considering the urban economic development can prevent ecological degradation to a certain extent and can better achieve the development goals of new urbanization and ecological civilization construction required by the Chinese government.

In the context of the total decline in the three scenarios simulated in the WMA in 2030, according to

Table 8. Changes in the value of individual ecological services in WMA under different scenarios (billion yuan).

Type	2020 (Current Situation)	2030 (ND)	Variation	2030 (EFD)	Variation	2030 (EECD)	Variation
Food production	20.70	20.00	−0.70	20.37	−0.33	20.20	−0.50
Raw material	33.03	31.99	−1.04	32.72	−0.31	32.25	−0.78
Gas regulation	51.09	49.48	−1.61	50.59	−0.50	49.88	−1.21
Climate regulation	56.03	53.94	−2.09	55.22	−0.81	54.55	−1.48
Water flow regulation	84.33	77.64	−6.69	80.60	−3.73	80.35	−3.98
Waste treatment	66.63	61.41	−5.22	63.63	−3.00	63.53	−3.10
Soil maintenance	61.26	59.39	−1.87	60.62	−0.64	59.85	−1.41
Biodiversity protection	63.47	60.85	−2.62	62.36	−1.11	61.66	−1.81
Aesthetic view supply	29.69	27.82	−1.87	28.74	−0.95	28.52	−1.17

, from the perspective of individual ESVs, in the ND scenario, it is mainly water regulation and waste disposal that have the largest decline, i.e., decreases of 669 million yuan and 552 million yuan, respectively. Under the EECD scenario, the change trends of various types of ESV are the same, with the smallest declines in food production and raw material production, which have dropped by RMB 50 million and 78 million, respectively, and with the largest declines in respect of water flow regulation and waste treatment, at 398 million and 310 million, respectively. Under the three scenarios, the order of each type of ESV is as follows: water flow regulation > waste treatment > maintenance of biodiversity > soil conservation > climate regulation > gas regulation > raw material production > aesthetic view supply > food production.

5. Discussion

LULC change is the main factor influencing the response of ESV. In this study, we used the PLUS-GMOP coupling model and land use data for WMA in 2000, 2010, and 2020 to simulate the land use situation of WMA in 2030 under various scenarios, and further measured the value of ecosystem services in WMA.

From the perspective of the characteristics of LULC change, the LULC change in WMA is a complex result of multiple factors. Farmland and water areas are the primary sources of transfer of built-up land, forest land, and grassland. The continuous reduction of farmland in the WMA is closely related to the government’s strict implementation of the policy of “returning farmland to forests and grasslands” and urban expansion under the background of new urbanization. This is consistent with research that concluded that the Jianghan Plain in the west of the WMA is a crucial food production base in China [30,31], thus it is urgent to protect the farmland. Forests and waters are of great significance to ecosystem services. Wuhan is known as the city of a thousand lakes. The northern and southern parts of the WMA also have two principal mountain systems south of China, namely the Dabie Mountains and the Mufu Mountains. WMA is not only a critical water source protection area but also a dense biological source area of ecological corridors, and is worthy of vigorous protection [32].

As shown by the changing trend of land use simulation results, the simulation results of various types of land are the same as the development trends under this scenario. The PLUS model has good accuracy in the simulation of different land types. By comparing the simulation results for 2020 and actual data, the simulation accuracies of the three types of farmland, forest land, and construction land that account for the most significant area are 0.896, 0.899, 0.896, and 0.912, guaranteeing the feasibility of future multi-scenario simulation accuracy. In the three scenarios, the built-up land shows different degrees of expansion. Among them, in the ND scenario, the built-up land has the fastest growth, increasing by 1940.9 km² in ten years, while the farmland has decreased by 872.47 km². The contradiction between the expansion of built-up land and environmental protection is a problem that must be addressed in the future development of WMA.

There is a certain degree of subjectivity in the study of the neighborhood weights set by the model; in the selection of land use driving factors, data availability are mainly considered, thus it is impossible to identify all the factors driving land use change comprehensively. Moreover, the Markov chain model for the simulation of total land use can only simulate the future according to the change law of existing land use types. Future policy factors cannot be taken into account [33]. In follow-up research, the accuracy of the regional ecosystem service value scale can be improved, and the ecosystem service value can be predicted with high-precision biological and socioeconomic data combined with land use data.

6. Conclusions

In this research, we proposed a coupled GMOP and PLUS model method for the multi-scenario simulation of LULC structure and ESV. Through the transformation of land use types in WMA from 2000 to 2020 and the changing trend of ESV, it was concluded that the changes in land use types in WMA in the past two decades are mainly reflected by the reduction in farmland and water area and the increase in built-up land. This shows that there is still a problem in respect of urban expansion and ecological protection. The ESV has experienced a trend of rising first and then falling, increasing by 3.42 billion yuan from 2000 to 2010 and decreasing by 4.61 billion yuan from 2010 to 2020. The overall decline of 11.86 billion yuan is also the inevitable result of the rapid growth of built-up land and the reduction in ecological lands such as farmland, water area, and forest land. In addition, we predicted the land use structure and ESV in 2030. It was found that the future ESV shows a downward trend under the three preset scenarios. Among them, the ESV is the highest in the EFD scenario, and the ND scenario presents the highest economic value. The ESV and economic benefit are not the highest in the EECD scenario. Compared with the EFD scenario, the value is slightly decreased. Still, the economic value is significantly increased, which represents a problem of coordination between economic development and ecological protection that must be considered in the development of such large urban areas. Therefore, these three scenarios can be used as research references for the future development of large metropolitan areas such as WMA, which are rapidly expanding in urban construction. Effective protection of the ecological environment and improving the service function of the ecosystem are the top priorities.

In this study, we proposed the use of coupled PLUS-GMOP models to simulate and predict multiple scenarios of land use and ecological service values in a large metropolitan area such as the Wuhan metropolitan area, innovatively solving the problem of spatial distribution that is difficult to achieve with traditional system dynamics models, and taking advantage of multi-objective algorithms to improve the referenceability of the simulation result. The research results show that the multi-scenario simulation method using the coupled PLUS-GMOP model can be used to scientifically intervene in the rational use of land resources and the maintenance of the ecological environment in large urban areas. It also provides methods and paths for sustainable urban development and spatial planning

In future land use structure simulation and ESV research in large urban areas, the microscopic reasons for changes in ESV are worth exploring. In addition, future development scenarios are unpredictable. The three scenarios proposed in this paper will not represent all possibilities. Making full use of machine learning algorithms, submitting intelligent adaptive simulation scenarios, and updating feedback in real-time with the city’s construction to assist in the urban planning of each region are of great significance for improving the science of planning and the effectiveness of the model.