A Dynamic Perspective on the Coevolution of Urbanization and Ecosystem Services in the Yangtze River Delta: Patterns, Processes, and Implications

1. Introduction

Urbanization is a global socio-economic development trend that has profound impacts on regional ecological environments. Urbanization is a multi-dimensional concept that includes demographic urbanization, economic urbanization (primarily the urbanization of industrial structures), spatial urbanization, and social and cultural urbanization (including the urbanization of living concepts and lifestyles). The urbanization process affects regional ecosystem services in complex ways by altering land use patterns, disrupting biodiversity, and influence the cycling of matter and energy [1]. Ecosystem services refer to the various benefits that humans derive from ecosystems, including provisioning services, regulating services, cultural services, and supporting services [2]. With the rapid advancement of urbanization, the ecosystem services in urban areas are increasingly challenged [3]. Studies have shown that urban expansion has encroached on large amounts of farmland, leading to a decline in provisioning services such as food production [4,5]; the increase in impervious surfaces in cities has weakened regulating services such as flood regulation [6]; the reduction of urban green spaces has diminished cultural services such as recreation and aesthetics [7]; and increased urban pollution has harmed supporting services such as soil formation and nutrient cycling [8]. Therefore, clarifying the relationship between urbanization and ecosystem services is crucial for promoting sustainable urban development.

As global urbanization accelerates, the impact of urbanization on regional ecosystem services has increasingly become a frontier and hotspot in ecological and sustainability research. Urbanization is a multi-faceted concept, and American geographer Friedmann divided the urbanization process into “Urbanization I” and “Urbanization II” [9]. He suggested that Urbanization I includes demographic urbanization, economic urbanization, and land urbanization, reflecting the explicit characteristics of urbanization; Urbanization II mainly refers to changes in people’s living concepts and lifestyles, reflecting the implicit characteristics of urbanization. Currently, scholars mainly explore the complex relationship between urbanization and ecosystem services from the following perspectives: (1) Using the ecosystem service value (ESV) equivalent factor method proposed by Costanza et al. [10] to quantitatively assess the impact of land use changes during Urbanization I on ESV [11,12,13]; (2) Analyzing the trade-offs and synergies among different ecosystem services in the context of urbanization to clarify the pathways for optimizing the configuration of multifunctional urban landscapes [14,15]; (3) Unveiling the multi-scale effects of urbanization on ecosystem services and exploring the spatial differentiation of ecosystem services within cities and across urban-rural gradients [16,17,18]. In terms of research methods, most studies have adopted spatial analysis and spatial statistics, using software such as ArcGIS, FRAGSTATS, and GeoDa to depict the impact of urbanization on the spatial patterns of ecosystem services from the perspectives of landscape patterns and spatial autocorrelation. Among these, the global spatial autocorrelation (Global Moran’s I) and local spatial autocorrelation (Anselin Local Moran’s I) methods in GeoDa effectively reveal the spatial aggregation characteristics and spatial heterogeneity of ecosystem services, gaining increasing attention from scholars [19,20,21]. Despite the substantial achievements in exploring the impact of Urbanization I on ecosystem services, there is insufficient attention to the implicit characteristics of Urbanization II. Moreover, most studies have focused only on the static relationship between urbanization and ecosystem services at a single point in time, lacking a systematic depiction of the interaction between the dynamic evolution of urbanization and the spatiotemporal changes in ecosystem services.

The Yangtze River Delta (YRD) region is one of the most urbanized and economically developed areas in China, facing significant challenges in the coordinated development of urbanization and ecosystem services. In recent years, scholars have conducted a series of studies on urbanization and ecosystem services in the YRD region. For instance, Li et al. [22] found that the coupling coordination degree between urbanization and ESV in the YRD region showed a slow upward trend during the study period, indicating an initial state of coordination. Mao and Niu [23] indicated that Jiangsu Province is a major provider of ecosystem services, while the YRD region as a whole is an ecological deficit area that needs to pay substantial ecological compensation. Yang et al. [24] demonstrated that the impact of urbanization on ecosystem services in the YRD region exhibits significant spatial heterogeneity, particularly evident along the urban-rural gradient. However, these studies also overlook the implicit characteristics of urbanization and lack a portrayal of the dynamic evolution process of urbanization, making it difficult to comprehensively understand the mechanisms through which urbanization affects ecosystem services.

