Spatial–Temporal Impacts of Urban Sprawl on Ecosystem Services: Implications for Urban Planning in the Process of Rapid Urbanization

Rapid urbanization drives land cover change, affecting urban ecosystems and inducing serious environmental issues. The study region of Changchun, China was divided into three urbanization categories according to different urbanization levels and the characteristics of urban sprawl and changes and relationships between typical ecosystem services (ESs) under rapid urbanization were analysed. The results showed that Changchun has undergone considerable urban expansion since 2000, which has significantly impacted all ESs in terms of spatial and temporal heterogeneity. Habitat suitability and crop yield have relatively stronger service capacity in the study area. Since the expansion of large-scale infrastructures, the mean ES values of developed urban areas are the lowest among the three zones, except for water retention and sandstorm prevention in 2015, when the balance between all services decreased. Over the past 16 years, habitat suitability in developing urban areas has decreased to a large extent due to urban sprawl. Because of the improvement in agricultural science and technology, crop yield in three regions increased, while the area of cropland reduced from 1720 km2 to 1560 km2 (9.3%). Synergies between habitat suitability and carbon storage and habitat suitability and soil retention were detected in three areas. A trade-off between habitat suitability and water retention was detected in three areas. The interactions between crop yield and carbon storage, habitat suitability, and soil retention were more complex in this study region. In addition to water retention, urbanization index has a negative correlation with ESs. According to the results, some suggestions to alleviate ES loss during the process of rapid urbanization were proposed, which may guide scientific urban planning for sustainable urban development.


Introduction
Urbanization is a global universal process and an inevitable phase towards modernization for the development of most countries [1]. According to a report released by the United Nations, the rate of global urbanization increased from 30% to 55% between 1950 and 2018 and will continue to increase to 68% by 2050 [2]. In China, the urbanization rate, which was only 10.64% in 1949, exceeded 50% for the first time in 2011 and reached 59.58% in 2018 [3]. There are many factors influencing urban expansion, including rural land reform, large-scale industrialization, incentive policy, migration, economic reform, and market-led urbanization, etc [3]. With the rapid development of urbanization, rural populations moved into cities and secondary and tertiary industries continued to gather in cities. China's seventh national census shows that 902 million people live in the city in 2021, reaching 63.89% of the total population [4]. The population agglomeration in urban areas in particular has led to a surge in demand for urban land resources. Statistics about urban As a famous old industrial base in China, urbanization of Changchun was once higher than the national average level. With the emergence of abuses in the development of resource-based cities, the economic development of Changchun has stagnated. The proposal and implementation of the Strategy of Revitalizing Northeast China in 2003 propelled Changchun into an economic boom and rapid urbanization, which has inevitably accelerated land cover changes, soil resource loss, the heat island effect, and other environmental issues [21].

Data Source and Processing
In this study, we collected and used multi-source geographic data to quantitatively assess ES changes during urban sprawl. The main data sources and their descriptions are as the follows (Table 1): As a famous old industrial base in China, urbanization of Changchun was once higher than the national average level. With the emergence of abuses in the development of resource-based cities, the economic development of Changchun has stagnated. The proposal and implementation of the Strategy of Revitalizing Northeast China in 2003 propelled Changchun into an economic boom and rapid urbanization, which has inevitably accelerated land cover changes, soil resource loss, the heat island effect, and other environmental issues [21].

