Analysis of Temporal and Spatial Dynamics of Ecosystem Services and Trade-Offs/Synergies during Urbanization in the Loess Plateau, China

: As a typical ecological fragile zone and an area with a high intensity of human activities, the Loess Plateau (LP) of China has signiﬁcantly altered its ecosystem and the corresponding services under the inﬂuence of urbanization processes. However, most existing studies focus on the spatial and temporal changes of ecosystem services (ESs) and their interrelationships under the inﬂuence of ecological restoration works in the LP, leaving limited research on the impacts of urbanization on ESs. Therefore, this study constructed a research framework for exploring the spatio-temporal dynamics and interactions of ESs under the inﬂuence of urbanization based on time series data from 2000 to 2020. The results showed that: (1) based on the comprehensive urbanization level (CUL), developed and developing areas accounted for 5.63% of the total area; (2) for the whole LP, all ESs except Habitat Quality (HQ) showed an increased trend. HQ showed a trade-off with the other services, while there was a clear synergy between the other three types of services; (3) in terms of processes of urbanization, Carbon Sequestration, Water Yield and HQ gradually decreased with increased levels of urbanization, and Soil Conservation increased the least in developing areas. The trade-off between HQ and the other three services decreased with increasing urbanization, while the synergy between the other three services strengthened as urbanization deepened. These ﬁndings suggest that urbanization signiﬁcantly impacts ESs. It is necessary to implement appropriate measures (e.g., sponge city construction, urban green space, etc.) to address the impacts of urbanization on ESs.


Introduction
The current acceleration of urbanization has brought great benefits to the world's economy and society, and urbanization has become one of the most prominent features of social development [1,2].China, as the world's largest developing country, has experienced a significant increase in its urbanization rate, from 36.09% in 2000 to 63.89% in 2020 [3].However, urbanization has also resulted in ecological degradation and fragmentation of landscapes in urban areas.These changes have altered ecosystem processes and consequently impacted the ecosystem services (ESs) that can be provided to humans [4].For instance, land use changes due to urbanization have impacted carbon storage, which in turn affects the climate regulation function of ecosystems [5].Urbanization also leads to the reduction in services such as food supply, habitat quality, and soil water storage, contributing to ecological degradation [6].Additionally, high-intensity human activities result in the shrinking of vegetation areas and soil erosion, which further affect regional sustainable development.Over the past 50 years, approximately 60% of ecosystem services have been degraded or unsustainably utilized due to human activities [7].These impacts are not only evident in changes in the spatial and temporal dynamics of individual Land 2023, 12, 2136 2 of 19 services, but also in the spatial relationships between trade-offs and synergies among services [8,9].Different stages of urbanization may exhibit significant spatial and temporal variations in trade-offs and synergies between the same combinations of ecosystem services [10].Understanding these dynamics during rapid urbanization is crucial for effective ecological management [11].Therefore, studying the spatial and temporal changes in ecosystem services and their trade-offs and synergies during urbanization is of practical importance [12,13].
In recent years, the study of ESs trade-offs and synergies has become an important topic in ESs research [14].The blind pursuit of certain ESs and the lack of comprehensive considerations in ecological management will always lead to the degradation of other ESs, thus raising the issue of trade-offs and synergies among ESs [15].Neglecting tradeoffs and synergies between different ESs may reduce the total value of regional ESs and is thought to affect the stability and sustainability of regional ecosystem structures and functioning.Currently, existing studies on ESs in the urbanization process are mostly from two perspectives.One is a comparative study on the valuation of ecosystem services at the scale of administrative units [16,17].The valuation of ESs converts different types of physical quantities into a uniform amount of value for further spatial analysis and comparison.However, studies on the valuation of ESs have ignored dynamic factors such as changes in ecological conditions, spatial heterogeneity and changes in human activities [18].Therefore, considering the impact of dynamically changing factors on the value of ESs can improve the accuracy of assessment and comparison.The second is the simulation and assessment of ESs at a large spatial scale, which is mostly from a regional perspective.Few studies have analyzed the process of urbanization by dividing the region into stages or subregions.Many studies neglect the varying effects of different urbanization stages on ESs; consequently, they do not formulate stage-specific management strategies [19].In summary, it is necessary to complement studies of the impact of different stages of urbanization on ESs by revealing spatial heterogeneity in the dynamics of ESs as well as trade-offs and synergies at long time series scales.
The Loess Plateau (LP) of China is considered to be one of the world's most seriously eroded regions with fragile ecological environments [20,21].China has always attached great importance to the ecological management of the LP, adopting a series of ecological projects such as Grain for Green.These ecological projects have achieved a great deal of success, which led to a great change in the ecological quality of the original degraded ecosystem [22].In order to better consolidate the existing restoration achievements and help the regional ecosystem to develop sustainably in the future, a correct understanding of the characteristics of ESs and their relationships on the LP is particularly important for enhancing ecological conservation in the region [23].Moreover, as urbanization in the LP has accelerated in recent years, exploring ESs interrelationships during urbanization in the region is of great significance for promoting sustainable development in the area.
This study aims to explore the spatial and temporal heterogeneity of key ESs in urbanization processes in the LP and to reveal the spatial differences in trade-offs/synergies between them.Based on the long time series data from 2000 to 2020, a comprehensive urbanization level (CUL) index is constructed to characterize the urbanization development.ESs were estimated using models such as InVEST, RUSLE and CASA.The specific objectives of this study were to: (1) identify and divide different urbanization zones in the LP region by using CUL; (2) explore spatial and temporal changing trends and trade-offs/synergies of ESs in the LP; (3) analyze changes in ESs and their trade-offs/synergies under different urbanization processes.The findings of this research are expected to offer valuable decision-supporting information for effective urban planning and ecological management in the region.

