Assessing spatiotemporal variations and predicting changes in ecosystem service values in the Guangdong–Hong Kong–Macao Greater Bay Area

ABSTRACT Rapid economic development and interference by human activities in rapid urbanization regions have caused great land use/land cover change (LUCC), which significantly affects ecosystem functions and services. It is crucial to assess the spatiotemporal evolution of ecosystem service value (ESV) in such regions, especially for the Guangdong–Hong Kong–Macao Greater Bay Area (GBA) in China. In this study, we investigated and predicted the effect of LUCC on the ESV in the GBA from 1990 to 2030 using the latest annual 30 m LUCC database of China, the future land use simulation (FLUS) model, and ecosystem service evaluation approaches. The study period was divided into the historical period (1990–2015) and the forecast period (2015–2030). The results showed that forest and cropland were the dominant land-use types (>77% of the GBA), and the expansion of built-up land (3822.4 km2) was the clearest process during 1990–2015. The reduction of cropland and forest contributed the most to the decrease in the total ESV. Moreover, the results confirmed that the FLUS model is effective at predicting future LUCC in the GBA. The ESV was predicted to decrease to 4962.23 × 100 million yuan in the 2030s under the current development mode if regional forest and waterbody reductions are not constrained. This study provides a reference for promoting the rational use of land resources and ecological construction in the GBA and can help to promote ecological planning and environmental protection.


Introduction
Land use/land cover change (LUCC) is the direct response of the regional natural environment to human economic activities (Lambin and Geist 2008;Pielke 2005). Serving as a bridge between human activities and the Earth's surface, LUCC plays an crucial role in global change (Ojima, Galvin, and Turner 1994;Smith et al. 2016;Wilbanks and Kates 1999). Human activities can directly or indirectly cause regional environmental change, which produce different land-use patterns (Veldkamp and Lambin 2001;Wu 2013). LUCC is also an indicator of human activities that modify the Earth's landscape and societal development (Liu et al. 2010). With the increase in global population, the demand for global land resources is increasing, and the health and productivity of the land are deteriorating (Cowie et al. 2018;Ramankutty et al. 2006;Valipour, Bateni, and Jun 2021). Moreover, land-use change, as the clearest manifestation of global environmental change, can indirectly affect surface material cycles and ecological processes such as climate change, biodiversity, biogeochemical cycles, and sustainable use of resources (Imhoff et al. 2004;Liang, Hashimoto, and Liu 2021;Ojima, Galvin, and Turner 1994;Song and Deng 2017;Valipour, Bateni, and Jun 2021). The ability of an ecosystem to provide services for human well-being is directly linked to the state of the ecosystem (its structure and processes) (Costanza 2000). The increased pressure on the ecosystem or change in land-use types can influence the ecosystem service (ES) supply or trade-offs between different services (Bateman et al. 2013;Ma et al. 2020).
Since Costanza et al. (1997) introduced the concept of ecosystem service value (ESV) in 1997, many studies have contributed to related theories, methods, and systems (Braat and de Groot 2012;Daily et al. 2000;Shoyama et al. 2017). Because the specific ecosystem situation varies across regions, direct use may lead to overestimation or underestimation of the results (Bryant et al. 2018;Harrison et al. 2018). Therefore, many scholars have tried to improve ES classification, service value, and service value evaluation methods (Costanza et al. 2014;Luck et al. 2009;Martínez-Harms and Balvanera 2012;Roche and Campagne 2019). Urban expansion is a land-use change process that transforms non-urban land into urban land, which directly drives LUCC and changes the ecosystem (Lawler et al. 2014;Li et al. 2016). Many scholars have discussed how urban expansion affects ESV at different spatial scales Kertész, Nagy, and Balázs 2019;Shiferaw et al. 2019). However, existing research on the relationship between urban expansion and ESV has mainly focused on the correlation between the two, and it is difficult to determine the detailed quantity and spatial distribution of the ESV caused by urban expansion. In addition, ignoring the impact of urban agglomeration LUCC on ESV under multiple scenarios will limit the actual reference value of the results in government land-use planning. Moreover, the temporal scale and spatial details of most of these studies are insufficient to illustrate internal differences.
The Guangdong-Hong Kong-Macao Greater Bay Area (GBA) has been at the forefront of China's reform and opening up . The economy has developed rapidly in recent years. With the release of "the outline development plan for the GBA" in February 2019, which aims to develop the GBA into "a role of high-quality model development," the GBA ushered in new development opportunities . The rapid urban expansion in the GBA has accelerated the transformation of rural landscapes into urban landscapes, thereby resulting in severe damage to the landscape pattern and urban ecology, which will inevitably cause significant changes in ESV (Peng et al. 2015).
The limitations of ESV-related studies, especially in spatiotemporal analysis and predictions, are the primary gaps in previous studies. To fill these gaps, the future land use simulation (FLUS) model and ES evaluation approaches are introduced in this study to simulate LUCC and the corresponding ESV in the GBA. Specifically, this study aimed to examine (1) the spatiotemporal LUCC dynamics in the GBA; (2) the resulting changes in ESV; and (3) the prediction of the future LUCC and ESV changes in 2030. This study can provide a reference for promoting the rational use of land resources and ecological construction in the GBA and help to promote ecological civilization planning and environmental protection in the GBA.

