Impacts of Historical Land Use Changes on Ecosystem Services in Guangdong Province, China

: Assessing land use change and its impacts on ecosystem services is of great signiﬁcance for optimizing land use management and enhancing ecosystem sustainability. This study explores land use changes and their impacts on ﬁve typical ecosystem services, namely grain production (GP), water yield (WY), soil conservation (SC), habitat quality (HQ), and carbon sequestration (CS), during 1990–2020 using the Integrated Valuation of Ecosystem Services and Trade-offs (InVEST) model in Guangdong province, which has experienced substantial land use change. During the study period, cultivated land, forest land, grassland, water areas, built-up land, and unused land correspond-ingly had changed by − 10.7%, − 1.9%, − 5.1%, 13.7%, 97.9%, and − 38.8%. For ecosystem services, the GP, SC, and HQ averagely decreased by − 8.66% ( − 12.3 t · km − 2 ), − 0.02% ( − 2 t · km − 2 ), and − 2.74% ( − 0.02), respectively, while WY and CS increased by 3.10% (22 mm) and 20.70% (515 t · km − 2 ), respectively. Land use changes that had the greatest average negative impacts on GP, WY, SC, HQ, and CS were cultivated land to built-up land ( − 150.9 t · km − 2 ), unused land to water areas ( − 1072 mm), grassland to unused land ( − 10,166 t · km − 2 ), forest land to built-up land ( − 0.65), and forest land to water areas ( − 2974 t · km − 2 ) respectively, and that had the greatest average positive impacts were grassland to cultivated land (78.8 t · km − 2 ), water areas to built-up land (943 mm), unused land to forest land (3552 t · km − 2 ), built-up land to forest land (0.40), and water areas to forest land (3338 t · km − 2 ), respectively. The results indicated that land use and its changes had a signiﬁcant impact on ecosystem services.


Introduction
As the benefits that humans derive from the structures, functions, and processes of ecosystems, ecosystem services play a critical role in achieving sustainable development goals [1][2][3].The terrestrial ecosystem is the closest ecosystem to human activities, providing ecosystem services such as fresh water, soil, and biological resources.It is also the spatial basis and material guarantee for sustainable development of human society [4].With rapid urbanization, the terrestrial ecosystem is increasingly modified, which has both negative and positive effects on the ability of the terrestrial ecosystem to provide services [5].The specific investigation of such effects will have a positive impact on the rational formulation of regional development plans, as well as the achievement of sustainable social development goals.
Land use changes can directly reflect the interaction between human activities and the terrestrial ecosystem in the process of socioeconomic development, especially in the context of the trade-off between rapid urban expansion and sustainable development in recent years [6].On the one hand, human destruction and interference with the natural environment have become more and more obvious [7].For example, human beings have expanded many cities to increase living space.The means include occupying the natural space and destroying the continuity of the natural environment, which makes them fragmented [8].At the same time, the transitional demand for resources in the ecosystem exceeds its ability to recover (such as excessive cutting of trees and excessive fishing), leading to irreversible damage to some ecosystems [9,10].However, on the other hand, appropriate policies and actions promote the restoration and protection of the natural environment [11].For example, China's Three Northern Protected Forests Project has provided a significant contribution to improving the environment in northern China [12].The United Nations Sustainable Development Agenda has also become an important guide for countries to balance environmental protection with development [13].Due to specific environments, species, resource types, and social values, land use types play different roles in ecosystems, with indispensable meanings in exploring regional ecosystem services [14][15][16].
Research on the relationship between land use and ecosystem services has been a hot topic in related fields.Many scholars have done relevant research.For example, Xie et al., used the improved unit-value-based method to evaluate the ecosystem service values of arid inland watershed, with the conclusion that grassland has the greatest ecosystem service value among all land use/cover types in arid inland watershed [17].Muhammad et al. assessed the dynamics of land use/cover change and associated ecosystem service values of coastal Bangladesh during 1999-2019 by analyzing historical Landsat land use/cover images and economic valuation techniques [18].Morshed et al. built the model of future ecosystem service value with land use/cover dynamics by using machine learning-based artificial neural network model for Jashore city, Bangladesh [19].Most of this research focused on the calculation and change analyses of ecosystem services under land use change background, while fewer studies looked into the ecosystem services induced by different land use conversions.Various complex mathematical and spatial models are used to quantify the impacts of land use on ecosystem services.Research interests include analyzing the history, linking the present, and predicting the future of ecosystem services.For example, the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) model spatializes the quantitative valuation of ecosystem services by simulating changes in the quantity and value of ecosystem goods under land cover scenarios [20].Cellular automata (CA) can simulate spatial patterns of land use changes, from the microscopic level to the macroscopic, by combining with GIS to establish interlinkages between land use categories and drivers [21].Multi-criteria decision analysis addresses the impacts of potential land use conversions on ecosystem services [22].The Land Change Modeler (LCM), as a mean of predicting likely future ecosystem conditions, enables good prediction of future land use uncertainty [23].The GLOBIO model is widely used to assess anthropogenic changes in biodiversity [24].Among the models, the InVEST model is widely applied all over the world to assess the impacts of land use change on ecosystem services [25], which is appropriate for reliably assessing multiple ecosystem services.
Short-term land use changes are highly susceptible to policy influence, and the benefits of land use changes can greatly influence policy makers' decision-making [26].Therefore, exploring the impacts of land use changes on ecosystem services can provide a clear reference for policy makers, which promotes the achievement of sustainable development goals.In view of this, this study applied the InVEST model to assess how land use changes impacted on ecosystem services from 1990 to 2020 in Guadong province.Considering the significance of ecosystem services and data availability, this study selected grain production (GP), water yield (WY), soil conservation (SC), habitat quality (HQ), and carbon sequestration (CS) services as typical representatives of ecosystem services.Guangdong province was chosen as an example since it is one of the most rapidly developing provinces in China that has experienced rapid urbanization, and most of its natural lands had been changed and developed to satisfy the needs of socioeconomic development [27].This study is expected to help policy makers to better understand the impacts of land use changes on ecosystem services, and to provide references for research in related fields.

