Multiscale Characteristics and Drivers of the Bundles of Ecosystem Service Budgets in the Su-Xi-Chang Region, China

Managing ecosystem services (ESs) to meet human needs is critical to achieving sustainable development in rapidly urbanizing regions. Identifying ES budget bundles and analyzing their drivers at a multiscale level can facilitate management decision-making; however, further research is required in areas undergoing rapid urbanization. This study quantified the supply, demand, and budgets of six typical ESs at the county, township, and village scales in the Su-Xi-Chang region in 2020. Additionally, the influence of natural environmental and socioeconomic factors on ES budget bundles was investigated based on K-means cluster analysis and the Geodetector model. The results showed that ESs on all three scales showed a mismatch between supply and demand. The similarity in the spatial pattern of supply, demand, and budgets of ESs at the township and village scales was higher than that at the township and county scales. The location and area of surplus, balance, and deficit varied with scale. We found that population density and the proportion of impervious surfaces are the main factors influencing the formation of the ES budget bundles at different scales. In addition, the diversity and degree of interpretation of drivers varied with scale. We believe that focusing on the overall situation on a large scale and implementing precise management on a small scale can make management decisions more effective. This study can provide a scientific basis for the sustainable utilization of ESs in the Su-Xi-Chang region, and the research results and methods can provide a reference for similar studies in other rapidly urbanizing areas in the world.


Introduction
Ecosystem services (ESs) refer to the various benefits humans directly or indirectly derive from ecosystems [1]. Trees, grasslands, and water systems in cities provide residents with multiple important ESs. By 2050, 68% of the world's population is expected to live in cities [2]. Urban ecosystem services will be more closely linked to human wellbeing [3]. However, with the acceleration of urbanization, 63% of global ESs have degraded sharply [4]. The reduction in ESs will constrain economic development and affect human well-being.
The supply and demand of ES linking natural ecosystems with socioeconomic systems have been a research hotspot in the past decade [5][6][7]. The supply of ESs is the products and services that the ecosystem provides to humans [8]. The demand for ESs is the preferential consumption and use of products and services provided by ecosystems by humans [9]. The evaluation methods of supply are mainly the model simulation method [10] and the value evaluation method [1,11]. Due to easy access to input parameters and accurate evaluation results, the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) developed by the Natural Capital Project has become one of the most commonly used models for evaluating ES supply. It can be used to quantitatively assess various ESs, such as flood mitigation, urban cooling, water yield, and carbon sequestration. The assessment methods of demand are mainly based on the LULC matrix method [12], questionnaire survey