Study area
The GBA is a typical example of rapid urbanization in China. It is located in China's southern coastal area (21°32ʹ-24°26ʹN, 111°20ʹ-115°24ʹE) and encompasses two special administrative regions (Hong Kong and Macao) and nine cities in Guangdong Province (Guangzhou, Shenzhen, Foshan, Jiangmen, Zhongshan, Huizhou, Dongguan, Zhaoqing, and Zhuhai) (Figure 1). With a population of approximately 70 million people and a total area of approximately 56,000 km 2 , the GBA is one of China's most economically developed and open areas (Fig. S1). It has a diverse range of high-tech industries, manufacturing plants, overseas enterprises, financial companies, and educational resources, accounting for 12.6% of China's gross domestic product (GDP) and 0.6% of the country's population in 2018 (China National Bureau of Statistics 2019). The GBA is one of the four major bay areas in the world (Zhou et al. 2018). The Chinese government places a high priority on GBA construction and development (The Central Government of China 2019). Many policies and plans to simulate regional development have been announced (Chen et al. 2020c;Zhang et al. 2020).

Methodology and Data Sources
Detailed methodology and data sources are shown in Figure 2. Annual LUCC data sources in the GBA were collected, as described in Section 3.1. The proposed calculation approach of ESV, including different land-  use types and the corresponding modified equivalent factor, is explained in Section 3.2. The FLUS model used to simulate LUCC through 2030 is described in Section 3.3. The predicted LUCC with satisfactory accuracy was used for corresponding future ESV analysis. We will describe these contents in detail in the following sections.

LUCC data sources and change analysis
The long-term series LUCC data were obtained from the latest annual 30 m LUCC database of China, which was developed by Tsinghua University, China . The database spans from 1980 to 2015 with 30 m spatial resolution and an annual time step. The data integrated images from the Moderate Resolution Imaging Spectroradiometer, Advanced Very High-Resolution Radiometer, and Landsat data using the Breaks for Additive Seasonal and Trend algorithm ). More details can be found in previous studies (Tu et al. 2021;Xu et al. 2020). The dataset performed well in annual classification and changedetection accuracy (overall 75.61%), with annual average accuracy of 72.10%, 78.93%, and 91.89% for cropland, forest, and built-up land, respectively. It provides basic data for research on urban planning and ecological assessment (Tu et al. 2021;Zhao et al. 2020).
The land-use transition matrix was used to quantify the land-use change (Ferreira Filho and Horridge 2014). The land-use transition matrixes were defined by comparing a successive set of images in different periods. In this study, we used ArcGIS 10.5 software to calculate transition matrixes of different periods, including 1990-2000, 2000-2010, and 2010-2015. As shown in Table 1, the transition matrix displays the types of land use at periods P 1 and P 2 . The symbol S ij represents the area of land transferred from type i to type j. The diagonal elements (i.e. S ii ) indicate the area of land retained in period P 2 . Elements excluded from the diagonal are land-use transitions from type i to type j. The area of land in different types (S i+ and S +j ) is the sum of S ij over all of j and i, respectively. A further computation based on the matrix table was created to calculate the gains (S +i -S jj ) and losses (S i+ -S ii ). The gain of one land-use type is equivalent to the transition from other land-use types between the study periods, whereas the loss represents the transition to other land-use types between periods.