Study Area
Guangdong province (20 • 13 N-25 • 31 N, 109 • 39 E-117 • 19 E), located at the southern part of mainland China, consists of four geographic and economic regions (Pearl River Delta, east Guangdong, west Guangdong, and north Guangdong), and has 21 cities.It is bordered by Hong Kong, Macau, Fujian, Guangxi, Jiangxi, and Hunan province, and is separated from Hainan province by the sea (Figure 1).The area of Guangdong province accounts for around 1.87% of China's total land area, and its average elevation is 198 m.The northern part is mostly hilly, while the southern part is mostly plain and tableland.The climate type of Guangdong province is mainly a subtropical monsoon climate, and the forest type is mainly a subtropical evergreen broad-leaved forest.As the southern gate of China, Guangdong province is the earliest birthplace of the Maritime Silk Road, located in the South China Sea shipping hub.Since the reform and opening-up of policies in 1978, Guangdong province has experienced rapid urbanization process and economic development, and it has become the most developed province in China [28].In 2021, Guangdong had the highest urbanization rate and gross domestic product (GDP), with its urbanization rate reaching 72.7% and its GDP reaching CNY 12.4 trillion, ranking the first in China.For more than 30 years, Guangdong province has been the most rapidly developing region in China, leading to its status as the province with the most typical land use change characteristics.

Data
The assessment of ecosystem services was based on multisource data, including land use data, meteorological data, socioeconomic data, soil property data, and topography data (DEM) etc., during the period 1990-2020.The details of the data required for this study are shown in Table 1.In order to leave land use change as the sole driver affecting ecosystem services changes, this study calculated the multi-year average meteorological indicators to keep climate constant from 1990-2020.Then, Kriging spatial interpolation was used to convert the meteorological data into raster data with a resolution of 1 km.The other raster data were also converted to a uniform resolution of 1 km.Data related to grain production were mainly from the statistical yearbooks of Guangdong province.

Methods
Guangdong province experienced rapid urbanization during 1990-2020, with builtup land expansion occupying the cultivated land, forest land, and water bodies, which had a great impact on the ecosystems.Referring to related studies and the principles of importance and data availability [14,30,31], five key ecosystem services (Table 2), including GP, WY, SC, HQ, and CS, were selected to analyze their spatiotemporal dynamics and responses to land use changes in Guangdong province during 1990-2020.Details about the methods used are as follows.
Table 2. Descriptions of the selected ecosystem services.

Ecosystem Services Description
Grain production (GP) GP is one of the most basic supply services of ecosystem and provides material guarantee for human survival and development [32].

