Identifying and Mapping the Responses of Ecosystem Services to Land Use Change in Rapidly Urbanizing Regions: A Case Study in Foshan City, China

: Rapid urbanization has degraded some important ecosystem services and threatens so-cioeconomic sustainability. Although many studies have focused on the effect of urbanization on ecosystem services, the effect and its threshold have not been well-identiﬁed spatially. In this study, we propose a research framework by integrating nighttime light data, the InVEST (Integrated Valuation of Environmental Service and Tradeoffs) model, and a spatial response index to characterize the response of ecosystem services to rapid urbanization. We considered Foshan City as a case study to explore the effect of rapid urbanization on ecosystem services during 2000–2018. Our results showed that rapid urbanization resulted in a 49.13% reduction in agricultural production and a 10.13% reduction in habitat quality. The spatial response index of agricultural production, habitat quality, soil retention, water yield, and carbon sequestration were 14.25%, 2.94%, 0.04%, 0.78%, and 0.07%, respectively. We found that developing urban areas had the highest spatial response index, indicating that this area was the crucial area for future land management. We consider that our research framework can help identify the key areas affected by rapid urbanization. Visualizing the spatial response index and extracting the threshold for different levels of urbanization will be conducive to sustainable urban management and planning.


Introduction
Urbanization is one of the most prominent features of contemporary global environmental change [1,2], and Chinese cities are undergoing rapid urbanization [3]. This is evident from the increase in the proportion of the urban population from 17.9% to 52.6% between 1978 and 2012 [4]. In addition, urbanization of the population has been accompanied by considerable land-use changes. These areas have suffered from environmental pollution, ecological deterioration, and economic loss [5,6]. Recently, the effect of rapid urbanization on ecosystems and resilient socio-ecological systems has attracted attention, with previous studies focusing on urban agriculture [1,7,8], urban forests [9,10], urban water [11,12], and habitat quality [13,14]; maintaining a healthy environment is relevant for sustainable urban construction and planning.
Complex urban socio-ecological systems result from the interaction between humans and nature. When planning sustainable cities, we need to move beyond city limits and consider the entirety of the human-dominated system, which depends on natural ecosystems [15]. Nighttime light data are closely linked to human-based systems and are typical signals of regional economic and social activities. Such data have been widely used to determine the extent of built-up land and spatio-temporal characteristics of urbanization dynamics [16,17]. These studies mainly focus on the relationship between nighttime light data and urban construction land expansion [18,19]. Meanwhile, with the expansion of the application range of nighttime light data, these data are applied to more research on the impact of urbanization, such as vegetation primary productivity [20], habitat quality [21], and other ecological effects [22]. Many previous studies use the digital numbers (DN) value of nighttime light data to identify urbanization areas as a main method. However, how to identify urban developed areas, developing areas, and rural areas through threshold division is a still problem worth discussing. Therefore, this study will use population, economy, and land use data to link the nighttime light data with the urbanization levels, and then to identify different levels of urbanization areas.
ES are recognized as the benefits people obtain from ecosystems, which incorporate society and ecosystems and support the survival and development of a population's wellbeing [23][24][25]. Studies on the quantification of ES mainly include physically quantifying services and assessing the economic value of interest [24,26]. These quantitative assessments, modeling, and mapping technologies of ES have been widely used for regional sustainable development [27,28], especially in rapidly urbanized areas [25,29]. In the previous studies, land use and land cover change (LULC) data were widely used as basic data to evaluate the ES changes in the process of urbanization [30][31][32][33]. These case studies, such as by Lyu [30], Long [31], Peng [32], and Liu [33], have gradually linked land use, ecosystem services, and urbanization. However, how ecosystem services respond to land use changes due to urbanization is not yet clear. In the meantime, from the perspective of land ecological management, the demand for quantitative and spatial visualization research is increasing. Therefore, it is necessary to clarify the quantitative and spatial response relationship between changes in ecosystem services and land use changes. In our study, we defined a spatial response index to characterize the impact of rapid urbanization on ES, and a threshold was spatially visualized to assist with urban planning and management.
The present study focused on Foshan City in the Pearl River Delta region, which is a typical representative of rapid urbanization in China. Based on the statistics, the unprecedented rate of urbanization increased from 75.06% to 95.00% in 2000-2018 and has had a considerable effect on Foshan. As an important city in the Guangdong-Hong Kong-Macao Greater Bay Area, the rapid development of advanced manufacturing has also affected the ecological environment quality, agricultural production, and human settlement of Foshan. Our objective was to propose a research framework to integrate nighttime light data, the Integrated Valuation of Environmental Service and Tradeoffs (InVEST) model, and a spatial response index to address the effect of rapid urbanization on ES. Based on multi-source data, we established the relationship between the nighttime light data index and the urbanization level to identify the different levels of urbanization areas. We aimed to resolve the following specific research problems: (1) How can the relationship between the nighttime light data index and the urbanization level be established to identify urbanization areas? (2) Which kind of ES is most affected by rapid urbanization? Which area exhibits the most considerable changes in ES and responsiveness of ES to rapid urbanization? We quantified and mapped five vital ES, including agricultural production, carbon sequestration, water yield, soil retention, and habitat quality under three different urbanization levels (developed urban, developing urban, and rural areas) between 2000 and 2018. We considered that our results and maps can be used for ecosystem protection or restoration activities to improve sustainable urban planning and management.