Study Area
The Loess Plateau is located in the middle reaches of the Yellow River Basin, between 100 • 54~114 • 33 E, 33 • 43 ~41 • 16 N, with a total area of about 648,700 km 2 , spanning 45 cities and 341 counties in 7 provinces: Inner Mongolia, Qinghai, Ningxia, Gansu, Shaanxi, Henan and Shanxi (Figure 1), carrying a population of 118 million [24].The terrain in this area is high in the northwest and low in the southeast.The climate is mainly semi-humid and semi-arid, and the precipitation is extremely uneven in time and space distribution.The LP is one of the regions with the most concentrated conflicts between resource development and environmental protection in China [25].For a long time, the development of the LP has faced prominent problems such as serious soil erosion, sparse vegetation, high population pressure, and low productivity [26].As an important ecological security barrier area in China, it is necessary to clarify the changes and their spatial interrelationships in ESs during the urbanization process in the LP, to provide decision-making support for regional ecological management.

Study Area
The Loess Plateau is located in the middle reaches of the Yellow River Basin, between 100°54~114°33′ E, 33°43′~41°16′ N, with a total area of about 648,700 km 2 , spanning 45 cities and 341 counties in 7 provinces: Inner Mongolia, Qinghai, Ningxia, Gansu, Shaanxi, Henan and Shanxi (Figure 1), carrying a population of 118 million [24].The terrain in this area is high in the northwest and low in the southeast.The climate is mainly semi-humid and semi-arid, and the precipitation is extremely uneven in time and space distribution.The LP is one of the regions with the most concentrated conflicts between resource development and environmental protection in China [25].For a long time, the development of the LP has faced prominent problems such as serious soil erosion, sparse vegetation, high population pressure, and low productivity [26].As an important ecological security barrier area in China, it is necessary to clarify the changes and their spatial interrelationships in ESs during the urbanization process in the LP, to provide decision-making support for regional ecological management.

Data Sources
The data used in this study mainly include land use/land cover (LULC) data, precipitation, evapotranspiration, normalized difference vegetation index (NDVI), soil texture data, population density, and night lighting data from 2000-2020.The spatial

Data Sources
The data used in this study mainly include land use/land cover (LULC) data, precipitation, evapotranspiration, normalized difference vegetation index (NDVI), soil texture data, population density, and night lighting data from 2000-2020.The spatial resolution and spatial coordinate systems of all the data are standardized.In order to standardize the spatial resolution of the multi-source data, 1 km × 1 km is taken as the basic unit, and all spatial data are projected using the Albers projection.The data on the population density of the municipal districts of the seven provincial capitals on the LP come from the China Statistical Yearbook.In this study, a higher resolution of 30 m of land use data are used to measure urbanization zoning (Figure 1b).Data descriptions, applications, and sources are shown in Table 1.

Urbanization Level Quantification and Zoning
The process of urbanization refers to the transfer of population structures from rural populations to urban populations, and the transfer of industrial structures from agriculture to industry and services, accompanied by the expansion of urban areas [27,28].Therefore, from the three aspects of population urbanization, economic urbanization and spatial urbanization, the comprehensive urbanization level (CUL) index is constructed to characterize the urbanization development.Population density (P) was selected as an indicator of population urbanization as it reveals how quickly or slowly the population is growing.Night lighting data are crucial for assessing economic urbanization, serving as a key indicator of urban development and economic activity by reflecting the level of illumination, prosperity, and population density in urban areas.So, the night lighting index (N) was selected as an indicator of economic urbanization [29,30].The expansion of urban land represents urbanization in spatial terms.Therefore, we utilized the proportion of urban construction land area (U) as an indicator for spatial urbanization.Population migration, structural and quantitative changes are the most visible manifestations of population urbanization, which brings about spatial and economic development of cities [31].Utilizing the expert scoring method [32] and considering the urbanization development within the study area [33], we assigned weights of 0.4, 0.3, and 0.3 to population density, the night lighting index, and the proportion of urban construction area after normalization.Subsequently, we employed the weighted summation method to compute the comprehensive urbanization levels for the LP in both 2000 and 2020.The computational method can be described as follows: where CUL is comprehensive urbanization level, P is population density, N is night lighting index, and U is the proportion of urban construction land area.At the pixel scale, the study area was classified into developed urban areas, developing urban areas, and rural areas based on the CUL values of each pixel.Specifically, the capitals of the seven provinces in the LP were selected to estimate the thresholds of different urbanization levels, including Xi'an, Lanzhou, Xining, Yinchuan, Huhehaote, Zhengzhou and Taiyuan.The average population density (P) of the main urban area in 2000 is taken as the threshold of population urbanization, which is 0.0094 after normalization.The minimum value of night lighting (N) on the construction land of the main urban area in 2000 is taken as the threshold of economic urbanization, which is 0.3059 after normalization.The proportion of the area of the construction land of the main urban area in 2000 is taken as the threshold of spatial urbanization, which is 0.0813 after normalization.Based on the formula (1), the threshold for the integrated level of urbanization is 0.12 in this study.Study areas were classified as developed urban areas when the CUL values for both 2000 and 2020 were greater than 0.12, as developing urban areas when the CUL value was less than 0.12 in 2000 but greater than 0.12 in 2020, and as rural areas when the CUL values for both 2000 and 2020 were less than 0.12 (Table 2).The overall framework of this study is shown in Figure 2.
Land 2023, 12, x FOR PEER REVIEW 5 of 20 a key indicator of urban development and economic activity by reflecting the level of illumination, prosperity, and population density in urban areas.So, the night lighting index (N) was selected as an indicator of economic urbanization [29,30].The expansion of urban land represents urbanization in spatial terms.Therefore, we utilized the proportion of urban construction land area (U) as an indicator for spatial urbanization.Population migration, structural and quantitative changes are the most visible manifestations of population urbanization, which brings about spatial and economic development of cities [31].Utilizing the expert scoring method [32] and considering the urbanization development within the study area [33], we assigned weights of 0.4, 0.3, and 0.3 to population density, the night lighting index, and the proportion of urban construction area after normalization.Subsequently, we employed the weighted summation method to compute the comprehensive urbanization levels for the LP in both 2000 and 2020.The computational method can be described as follows: where  is comprehensive urbanization level,  is population density,  is night lighting index, and  is the proportion of urban construction land area.At the pixel scale, the study area was classified into developed urban areas, developing urban areas, and rural areas based on the CUL values of each pixel.Specifically, the capitals of the seven provinces in the LP were selected to estimate the thresholds of different urbanization levels, including Xi'an, Lanzhou, Xining, Yinchuan, Huhehaote, Zhengzhou and Taiyuan.The average population density (P) of the main urban area in 2000 is taken as the threshold of population urbanization, which is 0.0094 after normalization.The minimum value of night lighting (N) on the construction land of the main urban area in 2000 is taken as the threshold of economic urbanization, which is 0.3059 after normalization.The proportion of the area of the construction land of the main urban area in 2000 is taken as the threshold of spatial urbanization, which is 0.0813 after normalization.Based on the formula (1), the threshold for the integrated level of urbanization is 0.12 in this study.Study areas were classified as developed urban areas when the CUL values for both 2000 and 2020 were greater than 0.12, as developing urban areas when the CUL value was less than 0.12 in 2000 but greater than 0.12 in 2020, and as rural areas when the CUL values for both 2000 and 2020 were less than 0.12 (Table 2).The overall framework of this study is shown in Figure 2.