Water yield (WY)
The provision of freshwater is one of the ecosystem service functions that plays a crucial role in promoting biological survival and ensuring ecological security [33].
Soil conservation (SC) SC is one of the important ecosystem services, representing the ability of ecosystems to protect soil and control erosion [34].
Habitat quality (HQ) Biodiversity is closely linked to the production of ecosystem services, and HQ is an expression of the ability of ecosystems to provide permanent living conditions for organisms [35].

Carbon sequestration (CS)
CS is considered as one of the most critical ecosystem services for assessing the response of productive capacity and ecological resilience to climate change [36].

Calculation of Grain Production
GP is mainly limited by cultivated land area and cultivation conditions.The cultivated land area was obtained by using the grid image of the proportion of cultivated land.Cultivation conditions were characterized by grain yield per unit area of each city.GP service in the study area was calculated using the following formula [29]: where P xy is the total GP of the cell (x,y) in cultivated land in units of tons; A xy is the area of cultivated land in the cell (x,y); and Y g C is the yield per unit area for grains on cultivated land (t/km 2 ).

Calculation of Water Yield
The water yield module of the InVEST is based on GIS raster data, and its core algorithm calculates the water production for each raster cell in a watershed using the water balance method in combination with climate, topography, soil characteristics, and land use parameters.The WY is the amount of rainfall for each grid cell in the area minus the actual evapotranspiration without upstream runoff recharge, where climatic elements, topographic factors, soil characteristics, and land use types mainly influence the balance between rainfall and evapotranspiration.The water yield module is based on the Budyko curve and the average annual rainfall.The annual WY of raster cells of different land use types Y ij is calculated as follows: where Y ij is the annual WY of land use type j in raster cell i; AET ij is the actual annual average evapotranspiration of land use type j in grid cell i; and P i is the average annual rainfall of raster cell i.
Since the actual annual evapotranspiration cannot be obtained by direct measurement, it can be approximated by Budkyo curves for . The approximate calculation is performed.
is an approximation of the Budkyo curve, which is calculated as follows: where R ij is the Budkyo dimensionless drying index for land use type j in raster cell i, which is the ratio of reference evapotranspiration to rainfall; and ω i is the ratio of the corrected annual available water for vegetation to the expected rainfall.The Budkyo drying index, R ij , is calculated by the following formula: where K ij is the vegetation evapotranspiration coefficient of land use type j in raster cell i; and ETO i is the potential evapotranspiration of raster cell i, also called the reference crop evapotranspiration, reflecting the evapotranspiration capacity determined by the climatic conditions.ω i is defined as a non-physical parameter describing natural climate-soil properties, and it is calculated as follows: where Z is the Zhang coefficient, which is used as a parameter representing seasonal rainfall distribution and rainfall depth, and is determined by regional rainfall characteristics; AWC i is the effective soil water content, also known as the vegetation available water content, and its value is determined by the soil texture and effective soil depth.

Calculation of Soil Conservation
The SDR sediment delivery ratio module of the InVEST model uses the US general soil loss equation for soil erosion (USLE), considering parameters such as rainfall erosion force (R), soil erodibility (K), slope length and slope gradient (LS), vegetation cover and management measures (C), and soil and water conservation measures (P).Using the land use type as the assessment unit, the SC amount of each assessment unit is obtained by subtracting the actual soil erosion amount from the potential soil erosion amount.The difference between the potential soil erosion amount and the actual soil erosion amount is that the former does not take into consideration the effect of vegetation cover and management factor, nor the soil and water conservation measure factor on soil erosion sequestration.That is, they have the same expression, but to calculate the potential soil erosion amount, the soil and water conservation measure factor (P) and vegetation cover and management factor (C) are assigned to 1.The potential soil erosion and actual soil erosion are calculated as follows: where RKLS is the potential soil erosion; USLE is the actual soil erosion; R is the rainfall ero- LS is the slope length and slope gradient factor; C is the vegetation cover and management factor; and P is soil and water conservation measure factor.SC quantity is the amount of RKLS minus USLE: where SC is the annual soil retention, t•km −2 •a −1 .

Calculation of Habitat Quality
The HQ module of the InVEST model was selected for the assessment.The calculation of HQ includes two parts: grid external stress and grid habitat suitability.The grid external stress D xj is calculated as follows: where D xj is the level of habitat stress in raster x in land cover (or habitat type) j; w r is the weight of the stressor, indicating the relative destructive power of a stressor to all habitats; r y is used to determine whether raster y is the source of the threat factor r; i rxy is the coercive effect of the stressor rxy in raster y on the habitats in raster x; β x is the accessibility of raster x under social, legal, and other protection states; and S jr is the sensitivity of land cover j to the stressor r, with higher values indicating high accessibility.
where d xy is the linear distance between grid x and grid y, and d rmax is the maximum influence distance of stress factor r.After considering the grid habitat suitability, the HQ was calculated as follows: where, Q xj is the HQ of raster x in land cover j; H j is the habitat suitability of land cover j for the species, with values from 0 to 1, with 1 indicating the strongest habitat suitability; k is the half-saturation constant, and when 1 − 5, the value of k equals the value of D; and the default value of z is 2.5.

