Spatiotemporal Changes in Ecosystem Services along a Urban ‐ Rural ‐ Natural Gradient ： A Case Study of Xi ʹ an, China

: Urban areas are the areas that are most strongly affected by human activities, which presents many challenges to the ecosystem and human well ‐ being. Ecosystem services (ES) are a comprehensive indicator to measure the ecological effects of urbanization. To effectively identify and evaluate the impact of urbanization on ES, the spatial ‐ temporal pattern of ES should be considered. According to the level of urbanization, Xi ʹ an city is divided into four regions: the urban core area, the urban extended area, the rural area, and the ecological conservation area, then, five comprehensive ecosystem services (CES) are evaluated by In VEST model. The results showed the following: (1) There is an obvious spatial heterogeneity in the distribution of ES. The ecological conservation area is the hot spot of ES supply, and the low value is mostly distributed in the urban core area. (2) The CES in the urban extended area that has undergone the greatest change between 2000 and 2015, and the rates of change in the ecological conservation area are the smallest. (3) There is a significant correlation between urbanization and ES, and the correction between landscape urbanization and ES is the most significant. (4) The agglomeration relationship between urbanization and ES in different regions is not consistent. Regional division provides a new way to understand the interaction between urbanization and ES in time and space, so as to provide better guidance for policy makers in formulating sustainable development policies to alleviate the loss of ES caused by the process of urbanization.


Introduction
Ecosystem services (ES) refer to the benefits derived directly or indirectly by human beings from the ecosystem [1]. Human life and economic development depend on the goods and services provided by the ecosystem. ES are important indicators linking human well-being and the ecosystem [2,3]. There have been many studies on ES, especially the assessment and quantification of urban ecosystem services (UES), which have attracted increasing interest in recent years [4][5][6][7][8][9]. The difference between urban ecosystems and other ecosystems is that urban areas are human-dominated complex ecosystems [10,11]. Over the past few decades, many ES have degraded [12]. As people's dependence on ecosystem is becoming less obvious, the importance of ES is more likely to be ignored [3,13]. Therefore, it is necessary not only to understand how to maintain and improve ES, but also to integrate this knowledge into planning and decision making [10,14]. of 10,752 km 2 . The urbanization of Xi'an has experienced rapid development in recent years; GDP has increased from 73.385 billion CNY to 834.986 billion CNY between 2000 and 2018 [18]. With the rapid increase of the urban population, the urban area also experienced rapid and sustained expansion. Between 2000 and 2017, the built-up area increased from 477.64 km 2 to 1287.05 km 2 [40]. The sustainable development of cities results in a higher demand for the improvement of ES, at the same time, it also brings more pressure to the ecosystem. Therefore, it is necessary to study the impact of urbanization on ES. Xi'an was chosen for this study because it fully reflects the different levels of urbanization development, which is typical and representative.