Assessment of ESs
With the help of ArcGIS 10.8 software, this research quantitatively evaluates the four types of ecosystem services-carbon sequestration, soil conservation, water yield and habitat quality-in the LP from 2000 to 2020.
(1) Carbon Sequestration Net primary production (NPP) is the key driver of terrestrial carbon cycling, balancing CO 2 and O 2 and regulating global temperatures [34].NPP was quantified using the Carnegie-Ames-Stanford approach (CASA) model [35].In the CASA model, NPP is calculated based on the absorbed photosynthetically active radiation (APAR) multiplied by the light energy conversion.
where APAR(x, t) refers to photosynthetically active radiation (g C•m 2 •month −1 ) assimilated by vegetation when pixel x is in month t, and ε(x, t) denotes actual light energy utilization (g C•MJ −1 ) of plants when pixel x is in month t.
The value of APAR is determined by total solar radiation and the fraction of photosynthetically active radiation by vegetation.The formula is as follows: where the SOL(x, t) denotes the total solar radiation (MJ•m −2 ), and FPAR(x, t) denotes the fraction of the incident photosynthetic effective radiation assimilated by vegetation.The constant 0.5 indicates the proportion of effective solar radiation accounted for by total solar radiation.
(2) Soil Conservation Soil Conservation (SC) service refers to the important functions of ecosystems in reducing the loss of soil fertility [36,37].The revised universal soil loss equation (RUSLE) model was used to calculate the amount of soil conservation.
where SC,A p and A m are soil conservation, potential soil erosion, and actual soil erosion (t•hm −2 •yr −1 ), respectively; R is the rainfall erosivity factor (MJ•mm•hm −2 •h −1 •yr −1 ); K is the soil erodibility factor (t•hm 2 h hm −2 MJ −1 mm −1 ); L is the slope length factor; S is the slope steepness factor; C is the vegetation coverage factor; and P is the soil and water conservation measure factor.
(3) Water Yield Water Yield (WY) is defined as the capacity of an ecosystem to retain water, such as by intercepting precipitation, enhancing soil infiltration, inhibiting evaporation, mitigating surface runoff, and increasing precipitation [38,39].In this research, WY was predicted using the water yield module in the InVEST model [40].The specific principles can be found in the InVEST 2.3.0 practical manual [41], and the related equations are shown as follows: where Y xj is the average annual water yield on pixel x for land cover type j; P x is the average annual precipitation on pixel x; AET xj is the average annual actual evapotranspiration on pixel x for land cover type j; ω x is a dimensionless non-physical parameter; R xj is the Budyko aridity index on pixel x for land cover type j; k is a crop factor based on the ratio between evapotranspiration ET and reference evapotranspiration ET o for crops in various developmental phases; Z is the Zhang coefficient (seasonality factor), which should be adjusted based on the simulation results of water yield and the annual runoff observations in model validation; and AWC x is the volumetric plant-available water content on pixel x.
(4) Habitat Quality Habitat Quality (HQ) service refers to the ability of an ecosystem to provide appropriate conditions for the persistence of individuals and populations, and it depends on the proximity of the habitat to human land uses and the intensity of these land uses [42,43].Combined with relevant literature, this study takes national highways, provincial highways, high-speed railways, arable land, construction land and bare land as threat sources, and obtains the habitat quality index based on the habitat quality module of the InVEST model.The formula is as follows: where Q xj is the habitat quality of grid x in landscape type j; H j indicates the habitat suitability of landscape type j; D xj is the degradation degree of grid x in landscape type j; z is the normalization constant; and k is the half-saturation constant, the values are 2.5.