Quantification of Supply and Demand for ESs
Combining the importance to the lives of local residents, the feasibility of assessment methods and data, and referring to the study of [31,[54][55][56][57][58], we selected six typical ESs (crop production (CP), water retention (WR), PM2.5 reduction (PR), heat mitigation (HM), flood mitigation (FM), and landscape recreation (LR)) as research objects and simulated their supply and demand. The required data are shown in Table 1. The resolution of the raster data in this study is uniformly set to 30 m × 30 m, and the projected coordinate system is WGS_1984_UTM_Zone_51N. The LULC type is divided into six categories: cultivated land, forest, grassland, water area, built-up area, and unused land.

Quantification of Supply and Demand for ESs
Combining the importance to the lives of local residents, the feasibility of assessment methods and data, and referring to the study of [31,[54][55][56][57][58], we selected six typical ESs (crop production (CP), water retention (WR), PM2.5 reduction (PR), heat mitigation (HM), flood mitigation (FM), and landscape recreation (LR)) as research objects and simulated their supply and demand. The required data are shown in Table 1. The resolution of the raster data in this study is uniformly set to 30 m × 30 m, and the projected coordinate system is WGS_1984_UTM_Zone_51N. The LULC type is divided into six categories: cultivated land, forest, grassland, water area, built-up area, and unused land. Referring to the method of [60], the crop production was distributed to each pixel according to the ratio of the Vegetation Condition Index (VCI) of each cultivated land pixel. The calculation is as follows: CP i is the crop production of the ith cultivated land pixel, t; GP s is the total crop production in the Su-Xi-Chang region in 2020, GP s = 2.09 × 106 t; N is the total number of pixels of cultivated land in the study area; NDVI i is the annual NDVI value of the ith cultivated pixel; and NDVI max and NDVI min represent the maximum and minimum annual NDVI values of the cultivated land in the Su-Xi-Chang region.
(2) Demand We defined the demand for crop production (CP) by crop consumption, mainly including wheat and grains. The per capita annual crop consumption comes from the 2020 Jiangsu Statistical Yearbook. The calculation formula is: D ci is the crop demand on pixel i, t; P ci is the population on pixel i, person; and C a is the per capita annual crop consumption in the Su-Xi-Chang region, C a = 0.1221 t/person.

Water Retention
(1) Supply The supply of water retention (WR) is the difference between water yield and surface runoff [61]. Water yield was calculated using the InVEST Water Yield model. The model defines the water yield as the difference between precipitation and actual evapotranspiration. The plant's available water content was calculated using the formula from [62]. The Biophysical Table refers to the study of [63]. The surface runoff coefficient refers to Guidelines for Delineation of Ecological Protection Red Lines. The calculation is as follows: where Y ij is the annual water yield of LULC type j on pixel i, mm; AET ij is the average annual actual evapotranspiration on pixel i, mm; P i is the average annual precipitation on pixel i, mm; R ij is the Budyko Dryness Index on pixel i, dimensionless; w i is a non-physical parameter characterizing natural climate and soil properties, dimensionless; WR ij is the annual water retention of LULC type j on pixel i, mm; Runoff ij is the annual surface runoff on pixel i, mm; and C j is the average surface runoff coefficient of LULC type j. The specific coefficient values are shown in Table 2. (2) Demand This study expressed the demand for WR as the sum of agricultural, domestic, industrial, and ecological water consumption [24]. We assigned the average water consumption per unit area of cultivated land, industrial land, forest and grassland in each subdistrict to the cultivated land, industrial land, forest and grassland in each subdistrict to obtain the agricultural water consumption (D ai ), industrial water consumption(D ii ), and ecological water consumption(D ei ) on pixel i. The assignment procedure was conducted using the analysis tool in ArcGIS 10.8 software. The calculation is as follows: A a = CA subdistrict SA subdistrict (10) where D wi is the annual water consumption on pixel i, m 3 ; D ai , D di , D ii , and D ei are the agricultural water consumption, domestic water consumption, industrial water consumption, and ecological water consumption on pixel i, respectively, m 3 ; W c is the per capita annual domestic water consumption in the Su-Xi-Chang region in 2020, m 3 /person; P wi is the population on pixel i, persons; A a represents the average water consumption per unit of cultivated land; CA subdistrict represents the total water consumption of each subdistrict; SA subdistrict represents the cultivated land area of each subdistrict; and D ii and D ei are calculated in the same way as D ai .

PM 2.5 Reduction
(1) Supply In this study, the supply of PM 2.5 reduction (PR) was characterized by the amount of PM 2.5 particles adsorbed by the forest, which was calculated using the dry deposition model [64]. The calculation is as follows: where Q y is the total reduction in PM 2.5 particles by a forest in one year, g; D is the number of non-rainy days in a year; Q d is the total reduction in PM 2.5 particles by a forest in one day, g; F is the PM 2.5 dry deposition flux, g/(m 2 ·h); TCLA is the leaf surface area of different types of forest, m 2 ; T is the evaluation time, the number of hours per day, T = 24; R is the resuspension rate, R = 0.03; V d is the deposition velocity, which is related to the wind speed, V d = 0.0009 m/s; R and V d refer to the value of [65]; C h is the hourly average PM 2.5 concentration, µg/m 3 ; T c is the area of forest, T c = 900 m 2 ; and LAI is the leaf area index.
(2) Demand This study expressed the demand for PR in terms of PM 2.5 particles that exceed the WHO standard (affecting human health). The calculation formula is: (15) where D i is the demand for PR, µg; C a is the average annual PM 2.5 concentration, µg/m 3 ; PM 2.5permitted is the standard threshold of average annual PM 2.5 concentration stipulated by the World Health Organization [66], PM 2.5permitted = 10 µg/m 3 ; H is the height of the atmospheric boundary layer, H = 200 m [67]; and A is the pixel area, A = 900 m 2 .

Flood Mitigation
(1) Supply This study used the Urban Flood Risk Mitigation model of the InVEST model to calculate the supply of flood mitigation (FM). This module refers to the SCS-CN model and uses the curve number method to estimate R in . The calculation is as follows: where FM i is the supply of FM on pixel i, which is dimensionless, with a value between 0 and 1; P is the depth of rainfall for the design storm, mm; R pi is the runoff on pixel i, mm; I i is the potential retention on pixel i, mm; CN i is the runoff curve value on pixel i, dimensionless; and CN i refers to [68]; CN i takes a value between 0 and 100. When CN i is 100, it means that the rainfall is completely converted into runoff. The values of the coefficients for calculating FM are shown in Table 3. (2) Demand We refer to the research methods of [69,70] to quantify the demand of FM. The calculation is as follows: EconomicScore i (21) where FRSD i is the demand of FM on pixel i; PVI i is the population vulnerability index on pixel i; EVI i is the economic vulnerability index on pixel i, each LULC type has a corresponding score (Table 4); FRSD i , PVI i , anI EVI i are dimensionless; a and b are the weights of the population vulnerability index and the economic vulnerability index, respectively, a = b = 1/2; x a , x b and x c are Min-Max normalization coefficients representing the FR demand, PVI and EVI on pixel i, respectively. With this coefficient, FR demand, PVI and EVI can be normalized to a value between 0 and 1; Pop i is the average population density on pixel i; Old i is the population density of the elderly on pixel i; Child i is the population density of children on pixel i; and α, β, and γ are the weights of the average population density, the population density of the elderly, and the population density of children, respectively, α = β = γ = 1/3.
HM i = CC i , if i is part of a large green area or CC i ≥ CC parki CC parki , otherwise (25) where CC i is the cooling capacity index on pixel i. Its value is comprised between 0 and 1, 0 means no cooling capacity, 1 means maximum cooling capacity; shade i indicates the ability of trees to provide shade, set to 1 for trees taller than 2 m or 0 for trees below 2 m; ETI i is the evapotranspiration index on pixel i; albedo i represents the proportion of solar radiation reflected by the LULC type. Its value is between 0 and 1, where 0 represents the maximum absorption rate and 1 represents the maximum reflectivity. albedo i refers to InVEST recommended data value [72]. ET 0 is the potential evapotranspiration, mm; ET max is the maximum value of ET 0 in the study area; K c is the evapotranspiration coefficient; g i indicates whether pixel i can be regarded as a green space. A value of 1 means LULC can be regarded as a green space, 0 means LULC cannot be regarded as a green space. shade i and g i refer to the value of [73]; d(i, j) is the distance between pixel i and pixel j; d cool is the distance for which a green space has a cooling effect; j ∈ d radius from i is the set of pixels whose distance to i is less than d cool ; and HM i is the demand for HM on pixel i. The specific coefficient values are shown in Table 5. We used the urban heat island effect intensity to express the demand for HM [74]. The calculation formula is: where HMD i is the demand of HM on pixel i, which is dimensionless; x d is the Min-Max normalization coefficient representing the HM demand on pixel i. With this coefficient, the HM demand can be normalized to a value between 0 and 1, removing the units of supply and demand, and then we can further calculate the ESDR to measure the relationship between supply and demand of HM; T i is the land surface temperature on pixel i, • C; n is the total number of pixels of cultivated land; T cj is the land surface temperature of cultivated land on pixel j, • C;

Landscape Recreation
(1) Supply In this paper, the supply of landscape recreation (LR) was reflected by the green area, including forest and grassland [57]. The calculation formula is: (28) where S r is the supply of LR, m 2 /m 2 ; A greenspace,subdistrict is the green area of each jurisdiction, m 2 ; and A subdistrist is the total area of each jurisdiction, m 2 .
(2) Demand In this paper, the demand for LR was reflected by the product of population density and the per capita green area planned by the government. The calculation formula is: (29) where D r is the demand for LR, m 2 /m 2 ; P pop is the population density, person/m 2 . A guided greenspace is the per capita green area planned by the government, A guided greenspace = 13 m 2 /person.

Supply and Demand Relationship
This paper used the ecological supply-demand ratio (ESDR) to reflect the supply and demand relationship of ES, which may be surplus, balance, or deficit [75]. The calculation formula is: where S and D represent the supply and demand of ES, respectively; S max and D max represent the maximum supply and demand of ES, respectively. ESDR is expressed as a surplus (supply > demand) when it is completely greater than the range of 0. ESDR is expressed as a deficit(supply < demand) when it is completely less than the range of 0. ESDR is expressed as a balance (supply ≈ demand) when the range is close to 0. The higher the absolute value of ESDR, the greater the supply or demand for ESs.

Identifying ES Budget Bundles
We clustered the ESDR using k-means cluster analysis in Geoda spatial analysis software. The principle of clustering is that the same ES budget bundle has the greatest similarity, and the difference between different clusters is the largest. According to the elbow method [40], the optimal number of clusters was determined to be 6, and ES budget bundles at the county, township, and village scales were obtained. Before clustering analysis, the Min-Max normalization method was used to unify the ESDR values into the range of [−1, 1] to reduce the influence of unit and size among ESs. The calculation formula is: (31) where x new is the result after normalization; x max and x min are the maximum and minimum values of sample x, respectively.

Driver Analysis of ES Budget Bundles
Referring to the research of [44,47,[76][77][78] and combining it with the actual situation in the Su-Xi-Chang region, we select 14 driving factors that may affect the formation of ES budget bundles, including socioeconomic factors and natural environment factors. This paper employed the factor detector in Geodetector to identify the driving factors for the formation and spatial distribution of ES budget bundles [79]. The factor detector analyzes the degree of explanation of the independent variable X to the dependent variable Y. The calculation formula is: where h = 1, . . . , L is a layer of Y or X; N h and N are the number of units in the layer h and the whole area, respectively; σ 2 h and σ 2 are the variance of the layer h and the whole area, respectively; q is the degree of explanation of the independent variable X to the dependent variable Y-that is, X explains q% of Y. The value of q is between 0 and 1; 0 means that there is no correlation, and 1 means that Y is completely determined by X.

Crop Production
The supply of CP is spatially clustered on all three scales. The supply of CP at the township and village scales is similar to the distribution of cultivated land, and the highvalue areas are concentrated in the west and northeast of the study area (Figure 2(a2,a3)). As the scale becomes larger, the high-value areas of supply increase (Figure 2(a1)). The distribution of CP demand at the township and village scales are similar, with a clustered distribution, and randomly distributed at the county scale. The high-value areas of CP demand are concentrated in the core areas of Suzhou, Wuxi, and Changzhou, and gradually decrease to the periphery (Figure 2(b1-b3)). The ESDR of CP is clustered at the township and village scales. At the township and village scales, deficit areas are concentrated in the urban core and surrounding areas, while other areas are in surplus (Figure 2(c2,c3)). The ESDR of CP is randomly distributed at the county scale. The deficit areas increased at the county scale, the high-value areas are concentrated in the urban core areas, and the surplus is located in the western part of the study area (Figure 2(c1)). As the scale grew, CP in the study area's northeast, south, and northwest changed from surplus to balance. in the urban core and surrounding areas, while other areas are in surplus (Figure 2(c2,c3)). The ESDR of CP is randomly distributed at the county scale. The deficit areas increased at the county scale, the high-value areas are concentrated in the urban core areas, and the surplus is located in the western part of the study area (Figure 2(c1)). As the scale grew, CP in the study area's northeast, south, and northwest changed from surplus to balance. . Spatial distribution of CP (unit: t/900 m 2 ) at the county (a1,b1,c1), township (a2,b2,c2), and village (a3,b3,c3) scales; supply of CP (a1-a3); demand for CP (b1-b3); ESDR (c1-c3).

Water Retention
The supply of WR has a similar spatial distribution pattern at the three scales. The high-value supply areas are distributed in the western part of the study area, followed by the southern part. The areas with the smallest supply are distributed in the northeastern and central parts of the study area ( Figure 3(a1-a3)). The demand for WR is spatially clustered at the township and village scales. The high-value demand areas are distributed in the study area's central, southwestern, and northeastern parts (Figure 3(b2,b3)). The demand for WR is randomly distributed at the county scale. As the scale becomes larger, the areas of higher demand gradually increase (Figure 3(b1)). The ESDR of WR is clustered at the township and village scales and randomly distributed at the county scale. At the village scale, the deficit and high-demand areas are distributed similarly, and supply and demand are balanced in the remaining areas (Figure 3(c3)). At the township scale, the southern part of the study area gradually turned from balance at the village scale to surplus. The high-value areas of the deficit are concentrated in the northeastern part of the study area (Figure 3(c2)), that is, along the Yangtze River, where there are many factories with high water demands. The deficit areas increase at the county scale and are distributed in the study area's northeastern, western, and southern parts. High-value deficit areas are Figure 2. Spatial distribution of CP (unit: t/900 m 2 ) at the county (a1,b1,c1), township (a2,b2,c2), and village (a3,b3,c3) scales; supply of CP (a1-a3); demand for CP (b1-b3); ESDR (c1-c3).

Water Retention
The supply of WR has a similar spatial distribution pattern at the three scales. The high-value supply areas are distributed in the western part of the study area, followed by the southern part. The areas with the smallest supply are distributed in the northeastern and central parts of the study area ( Figure 3(a1-a3)). The demand for WR is spatially clustered at the township and village scales. The high-value demand areas are distributed in the study area's central, southwestern, and northeastern parts (Figure 3(b2,b3)). The demand for WR is randomly distributed at the county scale. As the scale becomes larger, the areas of higher demand gradually increase (Figure 3(b1)). The ESDR of WR is clustered at the township and village scales and randomly distributed at the county scale. At the village scale, the deficit and high-demand areas are distributed similarly, and supply and demand are balanced in the remaining areas (Figure 3(c3)). At the township scale, the southern part of the study area gradually turned from balance at the village scale to surplus. The high-value areas of the deficit are concentrated in the northeastern part of the study area (Figure 3(c2)), that is, along the Yangtze River, where there are many factories with high water demands. The deficit areas increase at the county scale and are distributed in the study area's northeastern, western, and southern parts. High-value deficit areas are concentrated in urban core areas. The areas where supply and demand are balanced are distributed in the central and southern parts of the study area (Figure 3(c1)).

PM2.5 Reduction
The supply and demand of PR have significant spatial agglomeration at the county, township, and village scales. The high-value supply areas on the three scales are concentrated in the west and south of the study area (Figure 4(a1-a3)). As the scale becomes larger, the high-value areas of supply increase gradually. The demand for PR gradually increases from the southeast to the northwest of the study area (Figure 4(b1-b3)). As the scale becomes larger, the high-value demand areas gradually decrease and are concentrated northwest of the study area. Although the southwest of the study area provides some supply of PR, the spatial pattern of ESDR and demand is similar because of the large demand. The PR of the entire study area is in deficit on all three scales. The high-value deficit areas are concentrated in the northwest of the study area, showing a downward trend from northwest to southeast (Figure 4(c1-c3)).

PM 2.5 Reduction
The supply and demand of PR have significant spatial agglomeration at the county, township, and village scales. The high-value supply areas on the three scales are concentrated in the west and south of the study area (Figure 4(a1-a3)). As the scale becomes larger, the high-value areas of supply increase gradually. The demand for PR gradually increases from the southeast to the northwest of the study area (Figure 4(b1-b3)). As the scale becomes larger, the high-value demand areas gradually decrease and are concentrated northwest of the study area. Although the southwest of the study area provides some supply of PR, the spatial pattern of ESDR and demand is similar because of the large demand. The PR of the entire study area is in deficit on all three scales. The high-value deficit areas are concentrated in the northwest of the study area, showing a downward trend from northwest to southeast (Figure 4(c1-c3)).

Flood Mitigation
The high-value areas of FM supply at the three scales are distributed in the west, south, and north of the study area, with spatial aggregation. The distribution of FM demand at the township and village scales is similar, with a clustered distribution ( Figure 5(b2,b3)). The demand for FM is randomly distributed at the county scale, with high-value areas distributed in the central and northeastern parts of the study area ( Figure 5(b1)). The ESDR of FM is clustered at the three scales. The high-deficit areas of FM are concentrated in the urban core ( Figure 5(c1-c3)). At the county scale, the western part of the study area is in surplus, the southern part is in balance, and the northeastern part is in deficit (Figure 5(c1)). There are many surplus and balance areas at the township scale ( Figure 5(c2)). At the village scale, the western and southeastern parts of the study area are in surplus, and other areas are in deficit ( Figure 5(c3)).

Flood Mitigation
The high-value areas of FM supply at the three scales are distributed in the west, south, and north of the study area, with spatial aggregation. The distribution of FM demand at the township and village scales is similar, with a clustered distribution ( Figure  5(b2,b3)). The demand for FM is randomly distributed at the county scale, with high-value areas distributed in the central and northeastern parts of the study area ( Figure 5(b1)). The ESDR of FM is clustered at the three scales. The high-deficit areas of FM are concentrated in the urban core ( Figure 5(c1-c3)). At the county scale, the western part of the study area is in surplus, the southern part is in balance, and the northeastern part is in deficit ( Figure  5(c1)). There are many surplus and balance areas at the township scale ( Figure 5(c2)). At the village scale, the western and southeastern parts of the study area are in surplus, and other areas are in deficit ( Figure 5(c3)).

Heat Mitigation
The supply of HM is clustered at the township and village scales and randomly distributed at the county scale. The high-value supply areas are mainly distributed in the west and south of the study area. The low-value areas are located in the core areas and Figure 5. Spatial distribution of FM (unit: dimensionless, the value is between 0-1) at the county (a1,b1,c1), township (a2,b2,c2), and village (a3,b3,c3) scales; supply of FM (a1-a3); demand for FM (b1-b3); ESDR (c1-c3).

Heat Mitigation
The supply of HM is clustered at the township and village scales and randomly distributed at the county scale. The high-value supply areas are mainly distributed in the west and south of the study area. The low-value areas are located in the core areas and the surrounding areas of Suzhou, Wuxi, and Changzhou ( Figure 6(a1-a3)). The demand for HM is spatially clustered at three scales. Demand in the northeast of the study area is higher than in the southwest of the study area. The high-value demand areas are concentrated in the urban core areas (Figure 6(b1-b3)). The surplus level of HM in the west and southwest of the study area is optimal at the three scales. At the village scale, deficit areas are concentrated in the urban core and surrounding areas (Figure 6(c3)). The spatial distribution characteristics of ESDR also change as the scale becomes larger. For example, at the township scale, the HM deficit area decreases (Figure 6(c2)); at the county scale, the deficit area is transformed into balance (Figure 6(c1)).

Landscape Recreation
The supply and demand of LR are randomly distributed at the county scale and clustered at the township and village scales. The supply of LR is mainly located in the west of the study area and along Taihu Lake. As the scale becomes larger, the high-value areas of supply increase gradually (Figure 7(a1-a3)). The high-value areas of demand for LR are concentrated in the urban core and surrounding areas (Figure 7(b1-b3)). At the county scale, the high-value areas of demand are concentrated in the core areas of Suzhou, Wuxi, and Changzhou (Figure 7(b1)). The ESDR of LR at the county scale is randomly distributed. Most areas at the county scale are in surplus and balance. The high-value areas of surplus are concentrated in Changzhou. The deficit areas are concentrated in the core areas of Suzhou and Wuxi (Figure 7(c1)). The value of ESDR also changes as the scale be- Figure 6. Spatial distribution of HM (unit: dimensionless, the value is between 0-1) at the county (a1,b1,c1), township (a2,b2,c2), and village (a3,b3,c3) scales; supply of HM (a1-a3); demand for HM (b1-b3); ESDR (c1-c3).

Landscape Recreation
The supply and demand of LR are randomly distributed at the county scale and clustered at the township and village scales. The supply of LR is mainly located in the west of the study area and along Taihu Lake. As the scale becomes larger, the high-value areas of supply increase gradually (Figure 7(a1-a3)). The high-value areas of demand for LR are concentrated in the urban core and surrounding areas (Figure 7(b1-b3)). At the county scale, the high-value areas of demand are concentrated in the core areas of Suzhou, Wuxi, and Changzhou (Figure 7(b1)). The ESDR of LR at the county scale is randomly distributed. Most areas at the county scale are in surplus and balance. The high-value areas of surplus are concentrated in Changzhou. The deficit areas are concentrated in the core areas of Suzhou and Wuxi (Figure 7(c1)). The value of ESDR also changes as the scale becomes smaller. For example, some areas at the township scale change from surplus to balance (Figure 7(c2)), and many areas at the village scale change from balance to deficit (Figure 7(c3)).

Multiscale Pattern Characteristics of ES Budget Bundles
Based on K-means clustering analysis, this study identified six ES budget bundles at the county (Figure 8), township (Figure 9), and village scales( Figure 10). Each bundle has similar ES characteristics [38]. This study found that ES budget bundles at the township and village scales have similar spatial patterns. However, compared with the township scale, the village scale is finer and more similar to the spatial distribution of LULC. However, the distribution of ES budget bundles at the township/village scale and the county scale are significantly different. The spatial distribution of ES budget bundles at the county scale is similar to that of municipal administrative districts, forming the spatial pattern of Suzhou, Wuxi, Changzhou, Taihu Lake and its surrounding areas, and urban core areas.

County Scale
Bundle 1 is mainly distributed in the north of Wuxi and Changzhou, accounting for 23.3% of the total study area. CP, WR, and LR are slightly higher than the mean value of the study area; PR, HM, and FM are slightly lower than the mean value of the study area. Bundle 2 is distributed in the relatively flat eastern region, concentrated in Suzhou, accounting for 32.44% of the entire study area. Bundle 2 has the highest level of PR and Figure 7. Spatial distribution of LR (unit: m 2 /900 m 2 ) at the county (a1,b1,c1), township (a2,b2,c2), and village (a3,b3,c3) scales; supply of LR (a1-a3); demand for LR (b1-b3); ESDR (c1-c3).

Multiscale Pattern Characteristics of ES Budget Bundles
Based on K-means clustering analysis, this study identified six ES budget bundles at the county (Figure 8), township (Figure 9), and village scales ( Figure 10). Each bundle has similar ES characteristics [38]. This study found that ES budget bundles at the township and village scales have similar spatial patterns. However, compared with the township scale, the village scale is finer and more similar to the spatial distribution of LULC. However, the distribution of ES budget bundles at the township/village scale and the county scale are significantly different. The spatial distribution of ES budget bundles at the county scale is similar to that of municipal administrative districts, forming the spatial pattern of Suzhou, Wuxi, Changzhou, Taihu Lake and its surrounding areas, and urban core areas.
ESs are all higher than the mean value, especially WR and HM, which are at the highest level in the study area. Bundle 6 is distributed in the northwest of the study area, located in the core area of Changzhou, accounting for 1.63% of the total study area. The six ESs in bundle 6 are all below the mean value, especially PR, which is at the lowest level in the study area.

Township Scale
Bundle 1 is distributed in the surrounding areas of the urban core areas of Suzhou, Wuxi, and Changzhou, accounting for 2.21% of the total study area. All ESs in bundle 1 have below-average values, except for WR, which is slightly above the mean value of the study area. Bundle 2 has the largest area at the three scales, accounting for 65.77% of the total study area. LR is at the average level in the study area, and the remaining five ESs are above the mean value. Bundle 3 is distributed in the relatively flat northern area, accounting for 19.07% of the total study area. Although WR and CP are above average, PR is at the lowest level in the study area. In addition, FM, HM, and LR are the average levels of the study area. Bundle 4 is mainly distributed in the southwest of the study area, accounting for 8.24% of the total study area. The high supply and low demand in the area where bundle 4 is located make all of the six ESs in bundle 4 at the highest level in the study area, especially LR, FM, and HM. Bundle 5 is concentrated in the urban core areas of Suzhou, Wuxi, and Changzhou, accounting for 0.19% of the total study area. The six ESs in this region are all below the average, and CP is at the lowest level in the study area. Bundle 6 is distributed in the study area's northern, northeastern, and central parts, accounting for 4.52% of the total study area. WR in bundle 6 is at a lower level; PR, CP, and HM are higher than the average of the study area; and LR and FM are at the average level of the study area.

Village Scale
The LULC types of bundle 1 are mainly impervious surface (57.7%) and cultivated land (17.6%) (the area proportion of LULC types in ES budget bundles at the village scale was obtained through the zoning statistical tool of ArcGIS 10.8), accounting for 19.93% of the total study area. The HM, FM, LR, and PR of bundle 1 are lower than the average value of the study area; WR is higher than the average value of the study area; CP is the average level of the study area. Bundle 2 is distributed in the central and eastern parts of the study area, accounting for 45.91% of the total study area. The LULC types of bundle 2 are mainly water (44.76%) and cultivated land (30.21%). All ESs are above-average values except for LR, which is slightly below the mean value of the study area. In particular, PR is at the optimal level in the study area. Bundle 3 is located in the west and northwest of the study area, accounting for 21.4% of the total study area. These areas are mainly covered by cultivated land (54.94%) and impervious surfaces (18.21%). The PR of bundle 3 is much lower than the average of the study area; the LR is slightly lower than the average of the study area; and the CP, HM, WR, and FM are slightly higher than the average. Bundle 4 is distributed in the forest (64.5%) in the western and central parts of the study area. The areas of cultivated land (17.22%) and impervious surface (2.28%) in bundle 4 were smaller and less disturbed by human activities. Thus, bundle 4 provides the highest LR, FM, HM, and more PR levels in the study area. However, its area only accounts for 10.19% of the study area. Bundle 5 is mainly concentrated in the urban core areas of Suzhou, Wuxi, and Changzhou, accounting for 1.61% of the study area. CP, HM, FM, and LR are all at the lowest levels in the study area, except for WR, which is slightly above average. Bundle 6

Driver Analysis of ES Budget Bundles
This study analyzed the extent to which different drivers explain the formation of ES budget bundles using the factor detector in Geodetector. Table 6 shows the results of the factor detection at the county, township, and village scales.

County Scale
Bundle 1 is mainly distributed in the north of Wuxi and Changzhou, accounting for 23.3% of the total study area. CP, WR, and LR are slightly higher than the mean value of the study area; PR, HM, and FM are slightly lower than the mean value of the study area. Bundle 2 is distributed in the relatively flat eastern region, concentrated in Suzhou, accounting for 32.44% of the entire study area. Bundle 2 has the highest level of PR and relatively high CP. FM and HM are slightly above the mean values. WR and LR are below the mean. Bundle 3 is distributed in the slightly higher western region, concentrated in the south of Wuxi and Changzhou, accounting for 25.55% of the total study area. All ESs are above the mean value of the study area, except for PR, which is slightly below average. FM and LR are at the highest level in the study area. Bundle 4 is distributed in the middle of the study area, accounting for 0.88% of the total study area. Bundle 4 is concentrated in the urban core areas of Wuxi and Suzhou. Due to the high demand, CP, WR, HM, LR, and FM are all lower than the average value of the study area, except for the PR, which is slightly higher than the average value. Bundle 5 is distributed in the Taihu Lake Scenic Area in the south of the study area, accounting for 16.2% of the total study area. The six ESs are all higher than the mean value, especially WR and HM, which are at the highest level in the study area. Bundle 6 is distributed in the northwest of the study area, located in the core area of Changzhou, accounting for 1.63% of the total study area. The six ESs in bundle 6 are all below the mean value, especially PR, which is at the lowest level in the study area.

Township Scale
Bundle 1 is distributed in the surrounding areas of the urban core areas of Suzhou, Wuxi, and Changzhou, accounting for 2.21% of the total study area. All ESs in bundle 1 have below-average values, except for WR, which is slightly above the mean value of the study area. Bundle 2 has the largest area at the three scales, accounting for 65.77% of the total study area. LR is at the average level in the study area, and the remaining five ESs are above the mean value. Bundle 3 is distributed in the relatively flat northern area, accounting for 19.07% of the total study area. Although WR and CP are above average, PR is at the lowest level in the study area. In addition, FM, HM, and LR are the average levels of the study area. Bundle 4 is mainly distributed in the southwest of the study area, accounting for 8.24% of the total study area. The high supply and low demand in the area where bundle 4 is located make all of the six ESs in bundle 4 at the highest level in the study area, especially LR, FM, and HM. Bundle 5 is concentrated in the urban core areas of Suzhou, Wuxi, and Changzhou, accounting for 0.19% of the total study area. The six ESs in this region are all below the average, and CP is at the lowest level in the study area. Bundle 6 is distributed in the study area's northern, northeastern, and central parts, accounting for 4.52% of the total study area. WR in bundle 6 is at a lower level; PR, CP, and HM are higher than the average of the study area; and LR and FM are at the average level of the study area.

Village Scale
The LULC types of bundle 1 are mainly impervious surface (57.7%) and cultivated land (17.6%) (the area proportion of LULC types in ES budget bundles at the village scale was obtained through the zoning statistical tool of ArcGIS 10.8), accounting for 19.93% of the total study area. The HM, FM, LR, and PR of bundle 1 are lower than the average value of the study area; WR is higher than the average value of the study area; CP is the average level of the study area. Bundle 2 is distributed in the central and eastern parts of the study area, accounting for 45.91% of the total study area. The LULC types of bundle 2 are mainly water (44.76%) and cultivated land (30.21%). All ESs are above-average values except for LR, which is slightly below the mean value of the study area. In particular, PR is at the optimal level in the study area. Bundle 3 is located in the west and northwest of the study area, accounting for 21.4% of the total study area. These areas are mainly covered by cultivated land (54.94%) and impervious surfaces (18.21%). The PR of bundle 3 is much lower than the average of the study area; the LR is slightly lower than the average of the study area; and the CP, HM, WR, and FM are slightly higher than the average. Bundle 4 is distributed in the forest (64.5%) in the western and central parts of the study area. The areas of cultivated land (17.22%) and impervious surface (2.28%) in bundle 4 were smaller and less disturbed by human activities. Thus, bundle 4 provides the highest LR, FM, HM, and more PR levels in the study area. However, its area only accounts for 10.19% of the study area. Bundle 5 is mainly concentrated in the urban core areas of Suzhou, Wuxi, and Changzhou, accounting for 1.61% of the study area. CP, HM, FM, and LR are all at the lowest levels in the study area, except for WR, which is slightly above average. Bundle 6 is distributed in the northeastern part of the study area, accounting for 0.96% of the total study area. LULC types of bundle 6 are mainly cultivated land (33.55%) and industrial land (26.85%). WR in bundle 6 has the lowest level in the study area; PR and CP have above-average values; and FM, HM, and LR have the average level of the study area.

Driver Analysis of ES Budget Bundles
This study analyzed the extent to which different drivers explain the formation of ES budget bundles using the factor detector in Geodetector. Table 6 shows the results of the factor detection at the county, township, and village scales.  The larger the q-value, the stronger the interpretation of X by Y. The q-value means that X explains 100 × q% of Y; the p-value indicates the degree of significance. p < 0.1 indicates significance, and it is indicated by **.

Discussion
This paper analyzed the supply, demand, and budgets of ESs and the drivers of ES budget bundles at the county, township, and village scales in the Su-Xi-Chang region. Unlike previous studies, this study focused on multiscale analysis and analyzed the impact of multiple drivers on ES budget bundles in rapidly urbanizing regions. The issues we attempted to address were threefold. (1) Theoretical implications: to determine the characteristics of the multiscale pattern of ES supply and demand in the Su-Xi-Chang region and to analyze the impact of the natural environment and socioeconomic factors on ES budget bundles. (2) Practical implications: to propose a multiscale decision-making process. (3) To discuss the limitations and future perspectives of this study.

Multiscale Pattern of Supply and Demand of ESs
This paper evaluated the supply and demand patterns of six typical ESs at the county, township, and village scales in the Su-Xi-Chang region and identifies six ES budget bundles. The results showed that the supply, demand, and budgets of ES are scale-dependent, which is consistent with the findings of [25]. On the one hand, the large scale masks the spatial heterogeneity and clustering characteristics of the supply and demand distribution of ESs at the small scale, similar to the findings of [58]. In this study, the distribution of supply and demand of ESs was more refined at the township and village scales than at the county scale. In addition, the demand for CP, WR, FM, and LR, as well as the supply of HM and LR, were randomly distributed at the county scale but clustered at the township and village scales. Therefore, it is unscientific to simply add, generalize, or extrapolate the information represented by a certain scale. On the other hand, the large scale highlights high-value areas of supply and demand. In this study, the high-value areas of CP supply were concentrated in cultivated land. The high-value areas of supply of regulating services (WR, PR, HM, and FM) and cultural service (LR) were concentrated in the forest and water areas. The high-value areas of demand were mainly concentrated in the urban core areas.
Due to the inconsistency in the quantity and location of the supply and demand of various ESs, all ESs at the three scales showed a mismatch between supply and demand. The position and area of surplus, balance, and deficit vary with scale. For example, the supply of HM at the county scale can meet the demand, but the township and village scales cannot ( Figure 5(c1-c3)). On the three scales, CP, PR, and FM all showed a deficit in the urban core areas, unable to meet their own needs through local supply. ES budget bundles also exhibited scale dependence. The spatial patterns and inter-service composition of ES budget bundles at the township and village scales were similar (Figures 7 and 8), while there were significant differences between the township/village scale and the county scale ( Figure 6). This is similar to the findings of a study in Quebec, Canada [25], which demonstrated that the spatial pattern and inter-service composition of ES supply bundles were similar at the two smaller scales (1 km and 3 km) but very different at a larger scale (9 km). This suggests that the formation of ESB at similar scales is similar and stable. When the scale becomes larger or smaller, the spatial pattern and the inter-service composition of ESB will be reconfigured. In summary, we believe that the county scale in the Su-Xi-Chang region is suitable for analyzing the overall situation of ESs, while the township and village scales include detailed information and are suitable for implementing precise management.

The Impact of the Natural Environment and Socioeconomics on ES Budget Bundles
Clarifying the reasons for the formation of ESB will help to deepen the understanding of the coupling between the natural environment and the social economy and inform managers of the underlying mechanisms that affect decision-making [41]. In this study, social factors (POP and IS) determined the formation of ES budget bundles at three scales (Table 6), which is similar to the findings of [42,76]. Economic factors (GDP) only played a significant role at the township scale, which is consistent with the findings of [3]. Natural environment factors (such as DEM, SR, and PRE) played a role at the township and village scales. Our results show that, compared with natural environmental factors, socioeconomic factors play a dominant role in forming ES budget bundles in the Su-Xi-Chang region.
We found that the drivers of ES budget bundles were also scale-dependent ( Table 6). The larger the scale, the stronger the explanation of the driving factors, similar to the results of a study in the Yangtze River Delta region [56]. The larger the scale, the more single the driving factor, and the smaller the scale, the more diverse the driving factor, which is consistent with the research results of [3]. While small scales are better for identifying drivers [80], multiscale analysis of drivers of ESBs in decision-making helps to capture the key at different scales [3].
Therefore, landscape planners and managers must comprehensively consider local natural conditions, typical ES, and scales when making decisions [41]. On the other hand, it is also necessary to focus on regulating human activities [45]. Increasing supply will undoubtedly improve the mismatch between supply and demand of ESs, but it is unrealistic to only increase supply without reducing demand to achieve the sustainable use of ESs [18]. It may be difficult to solve the mismatch between the supply and demand of various ESs at different scales [57], but it will bring about a mismatch or trade-off between the supply and demand of other services [56].

Multiscale Decision-Making Process
We combined multiscale ES budget bundles and drivers to inform decision-making processes toward the sustainable use of ESs. The multiscale decision-making process consists of four steps. The first step is to determine who will govern based on the scale. The regions of Suzhou, Kunshan (the county-level administrative unit under Suzhou), and Kunshan development zone correspond to the management scope of managers at the county scale, the township scale, and the village scale, respectively ( Figure 11, step 1). The second step is to determine where to implement the decision based on the distribution of ES budget bundles ( Figure 11, step 2). The county and village scales include bundle 1, bundle 2, bundle 4, and bundle 5, and the township scale includes bundle 1 and bundle 2. The third step is to determine the ESs that need to be managed based on the histogram ( Figure 11, step 3). The fourth step is to formulate specific management countermeasures based on drivers, local natural conditions, typical ESs, scale, and stakeholder preferences ( Figure 11, step 4) [81].
In the example shown in Figure 11, the area of bundle 2 is the largest at the county scale, so the situation in the area where bundle 2 is located can best represent the overall situation of Suzhou. The PR and CP in bundle 2 are better, but the WR and LR are slightly lower than the mean value of the study area. Therefore, decision-makers at the county scale need to improve the WR and LR of the area where bundle 2 is located. When refining from the county scale to the township scale, Kunshan includes bundle 2, where all ESs are higher than the average level of the study area, and bundle 1, where all ESs except WR are lower than the mean value. Then, the decision-makers at the township scale need to include the region where bundle 1 is located in the decision-making scope and take the improvement of HM and FM as the main task. From the township scale to the village scale, bundle 1 and bundle 5 are the main ES budget bundles of the Kunshan Development Zone. The ESs in bundle 5 have more below-average values than the ESs in bundle 1, especially CP, HM, and FM. Therefore, when the workload is heavy, decision-makers can prioritize bundle 5 at the village scale as the priority area of management. Combined with the driving factors, policy-makers need to focus on regulating human activities and comprehensively consider local natural conditions, typical ES, and scale [45]. At the same time, cultivated land should be used intensively, and the conversion of other lands to urban land should be controlled [3]. Strengthening the cross-regional flow of ES, giving play to the leading role of the advantageous bundle, and complementing the disadvantaged bundle can help achieve a win-win situation [46].
In conclusion, decision-makers at the county scale should pay attention to the overall situation of ES and find the prominent locations of the imbalance between supply and demand. In this way, the decision-making time is shortened, decision-making efficiency is improved, and overall decisions that meet the interests of most people are made. Decisionmakers at the township scale have clearer images and can make decisions that are more in line with local conditions, with decisions that are consistent with the higher-level scale.
Managers at the village scale should incorporate the opinions of stakeholders and help decision-makers at the larger scale fill in the small but important issues that are easily overlooked. The large scale highlights the focus of decision-making, the small scale makes up for the lack of precision in the large scale, and the combination of multiple scales can make ES management more effective [24,37]. The multiscale decision-making process we propose can be broadly applied to making management decisions in other rapidly urbanizing regions of the world. In conclusion, decision-makers at the county scale should pay attention to the overall situation of ES and find the prominent locations of the imbalance between supply and demand. In this way, the decision-making time is shortened, decision-making efficiency is improved, and overall decisions that meet the interests of most people are made. Decision-makers at the township scale have clearer images and can make decisions that are more in line with local conditions, with decisions that are consistent with the higher-level scale. Managers at the village scale should incorporate the opinions of stakeholders and help decision-makers at the larger scale fill in the small but important issues that are easily overlooked. The large scale highlights the focus of decision-making, the small scale makes up for the lack of precision in the large scale, and the combination of multiple scales can make ES management more effective [24,37]. The multiscale decision-making process we propose can be broadly applied to making management decisions in other rapidly urbanizing regions of the world.

Limitation
(1) In terms of supply and demand assessment, from the perspective of demand, for example, the demand for WR is the actual demand calculated according to the statistical Figure 11. Four steps of the multiscale decision-making process (taking the Su-Xi-Chang region as an example).

Limitation
(1) In terms of supply and demand assessment, from the perspective of demand, for example, the demand for WR is the actual demand calculated according to the statistical data. At the same time, the demand for LR was calculated based on the minimum threshold. Therefore, the actual demand for LR may be greater. From the perspective of supply, there is no distinction between actual supply and potential supply. The supply of ES evaluated by the InVEST model is the potential supply, such as FM, which belongs to the potential supply (maximum threshold). (2) It is necessary to map and analyze the supply and demand of ESs in the three administrative spatial units of county, township, and village. This helps to link production to consumption and also meets the decision-making needs of managers [76]. However, the range of influence of ESs is of many shapes and may be larger or smaller than the range of administrative boundaries [36]. In this study, most ESs were generated locally, but WR, for example, was delineated by watershed, and administrative boundaries may have split some ecological boundaries. Therefore, decision-makers need to work with surrounding areas to incorporate the actual management scope of ES into the decision-making process [82]. (3) We did not consider the ES flow, which may alleviate the supply-demand imbalance to a certain extent, which will be the direction of further research in the future.

Conclusions
This paper analyzed the supply, demand, and budgets of ESs and the drivers of ES budget bundles at the county, township, and village scales in the Su-Xi-Chang region. The results showed that the high-value areas of ES supply were mainly distributed in the western part of the study area. The high-value areas of ES demand were mainly concentrated in urban core areas. Due to the inconsistency in the quantity and location of the supply and demand of various ESs, ESs on the three scales all showed a mismatch between supply and demand. Among them, CP, PR, and FM showed deficits in the urban core area and could not meet their own needs through local supply. The location and area of ES supply and demand surpluses, balances, and deficits varied with scale. The spatial patterns of ES supply, ES demand, and ES budget bundles at the township and village scales were similar, while the township/village and county scales were different. We found that socioeconomic factors significantly impact ES budget bundles in the Su-Xi-Chang area more than natural environment factors. The top two drivers on the three scales were population density and the percentage of impervious surfaces. The diversity and degree of interpretation of influencing factors varied with scale. We believe that focusing on the big picture on a large scale, implementing precise management on a small scale, and focusing on regulating human activities can make management decisions more effective. This study can provide a scientific basis for the sustainable utilization of ecosystem services in the Su-Xi-Chang region, and the research results and methods can provide a reference for similar studies in other rapidly urbanizing areas in the world.  Data Availability Statement: All data generated or analyzed during this study are included in the published article.