ESV analysis
According to the United Nations, ecological services are grouped into four groups, namely provisioning services, regulating services, supporting services, and cultural services (Costanza et al. 1997). Based on the evaluation model (Costanza et al. 1997) and the real situation in China, Xie et al. (2003) proposed the equivalent factor method to evaluate China's ESV. Subsequently, Xie et al. (2015) further revised the method and established an equivalent table of value coefficients (VCs) for different ESVs in China. In the equivalent table, the equivalent coefficient value of cropland is 1, which is the basic reference for regional correction and other land-use types (Xie et al. 2015). To ensure that the coefficients were appropriate to calculate the ESV for different regions, we revised the table for the GBA via correction coefficients based on the main grain yield in the GBA, as follows (Eq. 1): where β represents the correction coefficient at time t and y and Y are the grain yields of the GBA and China at time t, respectively. In 2015, y and Y were 5.25 t/ha and 5.48 t/ha, respectively. Thus, β was calculated as 0.96.
Correspondingly, the equivalent values per unit area for different ecosystems in the GBA were calculated as follows (Eq. 2): where Y and A g are the total yield and area of the main grains in Guangdong Province in 2015, respectively, and P is the average price of the main grains in 2015 (1939.7 yuan/t), which was acquired from the Grain Net of South China (http://www.grainmarket.com.cn). The modified equivalent factor method based on unit area value was then adopted to calculate the ESV. The ESV per unit area of different ecosystems of the GBA in 2015 is listed in Table S1.
Finally, the ESV was calculated based on Eq.
(3), as follows: where S k is the area of k land use (ha) and VC k is the ESV per unit area of k land use (yuan/ha). We also calculated the coefficient of sensitivity (CS) based on the ESV and VC (Bian and Lu 2013;Das and Das 2019;Kreuter et al. 2001). The results showed that the CSs of ESVs were all less than 1 (CS<1), which indicated that ESVs were inelastic in response to the VC and the modified ESV coefficient was applicable to the study area.

FLUS model
The FLUS model is an improved future land-use change model that can explicitly simulate the longterm spatial trajectories of many LUCCs while accounting for human and natural environmental impacts, as well as solving complex land-use demand prediction and land-use allocation issues . It consists of a top-down system dynamics model and a down-top cellular automata (CA) model (Liang et al. 2018a). The artificial neural network (ANN) model is used to train non-linear relationships between historical land-use types and complex driving factors so that the probability of distribution can be calculated.
Land use data from 2005, as well as data on relevant physical, social, and economic dimensions, was used to run the FLUS model (Table 2). To represent the changes in distribution of different land-use types better, many driving factors were considered based on existing studies Zhe et al. 2020), including the digital elevation model, aspect, slope, river system, all level roads, and other data (e.g. protection zone) ( Table 2). These data were used to train the ANN model for estimating the probability of occurrence in urban and non-urban areas. Ten neurons in the input layer and twelve neurons in the hidden layers made up the backpropagation ANN model used in this study. Onefifth of the pixels were chosen at random across the GBA as the training datasets. To normalize the probability values to the range of [0, 1], we used the sigmoid function as the activation function for output layers (Liang et al. 2018b). During the training process, the learning rate and terminal conditions of the ANN model are self-adaptive. In the simulation module, considering the filtering results and processing efficiency, the 3 × 3 Moore neighborhood was used to run the FLUS model based on repeated testing (Liang et al. 2018b;Wang et al. 2021). For the first iterations, the initial inertia coefficients were set to 1 and would self-adapt during the CA iteration.