Spatial Relationship Analysis
Pearson correlation analysis was used to discern the spatial changing trends of indices and trade-offs/synergies between ecosystem services from 2000 to 2020.The calculating function is as follows: where r xy represents the correlation coefficient between factors x and y, ranging from −1 to 1 (r xy > 0, positive correlation; r xy < 0, negative correlation); x i , y i are the values of two elements, respectively; and x, y denote the average values of two elements.According to the correlation coefficient test rules, r > r 0.05 = ±0.4329 is defined as significant correlation.

Urbanization Zoning
The urbanization area in the LP exhibited obvious spatial heterogeneity, as depicted in Figure 3.The distribution of urbanization zones was dispersed throughout the region, extending from the provincial capital city to the surrounding areas, without a distinctly discernible boundary.The developed areas accounted for 2.49% of the regional area and were mainly distributed in the south of the study area in Xi'an City, Zhengzhou City and their adjacent regions.In the central portion of the study area, the developed areas were mainly distributed in Taiyuan and Linfen.Conversely, in the northern and western parts of the study area, urbanization was characterized by sporadic distribution.Developing areas tended to cluster around the completed urbanized regions, accounting for 3.14% of the total area.For instance, Xianyang and Baoji demonstrated a notable level of urbanization due to the influence and expansion of the urban centers like Xi'an.The rapid urbanization observed in both developed and developing cities was primarily attributed to factors such as efficient transportation networks, favorable policy support, and higher economic prosperity.In stark contrast, rural areas still accounted for the vast majority, covering a substantial 94.36% of the total area.This indicated that the level of urbanization in the majority of the regions remains relatively low, leading to significant regional disparities between urban and rural areas.
urbanization observed in both developed and developing cities was primarily attributed to factors such as efficient transportation networks, favorable policy support, and higher economic prosperity.In stark contrast, rural areas still accounted for the vast majority, covering a substantial 94.36% of the total area.This indicated that the level of urbanization in the majority of the regions remains relatively low, leading to significant regional disparities between urban and rural areas.The spatial trend of each service shows that the proportion of areas with a significant increase in CS is the largest, accounting for 58.90% of the total area.Conversely, areas with a significant decrease in CS account for 8.88% of the total area, mainly concentrated in the southern part of the study area.The spatial distribution of SC shows similar characteristics to CS. SC increased significantly, by 36.83%,mainly distributed in the central and northern parts of the study area, and decreased significantly by 8.29% less than CS, mainly distributed in the Guanzhong Urban Agglomeration.WY increased by 83.07%, covering a significant portion of the study area.The areas with the most significant increase (10.03% of the total area) are concentrated in the southwestern part of the study area, where the terrain is higher.The overall trend of changes in HQ indicate a significant decrease along the road network, with the largest proportion of significantly decreased areas accounting for 41.25% of the total area.These areas are mainly located in urban agglomerations and their surroundings.

Changing Trends of ESs
The trade-off and synergistic relationships of the four types of services show (Figure 5) that CS vs. SC, CS vs. WY, and SC vs. WY are predominantly synergistic, with synergistic area proportions of 84.84%, 88.45%, and 72.85%, respectively.CS vs. SC has the largest significant synergistic area proportion of 56.97%; CS vs. WY is the next largest (13.71%),SC vs. WY is the smallest (9.89%), and SC vs. WY has a more pronounced synergistic effect in the southwestern part of the study area.In contrast, all services vs. HQ predominantly exhibit trade-offs.SC vs. HQ has the largest proportion of significant trade-offs (26.35%), followed by CS vs. HQ (22.26%) and WY vs. HQ (11.89%).The significant trade-off areas are primarily centered around towns and road distribution areas.The spatial trend of each service shows that the proportion of areas with a significant increase in CS is the largest, accounting for 58.90% of the total area.Conversely, areas with a significant decrease in CS account for 8.88% of the total area, mainly concentrated in the southern part of the study area.The spatial distribution of SC shows similar characteristics (13.71%),SC vs. WY is the smallest (9.89%), and SC vs. WY has a more pronounced synergistic effect in the southwestern part of the study area.In contrast, all services vs. HQ predominantly exhibit trade-offs.SC vs. HQ has the largest proportion of significant trade-offs (26.35%), followed by CS vs. HQ (22.26%) and WY vs. HQ (11.89%).The significant trade-off areas are primarily centered around towns and road distribution areas.