Calculation of Carbon Sequestration
The carbon module of InVEST model was used to calculate CS.Its core algorithm applies the inventory method to assign a minimum estimated amount of carbon to each land use type for the four basic carbon pools (aboveground biomass (C above ), belowground biomass (C below ), soil (C soil ), and dead organic matter (C dead )), and measures the CS capacity of the region by counting the total estimated amount of carbon assigned to each land use type per unit area [37].The carbon pool values assigned for each land use type were obtained from the InVEST model.The carbon density (C i ) for land use type i is equal to the sum of aboveground, belowground, soil carbon, and dead carbon densities for land use type i, represented as follows.
The total carbon storage (C total ) is equal to the sum of carbon density for land use type i multiplied by the area (A i ) for land use type i, with n as the number of land use types.In this study, seven phases of land use data of Guangdong province based on remote sensing image interpretation in 1990,1995,2000,2005,2010,2015, and 2020 were used to explore the characteristics of dynamic changes of land use.Forest land and cultivated land are the main land use types in Guangdong province, which accounted for 59.87% and 23.81% of the total land area in 2020, respectively.Following were built-up land (7.57%) and water areas (4.40%).Concerning the spatial pattern of land use, forest land and grassland were mainly concentrated in Zhaoqing, Qingyuan, and Heyuan, which are located in the middle and upper reaches of Xijiang River, Suijiang River and Dongjiang River.Among them, forest land is mainly distributed in the low mountainous and hilly areas, while grassland is distributed in the transition zone between cultivated land and forest land.Cultivated land and built-up land, showing a trinuclear distribution pattern, are concentrated in the urban cluster in the Pearl River Delta (Guangzhou, Shenzhen, Dongguan, Zhongshan, Zhuhai, and Foshan), the southwestern part of Guangdong (Maoming and Zhanjiang), and the northeastern part of Guangdong (Shantou).Furthermore, there is extraordinarily little unused land (Figure 2).The quantitative structure and spatial patterns of land use in Guangdong province reflect its regional landscape pattern of "forests and grasslands in low hills, cultivated land, and built-up land in plains" with high exploitation rate.

Analysis of Land Use Changes in Guangdong Province
With the development of social economy over the past 30 years, the intensity of human production activities had continued to increase, resulting in significant changes in land use.The analysis about changes of areas of land use types can provide the general change trend of land use and its structure.In general, land use changes in Guangdong province from 1990 to 2020 mainly showed a significant expansion of built-up land and water areas, with the area of built-up land almost doubling in size, while the area of cultivated land, forest land, grassland, and unused land decreased in different degrees (Table 3).
Considering the development stages of periods, the changes of land use structure in Guangdong province had the following characteristics.Over the past 30 years of rapid economic development, Guangdong has rapidly expanded its cities.As shown in Figures 3 and 4, from 1990-2020, the area of built-up land has been continuously expanding, and the expansion rate has been stable at a high level.From 1990From -1995From , 2000From -2005From , 2005-2010, 2010-2015, and 2015-2020, the area expanded by 1875 km 2 (27.7%), 2251 km 2 (27.2%), 900 km 2 (8.5%), 784 km 2 (6.9%), and 1193 km 2 (9.8%), respectively (Figure 3).Only within the period 1995-2000 did the area of built-up land decrease, by 366 km 2 (−4.2%).In the past 30 years, the overall expansion rate of built-up land in Guangdong province has reached 97.9%, indicating that the urban area has been continuously ex-panding, and the intensity of human activities, as well as the impact on the local natural environment, has been continuously increasing.Water areas is another land use type that had expanded to a certain extent.The water areas increased by 941 km 2

Changes of Ecosystem Services in Guangdong
Province from 1990 to 2020

Changes of Grain Production
Grain is one of the essential elements for human survival and development, and GP serves the most important supply service of the ecosystem.On the one hand, grain output can reflect the comprehensive situation of the regional ecological environment.On the other hand, the change in GP can also reflect the human attitude towards the use of the ecosystem.GP services and changes in Guangdong province from 1990-2020 are shown in Figure 5a-c.