Markov chain model
The Markov chain model can predict dynamic variations with high prediction precision and accuracy. It has been widely used in the prediction of land-use changes (Arsanjani et al. 2013). The land-use demand in 2030 was calculated using the Markov chain model in this study. We collected the land-use data in 2010 and 2015 and used the Markov chain transfer matrix to evaluate the reciprocal transform relationship between different land-use types. Future land-use types in 2030 were projected using a land conversion probability matrix with a fiveyear step, which was calculated using Eq. (4). The results are presented in Section 4.3.2.
where n is the number of land-use types and P ij is the probability that land-use type i is converted to landuse type j (0≤ P ij ≤1).

Model validation
To validate the accuracy of the FLUS model, we compared the simulated LUCC with the remotely sensed LUCC in 2015 using four indicators, namely producer's accuracy (PA), user's accuracy (UA), Kappa, and the figure of merit (FoM) (Pontius et al. 2008). PA and UA are commonly used indicators for accuracy validation, and the calculation equations can be found in the study by Shao and Wu (2008). Kappa was calculated using Eq. (5), as follows: where P 0 denotes the simulation's correct proportion, P c denotes the model's correct proportion in the random case, and P p denotes the proportion of the correct simulation with ideal classification. FoM (Eq. (6)) is the ratio of the intersection of the observed and predicted change to the union of the observed and predicted change (Murayama 2012;Perica and Foufoula-Georgiou 1996). The range of FoM is 0% to 100%, thereby indicating that there is no overlap between observed and predicted change and a perfectly accurate prediction (Wang and Li 2011).
where A represents the error area due to observed change predicted as persistence, B represents the correct area due to observed change predicted as change, C indicates the error area due to observed change predicted as the wrong gaining category, and D indicates the error area due to observed persistence predicted as change (Hu, Li, and Lu 2018;Perica and Foufoula-Georgiou 1996).

Changes in LUCC
Based on the results of land-use changes in the GBA during 1990-2015 ( Figure 3 and Table 3), forest and cropland are the land-use types with the largest areas.
In 2015, they accounted for over 77% of the total area of the GBA (54% and 23%, respectively), whereas the proportions of wetland and unutilized land were approximately 1%.  (Table S2), cropland (2147.30 km 2 ) and forest (1015.18 km 2 ) were the main two contributors to the increased area of built-up land. Because the total amount of grassland is small, the overall change was not clear.

Spatiotemporal changes in ESV in the GBA
The ESV in the GBA in 1990GBA in , 2000GBA in , 2010GBA in , and 2015.80 billion yuan, respectively. In general, the ESV increased first and then decreased from 1990 to 2015, decreasing by 16.92 billion yuan from 1990 to 2015. From 1990 to 2000, the ESV first increased by 4.20 billion and then decreased after 2000 (Table 4).
Because of the continuous decline in forest and cropland areas, their ESVs have continued to decrease. From 1990 to 2015, the ESV of forest and cropland decreased by 13.23 and 2.15 billion yuan, respectively. The variation in grassland showed a slow descending trend with a decrease of 0.54 billion yuan during 1990-2015. The changes in wetland ESV followed a similar trend to that of grassland, decreasing from 3.33 billion yuan in 1990 to 1.54 billion yuan in 2015. The ESV of unutilized land was stable owing to the small land-use change. Waterbody was the only ecosystem for which the ESV increased, but the increase was far smaller than the decrease of other ecosystems. The total ESV decreased by 3%.
The ESV of different ecosystem functions in the GBA is shown in Table 5. In 2015, the main contributors to each ecosystem function were hydrological regulation (32%), climate regulation (23%), soil conservation (9%), and biodiversity conservation (9%). The sum of hydrological regulation and climate regulation accounted for over 50% of the contribution, thereby indicating that regulating services are the  greatest. The sum of regulating services and supporting services was much greater than that of provisioning services and cultural services. Significant decreases in food production, maintain nutrient cycle, and gas regulation were detected owing to substantially reduced areas of wetland and cropland. Only water supply increased with a rate of change of 24% from 1990 (5.02 billion yuan) to 2015 (6.22 billion yuan), which might have been attributed to the increase in waterbody and cropland areas.

