Abstract

With the rapid development of the economy, the land use/cover change (LUCC) in the Greater Bay Area (GBA) has undergone tremendous changes, which have had directly negative effects on ecosystem functions and services. The development of sustainable land use strategies to quantitatively evaluate ecosystem services is required. Based on multitemporal land use data (2005, 2010, 2015, and 2020), the equivalent coefficients table method was used to assess the ecosystem service values (ESVs), and the impact of LUCC on ecosystem services was analyzed. A future land use simulation (FLUS) model and multiscenario simulations were employed to predict land use change in 2030. Our results indicated that the loss of ESVs decreased by 14.29 billion yuan from 2005 to 2020. The spatial distribution of the high-value ecosystem services was concentrated around the peripheral area of the northern regions in the GBA, and those areas had less land use development and human activity. Compared with those in 2020, the total ESVs of the business-as-usual (BAU) scenario, socioeconomic development (SED) scenario, and cultivated protection priority (CPP) scenario in 2030 decreased, while they increased in the ecological protection priority (EPP) scenario. In the CPP scenario, regulating, supporting, provisioning, and cultural services increased slightly, but they decreased in the other scenarios. The patterns of LUCC were the main reasons for the decrease in ESVs, such as the loss of land with high ecological value. Additionally, a four-quadrant analysis is introduced to determine which land use simulation will be expected to be adopted by the government. The findings of this study provide valuable information for decision-making and policy development in the coastal zones and for the sustainable management of ecosystems.

1. Introduction

With the rapid expansion of industrialization and urbanization, the population has grown to unprecedented levels, and natural resources are constantly being exploited [1], leading to the increasing encroachment of natural ecosystems in the last century. To counteract the increasing global pressure on ecosystems and improve the supply of natural resources [2, 3], the concept of ecosystem services (ESs) has been proposed, and research on ecosystem services has gained widespread attention [4, 5]. The term “natural services” was first used academically by Westman [6] to measure “How much are nature’s services worth?” and, simultaneously, the term “ecosystem services” was officially introduced by Ehrlich and Mooney [7]. Two pioneering studies on ecosystem services and natural capital were published in 1997 [8, 9]. Since then, ecosystem services have sparked a discussion, and the evaluation of ESs has become a central issue in decision-making and land management [4].

Ecosystem services refer to the natural environmental conditions and functions that are formed and maintained by the ecosystem and developed for human survival. These services are the benefits obtained by human beings, directly or indirectly, from ecosystem functions [9, 10], and they are divided into four categories: provisioning, regulating, supporting, and cultural services [11, 12]. Seventeen types of ecosystem service coefficients for 16 biomes were first proposed by Costanza et al. [9, 13]. These new methods can be applied for evaluating ecosystem service value (ESV), which is defined as the estimation of the marginal value of ecosystem services [14]; uncertainties in these coefficients have been studied [15]. In 2014, the evaluation of global ESVs was updated by Costanza et al. [10] based on the same methods applied in 1997. Since then, most scientific studies have referred to or improved Costanza’s methods to estimate ESVs on different scales, such as the provincial scale, regional scale, national scale, and global scale [2, 16]. The Chinese ESV coefficient was revised by Xie et al. [12] to include 11 service types according to China’s ecological characteristics. This classification is adopted in the present study, and the new revised coefficient is used to estimate the ESV that refers to local characteristics.

The sustainable development of the social economy depends on the sustainable supply of ecosystem services [17]. However, the world has faced environmental pollution and the increasing exploitation of ecosystems [18]. Some severer cases of ecosystem destruction associated with land use systems and different functions in specific areas have reached or exceeded the ecosystem carrying capacity and even become irreversible. Land use change/cover (LUCC) is a crucial way in which human activities affect ecosystems [5, 19]. LUCC alters the structure and function of ecosystems, and the supply of ecosystem products and services influences ecosystem patterns and processes [20, 21]. High-intensity human activities, including urban expansion and land development, have enhanced the LUCC process worldwide [22, 23], which has impacted ecosystems and exacerbated the loss of ESV. The ecological benefit losses greatly impacted ecosystem services and caused imbalances [19]. Because there are time lags between land use change and ecosystem responses [24], land use legacies affect ecosystem service provisions. Therefore, an understanding of the impact of LUCC on ecosystem service helps decision-makers to minimize the negative consequences of LUCC or to relieve the ecological pressure through targeted management measures. Analysis of the past dynamics of LUCC [25] can help researchers anticipate the potential effects in future land use simulations and reveal tradeoffs [2]. Consequently, research on the effect of LUCC on ecosystem services has attracted increasing attention [1].