To narrow the research gaps, our study focuses on the Yangtze River Delta region, comprehensively considering both the explicit and implicit characteristics of urbanization. Adopting a dynamic perspective, our study investigates the spatiotemporal evolution patterns and interaction between the urbanization process and ecosystem services. The main research contents include: (1) Assessing the urbanization indices of 41 cities in the region between 2010 and 2020 and calculating the change rates of urbanization indices during this period; (2) Quantitatively evaluating the spatial pattern characteristics and change rates of ESV in the YRD region; (3) Conducting global and local spatial autocorrelation analyses to reveal the interaction between urbanization and ecosystem services from a spatiotemporal dynamic perspective. This study aims to expand the understanding of urbanization, enrich the research theories on the impact of urbanization on ecosystem services, and provide scientific basis for regional sustainable development. Our findings not only have significant implications for the YRD region but also offer insights for formulating sustainable development strategies in other rapidly urbanizing areas around the world.

2. Research Methods

2.1. Study Area

The YRD urban agglomeration region, including Shanghai Municipality, Jiangsu Province, Zhejiang Province, and Anhui Province, is one of the most economically developed and highly urbanized areas in China (Figure 1). This region is located along the eastern coast of China and encompasses 41 cities, including Shanghai, Nanjing, Suzhou, Wuxi, Changzhou, Hangzhou, Ningbo, Wenzhou, Hefei, and Wuhu, forming a large urban cluster (Table 1). Although the area only occupies 2.2% of China’s territory (approximately 358,000 square kilometers), it is home to about 16% of the country’s population, approximately 240 million people. In 2020, the GDP of this region accounted for 24.1% of the national total GDP (approximately 24.47 trillion yuan). Over the past forty years, the YRD region has experienced significant urbanization, with an annual growth rate of 1.2%, increasing from 20.6% in 1980 to 83.4% in 2020. The YRD region is not only the economic center of China but also a core area for technological innovation and international openness. It houses a large concentration of high-tech industries and world-class enterprises, serving as a vital engine driving China’s economic development.

2.2. Data Sources

The socio-economic data used in this study primarily originate from the China Statistical Yearbook, China City Statistical Yearbook, Jiangsu Statistical Yearbook, Zhejiang Statistical Yearbook, and Shanghai Statistical Yearbook during the years 2010 to 2020. These sources provide a comprehensive range of socio-economic and urbanization indicators. To ensure data consistency and comparability, a standardized method was employed for data processing.