Results validation
As indicated in Table 7, the verification results revealed that Kappa and FoM reached 0.79 and 0.07, respectively, which illustrated that the simulation  We also conducted a comparison of actual LUCC with FLUS-predicted LUCC in 2015 (Figure 4). The actual and simulated areas of cropland, forest, grassland, waterbody, built-up land, and unutilized land in 2015 were 12 508. 43 (simulated: 11 933.67), 29 785.22 (29 500.27), 1258.43 (1149.36), 3864.53 (3796.29), 139.73 (136.84), 7566.06 (8611.64), and 14.08 (13.78) km 2 , respectively (Table 2 and S3). The spatial distribution of different land use also suggested that the performance of the FLUS model was satisfactory and could be further used for subsequent analysis.

LUCC prediction
The driving factors, including land use, socioeconomic data, terrain, and others, were used for the implementation of the FLUS model. The simulated results mainly reflected the future LUCC via general development mode in the GBA. The predicted areas of each land type in the GBA and subordinate cities were analyzed using 2015 as the basis. Under this specific scenario, built-up land was still predicted to be the fastest-growing land-use type  (9282.22 km 2 ) in 2030. However, the varying decreases in cropland (12 492.38 km 2 ), forest (28 671.07 km 2 ), grassland (1157.82 km 2 ), waterbody (3384.76 km 2 ), wetland (134.50 km 2 ), and unutilized land (13.73 km 2 ) were predicted (Table S4).
Regarding the spatial distribution, Jiangmen occupies the largest area of cropland (2499 km 2 ). The areas with the largest area of forest, grassland, waterbody, and built-up land in 2030 were Zhaoqing (10,925 km 2 ), Jiangmen (283 km 2 ), Foshan (814 km 2 ), and Guangzhou (1634 km 2 ), respectively ( Figure 5). The predicted LUCC of the GBA in 2030 suggests that its spatial distribution pattern will be similar to that of 2015.

ESV prediction
Based on the FLUS-based LUCC and calculation framework of ESV (Sections 3.2 and 4.3.2), we predicted the ESV in 2030. The simulation results indicated that the ESV would decrease to 4962.23 × 100 million yuan in 2030 (Fig. S2). The greatest ESVs of different ESs were food production (93.28 × 100 million yuan), raw materials (121.35 × 100 million yuan), and water supply (51.44 × 100 million yuan) (Fig. S2). Compared with the ESV in 2015, the total ESV in the GBA in 2030 decreased by 285.78 × 100 million yuan owing to the reduction in ecological lands (forest, grassland, wetland, and waterbody). Among all ES types, hydrological regulation, climate regulation, and gas regulation were the three with the greatest predicted decrease in ESV in the future. Among the individual cities in 2030, Zhaoqing was predicted to contribute the most to the total ESV in the GBA in 2030 (Table 8). The corresponding contribution in descending order was Zhaoqing (1623.22 × 100 million yuan), Jiangmen (822.80 × 100 million yuan), Foshan (334.05 × 100 million yuan), Shenzhen (127.86 × 100 million yuan), Dongguan (141.56 × 100 million yuan), Guangzhou (518.02 × 100 million yuan), Hong Kong (100.10 × 100 million yuan), Zhuhai (104.02 × 100 million yuan), and Macao (1.55 × 100 million yuan). Compared with historical ESV in 1990-2015, clear changes in ESV were found in the regions with rapid urban expansion, such as Guangzhou, Shenzhen, and Zhuhai. The decreased ESV was also relevant to the reduction in regional forest and waterbody.

