Spatial–Temporal Interaction Relationship between Ecosystem Services and Urbanization of Urban Agglomerations in the Transitional Zone of Three Natural Regions

Clarifying the spatial interaction relationship between urbanization and multiple ecosystem services (ESs) is a prerequisite for reducing the impact of urbanization on the ecological environment and coordinating urbanization and ecological environmental protection. Urbanization is especially significant for ecologically fragile areas, where ecosystems are particularly sensitive to changes in urban patterns. This study considered the Lan–Xi (LX) urban agglomeration in three natural transitional regions using socio-economic, ecological environment, and other data, through a variety of methods, to supersede administrative boundaries and to explore the relationship between urbanization and ESs on a grid scale. The results revealed a significant negative spatial correlation between the levels of urbanization and comprehensive ESs, indicating that rapid urbanization has led to a decline in regional ESs. However, with the rapid urbanization trend from 2010 to 2018, the ESs in the LX region showed an upward trend because the implementation of ecological protection and restoration projects greatly offset the decline in ESs caused by urban expansion. We found a positive correlation between nutrient purification and the levels of urbanization among various ESs and four principal types of spatial–temporal interactions between ESs and urbanization levels. Among them, the high–high cluster areas occupied the smallest proportion, and the low–low cluster areas occupied the largest proportion.