Changing Trends of ESs in Urbanization Zoning
The heterogeneity of different services in the LP is obvious in different urbanization areas (Figure 6).Specifically, the proportion of the area with changes in CS shows a clear stepwise change in urbanization in developed, developing and rural areas.The proportion of CS significantly increased areas varies with urbanization level: 24.94% in developed, 31.65% in developing, and 60.92% in rural areas.Conversely, the proportion of significantly decreased areas rises with urbanization: from 7.96% in rural, to 19.72% in developing, and 26.98% in developed regions.No significant difference is observed between developed and rural areas for SC.However, its proportion of significantly increased area is the lowest (29.20%) and the proportion of significantly decreased area is the largest (20.98%) in developing areas.As for WY, there are no substantial differences among regions, with all having the majority of increased areas.Notably, rural areas exhibit the largest proportion of increased areas (94.25%), but the proportion of significantly increased areas is the smallest at 9.62%.In terms of HQ, the most significant decreased occurred in developed and developing areas, accounting for 82.69% and 78.86% of the total area with sharp declines, respectively.
developing, and 26.98% in developed regions.No significant difference is observed between developed and rural areas for SC.However, its proportion of significantly increased area is the lowest (29.20%) and the proportion of significantly decreased area is the largest (20.98%) in developing areas.As for WY, there are no substantial differences among regions, with all having the majority of increased areas.Notably, rural areas exhibit the largest proportion of increased areas (94.25%), but the proportion of significantly increased areas is the smallest at 9.62%.In terms of HQ, the most significant decreased occurred in developed and developing areas, accounting for 82.69% and 78.86% of the total area with sharp declines, respectively.The synergies and trade-offs of each ES also exhibit significant variations with the level of urbanization (Figure 7).CS vs. SC and CS vs. WY are dominated by synergism in different levels of urbanization, but the degree of synergism increases with decreasing urbanization.The proportion of synergistic area is the lowest in developed areas for CS The synergies and trade-offs of each ES also exhibit significant variations with the level of urbanization (Figure 7).CS vs. SC and CS vs. WY are dominated by synergism in different levels of urbanization, but the degree of synergism increases with decreasing urbanization.The proportion of synergistic area is the lowest in developed areas for CS vs. SC and CS vs. WY, constituting 68.36% and 72.52%, respectively.Among these, 27.26% and 7.59% are significant synergies.Rural areas exhibit the highest proportion of synergistic areas at 85.76% and 89.19% for both CS vs. SC and CS vs. WY, of which 58.73% and 14.02% represent significant synergies.Similarly, SC vs. WY is also dominated by synergistic relationships but shows the lowest proportion of synergies in developing areas (65.07%), while the highest proportion of significant synergies is observed in developed areas, at 32.37%.In contrast, SC vs. HQ, WY vs. HQ, and CS vs. HQ are all dominated by tradeoffs in different regions.Significant trade-offs of SC vs. HQ and WY vs. HQ decrease with decreasing urbanization, at 54.81% and 36.79%,25.08% and 35.62%, and 22.62% and 10.85% in developed, developing, and rural areas, respectively.CS vs. HQ displays an even distribution of overall trade-offs and synergies, but it has the highest proportion of trade-offs in developing regions (59.08%).
(65.07%), while the highest proportion of significant synergies is observed in developed areas, at 32.37%.In contrast, SC vs. HQ, WY vs. HQ, and CS vs. HQ are all dominated by trade-offs in different regions.Significant trade-offs of SC vs. HQ and WY vs. HQ decrease with decreasing urbanization, at 54.81% and 36.79%,25.08% and 35.62%, and 22.62% and 10.85% in developed, developing, and rural areas, respectively.CS vs. HQ displays an even distribution of overall trade-offs and synergies, but it has the highest proportion of trade-offs in developing regions (59.08%).