The spatial distribution of ESVs in response to LUCC has been studied [26]. Huang et al. [27] used the InVEST model to analyze the impact of LUCC on ESV. The ESV loss caused by LUCC exhibited significant spatial heterogeneity due to the spatial difference in land use [28]. In recent decades, a series of land use simulation models with good accuracy have been developed, such as the logistic-CA [29] model, cellular automata-Markov (CA-Markov) model [30], conversion of land use and its effects at small regional extent (CLUE-S) model [31], and artificial neural network-cellular automata (ANN-CA) model [32]. Subsequently, the future land use simulation (FLUS) model based on the traditional CA model was developed by Liu et al. [33] and was shown to have higher simulation accuracy than other models. The FLUS model has become a major tool in ES research and can be implemented on various spatial scales and regions to analyze urban growth boundary simulation [17, 19, 34], flood risk assessments [35], typical developed areas [5], mountain regions [36], and metropolitan regions [37]. Although socioeconomic and geographic conditions have been included as factors in the FLUS model [38], few studies have considered the background climate factor, which has significant effects on LUCC. The climate factor is therefore addressed in this study, and POIs are also selected.

The Greater Bay Area (GBA) is one of the major coastal areas in China and a crucial pivotal region connecting the development of the Belt and Road. The region has experienced complex land use changes due to rapid urbanization and the high-intensity development of human activities [39]. Therefore, it is a good case study for assessing the impacts of LUCC on ecosystem services over long periods for the following reasons. A large proportion of the population (https://www.stats.gov.cn/) greatly benefits from providing highly diverse ecosystem services typical of coastal regions. In addition, coastal ecosystems and their related ecosystem services are highly susceptible to LUCC and climate change. In this study, we focused on the GBA. The main objectives of this study were to assess the ESVs, predict the future distribution of LUCC using multiscenario simulation analysis based on the FLUS model, and estimate the variations in ESV resulting from the impacts of LUCC in the GBA.

The innovations of this study are as follows. First, we used a future land use simulation model to evaluate the ESVs. The advantages are that the FLUS model can process the complex competition and interactions among the different land use types. Second, we focused on assessing the effect of social and economic factors, but climate, facilities, and other environmental factors in land use prediction were also considered. Finally, we developed a four-quadrant analysis to determine which land use simulation is expected to be adopted by the government and, thus, provide scientific decision-making support for sustainable land use and ecosystem management in the GBA.

2. Materials and Methods

2.1. Study Area

The GBA is located at 111°12′E to 115°35′E and 21°25′N to 24°30′N (Figure 1). The GBA lies in the fragile and complicated ecological environment of the southeast coastal region in China. It extends over an area of 56000 km2 with various topographies, such as mountains (low mountains and hills), plains (alluvial plain and delta plain), and islands. The GBA has a subtropical monsoon climate with an annual average temperature of approximately 21.4–22.4°C, an annual average precipitation of 1808 mm, and an annual average sunshine hours of 2000 h (https://data.cma.cn/).

The GBA is composed of nine prefecture-level cities of Guangdong Province and the two special administrative regions (Hong Kong and Macao) in the Pearl River Delta Region. Like other coastal regions, the GBA is an important economic zone, with a relatively concentrated population totaling 86.40 million people and increasingly concentrated industries, with a total GDP of 11.3 trillion yuan in 2020 (https://www.stats.gov.cn/).

With the rapid development of urbanization and industrialization around the GBA, increasingly frequent human activities, rapid urban sprawl, and gradually increasing regional inequality, land use has changed significantly and led to excessive development and utilization activities [5, 39]. At the same time, the ecological environment has been put under increasing pressure, resulting in continuous declines in the ecosystem service functions. Thus, it is necessary to assess changes in ESVs and simulate ESVs in response to LUCC in the future.

2.2. Data Sources and Land Classification

We collected historical land use data (2005, 2010, 2015, and 2020) and driving factor data (Table 1), which were obtained from the Resources and Environment Science and Data Center, Chinese Academy of Sciences (30 m resolution) (https://www.resdc.cn/). The land use data were used for LUCC multiscenario simulation analysis and ESV assessment. In this study, wetlands were listed separately [40], so that the types of LUCC were reclassified into seven categories: farmland, forestland, grassland, water body, wetland, built-up land, and unused land (Table 2). The dynamic information (Table3) on land use over a 15-year period was calculated using ArcGIS software. Then, we obtained a transition matrix that represented the quantitative transition between different land use types.

According to the results of previous research and considering the availability of data [17, 19], the spatial driving factor dataset for LUCC multiscenario simulations was selected, as shown in Table 1. The digital elevation model (DEM), with a 12.5-meter resolution, served as the basis for data on terrain heights and the calculation of slopes and aspects. Data on the soil characteristics (e.g., clay content, silt content, and sand content), temperature, precipitation, gross domestic product (GDP), and population were provided by the Data Center for Resources and Environmental Sciences of the Chinese Academy of Sciences. Socioeconomic data, including town data, residential data, railway data, traffic data, water system data, and facilities data (POI), were collected from open-source data retrieved from OpenStreetMap (OSM). A unified coordinate transformation was performed with ArcGIS. Socioeconomic data and facility data were calculated by Euclidean distance analysis. The spatial reference of the WGS_1984 coordinate system was constructed for matching land use data and driving factor data.

Statistical datasets on ecosystem services, including grain output per unit area and total sown area of farm crops, were obtained from the National Bureau of Statistics, China Statistical Yearbook, Guangdong Statistical Yearbook, and the statistical yearbooks of various cities in the GBA. In addition, some policy documents were obtained from government reports for the GBA.

2.3. Ecosystem Service Valuation

The equivalent coefficients table method, which is one of three main ESV evaluation methods [41], is more intuitive and suitable for the assessment of ESVs because it has fewer data requirements on the regional and global scales. The values per unit area of ecosystem service in this study were revised in reference to previous studies [42] based on the method of regional corrections on cropland. The formula is defined as follows:where is the revision factor. and are the annual average food production of cropland in the GBA and China, respectively. represents the revision equivalent factor for land use type i, while represents the standard equivalent factor which is quoted in Xie et al.[12] for land use type i of farmland, forestland, grassland, water body, wetland, built-up land, and unused land.

Because the food production function of cropland determined the equivalent ESV coefficient, the economic value of the average natural food production of cropland per hectare per year was a critical indicator. Generally, natural food production is equal to 1/7 of the actual food production [43]. We calculated the annual average food production of cropland in the GBA and China in 2020 and obtained values of approximately 5500 kg/hm2 and 5600 kg/hm2, respectively, and then determined the standard equivalent factor. Based on the average net profit per unit area of the major foods (rice and soybeans), the average ESV of one equivalent value for the GBA was 2215.61 RMB yuan hm−2·yr−1. The unit area values of some land use/cover types were calculated based on this equivalent coefficients’ table, in which the service value coefficient of built-up land was given as 0 [19]. Thus, the value tables of the ecosystem services per unit area of different land use types were obtained (Table 4). The ESV, including each land use type per hectare, was calculated, based on the following equations:where , , and refer to the ESV of the type i ecosystems, the ability category of , and the total ESV, respectively. is the area (in hm2) of the type i ecosystem, and (yuan/hm2) denotes the value coefficient for land use type i and ecosystem service function type j.

2.4. Spatial Simulation of Land Use Change

LUCC simulation models are effective and reproducible tools for analyzing the causes and consequences of land use patterns assuming various scenarios [33]. The FLUS model was first applied to simulated LUCC in mainland China from 2000 to 2010 [33]. The FLUS model was developed by introducing a self-adaptive inertia and competition mechanism into the cellular automata (CA) model to predict land use changes [33, 44]; the Markov chain (MC) approach was adopted to predict the scale of land use and multilayer artificial neural networks (ANNs) were used to estimate the probability of occurrence of different land use types. In this study, we used GeoSOS-FLUS software to develop both a geographic simulation and optimization system (GeoSOS) and a FLUS model. GeoSOS-FLUS software is a powerful tool to predict multiple LUCC scenarios by coupling human and natural effects. It can be used to make LUCC simulations more convenient and efficient [36]. It is important for decision-makers to predict future land use changes and dynamic LUCCs for multiscenario simulations. The general structure of the FLUS model is illustrated in Figure 2.

2.4.1. Prediction Type of Land Use

The pixel’s future type was predicted using the Markov chain model. The transfer probability matrix of land use types [45] was computed from the previous state at time t to the land use state at time t + 1 [46]. In this study, ten years was selected as the initial forecast unit to obtain the land use change, and the land use type in 2020 was predicted. The Kappa coefficient and figure of merit (FOM) were used to evaluate modeling accuracy to better compare our simulation results in 2020 with those of previous studies in 2020. Then, the land use type in 2030 was predicted by inputting historical data into the FLUS model. We also predicted multiscenario simulations in 2030. The model stopped outputting the results when it achieved the iteration target; otherwise, the land use demand was not met, and the model iterated 300 times in the forecast year.

2.4.2. Land Use Change Simulation Using Cellular Automata (CA)

The CA simulation was implemented in two steps. Step 1: an artificial neural network (ANN) is used to train and estimate the probability of occurrence of each specific pixel of land use types. Step 2: a reliable self-adaptive inertia and competition mechanism is designed to identify the competition and interactions among different land use types. The change in land use type was obtained through the high transformation of the probability of occurrence during the CA iteration; for example, a specific pixel either retains the current land use type or is transferred into another land use type. A schematic framework of the CA model is presented in Figure 2.

Artificial neural networks are typically used to capture the nonlinear relationship between factors that are dependent on independent variables. It fits the complicated relationships between input data and predicted targets through a number of learning-recall iterations [47] and is widely used in studies of land use change. In this study, DEM, aspect, slope, soil type, temperature, precipitation, potential crop, NDVI, ecological protection, GDP, population, distance to town, distance to residential area, distance to railway, distance to traffic, distance to water, and POI were selected as potential driving factors of LUCC (Table 1). The land use transfer matrix was determined, and the probability of occurrence of land use types was estimated. The expected output and the current output were compared by randomly selecting the primary weights [48] and then predicting the futures land use type. In general, an ANN is composed of three layers: an input layer, a hidden layer, and an output layer network. Its formula is as follows:where is the probability of pixel conversion from the original land use type to the target land use type k at iteration time t; is an adaptive weight between the input layer and the hidden layer; is the excitation function from the hidden layer to the output layer; and is the signal received by in the hidden layer n from the pixel p at time t.

The ANN model was developed to establish the relationship between the probability of occurrence surface for a specific land use type and the given spatial factors [33]. Whether the pixels of land use type are converted to a specific land use type not only is consistent with the probability of occurrence but also is related to other variable components accounting for different development statuses over the prediction period. Therefore, the probability of occurrence is combined with the conversion cost, neighborhood condition, and the competition and interactions of different land use types to estimate the combined probability for each land use pixel.

2.4.3. Multiscenario Land Use Prediction

The CA model determines the pixel land use type according to the pixel s previous state at a microscale, the neighbor pixel state, and the land use transfer rules [17]. The neighborhood weight value in CA describes the neighbor interaction between different land use types [17, 49]. In this study, a 3 × 3 Moore neighborhood model was adopted, and the weight of the neighborhood of each land use type ranged from 0 to 1; the larger the weight is, the stronger the expansion ability of the land use type is [46]. The conversion cost is defined as the difficulty of converting from the current land use type to another land use type. The conversion cost values are in a binary matrix consisting of zeros and ones, where zero means that the transfer of the land use from the row to the column type is not allowed, while one means that the transfer is allowed [34].

With the LUCC characteristics and the current situation of social and economic development in the GBA, we considered four scenarios: business-as-usual (BAU), socioeconomic development (SED), cultivated protection priority (CPP), and ecological protection priority (EPP). The four scenarios on the land use transfer matrix in 2030 were calculated according to the transfer probability matrix of land use types from 2010 to 2020. GeoSOS-FLUS software was used to calculate the land use adaptability probability based on the land use data and driving factors of LUCC (e.g., DEM, aspect, slope, soil type, temperature, precipitation, potential crop, NDVI, ecological protection, GDP, population, distance to town, distance to residential area, distance to railway, distance to traffic, distance to water, and POI). Then, the spatial distribution data of LUCC for the BAU, SED, CPP, and EPP scenarios were obtained.

3. Results

3.1. Changes in Land Use/Cover
3.1.1. Temporal Analysis of LUCC

The proportion of land use/cover types is as shown in Table 3. Forestland was the predominant land use type in the GBA in 2005, 2010, 2015, and 2020, accounting for 54.86%, 54.39%, 53.85%, and 53.70%, respectively, and its proportion exhibited a continuous decline. Farmland accounted for 23.84%, 22.89%, 22.43%, and 21.88%, respectively, and these types also appeared to undergo a sustained slow decline in 2005, 2010, 2015, and 2020. The water area declined continuously. Moreover, the proportion of built-up land exhibited a continuously increasing trend, accounting for 11.59%, 13.34%, 14.27%, and 15.02%, respectively. Grassland and wetland tended to fluctuate from 2005 to 2020, initially decreasing, then increasing, and finally decreasing. The proportion of unused land was the smallest, and it decreased, accounting for only less than 0.03% of the GBA.

From 2005 to 2020, built-up land exhibited the largest change in the GBA, with a total net area increase of 1921.61 km2 (Table 3). The net area changes in farmland and forestland were larger, followed by those of water area and wetland. Finally, the change in net area of unused land was the smallest. The area of land use/cover types showed varying characteristics during different phases. From 2005 to 2010, the areas of farmland and forestland declined significantly, with decreases of 531.61 km2 and 260.79 km2, respectively, while built-up land had the most significant increase (957.66 km2). From 2010 to 2015, the forestland area showed the greatest decreases (306.96 km2). The area decrease of farmland and water areas was 254.32 km2 and 105.74 km2, respectively. The area of grassland increased (139.50 km2) and that of wetland slightly increased (2.68 km2). Built-up land exhibited a continuously increasing trend, with an area of 525.95 km2. From 2015 to 2020, the area of farmland decreased the most, followed by that of forestland, grassland, and wetland, while the water area increased.

The change rates of the land use/cover types showed varying characteristics during different periods (Table 3). In 2005–2020, farmland, forestland, and water area showed stages of slow decline, while unused land and wetland decreased more rapidly, accounting for 66.10% and 34.84%, respectively. However, built-up land increased the fastest (29.60%).

3.1.2. Spatial Patterns of Land Use/Cover Change

The spatial distribution of land use/cover types in the GBA is shown in Figure 3. The primary land use types were forestland and farmland in the GBA. Forestland was mainly distributed in the western, northwestern, and eastern regions. Farmland was found in the central, southwestern, and eastern regions, including Guangdong, Foshan, Zhongshan, Jiangmen, and Huizhou. Built-up land was concentrated in the central region, which had a high level of urbanization, such as Guangdong, Shenzhen, Hong Kong, and Macao. Grassland was scattered throughout the GBA. Water areas were mainly distributed in the northern region of the GBA.

3.1.3. Transformation Patterns of LUCC

The total conversion area of land use/cover types occurring in the GBA from 2005 to 2020 reached 6242.17 km2. The transformed land use types were mainly farmland, forestland, water area, and built-up land, which accounted for 39.78%, 21.69%, 18.67%, and 15.36% of the total transformed area, respectively. These land use types were mainly converted into built-up land, farmland, forestland, and water area.

To illustrate where the main conversion of land use occurred between 2005 and 2020, Figure 4 shows the spatial distribution of the main land use changes. Built-up land showed a significant increase of 45.00% and mainly occupied farmland and forestland, with area increases of 1502.2 km2 and 723.07 km2, respectively. Farmland declined by 2483.32 km2, and this area was mainly converted into built-up land, water area, and forestland. The water area was mainly changed into farmland and built-up land, with areas of 553.41 km2 and 492.31 km2, respectively. A small amount of wetland and grassland was transformed and converted to forestland, built-up land, and farmland. The proportion of unused land was the lowest and it was mainly transformed into wetland and forestland. Finally, the transformations of land use/cover types were caused by human activities, one of the most significant driving factors.

3.2. Changes in Ecosystem Service Values
3.2.1. Temporal Analysis of Ecosystem Service Values

With the high-quality development of industrialization and urbanization in the GBA, constructed land is constantly expanding in urban and rural regions, and the expansion of construction related to urbanization, industry, and transportation has resulted in the occupation of forestland, water area, and farmland, representing the primary drivers of the reductions in ESVs. The value of the ecosystem services in the GBA in 2005, 2010, 2015, and 2020 is shown in Table 5 and Figure 4. The total ESVs were approximately 558.14 billion yuan in 2005, 550.69 billion yuan in 2010, 544.48 billion yuan in 2015, and 543.85 billion yuan in 2020, exhibiting a trend of continuous decline. The value of the ecosystem services provided by different land types was different. Forestland contributed the most, accounting for more than 75% of the total ESV, followed by water area, which accounted for approximately 20%. Among land use/cover types, the ESV losses of forestland and water area were significant, with decreases of 8.95 billion yuan and 4.12 billion yuan, respectively, followed by those of farmland and wetland, which decreased by 92.60 million yuan and 59.40 million yuan, respectively. The ESV loss of unused land was the lowest, with a decreasing value of 0.05 million yuan. In contrast, the ESV of grassland increased by 30.05 million yuan, accounting for less than 2% of the total ESV.

3.2.2. Spatial Changes in Ecosystem Service Values

Figure 5 also shows the spatial distribution of the ESV change. The high-value areas of ecosystem services were concentrated in almost every peripheral of the northern regions in the GBA; those areas had less land use development and human activity and were mainly located in Zhaoqing, Huizhou, and Jiangmen. Conversely, the low-value areas of ecosystem services were distributed in the central regions of the GBA, including Guangdong, Dongguan, Zhongshan, Hong Kong, and Macao; those areas were subjected to intense human activity and frequent land use change, especially the expansion of constructed land.

To determine which type of land transition that had a major effect on the total ESV change, we produced a transition matrix of ESV change between 2005 and 2020. The increasing and decreasing changes in ESV are shown in Figure 6, and, in general, the loss of ESVs was greater than the increasing value of ecological services. The areas of decreased ESVs were distributed throughout the entire GBA, and the greatest losses were distributed in Foshan city, while the areas of increased ESVs were scattered in the central region of the GBA, such as Zhuhai city.

3.2.3. Changes in Different Ecosystem Service Values

From 2005 to 2020, the regulating services contributed the most, accounting for the highest proportion (>70%) of ecosystem services. They were followed by supporting services, accounting for approximately 20%, while provisioning services and cultural services accounted for the lowest proportion of the total ESV, approximately 10%.

Regulating services, supporting services, and cultural services showed a continuous declining trend, while provisioning services showed first a minor decrease and then a slight increase, as shown in Figure 6. From 2005 to 2020, the primary ecosystem service types, including regulating services, supporting services, provisioning services, and cultural services, experienced ESV losses of 10.80 billion yuan, 2.38 billion yuan, 56.60 million yuan, and 54.80 million yuan, respectively.

In terms of the ecosystem service subtypes, the regulation of water flows contributed most to the total ESV, followed by climate regulation, maintenance of soil, biodiversity conservation, air quality regulation, waste treatment, and aesthetic inspiration. In contrast, maintenance of the nutrient cycle contributed the least. The ESV of the water supply exhibited an increasing-decreasing-increasing trend. The regulation of water flows experienced an initial decrease and then increase. The ESVs of all other ecosystem service subtypes decreased during all periods (Figure 7).

3.3. Future Changes in LUCC and EVS
3.3.1. LUCC in 2030

The ROC results are shown in Table S1. The ROC values were larger than 0.74, indicating that there was a good fit for each land type. Furthermore, the Kappa coefficient and the overall accuracy of land use simulation under the BAU scenario were 88.29% and 92.58%, respectively, which were used to assess the accuracy of the simulation results by comparing the differences between ground-truth land use and simulated land use in 2020. The FOM was 0.2535, which is a high accuracy value. In this study, multiscenario simulations of land use were carried out, including BAU, SED, CPP, and EPP (Figure 8). The results showed that the change in land use types was consistent with the spatial regional distribution (e.g., the distribution of built-up land) without any large deviation in spatial location, but there were differences in local areas. Compared with 2020, the change area of each land use type was considerable in 2030, and the degree of overall land use fragmentation was more significant.

To identify the main conversion types, we extracted the transformation of land use types from 2020 to 2030 from the land transfer matrix based on the BAU, CPP, SED, and EPP scenarios (Table 6, Tables S2S5). The results showed that most of land use types underwent conversion, including farmland, built-up land, forestland, and water area.

In the BAU scenario, the areas of built-up land and grassland were 484.33 km2 and 80.75 km2, respectively, while the areas of farmland, water area, and forestland were 451.46 km2, 66.02 km2, and 48.50 km2, respectively. The increase in built-up land was mostly due to the occupation of farmland, and the change in land use was consistent with that represented in the land transfer matrix; therefore, built-up land increased sharply by 60.62%, while the farmland declined by 70.37%. Under the SED scenario, to continuously meet rapid economic and social development, the demand for land increased. There was a rapid development of urban sprawl, which expanded from the periphery of the original urban boundary, and built-up land significantly increased (total area 459.40 km2), while the area of grassland increased only slightly. However, the forestland rapidly decreased (361.71 km2) followed by decreases in water area, farmland, and wetland and a slight decline in unused land. The expansion of built-up land resulted from the conversion of forestland, farmland, and water area, and these results were consistent with the land transfer matrix. In the CPP scenario, farmland played an important role in maintaining national food security. Within the framework of the national security system, the goal is to control the least cultivated areas and protect the basic farmland preservation area. The area of farmland increased the most (64.58 km2) and that of built-up land increased only slightly (29.37 km2). Forestland decreased the most, with a total area of 154.08 km2. Forestland was converted to expand farmland, and grassland was converted to urban construction in the land use transfer matrix. In the EPP scenario, there was less pressure on ecological land, such as forestland, water area, and grassland. Compared with 2020, the area of farmland declined, with a total reduction of 451.42 km2, while other land use types, including forestland, grassland, water area, and wetland, showed increases in area of 47.09 km2, 34.51 km2, 19.70 km2, and 2.79 km2, respectively. The ecological land provided high-quality service functions, although the area of farmland decreased the most and the area of built-up land increased the most.

3.3.2. Changes in ESVs in 2030

The changes in ESVs in the GBA under multiple scenario simulations in 2030 are presented in Table 7. The ESVs showed a downward trend in the BAU, SED, and CPP scenarios, with a total decrease in value of 2.18 billion yuan, 1.85 billion yuan, and 6.90 billion yuan, respectively. The ESV showed an upward trend in the EPP scenario, with a total increase of 1.43 billion yuan. In terms of the type of ecosystem service, in the BAU, SED, and CPP scenarios, the values of the regulating services, supporting services, provisioning services, and cultural services decreased compared with those in 2020, while they increased slightly in the EPP scenario.

The ESVs in different simulated scenarios were significantly different due to different land use types (Tables S2S5). In the BAU scenario, the decline in the total ESV was attributed to the decrease in the ESVs of farmland, with a loss of 2.57 billion yuan, and the change in other land use types led to a small increase or decrease in ESVs. In the SED scenario, the loss of ESVs was derived from the increase in the area of built-up land and decrease in the area of farmland, with total values of 3.47 billion yuan and 3.27 billion yuan, respectively, followed by those of forestland, grassland, and wetland, with decreases in ESVs of 0.13 billion yuan, 5.9 million yuan, and 1.5 million yuan, respectively. However, an increase in ESVs resulted from a slight increase in water area, with a total value of 1.4 million yuan. In the CPP scenario, the decrease in farmland was remarkably restrained, and the loss in ESVs was also significant, which increased 1.03 billion yuan, and the other land use types (e. g., forestland and water area) were converted to grassland and wetland, which increased the values of ecological services. In the EPP scenario, the increasing ESVs originated from the contributions of water area, forestland, and grassland, with the total values of 0.96 billion yuan, 0.54 billion yuan, and 0.22 billion yuan, respectively. This is because the protection of the ecological environment led the ESVs to significantly increase due to limitations on development activities and the partial restoration of damaged habitats.

4. Discussion and Conclusions

4.1. Discussion
4.1.1. Impact of LUCC on ESV

The change in land use was caused by the interaction between human activities and the ecological environment and led to the change in the ecological environment. The changes in patterns of land use and land cover types were affected by factors such as population growth, urbanization, land use policy, biology, climate change, and soil [5052], which included the natural environment, local conditions, socioeconomics, and policy orientation, and were important driving factors of the geographical environment and socioeconomic factors that affected the quality and quantity of land use change. For instance, changes in land use types were triggered by population growth and urbanization [53], such as the conversion of agricultural land to settlements. Therefore, the impact of LUCCs on ecosystem services is very complex. These changes significantly impacted ecosystem services and led to a change in the provision of functions [54, 55], and the losses in ESVs caused by land use change were experienced on different regional scales [53]. Therefore, it has been suggested that land use patterns, land management, and land use planning also affected ESV change. It is important to quantify and evaluate the impact of land use/cover change on ESV and its changes, and this is necessary for the sustainable development of land resources and ecological environment protection.

In line with other studies, our results indicated that rapid urban sprawl resulted in a vast loss of ESV via land use changes in China, especially in coastal regions [17, 19]. The GBA is a typical representative of coastal regions with fragile ecological environments. Land use types tended to be out of balance, which had an impact on ESVs under the dual effects of climate change and human activities. In addition, in recent decades, intensive development activities, such as urban sprawl and the construction of urban land, have accelerated profound LUCCs in the GBA. These changes have significantly impacted the ESVs. Hence, the impact of LUCCs on ecosystem services is ultimately a result of the relationship between land use development and ecological protection. We propose a four-quadrant analysis method for determining which land use simulation will be expected to be adopted by the government. In the four-quadrant analysis, the rapid expansion of urban construction indicates land use development, while less pressure on ecological land indicates ecological protection. The distribution of the multiscenario simulations is presented according to land use conversion and ESV transfer change (Tables S2S5) (Figure 9).

The positive and negative effects on ESVs were revealed by multiscenario simulations in the four-quadrant analysis, and the impacts on the landscape in terms of spatial expansion were different in different ecological regions. Changes in ESVs resulting from the decrease in farmland and built-up land were similar in the BAU and SED scenarios, which showed high land use development and low ecological protection in the HL quadrant. One reason is that farmland was occupied by the expansion of urban construction, and the ESVs exhibited a downward trend. A second reason is that legacy effects impacted the change in ESVs [56, 57]. In our study, we investigated the changes in ecosystem services in response to LUCC. The changes in the ESVs declined in the BAU and SED scenarios due to the continuous expansion of urban land, and the loss of ecological services in the SED scenario was higher than that in the BAU scenario, as there was a high speed of urbanization and rapid urban sprawl in the future. We believe that socioeconomic changes are the dominant driver of changes in ESVs [19]. In the CPP scenario, the ESVs for LUCC decreased the least, which showed low land use development and high ecological protection in the LH quadrant. Measures should be taken to protect farmland. On the one hand, the transformation of land use change should ensure a higher-level, higher-quality food supply, and a more efficient and more sustainable food security system at the regional level. On the other hand, land use should prevent the conversion of farmland into nongrain crop production areas and stabilize food production. The EPP scenario was the same as the CPP scenario, which is shown in the LH quadrant. Ecological environment protection should be considered in future land use change initiatives. The greatest conversion would be to forestland, grassland, and water area, increasing the ESVs in the multiscenario simulations.

The spatial spillover effects of land use changes on ESVs are significant and differ among different regions with different economic levels. There have been frequent changes in land use in the highly urbanized and economically important regions; for example, farmland, forestland, and water area have been converted to built-up land, promoting urbanization, and leading to more significant spatial spillover effects. The spatial distribution of land use led to the spatial heterogeneity in the loss of ESV. Furthermore, the ecological effects of land use change were also ignored by land use management and land use planning in China and caused by ecological environment deterioration [58] and the loss of ESVs. Our objective was to determine how to allocate and develop land use without compromising ecological sustainability so that the synergistic effect of urban high-quality development and ecological civilization construction (shown in red in Figure 9) will be generated.

In response to population growth, intensification of human activities, government policies, and initiatives that aim to improve the economy of the country without considering the environmental consequences [59], we recommend that the spatial effects of ESV changes be considered in response to the impact of LUCC on ESVs in the GBA. With economic and social development stage transformation, the land use transition in this area should be carefully controlled and, facing the degradation of ecosystem services, we should provide reasonable guidance for land use development.

4.1.2. Management Implications

Future spatiotemporal changes in LUCC and ESV were predicted through multiscenario simulations, and the results showed potential effects on long timescales in the GBA. This study will help to understand critical tradeoffs between ecosystem services caused by LUCC and provide valuable information for decision-making and policy development.

The loss of ESV has significant spatial heterogeneity, caused by the spatial difference in land use. To break the “Matthew effect” of ecosystem services, the low value of ecosystem services should therefore be prioritized for protection. We suggest that the characteristics and synergistic effect of different driving factors should be considered to optimize the ecosystem and control ecological risks, especially for low ESVs. The impact of complex human activities on ecosystem services made us adopt differentiated and diversified regulating strategies. The government should comprehensively consider the ESV changes in land use and emphasize the spatially explicit extent of the ecological environment effect in the GBA.

In addition, land use patterns are a critical way to protect ecosystem services, which should coordinate with local geographic conditions and socioeconomic development in multiscenario simulations. With the principle of the symbiotic development of land use development and ecological protection, the GBA needs to improve the spatial agglomeration of urban constructed land, balance the occupation and supplementary use of cultivated land, and maintain the high-value ecosystem service of ecological spaces. These are in response to a strategy for ecological civilization construction to improve the harmony between nature and human activities, which led to the publication of the Development Plan for the Guangdong-Hong Kong-Macao Greater Bay Area in 2018. Recently, the local government attempted to build an ecological compensation mechanism, following the principle that developers are responsible for protecting the environment and users must compensate for damages they cause. The ecological compensation mechanism should be further improved through ESVs and ecological protection costs and will be an important route for exploring ecological protection work.

4.1.3. Limitations and Future Works

Several limitations of the study should be considered. The first limitation of this study is associated with the ESV valuation method. In the ecosystem service valuation process, value coefficients are assigned to each land use type; in this study, the value coefficients used were those of Xie et al. [12] and the adopted local modified ecosystem services coefficients. However, with the social and economic development, there should also be differences in the understanding of ecological service value. Furthermore, with the rapid urbanization process, the demand for land resources should also affect ecological service functions. Therefore, the revising factors of social economic and resource scarcity could be considered with the value coefficient of ecological services. The second limitation is that the selection of driving factors could affect LUCC. Logically, these representative factors should have an important influence on simulation accuracy and directly affect the probability of occurrence of land use types. Several factors were considered in the study, such as geographical factors, environmental factors, socioeconomic factors, and POI. However, there was a lack of government-related decision-making factors, which may have led to uncertainties in land use simulation. Consequently, the introduction of the decision-making factors of land use-related stakeholders will enhance the scientific and comprehensive spatial simulation through the bottom-up analysis of natural carrying capacity and the top-down analysis of territorial spatial structure order. Urban expansion can be characterized by low speed and compact and sustainable growth under decision-making intervention. In addition, the accuracy and confidence level of the analyzed results would be more significant in future simulation research.

In future studies, we will test the applicability of higher-resolution simulations and focus on comparative experiments between the simulation results among more comprehensive restricted factors. Driving factors have relatively different impacts on land use, and the weights of different types of driving factors must be assigned. Future work may include adding the driving factor weights to the FLUS mode. The FLUS model transition rules should also be considered, as those rules may change over a long period. We hope to address this challenge in future works.

4.2. Conclusions

The spatiotemporal evolution characteristics of the impacts of LUCC on ecosystem services in the GBA by 2030 were analyzed in the multiscenario simulations, and the following conclusions were drawn.

In 2005–2020, forestland and farmland were the predominant land use/cover types, and their proportions continued to decline. Forestland was mainly distributed in the western, northwestern, and eastern regions, and farmland was found in the central, southwestern, and eastern regions. The water area declined continuously in the northern region of the GBA. In contrast, the proportion of built-up land increased continuously with high urbanization and economic levels as a consequence of sacrificing farmland, forestland, and water area, and it was concentrated in the central region. The proportion of grassland and wetland tended to fluctuate, initially decreasing, then increasing, and finally decreasing.

During 2005–2020, ESVs continuously declined and decreased significantly by 14.29 billion yuan. The ESV per unit area showed a significant decline, which was largely associated with a decrease in farmland and an increase in built-up land. Forestland contributed the most to the total ESV (>75%), and those areas experienced great land use transitions and human activity. The high-value ecosystem services were concentrated around the peripheral area of the northern regions in the GBA.

Compared with 2020, the total ESV in the SED, BAU, and CPP scenarios in 2030 decreased, while it increased in the EPP scenario. The transformation of land use change was the main reason for the decrease in ESV. For example, most land use types will be converted, changing to or from farmland, built-up land, forestland, and water area.

Data Availability

The data used to support the findings of this study can be obtained from the corresponding author upon reasonable request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by the Social Science Federation of Fujian Province, under Grant FJ2019C034, and Huaqiao University, under Grant 16SKGC-QG11.

Supplementary Materials

Table S1. The ROC of logistic regression. Table S2. The change of land use and ecosystem service value in BUA scenario in 2030. Table S3. The change of land use and ecosystem service value in SED scenario in 2030. Table S4. The change of land use and ecosystem service value in CPP scenario in 2030. Table S5. The change of land use and ecosystem service value in EPP scenario in 2030. (Supplementary Materials)