Changes of Water Yield
The ability to provide sufficient freshwater is a key element in evaluating whether the regional ecological environment is excellent.WY depends on the regional rainfall and water harvesting capacity.Changes in land use patterns affect climate and water distribution by altering the water cycle (e.g., urban construction hardens the subsurface, resulting in less infiltration and more water being retained and collected), while the construction of reservoirs and other water storage facilities can enhance water collection capacity to alter regional WY.The WY and its change in Guangdong province from 1990-2020 are shown in Figure 5d-f.From the perspective of spatial distribution, the WY from 1990-2020 generally followed the distribution pattern of rainfall, showing a gradual decrease from coastal areas to inland areas, while the results also clearly indicated that the WY had a close relationship with the nature of the underlying surface.Since most of the area was an impermeable urban hardening surface, the WY in the central urban agglomeration area breaks the laws of the following rainfall distribution, and is significantly higher than that in the surrounding forest land and grassland areas.Moreover, since most of the area was cultivated land with high infiltration and high-water consumption intensity, the WY in Zhanjiang was significantly lower than that in other areas with the same rainfall conditions.Specifically, the three cities with the lowest average WY from 1990-2020 were Chaozhou (502 mm), Zhanjiang (504 mm), and Meizhou (555 mm), and the three cities with the highest values were Dongguan (1097 mm), Shenzhen (1075 mm), and Guangzhou (1037 mm).From a temporal perspective, the average WY increased slightly from 716 mm in 1990 to 738 mm in 2020, and all cities in Guangdong province except Shantou and Zhanjiang showed increase trend with different magnitudes.Among them, Shenzhen had the largest growth rate of 24.59%, followed by Dongguan and Zhongshan, with an increase rate of 24.29% and 20.62%, respectively.Shantou and Zhanjiang had a very low decline of 0.73% and 0.22%, respectively.The regions with rapid growth in WY were mainly concentrated in areas with high precipitation and rapid expansion of urban hardened subsurface, such as the Pearl River Delta urban agglomeration.The regions with unchanged or decreasing WY were mainly scattered in the surrounding areas of forest land and cultivated land, and the area was originally used for construction.

Changes of Soil Conservation
Soil is an important cornerstone of ecosystem service delivery, closely inferenced by climate change, food production, and vegetation growth.SC capacity is limited by the relative relationship between soil erodibility and environmental erodibility, and changes in land use types will cause changes in the relative relationship between this pair of influencing factors.SC and its change from 1990 to 2020 are shown in Figure 5g-i.From the perspective of spatial distribution, the SC clearly showed a high distribution pattern o in forest land and grassland, and a low distribution in built-up land.Specifically, urban hard substrates had almost no SC capacity, and the small amount of SC was mainly provided by urban green areas, such as parks and urban forests.Forest land and grassland benefited from the soilfixing capacity of plant roots, and had better SC capacity.Specifically, the three cities with the lowest average SC from 1990-2020 were Zhanjiang (462 t•km −2 ), Foshan (1735 t•km −2 ), and Zhongshan (2569 t•km −2 ), and the three cities with the highest average SC were Shanwei (17,713 t•km −2 ), Qingyuan (17,038 t•km −2 ), and Shaoguan (14,324 t•km −2 ).From a temporal perspective, the average SC in Guangdong province decreased slightly from 10,531 t•km −2 in 1990 to 10,529 t•km −2 in 2020, with a total of 9 cities showing a decreasing trend, and the remaining 12 cities showing an increasing trend.The changes were extremely small, concretely all less than 1.5%.Among them, Zhanjiang had the largest increase of 1.28%, followed by Dongguan and Zhuhai with 0.67% and 0.47%, respectively.The largest decreases were in Chaozhou and Shanwei, with 0.48% and 0.32%, respectively.Over the past 30 years, most of the regions in Guangdong province had experienced a slight decrease in SC, with some mountainous areas sporadically distributed with obvious decreasing areas.The areas with rising SC mainly occured in areas with rapid urban development.