Urbanization Impact on ESs
Understanding the spatial and temporal heterogeneity of ESs changes is a prerequisite for regional ecological conservation, which can provide suggestions for urban ecological management.Following extensive restoration efforts in the LP, the regional ecological environment and service capacity have efficiently improved.Carbon sequestration and soil conservation capacities have particularly increased.Simultaneously, the influence of human activities on habitat quality has become more pronounced (Figure 4).When differentiating service performance by urbanization process, the impact of human activities on diverse service types becomes more evident and prominent.For example, carbon sequestration service decreases with increasing urbanization due to the negative impacts of urbanization on vegetation through the conversion of forested and arable land into construction land [44].However, complete urbanization does not necessarily equate to the disappearance of vegetation.In this study, developed areas still have 24.94% of significantly increased area (Figure 6a), indicating that efforts have been made for urban vegetation protection and greening.For example, in the western part of Datong City, which is located in the northeastern part of the LP, the CS can still experience a significant increase in the completed urbanization area, while the CS in the developing areas decreased significantly (Figure 8(A1)).Similarly, in the Guanzhong urban agglomeration in the southern part of the LP, the southern part of the completed urbanization area of Xi'an also showed an increased trend of CS (Figure 8(B1)).Several studies have shown that with the improvement of urbanization, the government's attention to urban vegetation has led to an increase in green spaces and vegetation [45][46][47].Additionally, the urban heat island effect has contributed to the acceleration of vegetation growth [48].Regarding SC, land use changes significantly alter surface and soil properties.Rural areas with good vegetation coverage conditions exhibit better soil conservation capabilities.Instead, in developing regions, the proportion of areas with significantly decreased SC reaches 20.98%.As shown in Figure 8(A2,B2), SC in the developing areas of Datong city and the Guanzhong urban agglomeration experienced a decreased trend, especially in the developing areas of Xi'an city, which showed a significantly decreased trend in SC (Figure 8(B2)).The reasons for this phenomenon are measures such as ground consolidation and urban greening, as well as regional land fragmentation and complex surface conditions during urbanization (Figure 6b) [49].WY is strongly correlated with precipitation and underlay.The area of the LP with significantly increased precipitation is 8.05%, mainly in the southwestern part of the study area, which is consistent with the distribution of significantly increased WY (Figures 1c and 4f).Urban expansion as a result of land use changes also poses a threat to habitat quality [50][51][52].In developed areas in this study, an obvious habitat fragmentation caused by human activities was observed and resulted in a significant decline in HQ, where the portion of decreased areas is 82.69% (Figure 6d).Significant increases in HQ are observed in rural areas, such as the rural areas in the south of Xi'an city, while decreases are predominantly identified in urbanized regions (Figure 8(B4)).Trade-offs and synergies among the four services are also observed within distinct urbanization regions.CS vs. SC and CS vs. WY both exhibit synergistic relationships covering a considerable area (Figure 7a,b).SC service requires woodland and grassland with higher vegetation cover, where the roles of plant interception and protection of slopes make soil conservation function enhanced [53].Thus, CS vs. SC has increasing Trade-offs and synergies among the four services are also observed within distinct urbanization regions.CS vs. SC and CS vs. WY both exhibit synergistic relationships covering a considerable area (Figure 7a,b).SC service requires woodland and grassland with higher vegetation cover, where the roles of plant interception and protection of slopes make soil conservation function enhanced [53].Thus, CS vs. SC has increasing synergies as urbanization decreases.Both CS and WY are closely related to vegetation cover [54].Under urbanization, CS and WY are more pronounced in rural areas with better vegetation cover, leading to a stronger synergistic relationship between the two.This pattern is particularly evident in the two urban agglomerations mentioned (Figure 9(A1,A2,B1,B2)).CS vs. SC synergies are stronger in the rural areas of the Datong urban agglomeration, while in the Guanzhong urban agglomeration, they are more pronounced in the rural areas of the northwestern part compared to urban areas (Figure 9(A1,B1)).In the Datong and Guanzhong urban agglomerations, a small number of developing areas around the developed urban areas have a trade-off about CS vs. WY, but rural areas predominantly demonstrate synergistic relationships (Figure 9(A2,B2)).While both CS vs. SC and CS vs. WY relationships are predominantly synergistic, the reasons differ.CS vs. SC is due to ecological restoration measures that have improved the vegetation quality of the region, while CS vs. WY is due to the fact that both CS and WY have increased while precipitation has increased (Figure 1c).Therefore, the proportion of significant synergisms for CS vs. WY (13.71%) is much lower than CS vs. SC (56.97%).When considering the SC vs. WY relationship, both developed and rural areas exhibit similar patterns in terms of trade-offs and synergies.Developed areas display higher synergy ratios due to extensive impermeable surfaces in cities, resulting in reduced SC and WY capacity.For instance, in the urban areas of the Datong and Guanzhong urban agglomerations, synergistic areas are more widespread, with stronger synergies concentrated primarily in developed areas (Figure 9(A3,B3)).In contrast, rural areas exhibit greater synergy ratios due to superior vegetation cover and higher SC and WY capacity.However, in developing areas, the construction of green facilities lags behind the expansion of construction land driven by economic interests, and land fragmentation exacerbates this situation, contributing to a decrease in synergy proportions and an increase in trade-offs [55].Urbanization has a profound impact on HQ, leading to the most evident trade-offs with SC and WY.HQ depends on the distribution of arable land as well as forest and grassland [56].Accompanied by extensive human activities, habitats experienced erosion and fragmentation [51].In rural areas, where human activities are less intensive, the trade-offs between HQ and other services are lower.However, in the developed urban regions of the two agglomerations, trade-offs are the most pronounced for both HQ vs. SC and HQ vs. WY, decreasing as the level of urbanization decreases (Figures 7 and 9(A5,A6,B5,B6)).Similarly, urban greening can offset the negative impacts of excessive human activities, and the construction of urban green spaces and green facilities can improve habitat quality [50].The government has implemented a series of environmental management strategies, including the expansion of public green spaces and the establishment of urban green parks and buffer forests, effectively alleviating the trade-off of CS and HQ in developed regions.
Our findings provide a perspective on ESs trade-offs and synergies for urban ecological governance and planning for sustainable urban development.To enhance the quality of urban ecological space, we suggest increasing the landscape diversity of cities, which can be achieved through the implementation of various measures, such as establishing green ecological corridors, urban green parks, and green transport systems [56,57].By adopting these strategies, cities can improve key ESs, including CS, SC and HQ.Additionally, urban expansion should minimize intense soil disturbance in order to mitigate severe soil erosion, thereby preserving the regional soil and water conservation capacity [58].Furthermore, the expansion of impermeable urban surfaces leads to decreased evapotranspiration and precipitation infiltration, resulting in higher localized WY [59,60].Therefore, the improvement of hydraulic engineering and the acceleration of intra-urban storage spaces associated with sponge city construction are essential to enhance regional WY [58].To optimize the Land 2023, 12, 2136 15 of 19 overall ecosystem benefits of the region, it is crucial to increase synergies and alleviate trade-offs among services.From the perspective of urbanization zoning, areas that have completed urbanization should slow down the rate of urban expansion and continue to promote the construction of urban green facilities [61].Specifically, core cities such as Xi'an and Taiyuan should play a leading role in demonstration construction and constructing a spatial system of urban ES such as ecological corridors.In contrast, rural areas should continue to implement ecological management measures such as Grain for Green and the Natural Forest Protection Project to preserve HQ [62].Developing areas should not only focus on economic benefits; it is essential for these regions to strictly control the scale of urban land use and allocate more land for ecological purposes [32].Emphasis should be placed on developing green spaces and implementing ecological protection and restoration strategies to strike a balance between urbanization and ecological preservation.
construction of green facilities lags behind the expansion of construction land driven by economic interests, and land fragmentation exacerbates this situation, contributing to a decrease in synergy proportions and an increase in trade-offs [55].Urbanization has a profound impact on HQ, leading to the most evident trade-offs with SC and WY.HQ depends on the distribution of arable land as well as forest and grassland [56].Accompanied by extensive human activities, habitats experienced erosion and fragmentation [51].In rural areas, where human activities are less intensive, the trade-offs between HQ and other services are lower.However, in the developed urban regions of the two agglomerations, trade-offs are the most pronounced for both HQ vs. SC and HQ vs. WY, decreasing as the level of urbanization decreases (Figures 7 and 9(A5,A6,B5,B6)).Similarly, urban greening can offset the negative impacts of excessive human activities, and the construction of urban green spaces and green facilities can improve habitat quality [50].The government has implemented a series of environmental management strategies, including the expansion of public green spaces and the establishment of urban green parks and buffer forests, effectively alleviating the trade-off of CS and HQ in developed regions.Our findings provide a perspective on ESs trade-offs and synergies for urban ecological governance and planning for sustainable urban development.To enhance the quality of urban ecological space, we suggest increasing the landscape diversity of cities, which can be achieved through the implementation of various measures, such as establishing green ecological corridors, urban green parks, and green transport systems [56,57].By adopting these strategies, cities can improve key ESs, including CS, SC and HQ.Additionally, urban expansion should minimize intense soil disturbance in order to mitigate severe soil erosion, thereby preserving the regional soil and water conservation capacity [58].Furthermore, the expansion of impermeable urban surfaces leads to decreased evapotranspiration and precipitation infiltration, resulting in higher localized