Data Source and Processing
The 2000, 2005, 2010, and 2015 national land use/cover datasets (NLCD) were obtained from the Resource and Environment Data Cloud Platform with a 30 × 30 m resolution (http://www.resdc.cn/, accessed 21 November 2019). The impervious surfaces data are derived from the urban impervious surface datasets of Chinese cities established by Gong Peng et al. [41] based on Landsat imagery and nighttime light (NTL) data, which can reflect the urbanization process timely and effectively. The precipitation, evapotranspiration, and other meteorological data were obtained from the Chinese meteorological data network (http://data.cma.cn/). The spatial meteorological data were interpolated by ArcGIS 10.2. Soil datasets, including soil depth, root depth, and soil texture, soil saturated hydraulic conductivity (PAWC), derived from Shaanxi Soil and related literature [42][43][44]. Statistics on grain output were obtained from the Xi'an Statistical Yearbook (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). The digital elevation model (DEM) with a spatial resolution of 30 m was obtained from the Geospatial Data Cloud, and the watershed and sub-watershed layers were generated based on the DEM data. The spatial distribution of population and economic data came from the resource environment data cloud platform. The spatial distribution of population and GDP with a spatial resolution of 1 × 1 km was supplied by the Resource and Environmental Science Data Center of the Chinese Academy of Sciences (http://www.resdc.cn).

Classification of Urbanization
There are diverse definitions and metrics for urbanization, such as population size, impervious surfaces area, building structure, and economic level [45]. High population density and extensive impervious surfaces are the two most significant factors defining urban characteristics [45,46]. The previous method of embodying urbanization mainly relied on the statistical data of administrative units, but there was no accurate method for spatializing urbanization. The impervious surfaces integrated Landsat imagery and nighttime lighting data, which is effective for spatial processes that reflect urban sprawl [47,48], so impervious surface data was used to determine the urban core area.
We divided Xi'an into four regions: the urban core area with a high development level (between 2000 and 2015), the urban extended area (between 2000 and 2015), the rural area (between 2000 and 2015), and the ecological conservation area (Figure 2). The level of urbanization in different areas of Xi'an is significantly different. The urban core area which impervious surfaces is widely distributed and concentrated has high economic and urbanization levels. Although the urban extended area has a low degree of development, the urbanization growth rate is relatively high. The rural area has a lower level of urbanization and population density. The southern part is an ecological conservation area and is an important ecological security barrier in China. It has always been the area that is least affected by human activities.

Quantifying ES
The northern part of Xi'an is the Loess Plateau, and the southern part is the Qinling Mountains. The rich geographical conditions have caused its diverse ecosystem types. Soil erosion and water pollution are important environmental problems in Xi'an, so it is of great significance to evaluate the function of soil conservation and water purification of the Xi'an ecosystem [49]. Food supply is the most basic ES. The vast farmland provides an important food supply service for the development of Xi'an. Water and the atmosphere are the most important material foundations for human survival and development, and also the basic conditions for sustainable development. Therefore, this study selected five important ES in Xi'an for assessment, which included water yield (WY), soil conservation (SC), water purification (WP), carbon sequestration (CS), and food supply (FS), and they were representative to reflect the ecosystem of Xi'an. First, the ecosystem service assessment and trade-off model (InVEST) was used to assess and quantify WY, SC, WP, and CS [50], and food calories were used to quantify FS [51]. Finally, comprehensive ecosystem services (CES) were used to compare the level of ES in different regions.

Water Yield
Ecosystems provide a number of vital services for humans, and the InVEST model can be used to map and evaluate goods and services from nature. WY is an important ES. Urbanization has altered the composition of the underlying surface. The impervious surfaces have replaced the original land use type, which has changed the surface water circulation process and also affected the supply of water. The WY of the InVEST model is based on the water balance formula. Subtracting the actual evapotranspiration from precipitation is the annual WY of pixel. The formula is as follows: where AET X represents the annual actual evapotranspiration of pixel X, and P X represents the annual precipitation of pixel X.

Water Purification
Water purification is an important ecosystems service [52]. The InVEST model estimates WP by calculating the retention of nutrients (nitrogen and phosphorus) by vegetation and soil. Higher nitrogen and phosphorus exports reflect lower WP. Calculated as follows: where ALV X refers to the load value adjusted by pixel X, HSS X is the hydrological sensitivity value of pixel X, and pol X represents the output coefficient of pixel X.

Soil Conservation
Soil conservation is an important measure to protect resources and improve the ecological environment. SC is related to climate, soil, topography, and vegetation. The SC module in the InVEST model uses the modified universal soil loss equation (USLE) to estimate soil erosion and sediment transport [53].
where R X is the rainfall erosivity factor, K X refer to the soil erodibility factor, LS X is the slope length factor, C X is the crop management factor, and P X represents the soil retention factor.

Carbon Sequestration
Carbon sequestration is a cumulative reserve of carbon in terrestrial ecosystems. It consists of the following four elements: above-ground carbon storage, underground carbon storage, soil organic carbon storage, and carbon storage of dead organic matter. Regional CS is closely related to ecosystem productivity and climate regulation [54]. CS are important for climate regulation. The InVEST model uses land-use maps and four carbon pools to estimate CS in the landscape. CS data were obtained by reviewing the existing literature [55,56] and the 2014 IPCC Guidelines for National Greenhouse Gas Inventories.

Food Supply
Food supply is important for food security and urban sustainability. We counted the production of food crops, cash crops, fruits, honey, aquatic products, etc. The FS services are evaluated by calculating food calories, which translates food production into energy values. Then, the visual output is obtained according to the land use map. Calculated as follows: where ES represents the total FS calculated by calories(kcal), i refers to the crop category from one to n, E i represents the food calories by category, A i (hm 2 ) is the area of category i, Y i (kg/hm 2 ) is the yield of food i, and kcal i represents the calories per 100 g of food by category.

CES
The CES are calculated based on the value of multiple ES to compare the comprehensive level of ES. It is expressed as follows: where ES j represents the CES of year j, W i represents the weight given to each service, S ij represents the standardized value of each ES, and n represents the service category of the assessment. In order to add different ES, the value of each ES is standardized. Then, an analytic hierarchy process is used to determine the weight of each ES analytic hierarchy process (AHP) in a multi-standard decision-making method based on hierarchical decomposition to compare evaluation preferences. Finally, yaahp software was used to derive the weight of various ES. In order to evaluate the loss of ES caused by urbanization along the urban-rural gradient, we calculated the loss values of five ES and CES in each pixel. Where Vi represents the value of the ES for pixel I, the reduced value of CES was expressed as follows:

The Correlation between Urbanization and ES
ES are influenced by many factors, such as population and land use, wherein land use change is the most critical driver [12]. In order to evaluate the relationship between urbanization and ES, GDP was selected to represent the development level of economic urbanization. Population density was used to represent the scale and agglomeration of the urban population [57]. The landscape urbanization was represented by the comprehensive index of land use degree. The method of Liu et al. was used to classify the land use degree and obtain different classification indexes [58], i.e., barren land, forest, grassland, water bodies, agricultural land, urban, rural, and construction land. The higher the comprehensive index of the degree of land use, the greater the intensity of the impact of human activities. The calculation method is as follows: where L i is the comprehensive index of the land-use degree, A i represents the grading index of the land-use degree of level I, and C i represents the area proportion of land-use type i.
SPSS software was used to analyze the correlation between population, economy, landscape urbanization, and each ecosystem service in Xi'an between 2000 and 2015. Due to the size and unit of population density and GDP being different from others, the original values were logarithmically transformed before correlation analysis.

Bivariate Spatial Autocorrelation Analysis for Urbanization and ES
In order to understand the spatial aggregation relationship between urbanization and ES, the Moran's I index (LISA) was calculated by GeoDa software. Local autocorrelation analysis can reflect the correlation between urbanization and ES in different regions. The calculation method is as follows: where Y i is the value of pixel I, n is the number of pixels, and Wij is the spatial weight matrix. The clustering map generated by bivariate local Moran's I method represents different types of correlation.

Land Cover Change under Urbanization
In order to study the impact of urbanization on ES, the first step is to make clear the temporal-spatial change of urbanization, and that land use is the most obvious embodiment of urbanization. The land use patterns in different areas of Xi'an differ greatly. The core area of the city is mainly construction land with concentrated human activities. The cultivated land accounts for a large proportion in Xi'an, and most of them are distributed in the northern rural areas. Forest and grassland are mainly distributed in the Southern Qinling area ( Figure 3). Generally speaking, changes in land use type were more pronounced in the urban extended area, which was manifested in the increase of urban construction land, and the decrease of grassland and cultivated land, whereas 51.2% of the urban extended area was due to the decrease in agricultural land. The change of land use indicated that Xi'an has experienced rapid urban sprawl in the past fifteen years, and the urban expansion process has been concentrated in the periphery of the urban core area.  Overall, the WY during the study period continued to decrease. The high value of WY is mainly distributed in the southern ecological conservation area, with the highest value of 11.8 ×10 8 m 3 , followed by the rural areas (Table 1); meanwhile, the urban core area is higher than the surrounding areas ( Figure 4). Between 2000 and 2015, the WY in most areas decreased, particularly in rural areas, which decreased by 22.94%, while the urban extended area showed an increasing trend. Overall, the WY increased by 19.65%.

Water Purification
The output value of nitrogen and phosphorus in Xi'an increased initially and, then, decreased during the study period, and overall, it showed an increasing trend. This reflects that the ability of WP has been weakened. From a spatial perspective, the total nitrogen and phosphorus output value in the north is the highest, indicating that the WP is poor. The output values of nitrogen and phosphorus in the urban core area and expansion area are low. During the study period, the output of nitrogen and phosphorus in rural areas changed significantly. Between 2000 and 2015, the output of nitrogen and phosphorus decreased by 15.97% and 9.81%, respectively (Table 1), indicating the improvement of WP. The output value of nitrogen and phosphorus in the ecological conservation area did not change significantly.

Soil Conservation
Although the SC fluctuated between 2000 and 2015, the overall trend showed an upward trend. The amount of SC in the four periods were 0.19, 0.28, 0.27, and 0.30 ×10 8 tons, respectively (Table 1). Between 2000 and 2015, SC increased by 57.09%, indicating that SC has improved in recent years. The SC of the southern ecological conservation area and the forest park in the east is relatively high, which provided 73% of the SC of Xi'an in 2015, while the SC of the urban core area and urban extended area in the north and the middle is poor. The SC of rural areas and ecological conservation areas increased by 49.01% and 70.90%, while that of urban core areas decreased by 94.60%, but this had little impact on the overall SC of Xi'an.

Carbon Sequestration
Between 2000 and 2015, the CS in Xi'an increased by 5.21%. The high value of S was mainly distributed in the nature reserves and forest parks in Zhouzhi County, Huyi District, and the southern part of Chang'an District, 62% of CS were provided by the ecological conservation area in 2015. The low value was distributed in the middle and north of the study area, and the CS was the worst in the urban core area, the lowest value was just 490 kg/km 2 . Between 2000 and 2015, the decline of CS was mainly concentrated in the urban extended areas, while urban core areas and rural areas have improved in recent years.

Food Supply
The high value of FS is mainly distributed in the rural area with the highest proportion of cultivated land. The highest value is distributed in Chang'an District and Lintong District, and the lowest in Beilin District, Lianhu District, and Xincheng District. The proportion of FS in the ecological conservation area is relatively low. Between 2000 and 2015, FS in the urban core areas and the urban extended areas decreased by 84.69% and 55.83%, respectively, while in rural areas, FS increased by 22.73%. The total FS in Xi'an increased by 16.99%.

Comprehensive Ecosystem Services
The spatial distribution of CES showed great heterogeneity, in which the high value was distributed in the southwest of the ecological conservation area and the low value was distributed in the urban core area, and gradually expanded to the urban extended area in 2015 ( Figure 5). Between 2000 and 2015, the CES of the four regions showed a downward trend, of which the urban extended area changed the most significantly, with a decrease of 33.27%, followed by the rural area and the urban core areas, while the change in the ecological conservation area was relatively small, with a change rate of only 1.57%.

Relationship between Urbanization and CES
This study calculated the correlation between demographic urbanization, economic urbanization, landscape urbanization, and ES (Table 2). There was a significant spatial negative correlation between CS, FS, SC, and demographic urbanization, economic urbanization, and landscape urbanization (P < 0.01), that is to say, urbanization led to the decline of ES. The correlation with landscape urbanization is the highest, which indicated that the change of land use degree has a greater impact on ES. As WP is a negative index, it is positively correlated with demographic, economy, and landscape urbanization (P < 0.01). WY was positively correlated with demographic urbanization and economic urbanization, but the correlation between WY and landscape urbanization was not significant.    Figure 6 shows that the spatial correlation between comprehensive urbanization and ES was significantly different in different regions. (1) The spatial pattern of CS and SC was similar. High-high correlation areas mainly occurred in the Southern Qinling Mountains where the ecological environment was well protected, and the urban core area showed low-high correlation. (2) The correlation between the output of nitrogen and phosphorus and urbanization showed low-low correlation in the ecological conservation area, while high-high correlation was scattered in the urban core area. The reason for this phenomenon was that the output of nitrogen and phosphorus was a negative index, and the higher output of nitrogen and phosphorus in the urban core area means the lower WP in the built-up area. (3) FS in the urban core area was mainly low-high concentration, with low-low concentration in the southeast of Xi'an City; the correlation was not significant in the areas with wide distribution of farmland. (4) WY has the most apparent correlation between the ecological conservation area and the core area of the city, where the ecological conservation area showed evident high-low correlation, and the high-high correlation was distributed in the core area of the city. In general, the correlation between services and urbanization was most significant in the urban core area and ecological conservation area, but not significant in rural areas.

Spatiotemporal Changes of ES in Different Regions
There are significant differences in the urbanization of Xi'an in different regions, and therefore we have mapped the distribution and spatiotemporal changes of ES in different regions between 2000 and 2015. In recent years, many scholars have studied the ES in different urban areas. Consistent with the previous results, the ES show significant spatial heterogeneity [24,27,59]. SC, WY, and CS showed a similar distribution pattern, i.e., the highest value was distributed in the ecological conservation area, followed by the rural area, and the lowest value was distributed in the urban core area. In addition, the WY in the urban core area is also high. This is mainly due to the better vegetation coverage of the southern mountain areas and the eastern forest park, which are the hot areas of ES supply. Farmland accounts for a large proportion in rural areas, and as a semi-natural ecosystem, farmland can also provide higher ES [60,61]. However, it should be noted that the loss of nitrogen and phosphorus in farmland is the main non-point source pollution, therefore, the WP is relatively poor. The construction of impervious surfaces in urban core areas has changed the water cycle process inside the city, which has affected the climate and nutrient output of the city, therefore, its ability to provide ES is poor [34,62]. Due to the emphasis on the ecological environment and the improvement of green infrastructure, CS in core areas have been improved effectively between 2000 and 2015. A large number of natural surfaces in the urban extended area are transformed into impervious surfaces, which changed the hydrological cycle in the urban area, these changes have resulted in an increase in surface runoff and a weakened surface infiltration capacity [62]. WP, CS, and FS declined, which could be due to the occupation of farmland, forest, grassland, and other natural ecosystems in the urban extended area, while the construction land area increased. The change of ES in the ecological conservation area is relatively small, mainly due to the conversion of forest land and grassland to cultivated land, resulting in the weakening of WP and WY. WP and FS services in rural areas have improved. The loss of cultivated land has reduced the output of nitrogen and phosphorus pollution. Meanwhile, the improvement of agricultural technology has made up for the losses caused by urban expansion.
Consistent with the previous research results, the impact of urbanization on ES has gradually decreased from the urban extended area to the urban core area, then, to rural area and ecological conservation area [29,63]. In the urban extended area, the land use and human activity intensity showed the most evident changes, and a large amount of cultivated land was converted into urban construction land. Therefore, the impact of urban sprawl on ES is more significant than in other regions. The core area of the city has the highest level of urbanization and is also the region with the most extensive human activities. The expansion of construction land and the loss of natural vegetation have caused the ES to be most significantly affected by urban sprawl [26]. The cultivated land in rural areas accounts for a very high proportion, and the intensity of human activities is relatively low, so it is less affected by urbanization. The ecological conservation area is a hot area for providing ES. The protection and supervision measures of the ecological environment are relatively complete, which are minimally affected by the ecological environment and human activities. Understanding how urban sprawl affects ES can guide future urban planning, promote sustainable development, and improve human well-being.

The Application of the Relationship between Urbanization and ES in Planning
Our research has determined the correlation between urbanization and ES. Consistent with previous research results, the correlation between ES and demographic urbanization, economic urbanization and landscape urbanization is very significant [25,26]. In areas with a high urbanization level, the capacity of CS, FS, and SC are low, so these ES are negatively correlated with urbanization. As WP is a negative index, it is positively related to urbanization. In general, the higher the level of urbanization, the lower the ES. Because WY is high both in ecological conservation areas and urban core areas, they have a lower correlation with urbanization.
In addition, the negative correlation between ES and landscape urbanization is the most significant, which shows that urbanization, especially the improvement of land use intensity, is the significant influencing factor of ES. The direct impact of urbanization is to change the method of land use. The ecological land such as forest land and grassland has been transformed into construction land [57]. At the same time, the impact of population growth and economic development on ES is also reflected in the expansion of construction land. Generally speaking, landscape urbanization, that is, the improvement of land use intensity, is an important factor in ES [12,64]. With the expansion of cities, the impact of land use change on ES is also increasing. Therefore, we should not blindly pursue the development speed, but pay attention to the coordinated development of economic development and ecological protection.
Bivariate correlation analysis revealed the spatial clustering relationship between urbanization and ES. It is not difficult to find that there are significant differences in the clustering relationship between different regions. This is related to the difference in the level of urbanization between different areas of Xi'an. The urban core area is the area where human activities are most concentrated and most extensive, and also the area with a high demand for ES; however, the supply of ES is the lowest. The ecological conservation area with the lowest level of urbanization provides the main ES. This also confirms previous research results that there are significant differences in ES along the urban-rural gradient [28,34]. Therefore, the assessment of ES needs to consider the differences of ES in different regions.
In the process of urban construction, the urban core areas should continue to improve green infrastructure construction, improve urban green coverage, increase the use of permeable ground, and continue to promote the construction of sponge cities. Urban green space can provide many ES, such as purifying air and water, regulating climate, reducing soil erosion, and enhancing human mental health [65][66][67]. Therefore, scientific and systematic planning of urban green space should be carried out to ensure the alternate distribution of green space and buildings [68], and therefore maximize the role of urban green space ES. The urban extended area is the main area of urbanization, and the loss of ES in this area is the largest during the process of urbanization. In urban planning, we should adhere to the principle of intensive use of land resources, improve land use efficiency, control the expansion of land in rapid development areas, and ensure the construction of ecological land. The cultivated land in rural areas of Xi'an is widely distributed. Because farmland provides important FS, urban expansion should pay attention to the protection of basic farmland and upgrade agricultural infrastructure to improve crop production. At the same time, human beings have a negative impact on the environment in the process of agricultural production, therefore, the sustainable development of agriculture is still a key problem to be solved in the future [60]. The Qinling Mountains in the south of Xi'an is a hotspot area of ES and also an important natural ecological security barrier. It provides the main ES in Xi'an and it is also a key area of ecosystem protection. In the future, we should continue to implement ecological protection policies and provide better ES.
Urban development must coordinate economic development and ecological protection. Urban construction should not only focus on the expansion of scale, but the lack of sufficient ecological protection. Eco-city planning should adhere to the principle of ecological priority in order to mitigate the impact of urbanization on ES, with a view to achieving a win-win situation for economic development and ecological protection. In the future, a growing number of people will live in urban areas, and cities will develop at a faster rate. For cities to achieve sustainable development, ES are an important part of sustainability, so we must pay attention to ES in urban areas.

Conclusions
In this paper, we presented an idea to identify the impact of urban sprawl on the supply of ES in different regions. Taking Xi'an as the research area, firstly, the level of urbanization is divided according to the impervious surface and, then, the spatial and temporal distribution pattern of ES in different regions is analyzed. On this basis, the spatial interaction characteristics of urbanization and ES are evaluated by correlation analysis and cluster analysis. Our results show the following: (1) Due to the impact of urbanization level and land use patterns, there were significant differences in the distribution of ES. The high value of SC, CS, and WY is distributed in the ecological conservation area, and the lowest value is distributed in the urban core area. In addition, the WP of urban core areas is also high, and the high value of supply services such as FS is distributed in rural areas. (2) Due to different driving mechanisms, the impact of urbanization on ES in different regions was significantly different. The CS in the urban core area has been improved. The biggest change in the urban extended area and rural area is observed in WP and FS. The service changes in the ecological conservation area are not evident. Overall, the changes from the urban extended area to the urban core area, rural area, and ecological conservation area gradually decrease. (3) ES are highly correlated with demographic urbanization, economic urbanization, and landscape urbanization, and the change of landscape urbanization has the most significant impact on ES. (4) ES and urbanization have significant agglomeration in urban core areas and ecological conservation areas, but the agglomeration relationship is manifested in different ways, while the agglomeration relationship in rural areas is not significant. This result is closely related to the difference of natural conditions and socioeconomic levels in Xi'an. These findings require that the spatial heterogeneity of ES should be taken into account in policy making, and the gradient analysis allows a better description of ES supply. Ecosystem protection should not only focus on ecological conservation areas, but also on reasonable planning and management in the limited urban areas, to achieve sustainable development of cities.  Acknowledgments: The author sincerely thank Hongjuan Zhang and Zhicheng Zhang for providing helpful discussions about this work. At the same time, we are very grateful to Marina Wang editor and two anonymous reviewers for their constructive suggestions, which have greatly helped to improve the manuscript. This study was jointly supported by the National Forestry Public Welfare Industry Research Project, grant number 201304309.

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