Estimating the Past and Future Trajectory of LUCC on Wetland Ecosystem Service Values in the Yellow River Delta Region of China

Abstract

Land use/cover change (LUCC) can impact the provision of ecosystem service values (ESVs), particularly in wetland regions that are subject to frequent and unsustainable land conversions. Exploring the past and future trajectory of LUCC and its effects on ESV has a great significance for wetland management and habitat stability. This study tried to reveal the patterns and magnitude of LUCC on ESV under varying land development scenarios in the Yellow River Delta region, which is a typical region undergoing serious degradation in China. In this study, a combined approach utilizing equivalent coefficients of ecosystem services was employed to determine the ESV of the wetland in relation to the major land use types (LUTs). The Markov–FLUS model was then used to simulate LUTs across multiple scenarios in 2030 and to clarify the relationship of ESV between wetland and other LUTs. The results indicated that the wetland was severely degraded, with a loss in area of 6679.89 ha between 2000 and 2020. Cropland and water body were the main sources of diversion and turnover for the wetland, respectively. Despite the multiple scenario projections revealed, the wetland area exhibited a similar growth rate and a homogeneity in ESV under the natural development (ND), urban construction and development (UCD), and the ecological development (ED) scenarios. The ED scenario was deemed the optimal development strategy for the wetland ecosystem. Our research will improve the comprehension of land development decisions and promote sustainable development in estuarine wetland areas.

1. Introduction

The meaning of ecosystem services (ESs) is the ecological functions from which peoples can attain diverse benefits in various natural ecosystems [1,2,3]. Ecosystem service values (ESVs) apply economic rules to estimate ESs and transform them into practical applications [4,5]. Wetlands play a crucial role in contributing to the global stability of ecosystems. Globally, more than 23% of the annual ESV is provided by wetlands [6]. However, as human activities intensify and bring about the loss of biodiversity and unstable change in the global climate, numerous wetlands, especially estuarine wetlands, have been severely degraded. Half of the world’s wetlands have disappeared since 1900, leading to a decline in their value and function [7,8]. Land use/cover change (LUCC) has a momentous influence on the provision of ES. Therefore, it is essential to examine the history and future trajectory of LUCC and its effects on ESV in order to make informed decisions on land development and achieve sustainable development in wetland regions [9,10].

Numerous works from the literature have probed the consequences of the land use conversion on ES through the analysis of land use types (LUTs) and ES assessment [11,12,13]. Originally, the current and former LUTs served as the basis for investigating the effects of LUCC on ESV [14,15]. Research on LUCC and its impact on ESV have primarily focused on evaluating the role of different LUTs in ESV [16,17], examining the drivers of ESV change by analyzing various LUTs in combination with natural and economic factors [18,19], and exploring trends in ESV changes due to the shifts in different LUTs [20,21]. These studies can facilitate the detection of spatial heterogeneity and trends in ESV, while the results are less practical for anticipating and mitigating future issues.

The modelling and prediction of land utilization can provide a scientific basis for future changes in ESV in wetland ecosystems. More and more land utilization models have been proposed and applied in the study of ESV, such as the cellular automata model, CA-Markov model, CLUE-S model, and so on [22,23]. However, these models can only simulate land utilization from quantity or space, and cannot take into account the impact of both on land utilization. To address this issue, the Markov–FLUS model has been developed and employed in several studies to simulate future patterns of land usage [24,25,26]. This method, which combines a roulette wheel selection mechanism and a neural network algorithm with the cellular automata (CA) model, can be achieved by generating an adaptive atlas that integrates a range of natural and social factors, and the impact on land utilization can be considered both quantitatively and spatially. Its results will make spatial trajectories more explicit, allowing for the simulation of various LUTs in different scenarios.

Based on the LUTs in a spatially explicit simulation, the utilization of biophysical models [27] and the equivalent-coefficient method [28] can be employed to evaluate future ESVs. The lack of available data hinders the use of biophysical models to predict ESV under alternative scenarios due to the need to parameterize many variables [29]. In contrast with biophysical models, equivalent-coefficient approaches predict ESVs more easily, with less data requirements and acceptable results [30]. As a result, these methods were widely used to assess changes in ESV in response to LUCC [31,32]. However, some ESV assessments using the equivalent-coefficient method assumed identical ESV per acre for a particular LUT without giving thought to the heterogeneity of ES for a particular LUT in diverse zones [33,34,35]. Therefore, it is crucial to establish appropriate guideline systems that are suitable for the generation of spatially explicit ESVs per unit area and allow for the analysis of spatial disparities in ESVs [36,37]. The adjusted equivalent-coefficient method allows for ESV calculations under data-limited conditions. Additionally, the Markov–FLUS model can simulate land utilization in both quantitative and spatial perspectives, enabling the prediction and assessment of ESVs accurately. These methods make it possible to study the impact of historical and future trajectories of LUCC on ESVs.