Limitations and Prospects
This study reveals the spatial and temporal changes of the four important ESs and their trade-offs and synergistic relationships from 2000 to 2020 throughout the whole LP and different urbanization sub-regions, which can provide a basis for exploring the impacts of the urbanization process on the ESs.But there are still some limitations.Firstly, the selection of ESs in this study was constrained by data accessibility and the quantitative characteristics of ESs quantification.As a result, only CS, SC, WY and HQ were considered to characterize the ESs of the LP.Given the diverse range of ESs, future research should adopt a more comprehensive approach, taking into account the four aspects of provisioning, supporting, regulating, and providing cultural services [63].Secondly, the urbanization zoning method employed in this study encompassed economic urbanization, population urbanization, and spatial urbanization, combining all three aspects.However, due to the difficulty in obtaining data on social urbanization indicators, more indicators were not incorporated.Future research should consider the indicators of social urbanization [64,65] as well as consider multiple indicators of various types of urbanization, to enhance the accuracy and comprehensiveness of the urbanization zoning method.

Conclusions
This study established a research framework to investigate the spatio-temporal dynamics of ecosystem services (ESs) and their interactions during the rapid urbanization period from 2000 to 2020 in the LP.The analysis focused on the impact of urbanization stages on ESs, revealing variations in overall area and zoning perspectives.From the perspective of the whole region, CS and SC both showed an increasing trend, while WY closely correlated with precipitation patterns, with the largest increase in area.Meanwhile, HQ experienced a substantial decrease due to anthropogenic disturbances.CS vs. SC, CS vs. WY, and SC vs. WY were dominated by synergism.From the perspective of the zoning, the four types of ESs and their trade-off/synergistic relationships also differed significantly across urbanization sub-regions, indicating that urbanization had an important impact on these ESs-notably in CS, with a higher increase in rural regions compared to developed areas.SC witnessed the least increase in developing regions due to high land fragmentation, while HQ exhibited the most pronounced increase in rural areas.The degree of synergy between CS vs. SC, CS vs. WY decreased as urbanization increased.In SC vs. WY, since land fragmentation was more serious in the developing areas, the proportion of regional synergy was the lowest in urbanization in progress.Trade-offs between HQ and other services dominated various regions.In summary, rapid urbanization significantly affects the synergies and trade-offs among key ESs.Effective measures, including green infrastructure and sponge cities, are vital to mitigate negative urbanization impacts and sustain ESs.Our findings provide valuable insights for informed urban planning and ecological management, ensuring the continued provision of key ESs in rapidly urbanizing regions.

Figure 1 .
Figure 1.(a) Overview of the study area; (b) Spatial distribution of land-use types in Loess Plateau in 2020; (c) Changing trend in precipitation from 2000 to 2020.

Figure 1 .
Figure 1.(a) Overview of the study area; (b) Spatial distribution of land-use types in Loess Plateau in 2020; (c) Changing trend in precipitation from 2000 to 2020.

Figure 2 .
Figure 2. The overall framework.Note: CUL refers to comprehensive urbanization level.Figure 2. The overall framework.Note: CUL refers to comprehensive urbanization level.

Figure 2 .
Figure 2. The overall framework.Note: CUL refers to comprehensive urbanization level.Figure 2. The overall framework.Note: CUL refers to comprehensive urbanization level.

Figure 3 .
Figure 3. Urbanization zoning in the study area.

Figure 4
Figure 4 illustrates the average distribution and trend of key ESs in the LP over the past 21 years.Specifically, the average values of CS and WY exhibit a notable pattern of higher values in the southern regions and lower values in the northern areas.The southern part, particularly near the Qinling Mountains, show concentrated areas with higher average values of CS and WY.The areas with low values of SC are primarily concentrated in the northern LP and near the Guanzhong urban agglomeration.In contrast, higher values for SC are observed in the southern Qinling region.Furthermore, concerning the HQ, the areas with low values are mainly concentrated in the southern Guanzhong urban agglomeration.

Figure 3 .
Figure 3. Urbanization zoning in the study area.