Changes of Habit Quality
The continuous development of cities will change the regional land use pattern, and the reduction of forest land and grassland will reduce the suitable space for maintaining biodiversity, resulting in the deterioration of HQ.The HQ and its change in Guangdong province from 1990-2020 are shown in Figure 5j-l.From the perspective of spatial distribution, the HQ showed a distribution pattern of low in the east and high in the west and north, with the rapid development of core cities leading to the continuous increase in the destructive effect of human activities on the ecological environment, thereby leading to a decrease in HQ.Specifically, the three cities with the lowest average HQ in 1990-2020 were Dongguan (0.42), Zhongshan (0.47), and Foshan (0.48), while the three cities with the highest average HQ were Heyuan (0.86), Zhaoqing (0.83), and Meizhou (0.83).On the urban scale, the HQ in the core urban areas was much lower than that in the suburban areas.Taking Guangzhou as an example, we used the third ring of Guangzhou as the boundary between core urban areas and suburban areas.The average HQ of core urban areas in Guangzhou was 0.26, which was much lower than that of suburban areas at 0.89.From a temporal perspective, the average HQ in Guangdong province decreaseed slightly from 0.73 in 1990 to 0.71 in 2020.All cities in Guangdong except for Maoming, Zhanjiang, and Shantou showed different magnitudes of decrease.Among them, Dongguan had the largest decrease rate of 34.62%, followed by Shenzhen and Foshan with a decrease rate of 32.81% and 24.53%, respectively.Maoming, Zhanjiang, and Shantou had extremely low increases of less than 1%.The areas with dramatic deterioration in HQ were mainly concentrated in areas with rapid urban development and built-up land expansion, such as the Pearl River Delta urban agglomeration.At the same time, the areas with improved HQ were mainly scattered in forest land areas or clustered in cultivated land expansion areas (e.g., Zhanjiang).

Changes of Carbon Sequestration
Ecosystems regulate regional climate conditions by increasing or decreasing atmospheric greenhouse gases (e.g., carbon dioxide).Therefore, the CS capacity of regions has an important impact on the local ecological environment.Forest land and grassland are the most important carbon pools, and changes in land use types will significantly alter the regional CS service.The CS and its change in Guangdong province from 1990-2020 are shown in Figure 5m-o.From the perspective of spatial distribution, the CS generally showed high values in forest land, followed by grassland and cropland, and low values in built-up land and water areas.It indicated that vegetation and soil were the most important carbon pools in the ecosystem [38,39].Forest land had a strong CS capacity due to its large amount of vegetation and deep soil, while grassland and cultivated land also had a certain amount of vegetation and soil with strong CS capacity.Most of the built-up land were hard substrates with almost zero CS capacity, while water areas relies almost only on water to absorb carbon dioxide to provide a little CS.Thus, their CS values were much lower than those of forest land, grassland, and cultivated land.Specifically, the three cities with the lowest average CS from 1990-2020 were Zhongshan (1437 t•km −2 ), Foshan (1463 t•km −2 ), and Dongguan (1535 t•km −2 ), and the three cities with the highest average CS were Heyuan (2897 t•km −2 ), Meizhou (2864 t•km −2 ), and Zhaoqing (2827 t•km −2 ).From a temporal perspective, the average CS increased from 2486 t•km −2 in 1990 to 3001 t•km −2 in 2020, and all cities in Guangdong province except for Dongguan and Shenzhen showed an increasing trend.Among them, Meizhou had the largest growth rate of 23.44%, followed by Yunfu and Shaoguan, with values of 23.38% and 23.02%, respectively.The CS in Dongguan and Shenzhen had decreased by 6.94% and 4.26%, respectively.Over the past 30 years, the rise of CS in most regions of Guangdong was within the range of 100-1000 t•km −2 , and there were significant decreases (>1000 t•km −2 ) in new areas of urban expansion.However, in existing urban areas, CS may increase as well as decrease.