Eco-environmental effects of LUCC
LUCC is the most direct interaction between anthropogenic activities and the natural environment Song and Deng 2017) and a main cause of environmental changes. Over the past decades, rapid LUCC changes have occurred in the GBA (Fang et al. 2020;Gong et al. 2020a;Liu et al. 2020). Population, economic activity changes, and transportation development have been recognized as the most significant drivers of LUCC and the key link of the coupled human-environment system Zhang et al. 2020), especially after China's reform and opening up in the 1970s (Govada and Rodgers 2019) (Jiang, Ye, and Ma 2014). Regional industrial structure, consumption patterns, and dietary structure can also affect LUCC (Meng et al. 2021).
In this study, we investigated and predicted the effect of LUCC on the ESV of the GBA from 1990 to 2030 using the latest annual 30 m LUCC database of China, the FLUS model, and ES evaluation approaches. Analysis of the spatiotemporal evolution process and characteristics of LUCC and ESV of the GBA is an important basis for optimizing ecosystem management and improving ecosystem quality. Correspondingly, the estimation of regional ESV can help to understand the synergistic effects of ESs with ecological coupling models. The related contents can provide insights into the relationships and mechanisms among different ESs and provide references for regional ES management and socioeconomic development. From this study, we can conclude that the ESV in the GBA is mainly composed of the ESs of forest, wetland, and waterbody. However, owing to the needs of urban development, the expansion of built-up land has crowded out other land-use types, thereby resulting in inevitable declines in other natural and semi-natural ecosystem areas. Therefore, the total ESV in the GBA decreased by 201.20 × 100 million yuan from 1990 to 2015. Based on our simulated results, the ESV will decrease to 4962.23 × 100 million yuan in the 2030s.
Urbanization is one of the most important drivers of the variations in ecosystem types and ESVs in the GBA. From 1986 to 2017, the total GDP increased from 284,245 million yuan to 10,032,690 million yuan, and the number of permanent residents increased from 2395.21 × 10 4 to 6797.7 × 10 4 , an increase of 34 times and two-fold, respectively. Meanwhile, the expansion of cities and construction land continue to occupy cropland, forests, and wetlands, and human activities also affect the pattern of ecosystems and their potential functions. Guangdong had the fourth largest reclamation area. Jiang et al. (2021) found that over the past 30 years, the area of the GBA affected by reclamation activities was close to 600 km 2 . The types of reclamation and utilization have changed from cropland and aquaculture pond land to urban expansion. Excessive reclamation has seriously affected the ecological functions of coastal wetland purification and pollution, beach protection, and bank protection. The saltwater marshes and mangroves along the coast of the GBA have been destroyed, which in turn has affected the water conservation function. This is consistent with our finding that the amount of water supply (hydrological regulation) in coastal cities such as Shenzhen, Dongguan, and Zhongshan has decreased. The decrease in supporting services (soil conservation and biodiversity conservation) may be attributed to the expansion of urban land and ineffective implementation of forest protection and restoration policies.

Policy implications and suggestions
Over the past few decades, China's urban agglomeration areas have developed rapidly. Adhering to the concept of development as the main line, many phenomena such as land urbanization and population urbanization are common, which result in the loss of the overall ecological space and the deterioration of the ecological environment of the urban area, thereby making it difficult to form a sustainable urbanization development process (Fang and Yu 2017). Therefore, the latest measures of land and space planning and new urbanization planning have raised economic and social development and ecological environmental protection to the highest level (Chen, Liu, and Lu 2016).
Our results show that significant reduction and degradation of ESV have occurred during the rapid urbanization process in the GBA, which have caused significant damage to the ecosystem and serious environmental problems. For example, as a result of human disturbance, such as coastal zone development and wetland reclamation, natural wetlands have been sharply reduced and their functions have declined. In accordance with our results, from 1980 to 2015, the coastal reclamation area reached 467.33 km 2 and the natural coastline declined by 289.62 km 2 (18.71%) (Zhao et al. 2019). Natural mangroves in the GBA face the problems of shrinking areas, alien species invasion, and biodiversity loss . Moreover, landscape fragmentation has been greatly increased in the GBA and the area of the traditional reservoir pit in the GBA has sharply decreased by 47.33% since 1980 (Gong et al. 2020b). In 2018, the population density of Shenzhen was 6522 people/km 2 , which was the highest among all cities in China. Dongguan, Foshan, and Guangzhou are also among the top 10 cities with the highest population density in China (China National Bureau of Statistics 2019). With the economic development and urban expansion of the GBA, the population will continue to grow. It is expected that the population of the GBA will reach 150 million in 2022 and the conflict between increasing population and limited land will increase. Thus, governments need to establish more scientific and rational urban planning policies and regulations with a strong emphasis on protecting the natural ecosystems of the GBA.