The land use data were obtained from the dataset titled “The 30 m annual land cover datasets and its dynamics in China from 1990 to 2020” (available at https://zenodo.org/records/5210928#.YeqApch1Mvd). This dataset offers 30-meter annual land cover data for China from 1990 to 2020. The land use types include farmland, forest, shrubland, grassland, water bodies, ice and snow, bare land, impervious surfaces, and wetlands. Data processing was conducted using ArcGIS 10.8 software, with the reclassification of land use data carried out using the lookup table tool.

2.3. Measurement of Urbanization Level

2.3.1. Construction of a Comprehensive Evaluation Index System

Constructing an evaluation index system for urbanization is a complex and multi-dimensional process, aimed at comprehensively reflecting both the explicit and implicit characteristics of urbanization. Specifically, the evaluation index system for urbanization can be divided into four dimensions: demographic urbanization, landscape urbanization, economic urbanization, and social urbanization. According to the theory of American geographer Friedmann [9], the urbanization process can be categorized into two types: “Urbanization I” and “Urbanization II.” Urbanization I focuses on the explicit geographical and economic characteristics, directly reflecting the development levels of economic activities and infrastructure within the city. Urbanization II emphasizes the implicit social and cultural characteristics, illustrating changes in the lifestyles and concepts of urban residents. Therefore, based on the specific conditions of the study area and the availability of data, a comprehensive set of indicators was developed to construct the urbanization evaluation system (Table 2).

2.3.2. Entropy Method

To scientifically and rationally determine the weights of each indicator, this study employs the entropy method. The entropy method is based on the principles of information theory, first proposed by Shannon and introduced into the field of economics by Theil to address weight distribution issues in multivariable systems [25]. The advantage of this method lies in its ability to reduce subjective judgment influence by objectively assigning weights based on the distribution characteristics of the data, thereby enhancing the scientific validity of the research. The entropy method calculates the information entropy of each indicator and assigns weights objectively based on the variability of the data itself, thereby increasing the objectivity and scientific accuracy of the weight distribution. The calculation steps are as follows:

Indicator standardization:

$$r_{ij}=\frac{x_{ij}-x_{imin}}{x_{imax}-x_{imin}}$$

(1)

$$r_{ij}=\frac{x_{imax}-x_{ij}}{x_{imax}-x_{imin}}$$

(2)

Proportion of the indicator $j$ in the year $i$:

$$p_{ij=\frac{r_{ij}}{\sum _{i=1}^m{r}_{ij}}}$$

(3)

The information entropy of the indicator $j$:

$$e_j=-k\sum _{i=1}^m{p}_{ij}\ln p_{ij}$$

(4)

The information redundancy of the indicator $j$:

$$d_j=1-e_j$$

(5)

The weight of the indicator $j$:

$$w_j=\frac{d_j}{\sum _{i=1}^m{d}_j}$$

(6)

The comprehensive evaluation index for the year $i$, that is, the level of urbanization:

$$\sum _{j=1}^n{w}_j\times r_{ij}$$

(7)

2.4. Measurement of ESV

The measurement of ESV employs the equivalent factor method, which was initially proposed by Costanza et al. in their global assessment of ESV. Subsequently, Xie et al. localized this method based on the conditions in China, making it more applicable to various Chinese ecosystems [26]. This method quantifies the contribution of ecosystems to human well-being and has been widely used globally in environmental policy and resource management decisions. Based on comprehensive research and statistical data, Xie et al. suggested that the economic value of natural food production per hectare of farmland can serve as a good proxy, as this value is trackable through effective markets [26]. Therefore, in this study, the average net profit of China’s ecosystems (excluding the cost value of human inputs) is used as the standard equivalent factor, making the equivalent factor method suitable for large-scale ESV accounting. The standard equivalent factor is estimated to be $ 503.2 per hectare.

The equivalent factor is essentially the weighting coefficient for evaluating various ecosystem services. Using this coefficient, the unit area value of a particular ecosystem service can be estimated as the product of the corresponding equivalent factor and the standard equivalent factor. The equivalent factor method is applied to estimate the ESV of the YRD urban agglomeration, and the calculation process is as follows:

$$V_{ES}=\sum _{j=1}^n\sum _{i=1}^n{A}_i\times E_{i,j}$$

(8)

$V_{ES}$ is the total ESV; $E_{i,j}$ is the value coefficient of ecosystem service for type $j$ of ecosystem type $i$; $A_i$ is the area of ecosystem type $i$. Simultaneously, the value coefficient is adjusted to $ 503.2 using regional crop production value per area, NPP (Net Primary Productivity), and precipitation data. The value equivalents for different land use and ecosystem service functions within the study area are determined through statistical analysis.

2.5. Dynamic Measurement of Urbanization and ESV

Unlike most previous studies that only focus on the static relationship between urbanization and ecosystem services at a single time point, this study attempts to explore the dynamic relationship between the evolution of urbanization and changes in ecosystem services over a continuous period of 11 years. We employed the linear regression analysis method to quantify the dynamic changes in urbanization development levels and ESV in the YRD region from 2010 to 2020. The formula is as follows:

$$V_{ES}=\sum _{j=1}^n\sum _{i=1}^n{A}_i\times E_{i,j}$$

(9)

In this context, $Y$ represents the score of urbanization development level or ecosystem service value, $X$ is the year, ${\beta }_0$ is the intercept, ${\beta }_1$ is the slope, and $\phi$is the error term. By ${\beta }_1$ calculating the annual change rates of urbanization development levels and ESV, we assess their dynamic trends during the study period. Combining the rates of change in urbanization processes with the changes in ecosystem service values, we comprehensively evaluate the impact of urbanization on ecosystem services in the YRD region.

2.6. Correlation Analysis

To reveal the spatiotemporal dynamic relationship between urbanization and ESV, this paper explores the spatial correlation between urbanization and ESV based on geographical locations. Therefore, this study first attempts to clarify the positive or negative correlation between urbanization and ESV using data from 41 units within the study area. Then, using GeoDa software, global and local spatial autocorrelation are measured through Moran’s I or Local Moran’s I (LISA) index [27]. The formula is as follows:

$$Mora{n}^{\mathrm{\text{'}}}sI=\frac{{\sum }_{i=1}^n{\sum }_{j\ne 1}^{n}Wij(Yi-\overline{Y}\left)\right(Yj-\overline{Y})}{S^2\sum _{i=1}^n\sum _{j\ne 1}^{n}Wij}$$

(10)

$$LocalMora{n}^{\mathrm{\text{'}}}sI=\frac{Yi-\overline{Y}}{{S_i}^2}\times \sum _{i=1j\ne 1}^{n}Wij\times (Yi-\overline{Y})$$

(11)

$$S^2=\frac{1}{n}\sum _{i=1}^n(Yi-\overline{Y})$$

(12)

$$\overline{Y}=\frac{1}{n}\sum _{i=1}^n{Y}_i$$

(13)

$Y_i$ and $Y_j$ respectively represent the values for unit $i$, unit $j$, and the mean value. $n$ is the number of spatial units (in this study, the study area is divided into 41 evaluation units), and $Wij$ is the spatial weight matrix.

3. Results

3.1. Spatiotemporal Dynamics of Urbanization

The overall urbanization level in the YRD region generally showed an upward trend (Figure 2). From 2010 to 2020, the level of Urbanization I gradually increased, while the development level of Urbanization II, although lower than that of Urbanization I, has accelerated significantly, especially accelerating significantly after 2017. This has led to a significant overall enhancement in the level of urbanization in the YRD region.

In terms of spatial distribution, the urbanization level in the YRD region exhibited distinct spatial heterogeneity (Figure 3). In 2010, the Urbanization I level of the study area showed a pattern of being higher in the east and lower in the west (Figure 3a1). The city of Shanghai had the highest Urbanization I level, followed by Nanjing, Suzhou, and Wuxi in Jiangsu Province, as well as Hefei and Ma’anshan in Anhui Province. The southern part of Jiangsu Province had a higher Urbanization I level than the northern part; the northern part of Zhejiang Province had a higher Urbanization I level than the southern part; Anhui Province exhibited a pattern of higher urbanization levels in the central area compared to the periphery. By 2020, the overall Urbanization I level in the YRD region had increased (Figure 3a2). Shanghai remained the area with the highest Urbanization I level. Most of Jiangsu Province, a large part of Zhejiang Province, and the areas around the provincial capital Hefei in Anhui Province all had high Urbanization I level, indicating fast economic development in these regions. We calculated the change rate of the Urbanization I index for 41 cities in the YRD region over 11 consecutive years to represent the development speed of Urbanization I within the region (Figure 3a3). The dynamic changes in the Urbanization I index revealed that coastal cities developed faster than inland cities in the YRD region. Cities like Shanghai, Nantong, Suzhou, Wuxi, Changzhou, and Nanjing in the southern part of Jiangsu Province, Hangzhou and Ningbo in the northern part of Zhejiang Province, Wenzhou in the southern part of Zhejiang Province, as well as Anqing and Chuzhou in Anhui Province, are all experiencing rapid development. Some cities, including Lianyungang in Jiangsu Province, Shaoxing in Zhejiang Province, and Huainan and Wuhu in Anhui Province, developed at a slower pace. Notably, cities like Wuhu and Tongling in Anhui Province experienced a regression in the development of Urbanization I.

In 2010, the development level of Urbanization II exhibited a pattern of being higher in the east and lower in the west (Figure 3b1). Areas with higher development levels included Shanghai and the provincial capitals of the three provinces (Nanjing, Hangzhou, and Hefei), followed by Suzhou and Wuxi in Jiangsu Province, and Ningbo and Wenzhou in Zhejiang Province. The development level of “Urbanization II” in cities of Anhui Province, excluding Hefei, and in the western cities of Jiangsu Province is relatively low. By 2020, the Urbanization II development level of the entire study area had significantly improved (Figure 3b2). In addition to Shanghai and the provincial capitals, cities like Suzhou and Ningbo saw substantial improvements in Urbanization II development. Urbanization II levels were higher in Jiangsu and Zhejiang provinces, although Anhui Province saw some improvements in Urbanization II, it still lagged behind the other provinces. It is noteworthy that the level of Urbanization II in Hangzhou and Ningbo, cities in Zhejiang Province, ranks among the highest, differing from the distribution of Urbanization I levels. Over the 11 years, the dynamic changes in the Urbanization II index showed that the development speed of Urbanization II in coastal areas was faster than in inland areas (Figure 3b3). Urbanization II in Jiangsu Province showed a trend of rapid development in the south and slow development in the north. In Zhejiang Province, Lishui experienced a regression in Urbanization II development, as did Chizhou in Anhui Province.

In 2010, cities with higher comprehensive urbanization levels were primarily Shanghai and the provincial capitals of the three provinces (Figure 3c1). The composite urbanization levels across the three provinces display distinct patterns. Jiangsu Province exhibited a pattern of being higher in the south and lower in the north; Zhejiang Province showed a pattern of being higher in the east and lower in the west; Anhui Province demonstrated a pronounced primacy distribution development pattern, characterized by high levels in Hefei and low levels in the surrounding areas. By 2020, most cities in Zhejiang Province and all cities in Jiangsu Province had high levels of comprehensive urbanization, particularly those around Shanghai, reflecting the significant radiating and driving effect of Shanghai (Figure 3c2). The comprehensive urbanization levels of cities in Anhui Province had significantly improved compared to 2010, reflecting a shift from a primacy distribution model to a more balanced distribution model in the urban development of the province. Over the 11 years, the dynamic changes in the comprehensive urbanization index showed that Shanghai had the fastest development speed (Figure 3c3). In Jiangsu Province, the development speed was faster in the south and slower in the north; the areas in Zhejiang Province with faster development were mainly located along the coast and in the provincial capital Hangzhou; Anhui Province showed a trend of fast development in the provincial capital Hefei and slower development in the surrounding areas. Notably, the comprehensive urbanization index of Tongling in Anhui Province was negative, and the comprehensive urbanization indices of Wuhu, Ma’anshan, and Huainan were low, with these cities all being located around Hefei.

3.2. Temporal and Spatial Changes in ESV

In 2020, the total ESV in the YRD region was 316,395.47 million dollars, with regulating services accounting for the largest proportion, followed by provisioning services, supporting services, and cultural services (Figure 4). From the perspective of different provinces, Anhui Province had the highest total ESV, followed by Zhejiang Province, Jiangsu Province, and Shanghai Municipality. Between 2010 and 2020, the overall ESV in the YRD region showed a declining trend, decreasing by 18,274.84 million dollars. The reduction in regulating services is the most pronounced., which decreased by 16,446.85 million dollars.

The total ecosystem service value (TESV) and its four subtypes in the YRD region exhibit significant spatial heterogeneity, with similar patterns observed in both 2010 and 2020 (Figure 5a1, a2). The TESV in the YRD region shows a distribution pattern with higher values in the south and lower values in the north, with higher values concentrated around water bodies, which is related to the strong regulating services provided by these areas. Unlike TESV, the provisioning ecosystem service value (PESV) in the YRD region has a north-high-south-low distribution pattern, with high-value areas mainly located in Anhui Province and the eastern part of Jiangsu Province, while low-value areas are mainly found in Shanghai and Zhejiang Province (Figure 5b1,b2). The high-value areas for regulating ecosystem service value (RESV) are primarily distributed in waterfront cities, such as Anqing, Chizhou, and Tongling along the Yangtze River, Suzhou, Wuxi, and Huzhou around Taihu Lake, and Yangzhou, Huai’an, and Suqian around Hongze Lake and Gaoyou Lake (Figure 5c1,c2). The supporting ecosystem service value (SESV) exhibits a distribution pattern of higher in the south and lower in the north, with higher values in southern Zhejiang and Anhui Provinces, and lower values in northern Anhui and Jiangsu Provinces (Figure 5d1,d2). The cultural ecosystem service value (CESV) displays a similar pattern characteristic to that of supporting services (Figure 5-e1, e2).

From a static spatial distribution perspective, comparing ESV between the years 2010 and 2020 alone does not reveal significant differences. However, an analysis of the ESV change rates over 11 consecutive years highlights noticeable differences between cities. In Jiangsu Province and Shanghai, the TESV decreases the fastest, followed by Zhejiang Province; Anhui Province experiences the slowest decrease in TESV, with some cities even showing a slow increase (Figure 5-a3). The regions with the fastest decrease in PESV are located in southern Jiangsu and Shanghai, while the slowest decreases are observed in southern Anhui and Zhejiang Provinces (Figure 5-b3). In Zhejiang, Lishui and Wenzhou show a slow increasing trend in PESV. The dynamic change pattern of RESV is similar to that of TESV, with the fastest decreases found in Jiangsu Province and Shanghai, and the slowest decreases in Anhui Province; some cities in Anhui also exhibit a slow increase in RESV (Figure 5-c3). The areas where the SESV is decreasing rapidly are primarily located in Zhejiang Province, with the fastest decreases observed in Shaoxing, Jinhua, and Wenzhou; in contrast, northern Anhui and northwestern Jiangsu experience slower decreases in SESV (Figure 5-d3). The areas where the CESV is decreasing rapidly include Zhejiang, Shanghai, and southern Jiangsu (Figure 5-e3). Similar to SESV, the regions with slower decreases in CESV are located in northern Anhui and northwestern Jiangsu.

3.3. Spatial Correlation Between Urbanization and ESV

We calculated the global Moran’s I statistics for the change rates of three urbanization indices (including Urbanization I, Urbanization II, and Comprehensive Urbanization) in relation to the change rates of total TESV between 2010 and 2020. The significance of the correlations was tested using Monte Carlo simulations (999 permutations; pseudo p-value: 0.01). Then, we analyzed the local spatial correlation between the change rates of urbanization indices and the change rate of TESV, and generated a LISA (Local Indicators of Spatial Association) map based on the z-test (p ≤ 0.05) to display spatial heterogeneity. The results of the global spatial correlation analysis indicated that all bivariate Moran’s I statistics were negative, suggesting a significant negative correlation between urbanization and TESV in the YRD region (Figure 6). The negative impact of Urbanization I on TESV was greater than that of Urbanization II, and the impact of Comprehensive Urbanization on TESV fell between that of Urbanization I and Urbanization II.

The results of the local spatial autocorrelation analysis divided the study area into five types of regions: high-high correlation, high-low correlation, low-high correlation, low-low correlation, and areas with no significant difference. According to the correlation analysis between Comprehensive Urbanization and TESV, low-high correlation and high-low correlation clusters appeared in the western and eastern parts of the study area, respectively (Figure 7a). In the low-high correlation clusters, where urbanization was slow but TESV declined slowly, the sacrifice of urban development speed resulted in ecological protection. In the high-low correlation clusters, where urbanization developed rapidly but TESV also declined quickly, urban development came at the cost of ecological degradation. Notably, Hefei, the provincial capital of Anhui Province, presented a high-high correlation, indicating a good balance between urban development and ecological protection. Low-low correlation areas, characterized by slow urban development and poor ecological protection, appeared next to high-low correlation areas, including Taizhou in Jiangsu Province and Huzhou in Zhejiang Province. From the correlation analysis between Urbanization I and TESV, it was observed that low-high correlation and high-low correlation clusters appeared on both sides of the study area, while high-high correlation and low-low correlation clusters were located around the low-high and high-low correlation clusters (Figure 7b). The local spatial correlation pattern between Urbanization II and TESV was similar to that between Comprehensive Urbanization and TESV, with the difference being that Jiaxing in Zhejiang Province exhibited a low-low correlation (Figure 7c).

4. Discussion

4.1. Temporal and Spatial Evolution Characteristics of Urbanization and ESV

This study, based on the multi-dimensional connotation of urbanization, reveals the phased differences in urbanization in the YRD region from two perspectives: Urbanization I (population, economic, and land urbanization) and Urbanization II (socio-cultural urbanization). The study finds that from 2010 to 2020, the level of Urbanization I gradually increased, while Urbanization II, although developing at a lower level than Urbanization I, progressed more rapidly, especially accelerating significantly after 2017 (Figure 2). This result is consistent with existing research, indicating that the urbanization process in the YRD region is continuously accelerating, particularly with a significant increase in the urbanization of living concepts and lifestyles [28]. This outcome enriches previous research that focused on single dimensions of urbanization and provides a new perspective for understanding the intrinsic mechanisms of urbanization. By comparing the spatiotemporal patterns of Urbanization I and Urbanization II, significant differences between the two were have been identified. Urbanization I exhibited a “T-shaped” distribution pattern along the eastern coast and the Yangtze River line, while Urbanization II exhibited a spatial organization pattern radiating from provincial capital cities to surrounding areas (Figure 3a2,b2). This indicates that different types of urbanization processes possess distinct spatial diffusion pathways and mechanisms, leading to varied impacts on the patterns of regional development.

The spatiotemporal patterns of the TESV indicate a severe ecological environment situation in the YRD region. During the study period, the TESV showed a declining trend, with the RESV experiencing the most significant decrease. This is consistent with the findings of Zhang et al. [29]. The spatial distribution of ESV exhibits significant heterogeneity, with high-value areas mainly concentrated around water bodies (Figure 5). This is associated with the strong regulatory service capacity provided by these water bodies. Regarding the spatial differentiation of changes in ESV, the areas with faster declines are mainly located in regions with higher levels of urbanization, such as Jiangsu and Shanghai, while areas with slower declines are concentrated in relatively less urbanized regions like Anhui and southern Zhejiang. This confirms the conclusions of many studies that faster urbanization often leads to a more rapid decline in ecosystem services [29].

4.2. Impacts of Urbanization on ESV

Global spatial autocorrelation analysis results indicate a significant negative correlation between urbanization and ESV in the YRD region, with the negative impact of Urbanization I on ESV being greater than that of Urbanization II (Figure 6). This result validates previous research, demonstrating that urbanization exerts complex effects on regional ecosystem services through mechanisms such as altering land use patterns, disturbing biodiversity, and impacting material and energy cycles [30]. A study by Cumming et al., published in Nature, pointed out that urbanization and population growth have altered the supply and demand relationship of ecosystem services, weakening the direct feedback between ecosystems and society, potentially pushing socio-ecological systems towards a collapse trap [31]. This indicates that the impact of urbanization on ecosystem services is long-term and potentially irreversible.

Local spatial autocorrelation analysis further reveals the spatial heterogeneity of urbanization’s impact on ESV. We found that low-high clusters and high-low clusters are distributed in two opposing “C-shaped” patterns (Figure 7). The eastern part of the study area is dominated by high urbanization-low ESV, characterized by high-low clustering, while the western part is dominated by low urbanization-high ESV, characterized by low-high clustering. This suggests that while eastern cities have developed rapidly, they have paid the price of environmental degradation, whereas western cities have achieved a balance between urban development and ecological protection to some extent. A study by Peng et al. also found that the impact of urbanization on ecosystem services exhibits threshold effects and spatial heterogeneity, implying that different regions should adopt differentiated management strategies [5].

Notably, Hefei, the provincial capital city, exhibits characteristics of high urbanization-high ESV, identified as high-high clustering, indicating that it has achieved good coordination in both urban development and ecological protection, which is worth emulating by other cities. This may be attributed to Hefei’s recent significant efforts in ecological civilization construction and green development strategy. A study by Grilo et al. showed that urban green spaces not only provide important cultural services for urban residents but also play a crucial role in maintaining urban biodiversity and ecosystem functions [32]. Therefore, strengthening urban green space construction is an important approach to reconcile the conflict between urban development and ecological protection.

Overall, urbanization significantly affects urban ecosystem services through various pathways such as land use changes, environmental pollution, and biodiversity loss. However, these impacts vary greatly between different regions and cities. A study by Eigenbrod et al. predicts that by 2050, urbanization will become the primary driver of spatial distribution changes in ecosystem services [33]. This means that future urban planning and management must fully consider the spatiotemporal dynamic changes in ecosystem services and formulate differentiated sustainable development strategies accordingly. Only by coordinating the relationship between urban development and ecological protection can we achieve a positive interaction between urbanization and ecosystem services, ultimately moving towards sustainable development.

4.3. Implications for Sustainable Urban Development

This study reveals the complex spatiotemporal interaction between urbanization processes and ecosystem services in the YRD region, providing important insights for regional sustainable development. Based on the study results, future urban planning and management need to focus on the following aspects:

(1) Optimizing urban development spatial patterns: Our study found that both the level of urbanization and the spatial distribution of ESV in the YRD region exhibit significant heterogeneity, with a significant negative correlation between them. This indicates that if we only pursuit urbanization development speed, while ignoring the optimization of urban spatial structure, will lead to regional ecological environment degradation. In the future, scientifically delineating urban development boundaries and prohibiting development zones is necessary to form a spatial pattern that coordinates urban development and ecological protection. In areas with high urban development intensity, it is essential to strengthen the construction of ecological infrastructure such as urban green spaces and water systems to enhance urban ecosystem services.

(2) Promoting transformation of urban development modes: This study shows that different types of urbanization processes have different impacts on ecosystem services, with the negative impacts of Urbanization I, characterized by “hard indicators” such as population and land use, being greater than those of Urbanization II, characterized by “soft indicators” such as the level of social civilization. This insight underscores that the key to achieving sustainable urban development lies in shifting the mode of urban development from extensive expansion to enhancement of intrinsic quality. We should abandon the traditional development model that sacrifices the environment for economic growth and actively develop a low-carbon economy and circular economy [34]. Special attention should be paid to changing the living concepts and lifestyles of “new citizens” in the urbanization process, helping them integrate into urban life.

(3) Improving regional ecological compensation mechanisms: The local spatial autocorrelation analysis has revealed the patterns of spatial differentiation in the impact of urbanization on ecosystem services. Eastern cities have paid the price of environmental degradation, while western cities have achieved a balance between urban development and ecological protection to some extent. This indicates the existence of an “ecological debt” issue within the urbanization process in the YRD region. It is necessary to accelerate the establishment of an ecological compensation mechanism that covers the entire area, reasonably adjust the relationships between stakeholders, and promote coordinated regional development [35]. Based on the results of this study, cities within the two opposing “C-shaped” areas on the flanks of the YRD region can be selected for paired and complementary development.

4.4. Limitations and Prospects

This study, based on the expanded connotation of urbanization, reveals the spatiotemporal evolution characteristics and interaction of urbanization and ESV in the YRD region. The study conclusions provide insights not only for sustainable development in the YRD region but also have reference value for formulating sustainable development strategies in other rapidly urbanizing areas. However, there are some limitations in this study: (1) Due to the availability of data, the construction of the urbanization evaluation index system is not comprehensive enough to fully capture the connotation of urbanization; (2) The equivalent factor method used for ESV accounting still needs improvement, and the quantitative characterization of ecosystem services needs to be strengthened; (3) The research scale is relatively coarse, making it difficult to reveal the spatial differentiation patterns of ecosystem services within cities. In the future, we plan to further expand data sources, refine evaluation indicators, improve accounting methods, and conduct multi-scale research to deeply analyze the interaction between urbanization and ecosystem services, providing more scientific decision-making basis for regional sustainable development.