The wetland in the Yellow River Delta (YRD) is a unique estuarine ecosystem. Compared with other estuarine wetlands, due to the large amount of sediment carried by the Yellow River into the sea, a large number of shallows and wetlands will be formed at the confluence of the river and the sea, and the wetland changes are very frequent, which is exacerbated by human activities. In addition, considering the availability of data, the YRD was chosen as the study region. The study had three objectives: (1) to investigate the LUCC over the last 20 years in the study region and the pertinence between wetland and various LUTs; (2) to summarize the change pattern of ESVs under past LUCC conditions and determine wetland’s contribution to ESV; and (3) to analyze the influence of different development scenarios on the regional wetland ecosystem and overall LUTs. Exploring the history and future trajectory of LUCC and its effects on ESV will enhance the comprehension of estuarine wetland ecosystem development and regional ecological conservation.

2. Materials and Methods

2.1. Study Area

The Yellow River Delta (YRD) is located between 118°33′~119°20′ E and 37°35′~38°12′ N, on the northern (Bohai Sea) coast of Shandong Province. It is mainly situated in Dongying (Figure 1). The climate of the region belongs to the temperate monsoon, with an annual temperature and precipitation of 12.8 °C and 555.9 mm, respectively. In accordance with the Shandong Bureau of Statistics, the population of Dongying was 2,193,500 at the end of 2020. The region’s GDP grew from CNY 46.511 million in 2000 to CNY 298.199 million in 2020. The area is an important petrochemical base in China, with topography that is lower in the northeast than in the southwest [38,39].

The protection and rejuvenation of wetlands are crucial for maintaining the security and biodiversity of the watershed ecosystem. YRD hosts a unique estuarine wetland ecosystem, which is rare on a global scale. However, the shortage of ecological water and other issues have resulted in a significant risk of deterioration. Studying the variation in ESVs in this region is of great significance for safeguarding and rehabilitating the wetland [40,41].

2.2. Data Collection and Preprocessing