Study Area
Foshan is spread across an area of approximately 3792.65 km 2 in the south-central area of Guangdong Province, southern China (22 • 38 -23 • 34 N, 112 • 22 -113 • 23 E). The Xijiang River, Beijiang River, and tributaries of the Pearl River run through Foshan. The study area has a southern subtropical monsoon climate with an average annual temperature of approximately 24 • C and average annual precipitation of 1600-2000 mm [34]. Abundant hydrothermal conditions have contributed to a distinct "dike-pond" agriculture-aquaculture multiplex system to achieve renewable and sustainable development of ecological agriculture [35]. Administratively, Foshan contains five districts: Sanshui, Nanhai, Chancheng, Shunde, and Gaoming ( Figure 1). Since the implementation of the Chinese economic reform policy in 1979, Foshan has experienced unprecedented rapid urbanization and economic growth. The urbanization rate in this area has been far higher than that across the rest of China [36,37]. During 2000-2018, the permanent population increased from 5.34 million to 7.96 million, the rate of urbanization increased from 75.06% to 95.00%, and the gross domestic product increased from 105.04 billion to 993.59 billion yuan [38]. study area has a southern subtropical monsoon climate with an average annual temperature of approximately 24 °C and average annual precipitation of 1600-2000 mm [34]. Abundant hydrothermal conditions have contributed to a distinct "dike-pond" agriculture-aquaculture multiplex system to achieve renewable and sustainable development of ecological agriculture [35].
Administratively, Foshan contains five districts: Sanshui, Nanhai, Chancheng, Shunde, and Gaoming ( Figure 1). Since the implementation of the Chinese economic reform policy in 1979, Foshan has experienced unprecedented rapid urbanization and economic growth. The urbanization rate in this area has been far higher than that across the rest of China [36,37]. During 2000-2018, the permanent population increased from 5.34 million to 7.96 million, the rate of urbanization increased from 75.06% to 95.00%, and the gross domestic product increased from 105.04 billion to 993.59 billion yuan [38].

Research Framework
We developed a framework to integrate nighttime light data, InVEST, and the spatial response index of urban ES. This framework was used to address problems pertaining to the effects of rapid urbanization on vital regional ES and the responsiveness of ES to changes in LULC ( Figure 2). The spatial response index and spatial change patterns of ES can help policy makers identify key areas that suffer from increased urbanization.
First, we constructed a series of consistent nighttime light datasets by integrating the Defense Meteorological Satellite Program's Operational Line-scan System (DMSP-OLS) and the Suomi National Polar-orbiting Partnership Visible Infrared Imaging Radiometer Suite (NPP-VIIRS). The nighttime light data from 2000 and 2018 were used to classify the urbanization level of Foshan and to identify the rural, developing, and developed areas. Second, we used the InVEST and empirical models to assess the changes in agricultural production, carbon sequestration, water yield, soil retention, and habitat quality. We quantitatively evaluated the changes in LULC and ES from 2000 to 2018 and calculated the spatial response index under different urbanization levels. Finally, the various ES were normalized, reclassified, and overlaid to determine the total ES. The spatial distribution and zoning thresholds of the response of total ES to urbanization were mapped.

Research Framework
We developed a framework to integrate nighttime light data, InVEST, and the spatial response index of urban ES. This framework was used to address problems pertaining to the effects of rapid urbanization on vital regional ES and the responsiveness of ES to changes in LULC ( Figure 2). The spatial response index and spatial change patterns of ES can help policy makers identify key areas that suffer from increased urbanization.
First, we constructed a series of consistent nighttime light datasets by integrating the Defense Meteorological Satellite Program's Operational Line-scan System (DMSP-OLS) and the Suomi National Polar-orbiting Partnership Visible Infrared Imaging Radiometer Suite (NPP-VIIRS). The nighttime light data from 2000 and 2018 were used to classify the urbanization level of Foshan and to identify the rural, developing, and developed areas. Second, we used the InVEST and empirical models to assess the changes in agricultural production, carbon sequestration, water yield, soil retention, and habitat quality. We quantitatively evaluated the changes in LULC and ES from 2000 to 2018 and calculated the spatial response index under different urbanization levels. Finally, the various ES were normalized, reclassified, and overlaid to determine the total ES. The spatial distribution and zoning thresholds of the response of total ES to urbanization were mapped.

Data Collection
We used nighttime light data to quantify and classify the level of urbanization in Foshan. The DMSP-OLS light product in 2000 and NPP-VIIRS monthly product in 2018 were downloaded from the website of the Earth Observation Group. Owing to the differences between NPP-VIIRS and DMSP-OLS images, we constructed a series of consistent nighttime light datasets through integration and calibration. The inter-calibration method was referred to in the literature [39]. Nighttime light data from two sensors in the overlapping year of 2012 were used to construct the calibration function ( Figure 3). The calibrated DN ranged from 0 to 63, with higher values representing higher levels of urbanization [40,41].

Data Collection
We used nighttime light data to quantify and classify the level of urbanization in Foshan. The DMSP-OLS light product in 2000 and NPP-VIIRS monthly product in 2018 were downloaded from the website of the Earth Observation Group. Owing to the differences between NPP-VIIRS and DMSP-OLS images, we constructed a series of consistent nighttime light datasets through integration and calibration. The inter-calibration method was referred to in the literature [39]. Nighttime light data from two sensors in the overlapping year of 2012 were used to construct the calibration function ( Figure 3). The calibrated DN ranged from 0 to 63, with higher values representing higher levels of urbanization [40,41].