Introduction
Urbanization is one of the most important global social and economic phenomena and involves the transformation of population structure, economic structure, geographical space, and lifestyle between rural and urban areas [1]. The urbanization of various countries primarily involves two kinds of regional spatial processes. The first is the regional concentration process of population and non-agricultural activities in urban environments of different scales and the regional promotion process of transforming non-urban landscapes into urban landscapes. The second is the process of the geographical expansion of urban culture, urban lifestyle, and values in rural areas. In the modern world, the process of urbanization is greatly accelerated, and metropolitan areas are expanding. Metropolitan circles and urban agglomerations have become core areas of national development. Urbanization has made significant contributions to social progress and economic development but poses a serious threat to the health of non-urban ecosystems. Urban expansion has proven to be a driver of the damage to key health ecosystem services [2][3][4]. Since the 21st century, the ecological and environmental issues caused by global urbanization have become increasingly serious and received extensive attention from countries and related international organizations [5].
Since reform and opening up, urbanization in China has developed rapidly. The level of urbanization increased from 17.92% in 1978 to 60.6% in 2019, with an average annual growth rate of more than 1% [6]. The rapid development of urbanization has had a profound impact on the natural ecosystem and on human economic and social systems. Current issues include the rapidly increasing population [7,8], the acceleration of industrialization [8][9][10], and extensive economic growth, accompanied by the emergence of climate warming [11]; soil, air and water pollution [12]; loss of biodiversity [13]; disordered land use; unbalanced regional development; and other ecological environment and economic social issues. Rapid urbanization has seriously aggravated the pressure on the sustainable development of society and ecosystems [14], inevitably disrupting the structure and function of the original ecosystem, leading to the loss of ecosystem services (ESs) [15,16] and seriously restricting the process of healthy urbanization and the pace of aesthetic construction in China [17]. Reducing the impact of urbanization on the ecological environment and achieving sustainable development has become a global concern [18][19][20]. Understanding the interactive relationship between urbanization and ESs is a practical need to provide a scientific basis for coordinating regional ecological environmental protection and high-quality economic and social development for national strategic needs.
Interactions between urbanization and ESs are complex; urbanization is the primary driving force for ES change and processes [21,22], and ESs restrict the development of urbanization [23,24]. Coordinating urbanization and ESs is, therefore, at the forefront of Earth system science and sustainability research [25]. Many scholars worldwide have discussed the relationship between urbanization, global change, resources, and environment from the perspective of sustainable development on a global scale [26][27][28], exploring the interaction between urbanization and natural ecosystems and global change at different spatial scales [29][30][31]. For example, André Mascarenhas [32] explored the impact of population and urbanization development in the Lisbon Metropolitan Area on land occupation and ecosystem services, and the results show that each pattern can have both positive and negative effects on the supply of ecosystem services. Fernanda Coyado Reverte [33], taking the Taubaté Basin in Brazil as the research area, explored the impact of biodiversity in highly urbanized areas on ecosystem services. The results showed that anthropogenic action is the main factor that alters the availability of local services, with emphasis on the supply of water, soils, and mineral resources, and it potentially influences the quality of life of certain species. Akanksha Balha [34], predicting the impact of urbanization on water resources in the megacity of Delhi, proposed a road map to assess future water availability in various megacities across the globe. Robert Hoyer's [35] assessment of freshwater ecosystem services in the Tualatin and Yamhill basins under climate change and urbanization showed that the estimated relative changes in ES provision under different climate and urbanization scenarios are valuable for land management decisions. Chinese scholars have primarily focused on the resources and environmental effects of urbanization [36,37], the coupling mechanism between urbanization and ecological environment [38], and the relationship between urbanization and the ecological environment (Ma et al., 2018), guided by new national urbanization development. For example, Liu [38] proposed a framework for the coupling mechanism between urbanization and the ecological environment called the Coupling Rubik's Cube. Fang [39] reviewed the major progress made in urban agglomerations in China and constructed a theoretical framework and technical path for analyzing the interaction and coupling effects of urbanization and ecological environment in megacities. The research was concentrated in the Jingjinji urban agglomeration [40], the Yangtze River Delta urban agglomeration [41], the Pearl River Delta urban agglomeration [42], and other developed coastal urban agglomerations. However, each urban agglomeration has different conditions due to development history, resource and environmental endowments, and economic foundations. Large differences exist in the level of urbanization development, and urbanization and ecological environment issues are Sustainability 2020, 12, 10211 3 of 14 complex and diverse. Ecologically fragile areas, such as the transition zones of natural areas and the structure, function, and quality of the ecological environment of the ecosystem, are particularly sensitive to changes in urban patterns [43][44][45]. Clarifying the spatial interaction relationship between urbanization and multiple ESs is a prerequisite for scientifically understanding and coordinating urbanization and ecological environmental protection from the perspective of regional functions, realizing the healthy development of urbanization, and maximizing the benefits of ESs [46].
The transitional zones of natural regions are important regional ecosystems in China. The three natural regions, the Eastern Monsoon Region, the Northwest Arid and Semi-Arid Region, and the Qinghai-Tibet Alpine Region, are typical natural geographic regions divided according to factors such as topography, geological structure, and soil vegetation [47]. The topography, climate, hydrology, soil, vegetation, and other elements of these three natural regions have obvious transitional characteristics [48], making the relationship between urbanization and ESs specific and regional; temporal and spatial changes are also more complex. Here, the Lan-Xi (LX) urban agglomeration in the transitional zone of the natural region was used as an example to (1) quantitatively assess the level of development of key ecosystem services and urbanization development of the Lan-Xi urban agglomeration in 2010 and 2018 and spatially clarify the changing process of urbanization and ecosystem services; (2) reveal the spatial-temporal interaction relationship between urbanization and ESs in the LX urban agglomeration based on the grid scale, to explore the characteristics of urbanization and ecological environment evolution in the rapid urbanization areas in the upper reaches of the Yellow River; and (3) discuss suggestions for coordinating urban development and ecological protection in the future.

Study Area
The Lan-Xi (LX) urban agglomeration is in the upper reaches of the Yellow River, at the intersection of the Loess Plateau and the Qinghai-Tibet Plateau and the three natural areas. The special geographical location and natural conditions have formed a distinctive transitional geographical and ecological pattern from west to east. The LX urban agglomeration is an important inter-provincial urban agglomeration in Western China, with relatively dense populations and towns. In 2018, the total population of the region was 11,890,300; the population density was 122 people per km 2 , which was much higher than the average of the Gansu and Qinghai provinces and the western region; and the region primarily includes 39 counties, with a total area of 97,500 km 2 . According to the scale structure, the city is divided into central cities, regional central cities, node cities, and small and medium-sized cities ( Figure 1). The level of urbanization in the LX urban agglomeration is generally low, with large spatial differences. The urbanization rate increased from 27.83% in 2000 to 61.24% in 2018-slightly higher than the national average of 59.58%. However, with rapid economic development and the continuous acceleration of urbanization, resource and environmental constraints are also increasingly intensified.