Impact of Land Use Change on Ecosystem Services
In order to further analyze the impact of land use change on the five key ecosystem services, our study superimposed the spatial distribution of land use conversions during the period 1990-2020 with the spatial distribution of ecosystem service changes during the same time.Based on spatial statistical analysis tools, we obtained the statistical information of ecosystem service changes corresponding to different land use conversion types (Figure 6).GP increased in the areas converted from other land use types to cultivated land, while in the areas with cultivated land converted to other types, it showed a decreasing trend.In the area of conversion from cultivated land to built-up land, grain production was the most deteriorated (−150.9t•km −2 ), and in areas of conversion from grassland to cultivated land, GP was the most optimized (78.8 t•km −2 ).GP was directly restricted by the area of cultivated land and farming conditions.Over the past 30 years, the rapid urbanization in Guangdong province with the rapid decrease of cultivated land area has led to the general decline of regional grain output.Cultivated land which was converted from grassland, forest land, and water areas had better GP values.This is due to the fact that grassland, forest land, and water areas themselves have relatively deep soil and moisture conditions, providing a good foundation for cultivated land after conversion.
WY it mostly showed an overall decreasing trend in areas converted to cultivated land, forest land, and water areas, while showed an overall increasing trend in areas converted to built-up land and unused land.Among them, the WY values in the areas that were converted from unused land to water areas were the most deteriorated (−1072 mm), and in the areas that were converted from water areas to built-up land, it was the most optimized (943 mm).It is worth noting that the WY in our study is mainly used to evaluate the ability of WY based on natural precipitation, that is, the ability to retain precipitation in situ.WY in water areas is the lowest because natural precipitation leaves the site in a short time with runoff, or evaporates violently in lakes.In cultivated land and forest land, precipitation infiltrates into the soil rapidly and is used by plants, thus maintaining a low WY.In built-up land, hard substrates prevent the infiltration and loss of precipitation.In addition, WY can be effectively improved by the city's well-developed water diversion and storage projects, such as sewers and cisterns.
For SC, there was an overall decreasing trend in areas converted to cultivated land and unused land, and an increasing trend in areas converted to forest land, water areas, and built-up land.The conversion from grassland to unused land is the most deteriorated (−10,166 t•km −2 ), and the conversion from unused land to forest land is the most optimized (3552 t•km −2 ) among all the conversion types.Overall, as forest land, water areas, and built-up land had higher soil retention capacity, the conversion from other land use to them would help to increase SC.In addition, for the mutual conversion between two land use types, the degree of change of SC was different.For example, the conversion from unused land to grassland led to an average SC increase of 2672 t•km −2 , while the conversion from grassland to unused land led to an average SC decrease of 10,166 t•km −2 .The main reason may be that land use conversion occurred at different positions, where the SC capacity would be influenced by its topography, soil texture, and regional climate conditions.
HQ it showed an overall decreasing trend in areas converted to cultivated land, builtup land, and unused land, while showing an overall increasing trend in areas converted to forest land, grassland, and water areas.The HQ in the areas converted from forest land to built-up land was the most deteriorated (−0.65), and the areas which were converted from built-up land to forest land was the most optimized (0.40).Forest land has excellent natural conditions, where a large number of plants, animals, and microorganisms living in a stress-free environment, leading to a high HQ rating.However, in the process of urban development and expansion, a large amount of forest land was converted to built-up land, and the original natural environment was deteriorated, along with its ecosystem structure being destroyed.This resulted in a sharp decline in HQ evaluation.However, some of the built-up land, such as abandoned industrial and mining areas, has been restored to forest land, which can significantly improve the local HQ on the contrary.
CS decreased in the areas converted to water areas, built-up land, and unused land, while increased in the areas converted to cultivated land, forest land, and grassland.The CS in the areas converted from forest land to water areas was the most deteriorated (−2974 t•km −2 ), and most optimized in the areas converted from water areas to forest land (3338 t•km −2 ).Cultivated land and grassland also have better CS values for the same reason.In water areas, water is the main carrier of CS, but in conventional environments water has little capacity to sequester carbon.In built-up land, a large amount of concrete, asphalt, and other major constituent materials have little CS capacity.Therefore, both of them have very poor CS values.