Data Collection
We used nighttime light data to quantify and classify the level of urbanization in Foshan. The DMSP-OLS light product in 2000 and NPP-VIIRS monthly product in 2018 were downloaded from the website of the Earth Observation Group. Owing to the differences between NPP-VIIRS and DMSP-OLS images, we constructed a series of consistent nighttime light datasets through integration and calibration. The inter-calibration method was referred to in the literature [39]. Nighttime light data from two sensors in the overlapping year of 2012 were used to construct the calibration function ( Figure 3). The calibrated DN ranged from 0 to 63, with higher values representing higher levels of urbanization [40,41]. We used data from 2000 and 2018 because the identification of land-use changes and changes in ES requires a long-term process. The data from 2018 were part of the latest and most complete data series that we could obtain. This period was also the period that experienced the fastest urbanization in Foshan. We utilized nighttime light data with 30 arc-second grid spatial resolution (approximately 1 km). Gross domestic product (GDP) and population density data were also at 1 km spatial resolution. The other data category included LULC data at 300 m spatial resolution, the digital elevation model at 30 m spatial resolution, soil data at 1 km, NDVI data at 1 km, and meteorological data from an observation station, which were resampled to 300 m for LULC and ES analysis. The data category and sources are shown in Table 1.

Classification of Urbanization Levels
To establish the relationship between nighttime light data and the degree of urbanization development, we used the DN threshold to identify spatial regions under different levels of urbanization development over a period of time (2000-2018) [42]. Based on population density, land use intensity, and GDP, we constructed a comprehensive urban development index that characterized Foshan's urbanization level. Land use intensity refers to the percentage of construction land area to the total area of the region, reflecting the intensity of human activities [43]. The greater the proportion of construction land area, the greater the intensity of human activities and the higher the level of urbanization. Population density can reflect the degree of population concentration in a city. The greater the population density, the higher the level of urbanization. GDP reflects the economic level of the region. The higher the GDP, the higher the level of urbanization. Due to the obvious differences in population density and GDP per area in some developed areas, the two indicators fluctuate greatly. Therefore, with the method of logarithm in statistics, the weakened volatility characteristics are substituted into the calculation to facilitate subsequent analysis. Based on the results of comprehensive urbanization indicators, the Jenks natural breakpoint method was used to divide the urbanization level into three levels: high, medium, and low, and the spatial distribution of Foshan's urbanization levels in 2000 and 2018 were obtained ( Figure 4). The formula for calculating the comprehensive level of urbanization is as follows: where D i is the comprehensive level of urbanization in year i, x i is the land use intensity, y i is the population density, and z i is the GDP in Foshan.
of the region. The higher the GDP, the higher the level of urbanization. Due to the obvious differences in population density and GDP per area in some developed areas, the two indicators fluctuate greatly. Therefore, with the method of logarithm in statistics, the weakened volatility characteristics are substituted into the calculation to facilitate subsequent analysis. Based on the results of comprehensive urbanization indicators, the Jenks natural breakpoint method was used to divide the urbanization level into three levels: high, medium, and low, and the spatial distribution of Foshan's urbanization levels in 2000 and 2018 were obtained ( Figure 4). The formula for calculating the comprehensive level of urbanization is as follows: where Di is the comprehensive level of urbanization in year i, xi is the land use intensity, yi is the population density, and zi is the GDP in Foshan. Based on the results of urbanization levels in Foshan in 2000 and 2018, we extracted the average DN value of nighttime light data in three areas with different levels of urbanization (high, medium, and low) as the threshold. Based on the calibrated nighttime light data, recorded DN and calibrated DN were used as indicators to identify the urban sprawl in Foshan from 2000 (DN2000) to 2018 (DN2018). We confirmed the DN threshold and classified Foshan into three different levels of urbanization: developed urban, developing urban, and rural areas ( Table 2).   (Table 2).

ES Assessment and Mapping
We focused on five critical ES, namely, provisioning services (agricultural production), regulating services (carbon sequestration, water yield, and soil retention), and supporting services (habitat quality). The regulating and supporting services were assessed by InVEST, which is an assemblage model jointly developed by Stanford University, the World Wildlife Fund, and the Nature Conservancy under the auspices of the Natural Capital Project. ArcGIS version 10.2 and the InVEST model were integrated for the assessment. The model was spatially explicit, using maps as information sources and producing maps as outputs. The model can be used to analyze the effect of land use on ES based on given land-use maps and biophysical data for the region [44].