The Data Sharing Service System (https://data.casearth.cn/ (accessed on 2 March 2023)) furnished the land utilization data with 30 m spatial resolution, the precision of land utilization classification for different years reached 82.5%, and the study showed that its classification products can meet the needs of applications at a global or regional scale [42]. In accordance with the data classification system and the current situation in the study region, the data were processed and classified into seven categories: artificial surface, cropland, bare land, forest, water body, wetland, and grassland. Six types of driving factors’ data were used in this study: DEM, road, temperature, precipitation, slope, and railways. The Resource and Environment Science and Data Center (https://www.resdc.cn/ (accessed on 2 March 2023)) furnished the administrative division of Dongying, road, and railway data. The data of the road and railway with 1 km spatial resolution and distances were calculated to roads and railways by using Euclidean distances in ArcGIS 10.2. Geospatial Data Cloud (https://www.gscloud.cn/ (accessed on 30 March 2023)) furnished the DEM data with 30 m spatial resolution. The slope data were obtained by processing the DEM data with ArcGIS 10.2. The China Meteorological Science Data Service Centre (https://data.cma.cn/ (accessed on 4 April 2023)) furnished the temperature and precipitation data with 1 km spatial resolution. The six driving factors were normalized using the Fuzzy Membership tool in ArcGIS 10.2. The Statistical Yearbook of China, Shandong, and Dongying furnish the socioeconomic data. The National Farm Product Cost–Benefit Survey furnished market prices for different grains in the study for various years.

2.3. Methods

2.3.1. Extent Dynamic of Land Utilization

The intensity of transformation between various LUTs can be denoted as the integrated extent dynamic of land utilization [43]:

$$ID=\frac{{\displaystyle \sum _{\mathrm{i}=1}^n\left|U_{in-i}-U_{out-i}\right|}}{2{\displaystyle \sum _{i=1}^n{U}_{bi}}}\times \frac{1}{C}\times 100\%$$

(1)

where ID is the integrated extent dynamic of land utilization. U_in-i is the total of the other LUTs converted to LUT i over the course of the study. U_out-i is the total of LUT i converted to the other LUTs over the course of the study. U_bi is the area of LUT i at the initial phase of the study. C is the study period.

The LUTs’ change rate can be denoted as the single extent dynamic of land utilization [44]:

$$SD=\frac{G_b-G_a}{G_a}\times \frac{1}{H}\times 100\%$$

(2)

where SD is the single extent dynamic of land utilization. G_a is the area of a LUT at the initial phase of the study. G_b is the area of a LUT at the final phase of the study. H is the study period.

2.3.2. Land Use Transfer Matrix

The land use transfer matrix can describe the conversion relation between different LUTs and reveal information about transfer in and out of each LUTs [45]:

$$A_{ij}=\left[\begin{array}{cccc}{A}_{11}& A_{12}& \cdots & A_{1n}\\ A_{21}& A_{22}& \cdots & A_{2n}\\ ⋮& ⋮& & ⋮\\ A_{n1}& A_{n2}& \cdots & A_{nn}\end{array}\right]$$

(3)

where A_ij is the area of LUT i converted to LUT j over the course of the study (i, j = 1, 2, …, n). In this study, n is equal to 7.

2.3.3. Assessment of Ecosystem Service Values (ESVs)

For an accurate assessment of ESVs in the study region, the ecosystem service equivalent value was adjusted by socioeconomic data in the study region [46].

• Adjustment of ecosystem service equivalent factor

Applying the ecosystem service equivalent factor of the other study region to the YRD in this paper would lead to inaccurate results [47]. The ecosystem service equivalent factor was adjusted by using a method based on food production per unit of area:

$$\epsilon =\frac{N_d}{N_c}$$

(4)

$$E_{\mathrm{i}}{}_{\mathrm{j}}=\epsilon \times E_{oij}$$

(5)

where ε is the regulatory factor. N_d is the harvested food per unit area in the study region. N_c is the harvested food per unit area in China. E_ij is the ecosystem service equivalent factor of ES function j for LUT i after adjustment. E_oij is the ecosystem service equivalent factor of ES function j for LUT i before adjustment. The calculated ε is 1.06.

2.

The determination of ecosystem service value onset per hectare

The ESV per hectare is generally considered to be equal to one-seventh of the economic worth of grain harvested on a hectare of cropland at prevailing market prices. The main grain products in this study region are soybean, maize, and wheat [48]:

$$E_x=\frac{{\displaystyle {\sum }_{n=1}^m\frac{u_n{v}_n{s}_n}{A}}}{7}$$

(6)

where E_x is the ESV per hectare (CNY ha⁻¹). The n is the grain type. The u_n (CNY ton⁻¹) is the average market price of the nth grain. The v_n (ton ha⁻¹) is the average unit yield of the nth grain. The s_n (ha) is the average area planted of the nth grain. The A (ha) is the area planted of all grains. The calculated E_x of YRD is 1425.33 CNY ha⁻¹.

The coefficients on ecosystem service value were corrected using the method of food production per unit of area and the determination of ecosystem service value onset per hectare. Based on the ecosystem services equivalent factor table applicable to China [49] and equivalent-coefficient method, the regulatory factor, ESV per hectare, and ES value for different land utilization types applicable to the study were calculated. The adjusted coefficients are listed in Appendix A at the end of the manuscript, specifically in

Table A1. Original ecosystem service equivalent factor.

FES	SES	Cropland	Forest	Grassland	Water Body	Wetland	Artificial Surface	Bare Land
Support services	Food production	0.85	0.23	0.23	0.80	0.51	0	0.00
	Raw material	0.4	0.54	0.34	0.23	0.50	0	0
	Water supply	0.02	0.28	0.19	8.29	2.59	0	0.00
Regulation services	Air quality regulation	0.67	1.76	1.21	0.77	1.90	0	0.02
	Climate regulation	0.36	5.27	3.19	2.29	3.60	0	0.00
	Waste treatment	0.1	1.57	1.05	5.55	3.60	0	0.10
	Regulation of water flows	0.27	3.81	2.34	102.24	24.23	0	0.03
Supply services	Erosion prevention	1.03	2.14	1.47	0.93	2.31	0	0.02
	Maintenance of soil fertility	0.12	0.16	0.11	0.07	0.18	0	0.00
	Habitat services	0.13	1.95	1.34	2.55	7.87	0	0.02
Cultural services	Cultural services	0.06	0.86	0.59	1.89	4.73	0	0.01

,

Table A2. Adjusted ecosystem service equivalent factor.

FES	SES	Cropland	Forest	Grassland	Water Body	Wetland	Artificial Surface	Bare Land
Support services	Food production	0.90	0.25	0.25	0.85	0.54	0.00	0.00
	Raw material	0.42	0.57	0.36	0.24	0.53	0.00	0.00
	Water supply	0.02	0.29	0.20	8.79	2.75	0.00	0.00
Regulation services	Air quality regulation	0.71	1.87	1.28	0.82	2.01	0.00	0.02
	Climate regulation	0.38	5.58	3.38	2.43	3.82	0.00	0.00
	Waste treatment	0.11	1.66	1.12	5.88	3.82	0.00	0.11
	Regulation of water flows	0.29	4.04	2.48	108.37	25.68	0.00	0.03
Supply services	Erosion prevention	1.09	2.27	1.56	0.99	2.45	0.00	0.02
	Maintenance of soil fertility	0.13	0.17	0.12	0.07	0.19	0.00	0.00
	Habitat services	0.14	2.07	1.42	2.70	8.34	0.00	0.02
Cultural services	Cultural services	0.06	0.91	0.63	2.00	5.01	0.00	0.01

and

Table A3. ES value (CNY/ha⁻¹) for different land use types.

FES	SES	Cropland	Forest	Grassland	Water Body	Wetland	Artificial Surface	Bare Land
Support services	Food production	1284.22	352.53	352.53	1208.68	770.53	0.00	0.00
	Raw material	604.34	810.82	518.73	347.50	755.42	0.00	0.00
	Water supply	30.22	418.00	287.06	12,524.94	3913.10	0.00	0.00
Regulation services	Air quality regulation	1012.27	2659.10	1823.09	1163.35	2870.61	0.00	30.22
	Climate regulation	543.91	7957.14	4819.61	3459.85	5439.06	0.00	0.00
	Waste treatment	151.08	2367.00	1591.43	8385.22	5439.06	0.00	151.08
	Regulation of water flows	407.93	5756.34	3530.35	154,469.28	36,607.89	0.00	45.33
Supply services	Erosion prevention	1556.18	3238.25	2220.95	1405.09	3490.06	0.00	30.22
	Maintenance of soil fertility	181.30	246.77	171.23	105.76	271.95	0.00	0.00
	Habitat services	196.41	2951.19	2019.50	3852.67	11,890.39	0.00	30.22
Cultural services	Cultural services	90.65	1294.29	891.40	2855.51	7146.32	0.00	15.11

.

3.

Calculation of ecosystem service values

The ESV is calculated using the benefit transfer method [50]:

$$Q_{\mathrm{i}}={\displaystyle \sum _{j=1}^n{E}_{ij}}\times E_x$$

(7)

$$TESV={\displaystyle \sum _{i=1}^n{Q}_i}\times S_i$$

(8)

where Q_i (CNY ha⁻¹) is the ESV of one-hectare LUT i. E_ij is the ecosystem service equivalent factor of ES function j for LUT i after adjustment. E_x is the ESV per hectare (CNY ha⁻¹). TESV is the ESV in the region. The S_i (ha) is the area of LUT i.

2.3.4. Markov–FLUS Model

• Markov

The Markov model will automatically create the conversion probability matrix, the conversion area matrix, and the transition conversion of various LUTs [51]:

$$R(s+1)=P_i{}_j\times R(s)$$

(9)

$$P_{ij}=\left[\begin{array}{cccc}{P}_{11}& P_{12}& \cdots & P_{1n}\\ P_{21}& P_{22}& \cdots & P_{2n}\\ ⋮& ⋮& & ⋮\\ P_{n1}& P_{n2}& \cdots & P_{nn}\end{array}\right]$$

(10)

where R(s + 1) is the land cover condition at time s + 1. R(s) is the land cover condition at time s. The P_ij is the probability of conversion from the previous condition to the next condition.

2.

FLUS model

The FLUS model’s structure mainly consists of two parts [52].

(1) Probability calculation of land use distribution suitability based on neural network.

The Artificial Neural Network (ANN) is a feedforward neural network with a multilayer structure. Its effect is to train and predict the suitability probability of each LUT on a specific grid cell, which can be used to simulate the nonlinear complicated relationship between land use design and the drivers that are likely to influence their development.

(2) Adaptive cellular automaton and competition mechanism based on inertia.

The interaction of different LUTs is resolved by adaptive inertia and competition mechanisms. During the iterative procedure of the meta-automata, the combination probability of LUTs is estimated on each specific grid cells, and the ruling LUTs are assigned to the grid cell.

The Markov–FLUS model is a coupled model that works well for predicting future changes in LUTs.

2.3.5. Methods for Evaluating the Accuracy of Model Results

In this study, the kappa coefficient and figure of merit (FOM) were utilized to assess the simulation results. The kappa coefficient is commonly employed to assess the accuracy of classified or modeled data derived from remotely sensed images, with a focus on evaluating the data’s accuracy in terms of spatial location:

$$K=\frac{L_p-L_q}{1-L_q}$$

(11)

where K is the kappa coefficient. L_p is observation consistency, L_q is L_p is observational consistency. The calculated K is 0.8509, indicating a successful simulation, as it exceeded the threshold of 0.8 for good simulation [53].

The FOM method more effectively evaluated the metacells that underwent change, enhanced the accuracy assessment of the simulation, and proved to be more objective and realistic, particularly when the extent of land utilization change is relatively small compared to the total area of the study region.

$$FOM=\frac{R_2}{R_1+R_2+R_3+R_4}$$

(12)

where R₁ represents a region where change is observed, and the region remains unchanged during the simulation. R₂ represents a region where change is observed, and the simulation correctly categorizes the region. R₃ represents an area where change is observed, but the simulation assigns it to the wrong category. R₄ indicates a region where no change was observed, but the simulation incorrectly changes the region. The calculated FOM value is 51.26%, which was below the standard range of 1% to 59%, and the accuracy rate increases with the increase in the FOM value [54].

3. Results

3.1. Spatial–Temporal Changes in Land Use

The major LUT in YRD was cropland, which constituted the highest average area of 62.28% in YRD from 2000 to 2020, followed by water body (16.23%), wetland (14.35%), and artificial surface (6.77%). Accordingly, grassland, bare land, and forest accounted for only 0.37% of the entire area. During 2000 to 2020, the land use trend displayed a significant decrease in cropland, forest, and grassland, accompanied by a significant increase in water body, artificial surface, and bare land (

Table 1. Area and scale of different land use types in 2000, 2005, 2010, 2015, and 2020.

	2000		2005		2010		2015		2020
	Area /ha	Proportion/%	Area /ha	Proportion/%	Area /ha	Proportion/%	Area /ha	Proportion/%	Area /ha	Proportion/%
Artificial surface	39,649.27	5.17	47,000.02	6.13	52,139.22	6.80	59,369.23	7.74	61,447.37	8.01
Cropland	504,680.14	65.79	487,304.23	63.52	478,055.59	62.32	465,109.30	60.63	453,630.25	59.13
Forest	79.26	0.01	68.19	0.01	48.22	0.01	41.26	0.01	0.36	0.00
Grassland	2817.72	0.37	2808.93	0.37	2383.19	0.31	2251.59	0.29	1899.45	0.25
Bare land	250.71	0.03	298.07	0.04	343.99	0.04	394.37	0.05	482.35	0.06
Water body	101,338.46	13.21	121,427.69	15.83	127,864.81	16.67	134,012.87	17.47	138,029.81	17.99
Wetland	118,310.28	15.42	108,216.64	14.11	106,287.81	13.86	105,941.78	13.81	111,630.38	14.55

).

The noteworthy point is that the wetland, as the major LUT in YRD, was characterized by a trend of decrease, followed by an increase. From 2000 to 2015, the proportion of wetland decreased from 15.42% to 13.81%, and the area was reduced by 12,368.5 ha. The reduction in wetland area was predominantly in the northwestern part of the study region. Subsequently, between 2015 and 2020, the proportion of wetland increased to 14.55%, and the area increased by 5688.6 ha. The expansion of the wetland occurred in the estuary region (Figure A1), which has been designated as a National Ecological Reserve. The results indicated that the wetland, forest, and grassland are exhibiting a clear trend of degradation under normal development patterns, yet this deterioration can be slowed down or ameliorated through policy and human interventions.

The integrated extent dynamic of land utilization in the study region was 0.38% in 2000–2020. The forest, bare land, artificial surface, water body, and grassland exhibited a dramatic dynamic degree during the four periods of 2000–2005, 2005–2010, 2010–2015, and 2015–2020. However, there were relatively minor single extent dynamics of land utilization for these periods. The single extent dynamic of land utilization in the wetland was negative during the first three periods and positive across the last period. Alternately, positive and negative changes during the four periods resulted in a less pronounced single extent dynamic of land utilization in wetland across the entire study period (2000–2020) (Figure 2).

Throughout the entire period, the amount of cropland transferred out exceeded the amount of cropland transferred in, resulting in a reduction of 62,563.50 ha, or about 31.47%, 34.47%, and 33.11% of cropland converted to artificial surface, water body, and wetland, respectively (Figure 3). Artificial surface was dominated by the transferred in process during the period of 2000–2020, approximately 21,798.10 ha expansion from cropland (90.31%), wetland (5.75%), and water body (3.54%). During the same period, the water body experienced an expansion of about 51,101.45 ha, with 42.2% from cropland and 57.44% from wetland. Over the past twenty years, wetland conversion exhibited a relatively complex process involving both reduction and expansion. From 2000 to 2020, a total of 36,658.85 ha of wetland was transferred out, of which the proportions transferred to artificial surface, cropland, and water body are 3.42%, 16.15%, and 80.08%, respectively. At the same time, a total of 29,979.52 ha of wetland was transferred in, of which the proportions transferred in from artificial surface and water body are 69.10% and 29.87%, respectively. In general, the reduction in cropland and expansion of artificial surface and water body were the characteristics of the LUCC transformation in YRD. Moreover, the transformation of the wetland featured the complex conversions of the water body and cropland.

3.2. Spatial and Temporal Characteristics of ESV in YRD from 2000 to 2020

The total ESV in YRD increased from CNY 31.64 billion in 2000 to CNY 37.75 billion in 2020 (

Table 2. The ESV of the study area from 2000 to 2020.

		Cropland	Forest	Grassland	Water Body	Wetland	Artificial Surface	Bare Land	Total
2000	ESV/109 CNY	3.06	0.00220	0.051	19.23	9.30	0.00	0.00008	31.64
	Proportion/%	9.66	0.01000	0.160	60.78	29.39	0.00	0.00024	100.00
2005	ESV/109 CNY	2.95	0.00190	0.051	23.04	8.51	0.00	0.00009	34.56
	Proportion/%	8.54	0.01000	0.150	66.69	24.61	0.00	0.00026	100.00
2010	ESV/109 CNY	2.90	0.00140	0.043	24.27	8.35	0.00	0.00010	35.56
	Proportion/%	8.14	0.00400	0.120	68.24	23.49	0.00	0.00029	100.00
2015	ESV/109 CNY	2.81	0.00120	0.041	25.43	8.33	0.00	0.00012	36.62
	Proportion/%	7.70	0.00300	0.110	69.45	22.74	0.00	0.00033	100.00
2020	ESV/109 CNY	2.74	0.00001	0.035	26.19	8.77	0.00	0.00015	37.75
	Proportion/%	7.28	0.00003	0.090	69.39	23.24	0.00	0.00039	100.00

); there was a growth of 19.31%. The trends in the ESV can be divided into two periods: rapid growth (2000–2005) and slight growth (2005–2020). During the first period (2000–2005), the total ESV in the study region increased by about CNY 2.92 billion, there was a growth of 9.20%. In the second period (2005–2020), the total ESV increased by CNY 3.19 billion, with an average growth rate over three periods of 2.99%. During 2000 to 2020, the ESV of the wetland, cropland, grassland, and forest exhibited a slight reduction, decreasing by approximately CNY 0.53 billion, CNY 0.31 billion, CNY 16.74 million, and CNY 2.21 million, respectively. However, the ESV of the water body and bare land increased by CNY 6.96 billion and CNY 0.07 million, respectively. These findings suggested that the ESV exhibited faster growth during the initial phase, when land use changes occurred more frequently, but decelerated as the spatial arrangement of land use stabilized.

The ES function was divided into four first-level categories (FESs) and eleven second-level categories (SESs). In terms of the four FESs, the category of regulation services has the largest proportion of total ESV, exceeding approximately 75% during the 2000–2020 period. Consequently, the other three FESs, i.e., cultural services, support services, and supply services, offered smaller proportions of the total ESV, accounting for approximately 3.40%, 9.30%, and 9.50%, respectively. From the SES perspective, the largest ESC expansions were in the regulation of water flows, water supply, and waste treatment, amounting to 26.72%, 24.69%, and 16.64%, respectively. The expansion of these sectors was primarily linked to increases in the area of the water body. Meanwhile, the three SECs with the largest reduction were raw material, maintenance of soil fertility, and erosion prevention, with rates of −5.49%, −5.46%, and −3.97%, respectively. These decreases were mainly due to a decrease in cropland. Wetland contributed most to the three SESs of habitat service, cultural service, and climate regulation, with more than 67.29%, 63.94%, and 44.19%, respectively (Figure 4). These findings showed that the wetland had a crucial role in regional ecosystems.

Figure 5 shows the ESV spatial distribution from 2000 to 2020 in the study region. The ESV reduction region was mainly situated at the estuary of the Yellow River and the eastern coastal wetland region, and the ESV expansion region was predominantly present in the northern part of the study region, where the wetland was transformed into a water body from 2000 to 2020. Although the area of the cropland decreased significantly and most of it was transferred to the artificial surface, the ESV provided by its unit area was small. As a result, the ESV of the cropland and the artificial surface still kept below CNY 5 hundred thousand.

3.3. Multi-Scenario Land Use Simulation

The future projections showed that the cropland was still the largest LUT for YRD in 2030, followed by water body, wetland, and artificial surface. Most of the wetland and water body were distributed in the northern and eastern regions. Grassland, bare land, and forest combined account for less than 0.27% of the entire study region.

Regarding the temporal change in the LUCC in 2030 (

Table A4. Land use area in 2020 and multi-scenario.

	Artificial Surface	Cropland	Forest	Grassland	Bare Land	Water Body	Wetland
2020	61,447.37	453,630.25	0.36	1899.45	482.35	138,029.81	111,630.38
ND	70,376.58	431,099.46	0.27	1633.68	462.33	147,339.54	116,204.67
CP	61,578.27	455,051.79	0.27	1598.22	461.07	137,130.12	111,296.79
UCD	72,782.28	428,807.52	0.27	1631.16	462.24	147,288.51	116,144.55
ED	67,695.39	433,402.83	0.45	1638.99	460.71	147,511.89	116,406.27

), the wetland area increased significantly from 2020 to 2030 under the natural development (ND), urban construction and development (UCD), and ecological development (ED) scenarios, with an expansion of 4.10% for the ND scenario, 4.04% for the UCD scenario, and 4.28% for the ED scenario. In addition, the wetland area decreased by 0.3% under the cropland protection (CP) scenario at the same time. The results indicated that the cropland was the primary expansion source of wetland in this study area. The water body-changes trend was the same as for the wetland, and increased by 6.74% in the ND scenario, 6.71% in the UCD scenario, and 6.87% in the ED scenario; it also decreased by 0.65% in the CP scenario. Only in the CP scenario was the cropland area increased, growing by about 0.31%. The grassland and bare land decreased dramatically under all scenarios. In the context of the LUCC spatial change in 2030, the growth of the water body and wetland was mainly concentrated near the Yellow River estuary in the ND, UCD, and ED scenarios (Figure 6). The reduction in the cropland has caused the interface between cropland and wetland, as well as the water body shifting far away from the shoreline. The future expansion of the artificial surface was mainly concentrated in the central construction zone under the UCD scenario. The ED scenario was found to be the most suitable for regional development.

3.4. Characteristics of Changes in ESV under Multi-Scenarios

The total ESVs in 2030 in YRD under NP, UCD, CP, and ED scenarios were CNY 39.74 billion, CNY 37.56 billion, CNY 39.71 billion, and CNY 39.8 billion, respectively, with the change rate of 5.26%, 5.18%, −0.51%, and 5.42% compared with the total ESV in 2020, respectively. In the midst of them, the ND, UCD, and ED scenarios had a facilitating effect on the ESV improvement, and all of them had an increase of around 5%; correspondingly, the CP scenario showed a decrease in ESV by about −0.51% (

Table A5. ESV (10⁶ CNY) for different types of land use in different scenarios after simulation.

	Cropland	Forest	Grassland	Water Body	Wetland	Artificial Surface	Bare Land	Total
2020	2748.32	0.01	34.62	26,195.00	8773.52	0.00	0.15	37,751.62
ND	2611.82	0.01	29.78	27,961.78	9133.04	0.00	0.14	39,736.56
CP	2756.93	0.01	29.13	26,024.26	8747.31	0.00	0.14	37,557.77
UCD	2597.93	0.01	29.73	27,952.10	9128.31	0.00	0.14	39,708.22
ED	2625.77	0.01	29.87	27,994.49	9148.88	0.00	0.14	39,799.17

).

In the ND scenario, the total ESV increased by CNY 1.98 billion from 2020 to 2030, primarily due to the conversion of cropland to water body and wetland. The decrease in ESV resulting from the conversion of cropland to artificial surface, water body, and wetland was a total of CNY 0.14 billion. In the CP scenario, the total ESV experienced a reduction of CNY 0.19 billion from 2020 to 2030, primarily due to restrictions on cropland transfer, leading to a decrease in water body and wetland areas. The cropland ESV experienced an increase of CNY 8.61 million, while the overall decrease in ESVs of water body and wetland amounted to CNY 0.2 billion. In the UCD scenario, the total ESV experienced an increase of CNY 1.96 billion, mainly due to the same reasons as in the ND scenario. Even though the expansion of the artificial surface caused a decrease in ESVs, it was much smaller in comparison to the increase brought by water body and wetland. In the ED scenario, the total ESV showed the largest growth, reaching a total of CNY 2.04 billion from 2000 to 2030, with the reason for the increase being the same as in the ND scenario. The ESVs of water body and wetland also recorded the largest increase in the ED scenario, with CNY 1.8 billion and CNY 0.38 billion, respectively. The ESV of the forest increased by 25% for the first time, while cropland and grassland exhibited the smallest ESV decrease, with −4.46% and 13.71%, respectively. In the ED scenario, not only wetlands were well protected and developed, but other land uses were also reasonably developed.

Figure 7 illustrates the distribution of the ESV across varying scenarios in 2030. With the exception of the CP scenario, the central area of the study region experienced a significant increase in the ESV range of CNY 0.5 to 1 million for other scenarios. The distribution of ESVs under the EP scenario exhibited a significant reduction, and the ESVs shifted from the range of CNY 2–3 million to the range of CNY 1–2 million. The ESV supplied by the water body on a per-area basis was significantly greater than that of the wetland, which predominantly explained this alteration. Additionally, the transformation of the water body into wetland led to an ESV decrease in particular regions of this area. These findings demonstrated that the advancement of diverse ecological service functions depended on a well-planned allocation of varying land use categories in the region.

3.5. Hot Spot Assessment of ESV across Different Scenarios

To assess the ESV spatial clustering, the spatial statistical tool of ArcGIS was utilized to examine hot spots and cold spots within the study area. In 2020, ESV hot and cold spots were widespread and highly clustered, with hot spots predominantly located in the eastern and northern water body region, exhibiting a confidence level of nearly 99%. Meanwhile, the concentration of cold spots was in the southern region and western region. In comparison to 2020, the four scenarios that were modelled for 2030 were typically marked by a slight decrease in the range of hotspots’ distribution and a minor increase in cold spot regions (Figure 8).

Under the ND, UCD, and ED scenarios, a few small hot spots have emerged in the study region. The emergence of this phenomenon is the result of the increasing wetland area, which has led to a rise in ESV in this region. It was especially important to note that the 99% confidence level cold spots in the study region’s southern, western, and central sectors increased significantly under the UCD scenario, suggesting that the dispersion of artificial surface was mainly concentrated in this region, with a relatively minor impact on the wetland. The decrease in hot spots in the study region was caused by two factors: (1) the sources of wetland and water body have been in decline under the UCD scenario; and (2) some water bodies have been transformed into wetland, resulting in a decline in ESV and a weakening of aggregation.

4. Discussion

4.1. The Impact of Past and Future Trajectories of LUCC on ESVs in the YRD

During the study period, the wetland ESVs in the study area showed a decrease in total, consistent with the findings of previous research [55,56]. From an ESV contribution perspective, the water body was the highest contributor, followed by the wetland, together contributing nearly 90% of the ESV in the study region. Among them, the wetland alone accounts for over 22% of the ESVs. Wetland ecosystems can provide high ESVs [57] and serve several distinctive ecological functions, including ecological balance maintenance, biodiversity preservation, and water resource conservation [58]. The conversion of wetland to water body had increased regional ESVs. However, the wetland provided significantly higher ESVs than other LUTs in habitat and cultural services. The shift in wetland ESVs was closely related to wetland area changes [59]. In recent years, there was a fluctuation trend between the decrease and increase in the wetland area in YRD. The initial decrease in wetland area was primarily caused by the prolonged river cutoff time and expanded river cutoff range of the Yellow River, which resulted in a reduction in water flow from the river upstream. This resulted in wetland degradation [60]. In recent years, numerous wetland ecosystems, including intertidal zones, have been created due to the continuous expansion of aquaculture ponds and coastal salt production in the YRD. However, this has also resulted in a decline in the overall area of wetland ecosystems [61]. From the standpoint of land utilization transformation, cropland served as the principal source for the transfer in the wetland, and wetland predominantly transferred out to a water body. The amount of forest, grassland, and bare land within this study region was limited and had a relatively minor impact on alterations to the wetland. The artificial surface and wetland were spatially separated; thus, the artificial surface had a lower effect on the natural wetland environment.

In all simulation scenarios, the artificial surface increased at varying scales, which was in line with the findings of previous research [62]. The wetland area and ESVs decreased in the CP scenario due to the restricted alteration of cropland to other LUTs. In the ND, UCD, and ED scenarios, the wetland area exhibited a similar growth rate, and the increase in ESV showed a homogeneity in different scenarios. This provided additional evidence that the conversion of cropland in the study area significantly affected wetland. According to [63,64], a surge in population and economic growth has driven the rising need for artificial surfaces, further augmenting the expansion rate of the artificial surface in the UCD scenario. Under the ED scenario, the ecosystem was safeguarded and enhanced appropriately, leading to the greatest ESVs compared to other scenarios. This outcome was in line with the discoveries made by previous studies [65]. The distribution of ESV hot spots and cold spots revealed that a majority of ESV hot spots in the designated research zone were concentrated close to the water body. The distribution of hot spots in the wetland concentrated area was not pronounced, as the water body covered a larger area, and it can provide a higher ESV on a per-area basis compared to the wetland. However, when it comes to supply services and cultural services, wetlands played an irreplaceable role compared to other LUTs.

4.2. Policy Implications and Suggestions on YRD

Policy has a significant impact on land utilization. The results of the study showed that the area of aquaculture in YRD has expanded by 8.4 times in recent years [66], and the aquaculture area is mainly developed from wetland ecosystems such as mudflats, and the wetlands in the study area have been greatly damaged. It is important to note that policies such as ‘Opinions of Further Accelerating the Development of the Aquaculture Industry’ have played a significant role in this outcome. The implementation of policies such as the ‘National Wetland Conservation Action Plan’ has significantly mitigated the expansion of aquaculture in the study area, and some areas have begun to reduce aquaculture and restore wetlands that have been effectively protected and expanded [67].

Based on the obvious correlation between wetland ESV changes and other land utilization in the study area, wetland conservation policy in the future can be formulated in the following three aspects: (1) reduce the transfer of wetland to water bodies under the premise of securing water for domestic and production use; (2) since bare land can provide a negligible ESV per unit area, bare land will be selected for land utilization development instead of wetland, so as to reduce the degradation of wetland; and (3) aquaculture is a significant industry in the study area, occupying a large number of wetlands. To promote green development, it is necessary to implement relevant policies, transform the industry from high-speed to high-quality development, and reduce the number of aquaculture areas to restore wetland areas.

4.3. Limitations and Future Trends

The adjusted equivalence-coefficient method was used to estimate ESVs in YRD. Parameters such as grain production, the area of major grains, and the market price of major grains were selected to adjust the equivalence-coefficient method. Compared with other methods [68,69,70], this method requires fewer parameters and is easier to obtain. Therefore, it can calculate the ESV of global wetland ecosystems under the condition of limited data. However, this method is susceptible to inflation and the market economy [71]. The large fluctuations in market prices of major foodstuffs over a long time series lead to significant differences in the assessed ESVs of regions with similar land composition [72]. The accuracy and applicability of ESV assessments should be improved by incorporating various factors into the model.

The Markov–FLUS model can generate adaptive probabilities for different LUTs based on current land utilization data and six driving factors (DEM, road, temperature, precipitation, slope, and railways). It combines with an adaptive inertial competition mechanism to predict the distribution pattern of land utilization [73]. The unique mechanism can predict land use from both quantitative and spatial perspectives [74]. However, the accuracy of the model is affected by the combination of driving factors under different policies. For instance, the exclusion of ecosystem services in global policy frameworks is a reason for differences in wetland protection [75]. Additionally, the model’s accuracy decreases in regions other than the YRD due to varying driver choices under different policies. Outside the YRD region, the model’s accuracy decreased due to the use of different drivers resulting from policies. In the future, a global unified policy framework that combines multiple drivers could enhance the model’s predictive accuracy. This approach would be favorable for exploring the commonalities of wetland development across various regions.

5. Conclusions

This study calculated the ESV of major LUTs by combining the equivalent coefficients of ES after adjustments. The Markov–FLUS model was used to explore the changes in land utilization and ESV in the study area in 2030. The results indicate that wetland changes in the Yellow River Delta region were more frequent during the study period, with a trend of decrease followed by an increase. The ESV of wetland decreased by CNY 525 million. In total, 80.08% of the decrease in ESV in wetland was attributed to the conversion of wetlands into water body, while 69.1% of the increase in ESV in wetland was due to the conversion of cropland to wetland. The main LUTs affecting the change in wetland ESV in the study region were cropland and water body. The ED scenario was the most appropriate development mode for the wetland ecosystem. In the ND, UCD, and ED scenarios, the wetland area and ESV growth rate exceeded 4%. The wetland area and ESV exhibited a decline trend under the CP scenario, further suggesting that cropland is a major source of wetland transfer. This also suggests that the development of appropriate land policies is an effective means of achieving wetland conservation.