Discussion
Our study finds that land use and its change has had an extremely important influence on ecosystem services, and this influence was particularly evident in the changes of forest land and built-up land.In the context of rapid urbanization, built-up land changed the most dramatically.Followed by two distinct policies of deforestation and afforestation in different regions, it led to significant spatial heterogeneity in the change of forested land.The difference above had profound impacts on ecosystem services.The study of Liu et al. concerning the PRD urban agglomeration has similar findings [14].However, our study found that forest land in Guangdong Province has the best performance in the evaluation of ecosystem service values, while the study of Xie et al., shows that grassland is the best in the Aksu Basin [17].We believe that land use/cover contributes differently to the value of ecosystem services under natural conditions (e.g., climate, topography, etc.), which is closely related to the area of various land use/cover patterns in the region and the degree of interaction with other natural elements.It is important to note that our study only considered the "natural" component of ecosystem service value.In fact, humans' social, cultural, and emotional needs concerning the natural environment are also an important part of the value.The study of Tasser et al., included a questionnaire survey on the social value of ecosystem services, which is instructive for our subsequent research [40].
Ecosystem services respond differently to land use change, exhibiting local environmental influences.In this study, GP, SC, and HQ all tended to deteriorate overall over the 30-year period, deteriorating or optimizing locally (GP and HQ deteriorated particularly in urban built-up lands, while SC was optimized in urban built-up lands).WY and CS tended to be optimized overall over the 30-year period, with a significant trend of deterioration locally (WY decreased mainly in the surrounding areas of forest land and cultivated land, while CS significantly decreased in new areas of urban expansion).This kind of change in ecosystem services suggests that land use change does not simply change the landscape pattern of the terrestrial surface, but may also affect important natural cycles, such as the water cycle, the energy cycle, and the material cycle, thus changing the overall ecological environment, and causing an overall deterioration or optimization of ecosystem services.Hasan et al. integrated an analysis of studies related to land use and ecosystem services, and concluded that land use can greatly influence various elements of ecosystems, from biodiversity to climate, and, in this way, affect ecosystem services [41].Our study does not specifically analyze the effects of land use on various elements of ecosystems, but our results of changes in ecosystem services support the conclusions of Hasan et al.
It is worth noting that we did not make a qualitative judgment of "need to maintain or improve", since ecosystem services are defined as the "benefits" that humans can derive from ecosystems [42][43][44].This means that the value of ecosystem services should be comprehensively decided upon with consideration of human needs, and not just considered as a combination of natural attributes.Therefore, we cannot ignore the needs of economic and social development to pursue better ecosystem services, but we should make more reasonable development plans to achieve sustainable development while meeting our own economic and social development needs.
Our study area is Guangdong province, the most rapidly urbanized and developed province in China.As such, the quantitative results of the study are distinctly geographic in nature, and the conclusions drawn are not fully applicable to other regions.Moreover, we only considered the direct impacts of land use and its change on ecosystem services, and did not include socioeconomic factors in the analysis.Thus, we were not able to show the socio-economically driven land use changes and their impacts on ecosystem services.Future studies can introduce socioeconomic factors to explore the impacts of land use change on ecosystem services in a more integrated and comprehensive manner.

Conclusions
In this study, our aim was to assess changes of land use and its impacts on typical ecosystem services (GP, WY, SC, HQ and CS) in Guangdong province in 1990, 1995, 2000, 2005, 2010, 2015, and 2020.The results of the analysis indicated that there had been significant changes in land use.Specifically, cultivated land, forest land, grassland, water areas, built-up land, and unused land correspondingly changed at a rate of −10.7%, −1.9%, −5.1%, 13.7%, 97.9%, and −38.8% from 1990 to 2020, indicating that Guangdong province had experienced rapid urban area expansion.Responding to land use changes, the average GP, WY, SC, HQ, and CS had correspondingly changed by −8.66% (−12.3 t•km −2 ), 3.10% (22 mm), −0.02% (−2 t•km −2 ), −2.74% (−0.02), and 20.70% (515 t•km −2 ) from 1990 to 2020, showing different spatial heterogeneity characteristics.Land use changed dramatically in relation to the increasing of urbanization in Guangdong province, and affected the ecosystem services directly.Among them, the types of land use conversion with the average greatest negative impacts on GP, WY, SC, HQ, and CS were cultivated land to builtup land (−150.9t•km −2 ), unused land to water areas (−1072 mm), grassland to unused land (−10,166 t•km −2 ), forest land to built-up land (−0.65), and forest land to water areas (−2974 t•km −2 ), respectively, and with average greatest positive impacts were grassland to cultivated land (78.8 t•km −2 ), water areas to built-up land (943 mm), unused land to forest land (3552 t•km −2 ), built-up land to forest land (0.40), and water areas to forest land (3338 t•km −2 ), respectively.Therefore, it is evident that land serves as a cornerstone for the provision of ecosystem services that meet human requirements.In other words, land use patterns have unignorable impacts on ecosystem services.In future studies, we will add the social element of ecosystem services into the value evaluation, and the balance between the value of ecosystem services and the coordination of socioeconomic development will also be considered.

Figure 2 .
Figure 2. Land use of Guangdong province in 2020.

Figure 3 .
Figure 3. Area changes of each land use type in Guangdong, 1990-2020.

Figure 4 .
Figure 4. Area change rates of each land use type in Guangdong, 1990-2020.

Figure 6 .
Figure 6.Average change of the five key ecosystem services corresponding to different land use conversion types in Guangdong province, 1990-2020.

Table 1 .
Input data for ecosystem services assessment in Guangdong province.