Limitations and perspectives
Implementing the ES compensation system has emphasized ESV-related research Shiferaw et al. 2019;Zhang, Yushanjiang, and Jing 2019). Although this study tried to reveal the regional past-present-future spatiotemporal evolution patterns of ESV and provide necessary references for effective ecological planning and sustainable development, there are still some unavoidable uncertainties. As ecological services and ecological functions cannot be completely matched, there are insurmountable obstacles in the accurate calculation of ESV (Costanza et al. 2014;Shoyama et al. 2017). In our future research, various ES functions need to be refined. Moreover, a more accurate accounting system (classification and design) will help to accurately explain the variations in the ESV in the GBA. In this study, the FLUS model was used to predict the ecological landscapes and the spatial evolution of ESV. Considering the inherent principle of the FLUS model, we resampled the LUCC data. Conventional studies have only relied on transition probability to predict future evolution trends Song and Deng 2017). In comparison, this study has made progress in spatial simulation of ESV, and the corresponding simulation results are satisfactory. Quantifying the regional ESV is the first step; a deep understanding of ESV is the key. Theoretically, our understanding of ecosystem value is constantly changing. Many ESVs exist, but they have not been included in the existing ESV accounting system because a universal measurement method has yet to be established. As a result, we should conduct related explorations, further expand the scope of accounting, and explore the value measurement methods of the contribution of ecosystems to urban development. In terms of specific operational applications, we intend to develop a comprehensive platform that includes both accurate extraction of LUCC and comprehensive measurement of ESV and deploy the platform in a cloud computing environment, such as Google Earth Engine. It is expected that the developed platform will link contributions from government departments and scientific simulations via remotely sensed data.
ESV-related studies also need to consider the synergy and restrictions between economic development and the ecological environment, as well as the correlation between economic gains and ecological losses (Ma et al. 2020). Therefore, we need to analyze the coupling trend between economic development and environmental change under different development paths using dynamic change models and to establish a database of ecological service trade-off thresholds. Although not implemented in this study, we intend to develop a spatial allocation model by coupling the multiple LUCC dynamic simulation model to predict the spatial distribution of land use under more scenarios and depict different socioeconomic developments (population, GDP, and technological innovation) and natural climate changes (temperature and precipitation). Our study provides support to improve policymaking and understanding of the laws of LUCC.

Conclusion
This study analyzed the spatiotemporal LUCC dynamics and the resulting changes in ESV in the GBA from 1990 to 2015. Based on the results and other auxiliary data, the FLUS model was further used to simulate future ESV in the 2030s in the GBA.
Over the past decades, forest and cropland were the dominant land-use types, covering >77% of the GBA. From 1990 to 2015, the expansion of built-up land (3822.4 km 2 ) was the clearest process. The reduction in cropland (2606.74 km 2 ) and forest (983.93 km 2 ) determined the decrease in the total ESV in the GBA during 1990-2015 to a great extent. The sum of hydrological regulation and climate regulation accounted for over 50% of the contribution, thereby indicating that the specific service type of regulating services is the greatest contributor. In terms of historical ESV changes in different cities in the GBA, Zhaoqing and Macao had the maximum and minimum ESV, respectively. Moreover, the results suggest that the FLUS model is an effective method to predict FLUS results satisfactorily, with a Kappa coefficient and FoM of 0.79 and 0.07, respectively. Based on the FLUS-based LUCC, the simulation results indicate that the ESV will decrease to almost 4962.23 × 100 million yuan in the 2030s. The decrease in ESV is also relevant to the reduction of regional forest and waterbody.