Data Sources and Preparation
The principal data used in this study consisted of the following types and came from several sources ( Table 1): (1) Land use data for 2010 and 2018 were obtained from the Resource and Environment data center of the Chinese Academy of Sciences (Resolution 30m, http://www.resdc.CN). (2) Meteorological data were obtained from the China Meteorological Science Data Sharing Service Network; daily data such as temperature, precipitation, relative humidity, wind speed, and sunshine duration were the average of the three years before and after the study period. obtained from the Statistical Yearbook of Gansu Province and the Statistical Yearbook of Qinghai Province. After considering land use and the night light index, the population density and GDP density data were downscaled to a 100 × 100 m grid by interpolation. CGCS2000 coordinates were used for spatial data with a resolution of 100 m.

Data Sources and Preparation
The principal data used in this study consisted of the following types and came from several sources ( Table 1): (1) Land use data for 2010 and 2018 were obtained from the Resource and Environment data center of the Chinese Academy of Sciences (Resolution 30m, http://www.resdc.CN). (2) Meteorological data were obtained from the China Meteorological Science Data Sharing Service Network; daily data such as temperature, precipitation, relative humidity, wind speed, and sunshine duration were the average of the three years before and after the study period.
(3) Soil data were obtained from the World Soil Database (Resolution 1000m, http://www.fao.org/). (4) The digital elevation model (DEM)was downloaded from the United States Geological Survey (USGS) website (ASTER GDEM Resolution 30m, https://www.usgs.gov/). (5) The economic and social statistical data for 2010 and 2018 were obtained from the Statistical Yearbook of Gansu Province and the Statistical Yearbook of Qinghai Province. After considering land use and the night light index, the population density and GDP density data were downscaled to a 100 × 100 m grid by interpolation. CGCS2000 coordinates were used for spatial data with a resolution of 100 m. When the threshold is 1500, the generated river network is more consistent with the actual water system

Analysis of Spatial-Temporal Interaction Relationships
This study used exploratory spatial data analysis methods to analyze the spatial relationship between urbanization and ESs. The spatial association pattern was measured and tested using the global bivariate Moran's I index and local bivariate Moran's I index (Local indicators of spatial association, LISA) [49]. The global bivariate Moran's I index was used to test for spatial correlation between ESs and urbanization, and the degree of correlation. The local bivariate Moran's I was used to explore the spatial-temporal relationship between ESs and urbanization based on a grid scale. The global Moran's I index values range from −1 to 1. Values greater than 0 indicate that the data are spatially positively correlated, and larger values are more significantly correlated. Values less than Sustainability 2020, 12, 10211 5 of 14 0 indicate that the data are negatively correlated, and larger values indicate greater spatial variability. A value of 0 indicates random space. The calculation formula is as follows [50]: where I eu is the global bivariate Moran's I index, I eu is the local bivariate Moran's I index, W ij is an n × n spatial weight matrix, z e i is the standardized value of ESs of grid i, and z u j is the standardized value of the urbanization level of grid j.

Accounting for ESs
Current ES accounting methods primarily include the value equivalent, the NPP (Net Primary Productivity) quantitative index, and the model calculation methods. In comparison, the InVEST model overcomes the unintuitive problem caused by the abstract expression of text and can visualize the evaluation result space. Therefore, the InVEST model-comprising a set of models-is recognized by most scholars and is widely used to evaluate the ESs of the LX urban agglomeration. Four sub-modules of water yield, sediment delivery ratio, nutrient delivery ratio, and carbon storage and sequestration were chosen for this study to evaluate various ESs. Owing to space limitations, the calculation method of each module is shown in reference [51].

Measurement of Urbanization Level
Urbanization is a complex system that includes historical processes of social productivity development, technological progress, and industrial structure adjustment [52][53][54]. Urbanization reflects changes in multiple systems, such as regional population, economy, society, and environment [55]. Based on the research results of the level of recent urbanization development and the development characteristics and data availability of the study area, this study chose population density, economic density, and proportion of construction land to comprehensively measure the level of urbanization from the perspective of population, economic, and land urbanization. On this basis, the range standardization method was used to standardize the abovementioned indicators on the grid scale, and the comprehensive urbanization level index was obtained by equal weight weighted calculation. The formula is as follows: where X ij is the original value of the urbanization index i and grid j, and X ij is the normalized value of X ij . X imax and are X imin the maximum and minimum values of the urbanization index i in the study area, respectively.

Spatial-Temporal Characteristics of ESs
From 2010 to 2018, all types of ESs in the LX urban agglomeration showed an upward trend ( Table 2). Among them, sediment retention and amount of carbon stored experienced the largest (151.56%) and smallest (5.82%) increases, respectively. Water yield and nutrient purification increased by more than 15%. The overall distribution pattern of the four types of ESs in the LX urban agglomeration in 2018 was consistent with that of 2010, but the spatial distribution of various services was significantly different (Figure 2). Of these, water yield was higher in the southwest and lower in the east. High-value areas were primarily concentrated in the Qinghai-Tibet Plateau in the southwestern portion of the LX urban agglomeration, such as in Guinan, Gonghe, Jianzha, Tongren, and other counties. The region has good natural background conditions, high vegetation coverage, and abundant rainfall. The maximum water yields per unit area in 2010 and 2018 were 487.28 and 828.41 mm/ha, respectively. The overall soil retention was higher in the central part and lower in the east and west. In comparison, the amount of sediment retention per unit area was much greater in 2018 than in 2010; this was closely related to the continued advancement of returning farmland to forest, returning grazing to grassland, and other ecological protection and restoration projects. Nutrient purification was consistent with the farmland distribution pattern, trending high in the southeast and low in other regions. In the southeast (i.e., Weiyuan County, Longxi County, and Lintao County), the proportion of farmland was relatively large, the natural background conditions were ordinary, and agricultural non-point-source pollution was relatively large, resulting in a large amount of nutrient purification of 17.67 kg/ha/year from 2010 to 2018. The spatial distribution pattern of high-value carbon stored areas was consistent with the distribution of forestland, and the maximum carbon stored per unit area was 430 t/ha.
The ESs of central regional towns were the least affected by rapid urbanization, and the ESs of node towns were relatively low. In comparison, the levels of urbanization in small and medium towns were not high, but the ESs were relatively high.

Spatial-Temporal Characteristics of Urbanization Level
The spatial distribution of urbanization levels of the LX urban agglomeration in 2018 and 2010 was basically the same, whereas the spatial distribution of urbanization levels within the region was significantly different (Figure 3). The central cities of Lanzhou City (including Chengguan District, Qilihe District, and Anning District) and Xining City (including Chengzhong District, Chengdong District, Chengxi District, and Chengbei District) had the highest levels of urbanization. The urbanization level index increased from 0.31 to 0.39 from 2010 to 2018 (an increase of 25.81%). The urbanization level index of regional central cities increased from 0.021 to 0.023 from 2010 to 2018 (an increase of 9.52%). Among the node cities, the urbanization level of Lanzhou New District and Linxia City was relatively high. The urbanization level of the Lanzhou area increased from 0.02 to 0.09 from 2010 to 2018 (an increase of 350%), which was consistent with the establishment and development of the Lanzhou New District. In 2012, the State Council approved the establishment of the Lanzhou New District. Since then, the Lanzhou New District has used its location advantages and resource conditions to accelerate the process of new urbanization. Linxia City is affiliated with the Linxia Hui Autonomous Prefecture, which enjoys policy advantages in ethnic areas and relatively high levels of urbanization. The urbanization level indexes increased from 0.27 to 0.29 from 2010 to 2018. The levels of urbanization in other small and medium-sized cities and towns were relatively low due to the combined effects of factors such as topography, landform, and location.

Spatial-Temporal Characteristics of Urbanization Level
The spatial distribution of urbanization levels of the LX urban agglomeration in 2018 and 2010 was basically the same, whereas the spatial distribution of urbanization levels within the region was significantly different (Figure 3). The central cities of Lanzhou City (including Chengguan District, Qilihe District, and Anning District) and Xining City (including Chengzhong District, Chengdong District, Chengxi District, and Chengbei District) had the highest levels of urbanization. The urbanization level index increased from 0.31 to 0.39 from 2010 to 2018 (an increase of 25.81%). The urbanization level index of regional central cities increased from 0.021 to 0.023 from 2010 to 2018 (an increase of 9.52%). Among the node cities, the urbanization level of Lanzhou New District and Linxia City was relatively high. The urbanization level of the Lanzhou area increased from 0.02 to 0.09 from 2010 to 2018 (an increase of 350%), which was consistent with the establishment and development of the Lanzhou New District. In 2012, the State Council approved the establishment of the Lanzhou New District. Since then, the Lanzhou New District has used its location advantages and resource conditions to accelerate the process of new urbanization. Linxia City is affiliated with the Linxia Hui Autonomous Prefecture, which enjoys policy advantages in ethnic areas and relatively high levels of urbanization. The urbanization level indexes increased from 0.27 to 0.29 from 2010 to 2018. The levels of urbanization in other small and medium-sized cities and towns were relatively low due to the combined effects of factors such as topography, landform, and location.

Spatial-Temporal Interaction Relationship between Urbanization and ESs
GeoDa software was used to measure the global bivariate Moran's I index between ESs and the levels of urbanization in 2010 and 2018 based on the grid scale (Table 3). The global bivariate Moran's I index of comprehensive ESs and the urbanization level of the LX urban agglomeration in 2010 and 2018 were −0.082 and −0.90, respectively. During the study period, there was a significant negative spatial correlation between the level of urbanization and comprehensive ESs, indicating that rapid urbanization decreased ESs in the study area. All ESs services-except for nutrient purification, which was positively correlated-had significant negative correlations with the level of urbanization.

Spatial-Temporal Interaction Relationship between Urbanization and ESs
GeoDa software was used to measure the global bivariate Moran's I index between ESs and the levels of urbanization in 2010 and 2018 based on the grid scale (Table 3) 2010 and 2018 were −0.082 and −0.90, respectively. During the study period, there was a significant negative spatial correlation between the level of urbanization and comprehensive ESs, indicating that rapid urbanization decreased ESs in the study area. All ESs services-except for nutrient purification, which was positively correlated-had significant negative correlations with the level of urbanization. The impact of urbanization levels in the study area on water yield and sediment retention were relatively stable, but the impact on carbon storage increased. The Moran's I index for the levels of urbanization and carbon stored in 2010 and 2018 was −0.056 and −0.1, respectively. Overall, there was a negative correlation between ESs and the level of urbanization, and the degrees of negative correlation between various ESs and the levels of urbanization were heterogeneous. Through local bivariate Moran's I index analysis, the spatial correlation between ESs and the level of urbanization could be divided into four types: high-high cluster areas, high-low cluster areas, low-low cluster areas, and low-high cluster areas. The spatial distribution characteristics of the local correlations between various ESs and the level of urbanization were similar (Figure 4). Comparing the LISA aggregation maps of the comprehensive ESs and the level of urbanization in 2010 and 2018 revealed consistent spatial distribution characteristics. Specifically: (1) High-high cluster areas are primarily distributed in Xining City, Datong County, Huangzhong County, Huangyuan County, Huzhu County in the Qinghai portion of the LX urban agglomeration, and Linxia City, Jishishan County, Weiyuan, and Longxi County in the Gansu region, all of which have relatively high urbanization levels and good ESs. In 2010 and 2018, the high-high areas accounted for 3.21% and 2.14% of the total area of the study area, respectively. (2) The high-low cluster areas are primarily concentrated in the northern portion of Gansu in the LX urban agglomeration, dominated by Baiyin City and belonging to the transition zone from the Tengger Desert and the Qilian Mountains to the Loess Plateau, in which the natural vegetation coverage is relatively high, with natural forests, shrubs, and meadows, and the levels of urbanization are relatively low. (3) The area of low-low cluster areas gradually decreased from 26.43% to 25.92% from 2010 to 2018 and was primarily distributed in the northern and southern parts of Qinghai Province in the LX urban agglomeration, which had few ESs and low levels of urbanization. (4) Low-high cluster areas are primarily distributed in central Lanzhou City, the regional central city of Baiyin, the area surrounding Lanzhou New District, and (to a lesser extent) in Xining City. In 2010 and 2018, the areas of low-low cluster areas accounted for 5.27% and 5.18% of the total area of the study area, respectively. The population density and GDP density in this area are relatively large, the proportion of urban construction land is high, and the overall urbanization level is high. However, rapid urbanization stresses and destroys ecological resources more seriously, resulting in low ES supply capacity.
primarily distributed in central Lanzhou City, the regional central city of Baiyin, the area surrounding Lanzhou New District, and (to a lesser extent) in Xining City. In 2010 and 2018, the areas of low-low cluster areas accounted for 5.27% and 5.18% of the total area of the study area, respectively. The population density and GDP density in this area are relatively large, the proportion of urban construction land is high, and the overall urbanization level is high. However, rapid urbanization stresses and destroys ecological resources more seriously, resulting in low ES supply capacity.

The ESs of the LX Urban Agglomeration Gradually Improved
Overall, ESs in the LX urban agglomeration grew, although the region experienced rapid and intense urbanization from 2010 to 2018. The improved ESs were consistent with those from the China Ecosystem Assessment at the national scale [56] and from other studies [57]. Such improvement, similar to that at the national scale, may be primarily attributed to national conservation policies started in 2000, which effectively alleviated ecological deterioration [54,56].

The Regional Spatial Heterogeneity of Urbanization Level was More Significant in the LX Urban Agglomeration
The overall urbanization level of the LX urban agglomeration gradually improved, but the spatial difference was more significant. From 2010 to 2018, the urbanization levels of central cities

The ESs of the LX Urban Agglomeration Gradually Improved
Overall, ESs in the LX urban agglomeration grew, although the region experienced rapid and intense urbanization from 2010 to 2018. The improved ESs were consistent with those from the China Ecosystem Assessment at the national scale [56] and from other studies [57]. Such improvement, similar to that at the national scale, may be primarily attributed to national conservation policies started in 2000, which effectively alleviated ecological deterioration [54,56].

The Regional Spatial Heterogeneity of Urbanization Level was More Significant in the LX Urban Agglomeration
The overall urbanization level of the LX urban agglomeration gradually improved, but the spatial difference was more significant. From 2010 to 2018, the urbanization levels of central cities and regional central cities in the LX urban agglomeration increased by 25.81% and 9.52%, respectively. The urbanization level of node cities and small and medium-sized cities has not increased significantly. The spatial characteristics of urbanization levels do not show a spatial weakening pattern from core to periphery, unlike the results for the Jingjiji urban agglomeration and the middle Yangtze River urban agglomeration in Eastern China obtained by other researchers [58,59]. This is because the LX urban agglomeration is located in the transition zone between the Qinghai Tibet Plateau and the Loess Plateau, and the development of the urban agglomeration is supported and constrained by the complex and diverse topography. The western portion of the LX urban agglomeration is affected by the uplift of the Qinghai-Tibet Plateau and has formed a beaded river valley landform with numerous staggered canyons and basins. The western portion is the Loess Plateau, which contains abundant loess deposits and is affected by the rising and falling of the Qilian Mountain fault block and the erosion of rivers, forming discontinuous loess platforms and multi-level river terraces. Therefore, constrained by rivers, canyons, basins, and river terraces, to a certain extent, the radiative driving effect of the central cities of the LX urban agglomeration has been hindered. Affected by the policy, the urbanization level of central cities has increased significantly, and the urbanization of regional central cities, node cities, and small and medium-sized cities has slowly developed; additionally, the urbanization level of the LX urban agglomeration has significant spatial and regional differences.

The Relationship between Urbanization and ESs
The relationship between urbanization and ecosystem services has received widespread attention in recent years [58]. Many studies have shown that urbanization causes a decline in the ecological environment [49] and ESs [14,60] and affects the landscape patterns of urban space, which then influences ESs [61]. The results of this study also confirm this. There was a significant negative spatial correlation between the urbanization level of the LX urban agglomeration and the integrated ESs. The global bivariate Moran's I index values in 2010 and 2018 were −0.082 and −0.90, respectively. The increase in absolute value indicates that the spatial difference between the two gradually increased. The rapid progress of urbanization in the LX urban agglomeration has brought rapid economic and population growth and sharply increased demand for housing, transportation, and public service facilities in urbanized areas, and the rapid development of industrial structure and non-agricultural development has led to changes in land use. The conversion of forest land, grassland, and farmland with many ESs into construction land with few ESs has seriously affected the supply capacity of ESs in the region. The impact of urbanization depends on the type of ES. Specifically, the level of urbanization has a significant negative spatial correlation with water yield, sediment retention, and carbon storage but a positive spatial correlation with nutrient purification, which indicates that the improvement of the urbanization level of the LX urban agglomeration has promoted an increase in nutrient purification in the region; this is contrary to the conclusions obtained by other scholars in the central and eastern urban agglomeration [50,54] because the InVEST nutrient purification module is primarily used to evaluate the purification function of nitrogen and phosphorus nutrients in the ecosystem. The amount of nutrient purification depends on the pollutant load and the retention rate of nitrogen and phosphorus. Point-source pollutants, such as industrial wastewater and sewage treatment plants, and agricultural non-point-source pollutants are the primary sources of nitrogen and phosphorus. On the one hand, the amount of urban sewage discharge increased significantly during the process of urbanization, resulting in a significant increase in pollutant load. On the other hand, with increased levels of urbanization, awareness among the residents regarding ecological and environmental protection and the technical means of ecological environmental governance have improved. This partially promotes the maintenance of nitrogen and phosphorus by plants and microorganisms.

Suggestions for Coordinating Urban Development and Ecological Protection
As the primary area of ecological protection and high-quality development in the transition zone of three natural regions, the LX urban agglomeration has a complex relationship with urbanization and the ecological environment [25]. Therefore, from the perspective of system science, dialectical thinking should be used to understand the spatial interaction relationship between urbanization and ecosystems to coordinate ecological protection and urban development [17]. High-high cluster areas of ESs and urbanization should continue to strengthen the monitoring and protection of the ecological environment, strictly prevent the destruction of the ecological source of high ecological service functions, improve the intensive and economical utilization of urban construction land, accelerate industrial transformation, and upgrade in a green and high-quality direction. High-low cluster areas should rationally plan urban construction land, improve land use efficiency, and develop new industries that do not rely on natural resources. Low-low cluster areas are primarily distributed in high-altitude areas. More ecological restoration projects should be planned to accelerate the ecological transformation of traditional agriculture and animal husbandry. Low-high cluster areas should strictly control the scale of urban land use and layout ecological restoration projects according to the differentiated ecological environment to avoid the destruction of the original ecological environment.

Limitations and Future Perspectives
This study superseded administrative boundaries and explored the spatial-temporal distribution pattern of ESs and urbanization in the LX urban agglomeration in the transition zone of three natural regions based on grid and regional scales to reveal the temporal and spatial interaction relationship between ESs and urbanization. However, due to research data limitations, the urbanization development level measured by population, economy, and land using 2010 and 2018 as the time node may not be comprehensive. Future research should fully utilize the advantages of modern information technology, such as big data, and consider cultural, social, and other factors to build a long-term database and to clarify the dynamic evolution of the spatial interaction between ESs and urbanization to allow more in-depth research. Additionally, we only considered four types of ESs and excluded others, such as cultural services. The beautiful scenery and invaluable cultural heritage of the study area form an important national and international tourist destination in China, and the development of tourism influences ecosystem services [62]. Cultural ESs needs to be considered in future research.

Conclusions
(1) From 2010 to 2018, various ecosystem services in the LX urban agglomeration showed an overall upward trend. The overall distribution pattern of various ecosystem services is consistent, but the spatial differences are obvious. The central and regional center cities have the fewest ecosystem services, the node cities have relatively few ecosystem services, and the small and medium-sized cities have a relatively high number of ecosystem services. (2) In 2018 and 2010, the spatial distribution pattern of urbanization level in the LX urban agglomeration was basically the same, and the regional spatial differences were significant. The urbanization level is the highest in the central cities of Lanzhou and Xining, followed by the regional central cities.
(3) The relationship between urbanization and ESs is complex. There is a negative correlation between urbanization and ESs in the LX urban agglomeration, and the degree of negative correlation between urbanization and ESs is heterogeneric. The spatial-temporal relationships between urbanization and water yield, sediment retention, and carbon storage are significantly negatively correlated, but the relationship between urbanization and nutrient purification is significantly positively correlated. Future effective measures should be taken to promote the rational population and economic development, to promote intensive and economical land use, and to achieve the sustainable development of cities.
Author Contributions: Conceptualization, X.P. and P.S.; methodology, X.P. and N.W.; software, X.P.; formal analysis, X.P. and N.W.; writing-original draft preparation, X.P.; writing-review and editing, X.P. and N.W.; funding acquisition, P.S. All authors have read and agreed to the published version of the manuscript.