Figure 4
Figure 4 illustrates the average distribution and trend of key ESs in the LP over the past 21 years.Specifically, the average values of CS and WY exhibit a notable pattern of higher values in the southern regions and lower values in the northern areas.The southern part, particularly near the Qinling Mountains, show concentrated areas with higher average values of CS and WY.The areas with low values of SC are primarily concentrated in the northern LP and near the Guanzhong urban agglomeration.In contrast, higher values for SC are observed in the southern Qinling region.Furthermore, concerning the HQ, the areas with low values are mainly concentrated in the southern Guanzhong urban agglomeration.The spatial trend of each service shows that the proportion of areas with a significant increase in CS is the largest, accounting for 58.90% of the total area.Conversely, areas with a significant decrease in CS account for 8.88% of the total area, mainly concentrated in the southern part of the study area.The spatial distribution of SC shows similar characteristics to CS. SC increased significantly, by 36.83%,mainly distributed in the central and northern parts of the study area, and decreased significantly by 8.29% less than CS, mainly distributed in the Guanzhong Urban Agglomeration.WY increased by 83.07%, covering a significant portion of the study area.The areas with the most significant increase (10.03% of the total area) are concentrated in the southwestern part of the study area, where the terrain is higher.The overall trend of changes in HQ indicate a significant decrease along the road network, with the largest proportion of significantly decreased areas accounting for 41.25% of the total area.These areas are mainly located in urban agglomerations and their surroundings.The trade-off and synergistic relationships of the four types of services show (Figure5) that CS vs. SC, CS vs. WY, and SC vs. WY are predominantly synergistic, with synergistic area proportions of 84.84%, 88.45%, and 72.85%, respectively.CS vs. SC has the largest significant synergistic area proportion of 56.97%; CS vs. WY is the next largest (13.71%),SC vs. WY is the smallest (9.89%), and SC vs. WY has a more pronounced synergistic effect

Land 2023 , 20 Figure 4 .
Figure 4. Spatial distribution of mean values and trends for the four ESs from 2000 to 2020.(a) the mean values of CS; (b) the changing trends of CS; (c) the mean values of SC; (d) the changing trends of SC; (e) the mean values of WY; (f) the changing trends of WY; (g) the mean values of HQ; (h) the changing trends of HQ.

Figure 4 .
Figure 4. Spatial distribution of mean values and trends for the four ESs from 2000 to 2020.(a) the mean values of CS; (b) the changing trends of CS; (c) the mean values of SC; (d) the changing trends of SC; (e) the mean values of WY; (f) the changing trends of WY; (g) the mean values of HQ; (h) the changing trends of HQ.

Figure 5 .
Figure 5. Spatial distribution of trade-off/synergy for the four ESs from 2000 to 2020.(a) tradeoff/synergy of CS vs. SC; (b) trade-off/synergy of CS vs. WY; (c) trade-off/synergy of CS vs. HQ; (d) trade-off/synergy of SC vs. WY; (e) trade-off/synergy of SC vs. HQ; (f) trade-off/synergy of WY vs. HQ.** means significant, * means normal.

Figure 5 .
Figure 5. Spatial distribution of trade-off/synergy for the four ESs from 2000 to 2020.(a) tradeoff/synergy of CS vs. SC; (b) trade-off/synergy of CS vs. WY; (c) trade-off/synergy of CS vs. HQ; (d) trade-off/synergy of SC vs. WY; (e) trade-off/synergy of SC vs. HQ; (f) trade-off/synergy of WY vs. HQ.** means significant, * means normal.

Figure 6 .
Figure 6.Histogram of changing trends in ESs in different urbanization zoning.(a) changing trends of CS; (b) changing trends of SC; (c) changing trends of WY; (d) changing trends of HQ.

Figure 6 .
Figure 6.Histogram of changing trends in ESs in different urbanization zoning.(a) changing trends of CS; (b) changing trends of SC; (c) changing trends of WY; (d) changing trends of HQ.

Figure 7 .
Figure 7. Histogram of trade-off/synergy in ESs in different urbanization zoning.(a) tradeoff/synergy of CS vs. SC in different urbanization zoning; (b) trade-off/synergy of CS vs. WY in different urbanization zoning; (c) trade-off/synergy of CS vs. HQ in different urbanization zoning; (d) trade-off/synergy of SC vs. WY in different urbanization zoning; (e) trade-off/synergy of SC vs. HQ in different urbanization zoning; (f) trade-off/synergy of WY vs. HQ in different urbanization zoning.** means significant, * means normal.

Figure 7 .
Figure 7. Histogram of trade-off/synergy in ESs in different urbanization zoning.(a) tradeoff/synergy of CS vs. SC in different urbanization zoning; (b) trade-off/synergy of CS vs. WY in different urbanization zoning; (c) trade-off/synergy of CS vs. HQ in different urbanization zoning; (d) trade-off/synergy of SC vs. WY in different urbanization zoning; (e) trade-off/synergy of SC vs. HQ in different urbanization zoning; (f) trade-off/synergy of WY vs. HQ in different urbanization zoning.** means significant, * means normal.

Figure 8 .
Figure 8. Spatial trends for the four ESs from 2000 to 2020 around Datong city and Guanzhong urban agglomerations.(A) is the area around Datong city, (B) is the area around Guanzhong urban agglomerations, (A1-A4) are the trends of CS, SC, WY, and HQ around Datong city extracted from Figure 4, (B1-B4) are the trends of CS, SC, WY, and HQ around Guanzhong urban agglomeration extracted from Figure 4. I** means significant increase, I* means increase, D** means significant decrease, D* means decrease.

Figure 8 .
Figure 8. Spatial trends for the four ESs from 2000 to 2020 around Datong city and Guanzhong urban agglomerations.(A) is the area around Datong city, (B) is the area around Guanzhong urban agglomerations, (A1-A4) are the trends of CS, SC, WY, and HQ around Datong city extracted from Figure 4, (B1-B4) are the trends of CS, SC, WY, and HQ around Guanzhong urban agglomeration extracted from Figure 4. I** means significant increase, I* means increase, D** means significant decrease, D* means decrease.

Table 1 .
Data sources used in this study.

Table 2 .
Classification of urbanization areas in the Loess Plateau.