Agricultural Production
Agricultural production was used as an indicator of provisioning services. We defined the agricultural production value (based on the 2018 value of the USD) to include the total production value of agriculture, forestry, animal husbandry, and fisheries. Based on the LULC data and the Statistical Yearbook in Foshan (http://www.foshan.gov.cn/gzjg/ stjj/tjsj_1110964/lssj/index.html accessed on 8 May 2021), we included the value of food production on the grid. For different types of total output value, agriculture corresponded to cropland, forestry corresponded to woodland, animal husbandry corresponded to grassland, and fisheries corresponded to water bodies. The formula is as follows: where G i is the total output value of food type i for a grid, A i is the area of land-use type corresponding to food type i under the 300-m grid, F i is the total output value of food type i in Foshan, and S i is the land-use type area corresponding to food type i.

Carbon Sequestration
Carbon sequestration is an important indicator of the terrestrial carbon cycle and regional regulating services [45] and indicates the transfer of atmospheric CO 2 into different carbon pools [46]. As used in the InVEST model, carbon sequestration aggregates the amount of carbon stored in four carbon pools (aboveground biomass, belowground biomass, soil, and dead organic matter) according to the land-use maps and carbon density. The formula is as follows: where C t is the total carbon storage (MgC), C above is the carbon storage (MgC) in aboveground biomass, C below is the carbon storage (MgC) in belowground biomass, C soil is the carbon storage (MgC) in soil, and C dead is the carbon storage (MgC) in dead organic matter. The other carbon density parameters refer to previous studies [29,47,48].

Water Yield
The water yield in the InVEST model was defined as the amount of water runoff from the landscape. The water yield was calculated based on the principle of the Budyko curve and annual average precipitation [49]. The formula is as follows [50]: where Y xj is the annual water yield (mm) per pixel x of LULC j, AET xj is the annual actual evapotranspiration (mm) per pixel x of LULC j, and P x is the annual average precipitation (mm) per pixel x. These parameters were derived from meteorological observation station data (Table 1).

Soil Retention
The soil retention service was calculated based on the revised universal soil loss equation (RUSLE) [51]. The input data included geomorphology, climate, vegetation cover and management, and support practice. The module considers the interception capabilities of the ground. The formulae are as follows [50]: Remote Sens. 2021, 13, 4374 where USLE i is the amount of actual soil erosion on pixel i; RKLS i is the amount of potential soil erosion on pixel i; SD i is the amount of soil retention (t) on pixel i; R i is the rainfall erosivity (M J · mm (ha · h) −1 ), which is calculated by using the Wischmeier formula [52] based on the average monthly and annual precipitation; K i is the soil erodibility (t· ha · h (M J · ha · mm) −1 ), which is calculated by the erosion-productivity impact calculator model [53]; and LS i is the slope length-gradient factor.

Habitat Quality
Habitat quality is defined as a supporting ES that can provide conditions appropriate for individual and population persistence [54]. In this study, we integrated the habitat quality results of the InVEST model with NDVI and PM2.5 data (Table 1) to calculate the habitat quality in Foshan. These data and results were overlaid to get a comprehensive habitat quality map. The InVEST habitat quality model was based on a coarse-filter approach, which combined the information on LULC and threats to biodiversity [50]. These threats were considered as the sources of degradation, especially for human-modified LULC types. The sensitivity of each habitat type to degradation was based on general principles of landscape ecology and conservation biology [55,56]. A habitat quality landscape score was translated into a habitat quality index using a half-saturation function, which is simply the aggregate of all grid cell-level scores. Each LULC type was assigned a habitat quality index, and the formula is as follows [50]: where Q xj is the habitat quality in grid cell x, that is, in LULC type j; Q xj can never be >1; H j represents habitat suitability; k is the half-saturation constant and was set as 0.5; and D xj is the total threat level in grid cell x with LULC type j, and the related parameters refer to previous literature [29]. The formula is as follows [50]: where R is the number of threats, r is the threat, y indexes all grid cells on r's raster map, Y r is the set of grid cells on r's raster map, W r is the threat weight and indicates the relative destructiveness of a degradation source to all habitats, β x indicates the level of accessibility in grid cell x, where 1 indicates complete accessibility, and S jr is the sensitivity of habitat type j to threat r, where values closer to 1 indicate greater sensitivity. The effect of threat r that originates in grid cell y, r y , on the habitat in grid cell x is given by i rxy , and the related parameters refer to previous literature [29]. The formula is as follows [50]: where d xy is the linear distance between grid cells x and y, and d r max is the maximum effective distance of threat r's reach across space. If i rxy > 0, then grid cell x is in degradation source r y 's disturbance zone.
Based on historical land-use changes in Foshan, we selected cropland, built-up land, and unused land as degradation sources. The input parameters were determined from the literature [50,57,58]. The required inputs included the LULC map, threat factor layers (cropland, built-up land, and unused land), the accessibility to threats, the weight of the threat factor (Table 3), and the sensitivity of land types to each threat. These historical data refer to previous literature [29] to establish threat source parameters. We assumed that these empirical relationships were constant during long-term land-use changes. Habitat quality scores in the maps should be interpreted as relative scores of 0 to 1, with a perfectly suitable habitat scored as 1 and a non-habitat scored as 0.

Spatial Response of Ecosystem Services to Rapid Urbanization
We defined a spatial response index to characterize the impact of rapid urbanization on ecosystem services. The spatial response index is to measure the percentage change in ES due to the percentage change in LULC. The formula was based on the land-use dynamic index and elasticity for ES values as follows [6]: where R is the spatial response index of changes in ES in response to LULC, LCD represents the land change dynamics in a period, ∆LUT i is the converted area of LULC type i, LUT i is the area of LULC type i, ES j2018 represents the total amount of ES type j in 2018, and ES j2000 represents the total amount of ES type j in 2000. For the threshold identification and mapping of the spatial response index, we first normalized the results of all ES to compare the spatio-temporal heterogeneity of effects across ES. Second, we reclassified and recoded the value of each ES based on the natural breaking point of the data. For example, the corresponding codes for high, medium, and low carbon storage were 1, 2, and 3, respectively. The codes corresponding to the high, middle, and low values of agricultural production were 10, 20, and 30, respectively. The codes corresponding to the high, middle, and low values of habitat quality were 1000, 2000, and 3000, respectively. Third, we assumed that the weight of each ES was the same and spatially superimposed the results. Then, we reclassified the coded results again to obtain the total ES. Finally, to determine the threshold of spatial continuity, we converted the results into point data and performed a spatial kriging interpolation. The value of the spatial responsive index was derived from the results of spatially interpolated total ES, ranging from 1 to 3 (a higher value reflects greater spatial responsiveness on the grid).

Sensitivity Analysis
We used the Morris method [59] to test the sensitivity of the InVEST model to the variations in precipitation, actual evapotranspiration, volumetric plant available water content, vegetation rooting depth, rainfall erosivity, soil erodibility, attributes of threat sources, habitat suitability, and sensitivity. When other parameters remained unchanged, an original value of ± 10% was taken to disturb a certain parameter. Then, the changed parameters were input into the InVEST model and the results were compared. Sensitivity to agricultural production and carbon sequestration were based on empirical formulas with a fixed proportion. The model was run independently for each of these variations.

Urbanization Levels and LULC in Foshan
The different urbanization levels in Foshan are shown in Figure 5. with a fixed proportion. The model was run independently for each of these variations.

Urbanization Levels and LULC in Foshan
The different urbanization levels in Foshan are shown in Figure 5.   Cropland comprised the largest area in Foshan in 2000, accounting for 48.43% of the total area. In 2018, built-up land was the largest area, accounting for 35.99% of Foshan ( Figure 6). Urbanized land increased in area, whereas all other land-use types decreased in area. Regarding the transformation of land-use types, cropland and grassland decreased the most. The area of cropland, woodland, grassland, and water bodies that were transformed to built-up land was 523.53, 56.52, 175.14, and 93.33 km 2 , respectively. For different urbanization levels, cropland decreased by 66.97%, woodland by 53.34%, and grassland by 87.75%, and built-up land increased by 86.26% in the developed urban area. The developing urban area had the largest increase in built-up land, by 415.77%. Such a rapid expansion in built-up land was also observed in rural areas, with a 286.08% increase.

Agricultural Production
The changes in agricultural production services in Foshan from 2000 to 2018 are shown in Figure 7. The results show that the overall capacity of provisioning services dropped during this period. Areas with higher agricultural production were mainly distributed along the rivers and north-central rural areas. The gross output of agricultural production decreased by 1.20 million USD (based on the 2018 value of the USD) during this period ( Table 4). The decline in agricultural production was mainly concentrated in developing urban areas, reduced by 963.1 × 10 3 USD. However, the area with the largest relative decrease in agricultural production was located in developed urban areas, a decrease of 77%. The developed urban and rural areas reduced by 136.7 × 10 3 and 103.0 × 10 3 USD, respectively. Remote Sens. 2021, 13, x FOR PEER REVIEW 11 of 24

Agricultural Production
The changes in agricultural production services in Foshan from 2000 to 2018 are shown in Figure 7. The results show that the overall capacity of provisioning services dropped during this period. Areas with higher agricultural production were mainly distributed along the rivers and north-central rural areas. The gross output of agricultural production decreased by 1.20 million USD (based on the 2018 value of the USD) during this period ( Table 4). The decline in agricultural production was mainly concentrated in developing urban areas, reduced by 963.1 × 10 3 USD. However, the area with the largest relative decrease in agricultural production was located in developed urban areas, a decrease of 77%. The developed urban and rural areas reduced by 136.7 × 10 3 and 103.0 × 10 3 USD, respectively.

Carbon Sequestration
The changes in carbon sequestration in Foshan from 2000 to 2018 are shown in Figure 8. The changes in carbon sequestration were concentrated in developing urban areas. The areas of high carbon sequestration were mainly distributed in the southwestern area covered with woodland. Carbon sequestration increased by 84.7 × 10 3 MgC during 2000-2018 ( Table 4). The rural areas had the highest carbon sequestration (17,850.6 × 10 3 MgC in 2000; 17,860.9 × 10 3 MgC in 2018). Regarding the overall rate of change in carbon sequestration, the developing and rural areas increased by 0.55% and 0.06%, respectively. In the developed area, the carbon sequestration decreased by 0.16%. The area where carbon sequestration had not changed accounted for 65% of the total area. The area with increased carbon sequestration was mainly in the developing and northern developed areas, accounting for 22% of the total area. The area with decreased sequestration was in the southern developed urban area, accounting for 13% of the total area.

Water Yield
The changes in water yield in Foshan from 2000 to 2018 are shown in Figure 9. The results show that water yield increased in developed and developing urban areas. Higher values for water yield were observed in the eastern developed urban area. The water yield increased by 1504.7 × 10 3 m 3 in Foshan from 2000 to 2018 (Table 4). In 2000, the developed, developing, and rural areas showed water yields of 556.51 × 10 4 m 3 , 2408.49 × 10 4 m 3 , and 2621.68 × 10 4 m 3 , respectively. In 2018, the total water yield in the developed and developing urban areas increased by 4.96% and 3.04%, respectively. The changes in water yield showed that the area with increased yield accounted for 28% of the total area.

Water Yield
The changes in water yield in Foshan from 2000 to 2018 are shown in Figure 9. The results show that water yield increased in developed and developing urban areas. Higher values for water yield were observed in the eastern developed urban area. The water yield increased by 1504.7 × 10 3 m 3 in Foshan from 2000 to 2018 (Table 4). In 2000, the developed, developing, and rural areas showed water yields of 556.51 × 10 4 m 3 , 2408.49 × 10 4 m 3 , and 2621.68 × 10 4 m 3 , respectively. In 2018, the total water yield in the developed and developing urban areas increased by 4.96% and 3.04%, respectively. The changes in water yield showed that the area with increased yield accounted for 28% of the total area. values for water yield were observed in the eastern developed urban area. The water yield increased by 1504.7 × 10 3 m 3 in Foshan from 2000 to 2018 (Table 4). In 2000, the developed, developing, and rural areas showed water yields of 556.51 × 10 4 m 3 , 2408.49 × 10 4 m 3 , and 2621.68 × 10 4 m 3 , respectively. In 2018, the total water yield in the developed and developing urban areas increased by 4.96% and 3.04%, respectively. The changes in water yield showed that the area with increased yield accounted for 28% of the total area.

Soil Retention
The changes in soil retention in Foshan from 2000 to 2018 are shown in Figure 10. The results show a small change (0.15% decrease) in soil retention for the study area from 2000

Soil Retention
The changes in soil retention in Foshan from 2000 to 2018 are shown in Figure 10. The results show a small change (0.15% decrease) in soil retention for the study area from 2000 to 2018 (Table 4). The decreased area was mainly distributed in developed and developing urban areas. Areas with higher values of soil retention were toward the southwest of the rural areas. The total amount of soil retention in Foshan was 139.07 × 10 4 t in 2000 and 138.86 × 10 4 t in 2018. The rural areas accounted for the largest share of soil retention. The developed and developing urban areas showed a 4.01% and 3.64% decrease in soil retention, respectively. Areas with decreased soil retention accounted for 26% of the total area.

Habitat Quality
The changes in the habitat quality index in Foshan from 2000 to 2018 are shown in Figure 11. The results show a trend of degraded habitat quality, especially in developing and developed urban areas. Areas with higher habitat quality indices were mainly found in the southwest of the study area. The lower value areas were distributed to the east, gradually expanding to the northwest during the study period. The habitat quality index apparently declined by 17.79% in developed urban areas and 21.63% in developing urban areas. However, the habitat quality in rural areas increased by 0.14% (Table 4). Areas with decreased habitat quality accounted for 48% of the total area.

Habitat Quality
The changes in the habitat quality index in Foshan from 2000 to 2018 are shown in Figure 11. The results show a trend of degraded habitat quality, especially in developing and developed urban areas. Areas with higher habitat quality indices were mainly found in the southwest of the study area. The lower value areas were distributed to the east, gradually expanding to the northwest during the study period. The habitat quality index apparently declined by 17.79% in developed urban areas and 21.63% in developing urban areas. However, the habitat quality in rural areas increased by 0.14% (Table 4). Areas with decreased habitat quality accounted for 48% of the total area.

Spatial Response of ES Change with Respect to LULC
In Foshan, the spatial response of agricultural production (16.23%) was higher than that of other ES, followed by soil retention (5.75%), habitat quality (3.47%), water yield (2.38%), and carbon sequestration (0.05%).
In areas with different urbanization levels, the spatial response index differed for different ES (Figure 12). The sum of spatial response for the total ES in rural, developing, and developed areas was 14.15%, 22.50%, and 33.44%, respectively. The spatial response index of agricultural production was 4.05% in rural areas, 13.21% in developing urban areas, and 20.53% in developed urban areas. The spatial response of carbon sequestration in developing areas was higher than that in developed and rural areas. The spatial response of water yield changes was 4.24% in rural areas, 2.95% in developing urban areas, and 2.16% in developed urban areas. The highest soil retention spatial response was observed in rural areas (6.02%). The spatial response index of habitat quality was 0.09% in rural areas, 4.50% in developing urban areas, and 4.69% in developed urban areas.

Spatial Response of ES Change with Respect to LULC
In Foshan, the spatial response of agricultural production (16.23%) was higher than that of other ES, followed by soil retention (5.75%), habitat quality (3.47%), water yield (2.38%), and carbon sequestration (0.05%).
In areas with different urbanization levels, the spatial response index differed for different ES (Figure 12). The sum of spatial response for the total ES in rural, developing, and developed areas was 14.15%, 22.50%, and 33.44%, respectively. The spatial response index of agricultural production was 4.05% in rural areas, 13.21% in developing urban areas, and 20.53% in developed urban areas. The spatial response of carbon sequestration in developing areas was higher than that in developed and rural areas. The spatial response of water yield changes was 4.24% in rural areas, 2.95% in developing urban areas, and 2.16% in developed urban areas. The highest soil retention spatial response was observed in rural areas (6.02%). The spatial response index of habitat quality was 0.09% in rural areas, 4.50% in developing urban areas, and 4.69% in developed urban areas.
The spatial distribution of spatial response index derived from the total ES with respect to LULC is shown in Figure 13. The areas with higher values were more consistent with the distribution of water bodies. Four high-value accumulated areas were mainly distributed in the developing urban area. The spatial response index of the rural areas was lower. The average spatial response index in rural, developing, and developed urban areas was 0.39, 0.71, and 0.45, respectively. index of agricultural production was 4.05% in rural areas, 13.21% in developing urban areas, and 20.53% in developed urban areas. The spatial response of carbon sequestration in developing areas was higher than that in developed and rural areas. The spatial response of water yield changes was 4.24% in rural areas, 2.95% in developing urban areas, and 2.16% in developed urban areas. The highest soil retention spatial response was observed in rural areas (6.02%). The spatial response index of habitat quality was 0.09% in rural areas, 4.50% in developing urban areas, and 4.69% in developed urban areas. The spatial distribution of spatial response index derived from the total ES with respect to LULC is shown in Figure 13. The areas with higher values were more consistent with the distribution of water bodies. Four high-value accumulated areas were mainly distributed in the developing urban area. The spatial response index of the rural areas was lower. The average spatial response index in rural, developing, and developed urban areas was 0.39, 0.71, and 0.45, respectively.

Influence of Urbanization on ES in Foshan
Our results showed that rapid urbanization and expansion of built-up land have led to a significant reduction in agricultural production. This reduction in provisioning services was also found in another study that focused on the changes in ES in Foshan [60]. The main reason for the decrease in agricultural production was a severe labor shortage in rural communities [3]. Many people have migrated to urban areas to engage in industrial and commercial activities in search of higher wages. Labor shortages and industrialization have hastened the considerable changes in traditional agricultural practices. This

Influence of Urbanization on ES in Foshan
Our results showed that rapid urbanization and expansion of built-up land have led to a significant reduction in agricultural production. This reduction in provisioning services was also found in another study that focused on the changes in ES in Foshan [60]. The main reason for the decrease in agricultural production was a severe labor shortage in rural communities [3]. Many people have migrated to urban areas to engage in industrial and commercial activities in search of higher wages. Labor shortages and industrialization have hastened the considerable changes in traditional agricultural practices. This is a common problem in the Pearl River Delta region; decreases in agricultural production have also occurred in other cities in this region [61][62][63].
Interestingly, our results showed an increase in carbon sequestration, which may conflict with other published work [44]. This is because built-up land generally includes a certain proportion of artificial urban greening, which could have a higher carbon density than a semi-artificial cropland ecosystem [29,64]. In addition, a decrease in agricultural production indicates that the region needed to import resources from other regions, generating higher transport costs and transport-related emissions. This feedback of carbon emissions also brings more challenges to carbon sequestration.
The decrease in habitat quality was mainly due to the increase in built-up land and human activities. We treated the human-modified land types (e.g., cropland, built-up land, and unused land) as the sources of degradation. The concentration of population (based on nighttime light data) and the increase in built-up land (based on habitat accessibility and sensitivity) have reduced the habitat quality, especially in developed and developing urban areas. In addition, the increase in these threat sources could cause edge effects, with potential changes in the biological and physical conditions that occur at patch boundaries and within adjacent patches [65]. The decrease in habitat quality has also been reported by other studies. For example, Xiao et al. [66] proposed an ecological environment quality index based on non-construction land, and they also observed a loss of biodiversity conservation value in Foshan. Habitat loss has also been observed in surrounding cities with similar urbanization processes, such as Guangzhou [67], Shenzhen [68], and Dongguan [69].
The water yield increased slightly because several considerable measures have been undertaken to repair and restore abandoned fishponds and to build new fishponds, particularly since 2000 [67]. The increase in water yield under urbanization is consistent with that reported by other studies [70,71]. However, some studies reported that urbanization would lead to reduced water yield [72,73]. This conflict with our findings may be due to precipitation differences and reservoir rivers that may be affected by urbanization. The increase in impermeable pavements in urban built-up land may also contribute to evapotranspiration and surface runoff [74]. However, despite a slight increase in water yield, impermeable pavements are harmful to underground water sources, preventing aquifer recharge. Lowering the underground water table by preventing recharge may result in serious water shortages and disasters.
In terms of spatial variation, our results showed that the changes in ES were mainly concentrated in developed and developing urban areas. These results provide answers to our second research problem regarding the areas that had the greatest changes in ES under different urbanization levels. First, there was a reduction in agricultural production across Foshan. The developed and developing urban areas had a higher rate of change in agricultural production than rural areas (−77.98% and −66.33%, respectively) ( Table  4). The main reason for this result was the rapid development of real estate and the remarkable changes that have taken place in farmers' housing and lifestyle, especially in the Sanshui District. Second, the decrease in the habitat quality index (−21.63%) was mainly distributed in developing urban areas. This is due to the large increase in industrial land and population, leading to an increase in threat sources' peri-urban area. In addition, the variation in other ES was also concentrated in developed and developing urban areas, such as north of Nanhai District. Industrial development and agglomeration, such as from metal industries, furniture manufacturing, and auto manufacturing, are important drivers of changes to land use and ES in this region. The spatial heterogeneity of ES also highlighted the changes in urban-rural inequality [75].
Based on the statistics, the unprecedented rate of urbanization from 75.06% to 95.00% from 2000 to 2018 has had a considerable effect on Foshan. Nighttime light data confirmed the process of expanding urbanization in Foshan. As the second-most urbanized city in China, Foshan is making great efforts to achieve China's National New-type Urbanization Plan targets [4]. The plan covers almost every conceivable aspect of urbanization, such as sustainable development, institutional arrangements, and implementation. It sets key indicators, some with numerical goals, for urbanization level, public services, infrastructure, and resources and environment, as per the guiding principle emphasizing a sustainable and people-centered approach. Although Foshan's indicators of urbanization have far exceeded national targets, our results show that rapid urbanization has an inevitable effect on the natural environment and ES, and these changes will feed back into the social system and become a key constraint of economic development [76,77].

Implications for Urbanization Development and Sustainable Planning
With the implementation of the development plan for the Guangdong-Hong Kong-Macao Greater Bay Area as part of the national strategy in 2019, Foshan plays an important role in the development of advanced manufacturing. To consider the ecological effects of urbanization, the local government will implement innovative solutions that balance economic growth and sustainable development [67]. Several sustainable measures and plans have been implemented in recent years, including "High-quality forest city construction planning in Foshan", "Natural ecological civilization construction planning", and the ecological restoration project for rivers in Foshan.
Our results have implications for urban development and sustainable planning. We believe that sustainable planning should consider the living environment of human settlements and the effects of the ecological environment. First, our results showed that rapid urbanization greatly affects ES in Foshan. Agricultural provisioning service transitions directly affect resource consumption in the city and are related to regional food provisioning security [1,78]. We recommend that urban development should not encroach on cropland and that priority should be given to increasing the intensive use of built-up areas. Our results showed that decreases in habitat quality are mainly distributed in developed urban areas. The threat sources have increased and intensified, altering habitat accessibility and sensitivity [67]. Thus, it is necessary to plan and construct urban ecological corridors, such as channel and urban green space connections, to improve landscape connectivity.
Second, areas with different levels of urbanization should respond to important ES changes and carry out corresponding ecological restoration measures. Our spatial response results showed the responsiveness of ES to changes in land use, indicating that the developed and developing urban areas were the most important areas for ES changes during 2000-2018. Agricultural production and habitat quality, which are closely related to human survival, are more sensitive to urbanization. However, water yield services did not pose a threat in our study. Rural areas may soon enter a critical state regarding the loss of provisioning services due to increased industrial demand in Foshan.
Although urbanization and industrialization have greatly changed the social economy, natural ecosystems and their services are still the cornerstone of social development. While paying attention to the changes in industry, economy, population, and land use, we should also understand the changes in the ecological environment and the spatial differences in supply and services. ES is important for land-use planning and urban sustainable development, which can help prevent the loss of natural ES from offsetting the benefits of urbanization to the people [3]. For example, in the construction and planning processes for new cities, the amount of urban forest and green space should be increased and integrated into the architecture of the urban landscape. River dredging projects should be completely integrated with regional water resource management strategies to improve ecological connectivity. We consider that our spatial response map can directly help with practical urban planning. The identification of spatial thresholds can help locate key areas and implement land-use improvement and preventive measures. In the future, challenges for managing natural resources in urban areas may involve improving land-use efficiency within a limited area. The protection and improvement of non-construction land (e.g., cropland, woodland, grassland, and water bodies) are effective ways to ensure ecological security.
Moreover, we should not only consider the change in ES provision and distribution but also pay attention to the improvement and protection of ES quality. Urbanization and land-use planning should consider the conservation of high-priority habitats and high-productivity cropland. Though some indicators and numerical targets have been proposed in China's National New-type Urbanization Plan, such as the proportion of "green" buildings in new constructions in cities or meeting air-quality targets [4,79], these are not enough to ensure urbanization that is more comfortable, safer, and more harmonious with natural development. More ecological indicators and ES, such as air quality, recreation and tourism, and renewable energy, related to human well-being should be utilized in urban planning to reduce risk from a natural ecosystem perspective [80,81].

Limitations and Caveats
We used ArcGIS to combine the InVEST model with remote sensing data to quantify the effect of urbanization on ES. Some limitations and uncertainties exist in the physical assessment of ES. First, the model framework employs a relatively simple functional approach to assess ES. For example, the calculation of agricultural production by monetary value is a compromise between data source and data spatialization. More detailed agriculturalrelated spatial statistics, such as crop types and cultivation methods, are currently not available in this study. However, the monetary value of different agricultural industry statistics can realize the spatialization of agricultural production services. Although there are certain fluctuations in the prices of agricultural products, our results characterize the changes in production capacity and production levels over a period of time to a certain extent. The estimation of soil retention involved the RUSLE formula, which is based on a statistical relationship from many plot-scale experiments [82]. However, in practice, the slope has nonlinear effects on vegetation restoration and soil retention, which may overestimate the ES [83]. Meanwhile, the methods for assessing carbon sequestration and water yield were also the simplifications of complex ecological processes related to the carbon-water cycle [84]. The inter-annual variation of rainfall will also have an impact on regional water yield. Future research will consider collecting and averaging rainfall data over many years to reduce the impact of the volatility of data changes. Some changes in ecological processes do not consider the changes in species, age, and environmental factors over time [85].
Second, we only assessed provisioning, regulating, and supporting services and did not include cultural services [86]. This was because our goal was to identify the objective effects of different levels of urbanization on natural ecosystems, whereas subjective factors and a lack of material benefits may be considered when quantifying cultural services [87].
Finally, the spatial resolution contributed to the uncertainties in our findings. Land use has been identified as one of the most important drivers of the changes in ES provision [88][89][90]. We used LULC maps with a spatial resolution of 300 m in 2000 and 2018. Although higherprecision land-use data can improve the accuracy of ES assessment, the data we used were helpful for achieving the research objectives [91]. We paid more attention to the spatio-temporal variation rather than the improvement in absolute accuracy.
In addition, we revised the formula to determine the spatial response of ES changes with respect to LULC [6]. The spatial response indicator was revised to include a measure of the percentage change in ES due to a percentage change in LULC. Song and Deng [6] used the land-use dynamic index to calculate spatial response. However, this index represents the annual rate of land-use change, whereas the changes in ES were within a specific research period. Therefore, we removed the research period from the original LCD formula and calculated the land-use change dynamics within a specific period. We presume this revision better helped us to understand the relationship between ES changes and land-use changes. In summary, although there were some uncertainties in the simulation, we have certainty in our results. These methods are ultimately a trade-off between the complexity of ecosystem processes and the operability of model simulations.

Conclusions
In the present study, we proposed a research framework by integrating nighttime light data, the InVEST model, and spatial response. Based on multi-source data, we established the relationship between the nighttime light data index and the urbanization levels, and identified different urbanization areas. A spatial response index was defined, and its threshold was identified to characterize the response of ES to urbanization. We quantified and mapped the changes in land use along with five vital ES under three different urbanization levels in Foshan City between 2000 and 2018. Our results explain the effect of rapid urbanization on regional vital ES.
First, our results show that rapid urbanization greatly affected agricultural production and habitat quality services. Agricultural production reduced by 49.13% during 2000-2018 in Foshan, followed by a 10.13% reduction in habitat quality. Rapid urbanization had little effect on regulating services in which carbon sequestration and water yield had slightly improved, whereas soil retention had slightly decreased. Second, the changes in ES were mainly concentrated in developing urban areas. The spatial response index of agricultural production, habitat quality, soil retention, water yield, and carbon sequestration were 14.25%, 2.94%, 0.04%, 0.78%, and 0.07%, respectively. The average values of spatial response index in rural, developing, and developed urban areas were 39%, 71%, and 45%, respectively. We found that developing urban areas had the highest spatial response index, which indicates that the ES in this area are more sensitive to land-use change.
We consider that our research framework can help identify the key urbanization areas by identifying the threshold of nighttime light data. We recommend that areas with different levels of urbanization should respond to important ES changes and implement corresponding ecological conservation and restoration measures. In developed and developing urban areas, it is necessary to plan and construct urban ecological corridors, such as urban green space connections. In rural areas, cropland protection is important. Our results can further increase our understanding of the effect of urbanization on ES. Visualizing the spatial response and extracting the spatial threshold will be conducive to sustainable urban management and planning.
Author Contributions: Z.W. designed the framework and wrote the manuscript; R.Z. used the software and analyzed the data; Z.Z. carried out the validation and drew the maps. All authors have read and agreed to the published version of the manuscript.