Data Source and Processing
In this study, we collected and used multi-source geographic data to quantitatively assess ES changes during urban sprawl. The main data sources and their descriptions are as the follows (Table 1) imagery data were downloaded from the website of NOAA/NGDC (https://ngdc.noaa. gov, accessed on 15 April 2021). Night-time light imagery data for different periods arrived from different sensors: data of the Operational Linescan System (OLS) arrived from the Defence Meteorological Satellite Program (DMSP) in 2000 and data from the Visible Infrared Imaging Radiometer (VIIRS) arrived from the NPOESS Preparatory Project (NPP) in 2015; therefore, resampling, mutual correction, and interannual fusion of the two periods of data were carried out according to the methods described in the literature to ensure data consistency [26,27]. Evapotranspiration products from a moderate resolution imaging spectroradiometer (MODIS) with a 500 m resolution were accessed from the NASA website and were used to calculate the water retention services. Land cover datasets in 2000 and 2015, interpreted from Landsat Thematic Mapper (TM) and Operational Land Imager (OLI) data were provided by the Northeast Institute of Geography and Agroecology, Chinese Academy of Sciences (IGA, CAS). Land cover types were divided into seven categories according to Wu [28]: cropland, grassland, forestland, wetland, built-up land, water body, and others. Based on the ground investigation data, the overall accuracy of the land cover classification data was higher than 90% for the two periods. DMSP/OLS or NPP/VIIRS sensors can detect the low-intensity light emitted by urban lights and even small-scale residential areas and traffic flow at night, which make urban areas different from dark rural areas. Therefore, night-time light imagery data, as an intuitive representation of human activities, have been widely used in modelling urban sprawl [26]. Generally, higher Digital Number (DN) values indicate higher levels of urbanization. According to previous studies, a DN value of 50 was used as a threshold for distinguishing urbanization areas from rural areas. Accordingly, the study area was divided into three areas according to different levels of urbanization [1]: developed urban, developing urban, and rural area. For developed urban area, the DN values of the whole study region were both greater than 50 in 2000 and 2015; for developing urban area, the DN value in 2000 was less than 50, but greater than 50 in 2015; for rural area, the DN values were both less than 50 in 2000 and 2015.

Estimation of ESs
In this study, six vital ESs, including one provision service (crop yield), four regulation services (carbon storage, sandstorm prevention, soil retention, and water retention), and Land 2021, 10, 1210 5 of 17 one supporting service (habitat suitability) were selected because they are closely related to the dwellers' health and security, and are sensitive to land cover changes.
The approaches used to quantify each service are briefly described in the following subsections:

Carbon Storage Service
Carbon storage is an important indicator of climate regulation [29]. The InVEST model was used to evaluate carbon stored in carbon pools, including above ground biomass, below ground biomass, soil, and dead organic matter [30]. The total carbon storage is the sum of those four carbon sources: where C i represents the total carbon storage of grid cell i.
Carbon storage was calculated based on the specific ecosystem classification [28]. Because the InVEST model was designed for worldwide application, some functions are general without considering regional heterogeneity [30]. For better results, the biomass carbon stocks per unit area of each land cover type were derived from the local research results and long-term observation data from field stations in northeast China [31].

Soil Retention Service
Soil retention is an important regulating ES, especially in the case of rapid urban sprawl and significant increases in impervious surfaces in the Black Soil Belt. With frequent human interference, regional soil erosion is becoming more and more serious, and the soil retention capacity of the study area has changed greatly. The Revised Universal Soil Loss Equation [32] was used to estimate the soil retention service of different ecosystems. Areas of land with high soil retention present less risk of erosion loss [33]. Soil retention is the difference between potential soil erosion and actual soil erosion; the formula is as follows [34]: where A is the annual soil loss (t km −2 a −1 ); A p is the potential annual soil loss (t km −2 a −1 ); A a is the actual annual soil loss (t km −2 a −1 ); R indicates the erosion factor of rainfall runoff (MJ mm km −2 ha −1 a −1 ), which is calculated according to the daily rainfall; K refers to the soil erosion factor (t h JM −1 mm −1 ); LS denotes the slope length and gradient, representing the terrain factor; C indicates the dimensionless vegetation cover factor, and P represents the soil retention measure. All factors were spatialized into a grid with a 30 m resolution in ArcGIS 10.3 software and grid calculations were performed to obtain the annual soil loss for each unit.

Water Retention Service
Water retention is one of the most critical services of the ecosystem, which closely connects the growing population with the ecosystem [34]. The amount of water retention for each pixel on the landscape was calculated using the water balance equation [35]. The model is expressed as follows [35]: where WR is the total amount of water retention, P i means annual precipitation (mm), R i indicates the surface runoff (mm), AET i represents annual actual evapotranspiration (mm), A i is area (km 2 ) of the ecosystems defined by land cover types, α is average surface runoff coefficient collected from the published research results in this area [36], i is the index of ecosystem types, and j denotes the total number of ecosystem types. Precipitation data were calculated using meteorological data from CMDSSS (http://data.cma.cn, accessed on 15 April 2021). The evapotranspiration data were calculated using MODIS products and meteorological records [35].

Habitat Suitability
Habitat suitability plays a crucial role in biological diffusion and biodiversity conservation, representing the ecosystems' ability to provide suitable conditions for the sustainable survival of individuals and populations [37]. The habitat suitability index was used to estimate the habitat and vegetation supporting ecosystems across the whole landscape, which considered the relative quality of the ecosystem type, the availability of the water source, the sensitivity to disturbing factors, and the abundance of food. The equation was as follows [38]: where HSI represents the total value of the comprehensive index of habitat suitability in the study year, w i is the weight for the ith factor assigned by experts, f i is the standardized value of the ith factor involved in the assessment of habitat suitability, and i is the number of evaluated factors. Spatialization and weighted summation for all factors were performed utilizing Ar-cGIS 10.3 software.

Sandstorm Prevention
Sandstorm prevention of the ecosystem was calculated using the revised wind erosion equation (RWEQ) [39]. The RWEQ is an empirical model for estimating long-term soil loss caused by wind erosion [34]. In this model, climate, topography, vegetation, and soil properties were considered to quantify wind erosion [40]. Potential and actual wind erosion intensity were estimated, and the difference between them was used as the sand fixation capacity of the ecosystem to assess the function of wind and sand fixation in the ecosystem. The formulas involved are as follows [39]: where ∆Q means the amount of sand fixation (t km −2 a −1 ), Q 0 refers to the potential soil erosion without vegetation cover (t km −2 a −1 ), and Q v indicates the actual soil erosion considering the current land cover and management conditions (t km −2 a −1 ). Q(x) represents the amount of soil transported by wind at point x, Q max indicates the maximum sand transport capacity by wind, and S describes critical field length. WF, EF, SCF, K , and C represent weather factor, soil erodible fraction, soil crust factor, surface roughness factor, and vegetation factor, respectively.

Crop Yield
Located in the typical Black Soil Belt of Northeast China, the food supply of the study area is a crucial service for regional food security. Crop yield was selected to estimate food provisioning services of Changchun. Crop yields data were obtained from the statistical yearbook of Changchun in 2001 and 2016.

Relationship between ESs
According to previous studies on the relationship between ESs [1,37,41], the changes of ES in divided urban area were analysed and compared to reveal the influences of urbanization on the relationships between ES changes.
There are three relationship types between ESs: synergy, trade-off, and neutral. Synergy represents a win-win or lose-lose situation of two ESs with a consistent trend of increases or losses; trade-off indicates a win-lose or lose-win situation in which one ES rises at the expense of the other; and neutral means that a change in one ES will not cause variation in another [1]. In this study, 754 random points (432 in rural area, 209 in developing area, and 113 in developed areas) for each ES layer, both in 2000 and 2015, were extracted using ArcGIS 10.3 software. Before comparing the 2000 and 2015 ecosystem services, all data were standardized using the following equations:  (Table 2).    The spatial distribution of carbon storage, soil retention, and habitat suitability services presented similar patterns with high values in eastern areas, where the slope gradient was relatively high, and the dominant land cover type was forestland ( Figure 3). The values range from 0 to 37.1 t, 0 to 4.2 × 10 3 t km −2 a −1 , and 0 to 60.2 for these three ecosystem services, respectively. The sandstorm prevention deceased from west to east, and the areas with high values of water retention were concentrated in the central parts ( Figure 3). From 2000 to 2015, the three ESs-carbon storage, soil retention, and water retention-increased by 2.01%, 124%, and 28.03% in the study area, while habitat suitability and sandstorm prevention services decreased by 3.55%, and 1.64%, respectively ( Figure 4). In the study area, the area of cropland was reduced from 1720.4 km 2 to 1560.4 km 2 by 9.3%, but the crop yield service increased from 668.5 × 10 6 t to 753.4 × 10 6 t by 12.87% ( Figure 4).  The spatial distribution of carbon storage, soil retention, and habitat suitability services presented similar patterns with high values in eastern areas, where the slope gradient was relatively high, and the dominant land cover type was forestland ( Figure 3). The values range from 0 to 37.1 t, 0 to 4.2 × 10 3 t km −2 a −1 , and 0 to 60.2 for these three ecosystem services, respectively. The sandstorm prevention deceased from west to east, and the areas with high values of water retention were concentrated in the central parts ( Figure 3). From 2000 to 2015, the three ESs-carbon storage, soil retention, and water retention-increased by 2.01%, 124%, and 28.03% in the study area, while habitat suitability and sandstorm prevention services decreased by 3.55%, and 1.64%, respectively ( Figure 4). In the study area, the area of cropland was reduced from 1720.4 km 2 to 1560.4 km 2 by 9.3%, but the crop yield service increased from 668.5 × 10 6 t to 753.4 × 10 6 t by 12.87% ( Figure 4).    The relative service capacities are shown using radar plots based on standardized values ( Figure 5). Each axis of the plot indicates a different ES; the outermost point on the axes means the highest level of ES, and the provisioning service decreases towards the centre. The symmetry of the plot represents relative balance of the ESs; consequently, the larger and more symmetrical they are, the higher the overall potential benefits of ESs [1]. It can be seen that six service capacities in the study area are extremely unbalanced. Habitat suitability and crop yield are relatively stronger than the others, while the service capacity of soil retention is the weakest.  The relative service capacities are shown using radar plots based on standardized values ( Figure 5). Each axis of the plot indicates a different ES; the outermost point on the axes means the highest level of ES, and the provisioning service decreases towards the centre. The symmetry of the plot represents relative balance of the ESs; consequently, the larger and more symmetrical they are, the higher the overall potential benefits of ESs [1]. It can be seen that six service capacities in the study area are extremely unbalanced. Habitat suitability and crop yield are relatively stronger than the others, while the service capacity of soil retention is the weakest.
centre. The symmetry of the plot represents relative balance of the ESs; consequen larger and more symmetrical they are, the higher the overall potential benefits of It can be seen that six service capacities in the study area are extremely unbalance itat suitability and crop yield are relatively stronger than the others, while the ser pacity of soil retention is the weakest.

Relationships among ES Changes
Interactions between ESs were analysed to evaluate the differential changes in the process of urbanization and to detect the urbanization problems presented in chun (Table 3). Habitat suitability significantly corresponded with carbon storage, tention, and water retention. Synergies between habitat suitability and carbon s and habitat suitability and soil retention were detected in the three areas, that is,

Relationships among ES Changes
Interactions between ESs were analysed to evaluate the differential changes in ESs in the process of urbanization and to detect the urbanization problems presented in Changchun (Table 3). Habitat suitability significantly corresponded with carbon storage, soil retention, and water retention. Synergies between habitat suitability and carbon storage, and habitat suitability and soil retention were detected in the three areas, that is, habitat suitability deteriorated with the decrease in carbon storage and soil retention. However, the correlation coefficient decreased from developed urban area to rural area in terms of carbon storage and increased from a developed urban area to a rural area for soil retention. A trade-off between habitat suitability and water retention was detected in the three zones, and the correlation coefficient of the rural area was the highest. The interactions between crop yield and carbon storage, habitat suitability, and soil retention were more complex. There was only a weak trade-off between crop yield and carbon storage in rural areas, while the correlation between the other two areas was not significant. For crop yield and habitat suitability, a higher trade-off in developed urban areas but a weaker trade-off in rural areas can be seen.
The correlation coefficients of sandstorm prevention with all of the other five ESs were insignificant in three areas. Insignificant results were also detected between water retention and soil retention, and water retention and crop yield.

ES Changes in Different Urbanization Areas
Average values of ESs in the three urban zones were calculated to compare the impact of rapid urban sprawl on ESs (Table 4). It can be seen that the mean values of soil retention, carbon density, habitat suitability, and crop yield in rural areas are significantly higher than those in the other two areas, both in 2000 and 2015. All of the mean ES values for urban areas are the lowest among the three areas, except for water retention and sandstorm prevention in 2015. For developed and developing urban areas, carbon density and habitat suitability decreased, while the other four ecological services improved slightly from 2000 to 2015. For rural area, only sandstorm prevention decreased from 6.0 × 10 3 t·km −2 to 4.44 × 10 3 t·km −2 during the study period.
We can also learn from radar maps that the relative balance of all of the services in three areas is quite different. For developed urban area (Figure 6a), habitat suitability and sandstorm prevention are stronger than other services and the balance between all services decreased from 2000 to 2015. For developing urban area (Figure 6b), habitat suitability is significantly higher than other ESs in 2000 but decreased to a large extent in 2015. For rural area (Figure 6c), crop yield increased significantly from 2000 to 2015, while other ESs hardly changed. The relative balance of all of the services in developing urban area is the most symmetrical among the three areas, but the overall potential ES benefit is not higher (Figure 6d).  We can also learn from radar maps that the relative balance of all of the services in three areas is quite different. For developed urban area (Figure 6a), habitat suitability and sandstorm prevention are stronger than other services and the balance between all services decreased from 2000 to 2015. For developing urban area (Figure 6b), habitat suitability is significantly higher than other ESs in 2000 but decreased to a large extent in 2015. For rural area (Figure 6c), crop yield increased significantly from 2000 to 2015, while other ESs hardly changed. The relative balance of all of the services in developing urban area is the most symmetrical among the three areas, but the overall potential ES benefit is not higher (Figure 6d).

Correlation Analysis between Urbanization and ESs
A total of 1500 sample points were randomly extracted, including 120 in developed urban areas, 221 in developing urban areas, and 1159 in rural areas. Taking nightlight

Correlation Analysis between Urbanization and ESs
A total of 1500 sample points were randomly extracted, including 120 in developed urban areas, 221 in developing urban areas, and 1159 in rural areas. Taking nightlight satellite measurements as a proxy for urbanization, correlation analysis were analysed between urbanization and ecosystem services in 2000 and 2015, respectively. Habitat suitability were negatively correlated with urbanization in all three areas ( Table 5). The relationship between water retention and urbanization was negative in developed urban areas, while it was positive in developing urban areas and rural areas. The relationship between crop yield and urbanization was negative in developed and developing urban areas, while it was insignificant in rural areas. Due to the extensiveness and continued growth of water-impervious surfaces at the loss of forest and crop land in developing urban areas, urbanization is negatively correlated with carbon density in 2015. There was no significant correlation between soil conservation, sandstorm prevention and urbanization in all areas.

Typical Impact of Urban Expansion on ESs in Changchun
Many studies have shown that changes in land cover driven by urbanization contribute greatly to spatial and temporal changes in ESs [12,33,41,42]. The influences brought about by urbanization not only reflected the urban-rural gradient changes, but also presented a spatial spill over effect [36]. In this study, we also found that rapid and extensive changes in land cover from rural to urban area in Changchun led to a remarkable influence on the spatial distribution of ESs. Some ESs, including carbon storage, soil retention, habitat suitability, and crop yield, increased with the distance from developed urban areas, which was consistent with the previous research results with regard to gradient change along the urban-rural gradient [1]. We found that water retention represented the opposite trend: the mean value in developed urban area was much higher than that in rural area. This is mainly because the evapotranspiration of impervious surface is small, which leads to larger water retention values.
Due to different environmental and land cover patterns in different areas, the dynamic trends of ESs are quite different. According to previous studies, an increased rate of water-impervious surfaces in developing urban areas is the highest compared to the other two types of areas; consequently, the greatest reduction in regulation services as well as food production services occurs in developing urban areas [1,38]. In our study, the mean value of carbon storage and habitat suitability declined, while other services increased in developing urban areas.
Crop yield increased in all three areas, which is consistent with the results obtained by Sun et al. in the Atlanta metropolitan area [16]. Although urbanization took up a large amount of cropland, especially in developing urban areas, the progress in agricultural science and technology and the improvement of farming management greatly increased crop yield per unit area, which made the overall level and average level of the crop yield in the three areas improve. In addition, the conversion of some dry farmland into paddy fields in rural areas also promoted the largest increase in crop yield among three areas. However, in the study of Qi et al., it was found that the food supply in this area increased first and then decreased [43]. Cropland protection cannot be ignored.
Urban sprawl declined the area of cropland, as well as grassland. Accordingly, the average amount of carbon storage in developing urban areas and developed urban areas reduced, while that in rural areas increased due to the implementation of the Grain to Green Program (the area of forest land increased by 13.90%, and water body increased by 128.95%), which offset the decrease in carbon storage due to the increase in impervious surface (Figure 7). It is worth mentioning that because of burgeoning spiritual demands, the three regions pay more attention to the environment and aesthetics. As a result, there was an increase in green infrastructure from 2000 to 2015, and the forestland area in the three regions increased by 13.90%, 88.9%, and 29.1%, compared with 2000, respectively. However, due to the rapid growth of artificial surface (the rate of increase was 26.55%, 91.9%, and 10.43% for rural area, developing urban area, and developed urban area, respectively, from 2000-2015), developing urban area showed the largest decline in their habitat suitability index from 40.38 to 36.17 (Figure 7). In rural area, though rural residential land expanded significantly due to urban sprawl, habitat suitability still improved a small amount because forestland and water body provide better habitat suitability services than cropland. Both water retention and soil retention services presented an inclining trend in all three areas and the largest increase in water retention occurred in developed urban area because of these areas saw the largest increase in forestland and bodies of water. Furthermore, the reduced evapotranspiration due to the large increase in impervious areas also resulted in the largest increase in water retention occurring in developed urban area. Since the increasing forestland as well as the high-density buildings reduced the wind speed and surface sand, the largest increase in sandstorm prevention services occurred in developed urban areas. Green Program (the area of forest land increased by 13.90%, and water body increased by 128.95%), which offset the decrease in carbon storage due to the increase in impervious surface ( Figure 7). It is worth mentioning that because of burgeoning spiritual demands, the three regions pay more attention to the environment and aesthetics. As a result, there was an increase in green infrastructure from 2000 to 2015, and the forestland area in the three regions increased by 13.90%, 88.9%, and 29.1%, compared with 2000, respectively. However, due to the rapid growth of artificial surface (the rate of increase was 26.55%, 91.9%, and 10.43% for rural area, developing urban area, and developed urban area, respectively, from 2000-2015), developing urban area showed the largest decline in their habitat suitability index from 40.38 to 36.17 (Figure 7). In rural area, though rural residential land expanded significantly due to urban sprawl, habitat suitability still improved a small amount because forestland and water body provide better habitat suitability services than cropland. Both water retention and soil retention services presented an inclining trend in all three areas and the largest increase in water retention occurred in developed urban area because of these areas saw the largest increase in forestland and bodies of water. Furthermore, the reduced evapotranspiration due to the large increase in impervious areas also resulted in the largest increase in water retention occurring in developed urban area. Since the increasing forestland as well as the high-density buildings reduced the wind speed and surface sand, the largest increase in sandstorm prevention services occurred in developed urban areas. Furthermore, correlation analyses showed that the urbanization index and habitant suitability, and crop yield had a negative correlation in Changchun, which are consistent with previous studies [9,41] and indicates that urbanization exerts significant influence on ecosystem services. In our study, water retention is significantly negatively correlated in Furthermore, correlation analyses showed that the urbanization index and habitant suitability, and crop yield had a negative correlation in Changchun, which are consistent with previous studies [9,41] and indicates that urbanization exerts significant influence on ecosystem services. In our study, water retention is significantly negatively correlated in urbanized areas, but positively correlated in developing urban areas and rural areas. However, the study of sun et al. in Guangdong-Hong Kong-Macao "Megacity" found that water retention was negatively correlated only in urban areas, but insignificant in rural areas [9]. This is because in addition to land use, other factors such as topography and meteorology, likely contributed to the relationship between urbanization and service capacity. In developed urban area, the terrain is so flat and high-density artificial surface caused large runoff and consequent reduction of water retention. While in developing urban area and rural area, water retention is more impacted by slope gradient and heterogeneity of evapotranspiration and precipitation. Carbon density only significantly negatively correlated with urbanization in developing areas, but not significantly in developed area and rural area. It shows that the large increase of urban land has a significant impact on carbon storage in the developing area than others.

Effects of Urbanization on the Relationships among ESs
Correlation analyses of multiple ESs could provide instructive information for the balanced development of ESs overall. The factors leading to changes in ES relationships varied spatiotemporally due to the different characteristics of the natural environment and the degree of human disturbance. This study mainly focused on the changes in ES relationships caused by land cover transformation.
Studies have shown that there are general trade-off relationships between provision services and regulation services, such as crop yield and carbon storage, water retention, and soil retention [44][45][46]. However, in our study, crop yield presented a more complex relationship with other ESs in different urbanization areas. Due to the large proportion of cropland (the proportion was 67.01% in 2000 and 53.08% in 2015), cropland plays a vital role as a carbon sink in rural area. Although there was an increase in forestland areas, carbon storage had barely changed due to the decrease in large area of cropland. However, impacted by advances in agricultural science and technology, and agricultural management factors, crop yield showed rapidly increasing trend, and became the most prominent ES in rural area (Figure 6), which caused a weak trade-off between crop yield and carbon storage in rural area. Large-scale expansion of urban areas leads to loss of large areas of suitable habitats, resulting in the largest reduction in habitat suitability in developing urban areas in this study. These two diverse trends resulted in a high trade-off between crop yield and habitat suitability in this region.
In addition, we found that habitat suitability significantly corresponded with carbon storage, soil retention, and water retention. In contrast, land cover changes driven by urbanization transformed most land surfaces into water-impervious surfaces, which raised water retention but reduced the soil retention and carbon storage. Therefore, habitat suitability has a trade-off relationship with water retention, but a synergistic relationship with soil retention and carbon storage.

Implications for Urban Planning
Due to rapid urban sprawl, coordinating development activities without destroying ecological sustainability has become a preliminary and vital issue in relation to urban management. Many studies have tried to solve this situation through delimiting ecological protection areas with high ESs [47]. However, the benefit of urbanization often promotes development activities in ecological protection zones, which lead to a failure to protect of these areas. The "Guiding Opinions on the Overall Delimitation and Implementation of Three Control Lines in National Territory Spatial Planning" issued by the Office of the State Council of China in 2019 emphasized the importance of the demarcation of the ecological protection red line, permanent basic cropland, and the urban development boundary [48]. Therefore, the contradiction between urban development and ecological protection can be mitigated by strict constraints on development activities in delineated areas. Taking ESs and urbanization as a whole system and analysing their interactions at local scale can provide more constructive suggestions for ecological protection in the process of urbanization. According to urbanization and its effect on ESs, our detailed recommendations for local urban development planning are as follows: First, Changchun, as the central city of the Harbin-Changchun Urban Agglomeration, is also the region where black soil is concentrated. The trade-off between food security and economic development is always the core problem of regional development. From 2000 to 2015, due to urbanization, 156.91 km 2 of cropland was transformed into urban built-up land, accounting for 8.2% of the total amount of cropland in 2000. Despite the improvement in agricultural science, technology, and management, crop yield showed an increasing trend in each of the three areas (Table 3). However, this benefit can only be maintained in a limited range. In the long run, the protection of cropland area is necessary. Therefore, urban sprawl should give full consideration to the red line of basic agriculture protection, as well as the sustainable use of black soil, so as to achieve a win-win situation of food security and economic development.
Second, we found that the increase in soil conservation services in the three areas and carbon storage in rural area due to Changchun's rapid urbanization benefited from the increase in forestland and water body with high ESs. Therefore, in urban landscape plan, we should pay more attention to the construction of green infrastructure, especially in developing urban areas where the area of water-impervious surfaces increases extensively, so as to reverse the phenomenon of the continuous loss of ESs in the process of urban sprawl.
Finally, it has been proven that the allocations and patterns of land cover have a significant impact on ESs. However, commonly, human activities changed land use types and tend to generate the most profit [41]. Therefore, it is critical that policy makers and planners should understand how to design the land use structure in order to promote sustainable urban development in a more reasonable and beneficial way.

Conclusions
Taking Changchun as study area, we studied the impact of rapid urbanization on changes in ESs in inland agricultural areas. We assessed six ESs and their interrelation in different urbanization areas. The urbanization changes and their impact on these ESs were analysed by integrating remote sensing and geospatial analysis methods.
Changchun experienced rapid urbanization between 2000 and 2015. At the same time, the ESs of Changchun changed significantly over the past 16 years, with increases in carbon storage, soil retention, water retention, and crop yield, and decreases in habitat suitability and sandstorm prevention service. Urbanization has a significant impact on ESs, and ESs change regularly in different urbanization areas. The relationships between ESs were not only simple trade-offs or synergies, but presented more complex issues, and the correlation characteristics were different in different urbanization areas. Habitant suitability and crop yield have a significantly negative relationship with urbanization. The increase in forestland and water with high ESs slowed down the decrease in ESs due to urbanization. In view of the impact of urbanization on ecosystem services in Changchun, we put forward three suggestions to avoid or minimize the ecological problems caused by urbanization. In this study, only the impact of land cover changes on ecosystem services was considered. The comprehensive impact of climate, terrain, and other factors on ESs should be considered in future work.  Data Availability Statement: The data are not publicly available for privacy reasons.

Conflicts of Interest:
The authors declare no conflict of interest.