Spatial scaling of land use/land cover and ecosystem services across urban hierarchical levels: patterns and relationships

Land use/land cover (LULC) patterns seriously affect ecosystem services (ESs), especially in highly developed urban agglomerations. Exploring how LULC and ESs change spatially across urban hierarchical levels and understanding the possible mechanisms can promote the sustainable planning of urban landscapes. By mapping the spatial patterns of LULC and ESs in the three largest urban agglomerations of China, this study aimed to (1) identify the scaling relations of LULC and ESs across different urban hierarchical levels, (2) explore the possible mechanisms of these two types of spatial scaling, and (3) examine how the scaling relations of ESs relate to LULC and the policy implications. Based on LULC, we used the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) model and other biophysical models to quantify ES indicators. Then, scalograms were used to quantify the scaling relations of LULC and ESs with respect to changing spatial extent. Developed land and cropland exhibited the most predictable responses with changing spatial extent. Compared to other ESs, provisioning services were the most predictable. The predictable scaling relations of ESs at different urban hierarchical levels fell into two general types: power laws at the city proper level and exponential relationships at the metropolitan region and urban agglomeration levels. The scaling relations of both LULC and ESs varied across urban hierarchical levels. The spatial scaling of ESs was closely related to LULC patterns. Integrating the scaling relations of ESs into land use planning can help decision-makers formulate multi-scale landscape conservation strategies.


Introduction
For the past few decades, many parts of the world have been undergoing dramatic and rapid urbanization (Stokes and Seto 2019). As the world's most populous country, China has also experienced dramatic urban expansion and socioeconomic development since the national reform and opening-up policies of 1978 Liu et al. 2020). Various studies have used population density (Normile 2016;Poku-Boansi 2021) or built-up land area (Seto et al. 2012;Xu et al. 2018) to represent or measure urbanization. The urbanization population rate in China increased from 17.9 in 1978 to 59.6% in 2018 (Wu et al. 2020a, b), and the urban land area expanded from 21,770 to 74,827 km 2 in the same period (Kuang 2020a). This rapid urbanization has caused dramatic changes in land use/land cover (LULC) compositions and configurations with high spatial heterogeneity (Lawler et al. 2014;Hasan et al. 2020). As the basis of assessing ESs, LULC changes have resulted in myriad impacts on natural resources and further threatened ESs from a local to global scale, such as water scarcity , climate change (Patra et al. 2018), soil erosion , and habitat loss (Swenson and Franklin 2000). Studying how LULC and ESs change spatially in rapidly urbanizing regions is important for better understanding and predicting the patterns and processes of urbanization across multiple scales, and the findings may further benefit urban sustainable planning and ecological protection (Batty 2008;Wu et al. 2011;Zhao et al. 2018).
Studies have explored the spatial characteristics of LULC patterns and ESs (Viglizzo et al. 2012;Hasan et al. 2020). However, we still lack further information about the spatial scale effects of LULC and ESs, and their underlying mechanisms. The hierarchical scaling strategy provides an effective method of investigating spatial heterogeneity of complex urban systems over a range of scales (O'Neill et al. 1986;Wu 1999;Wu and David 2002). First, urban systems exhibit a hierarchically structured system wherein large urban regions consist of individual cities, which in turn consist of a smaller city proper (Ma et al. 2019). Thus, the hierarchical perspective can systematically characterize the impacts of landscape patterns on ecological processes in urban regions (Ma et al. 2016b;Bian et al. 2020;Zhang et al. 2020). Second, the scaling approach is a powerful tool that can quantify multiscale characteristics explicitly and describe how various urban attributes scale with city size (Bettencourt 2013;Lobo et al. 2020).
Until now, the scaling approach has predominantly been used to focus on two themes: The first theme is how urban metrics (e.g., physical, demographic, economic, or environmental attributes) change with city sizes (usually represented by the population size or urban area) (Brock 1999;Fuller and Gaston 2009;Zhao et al. 2018). Studies illustrated that most urban attributes follow approximate power-law scaling, with different scaling exponents (e.g., [ 1 for certain quantities reflecting wealth creation and innovation and \ 1 for certain material infrastructural quantities) (Bettencourt et al. 2007(Bettencourt et al. , 2010Ribeiro et al. 2020). The second theme is how landscape patterns (e.g., landscape shape) or ecological processes (e.g., the relations between impervious surface area and land surface temperature) change with spatial scale (e.g., grain size or extent). For example, Wu (2004) systematically investigated the responses of landscape metrics to changing scale and classified the responses into two categories: simple scaling functions and unpredictable behavior. In general, most of the existing studies have focused on landscape patterns (Kedron et al. 2018;Wu et al. 2004;Zhang et al. 2009) but neglected to explore the scaling relations of ecological processes and ES indicators. In the studies on landscape patterns, researchers rarely investigate the specific scaling relations for different LULC types. Moreover, it is still unclear how the scaling relations of different LULC types and multiple ESs change in rapidly urbanized regions with a hierarchical structure.
Urban agglomeration is a highly developed spatial form of integrated cities and exhibits a hierarchically structured system (Ma et al. 2019;Wang et al. 2019b). As urban complex systems, urban agglomerations suffer the most serious changes in terms of LULC and ESs (Zhou et al. 2018;Yu et al. 2019;Shen et al. 2020). Thus, we selected the three largest urban agglomerations in China, the Beijing-Tianjin-Hebei (BTH), the Yangtze River Delta (YRD), and the Pearl River Delta (PRD) areas (Du et al. 2018) to explore general ''rules of thumb'' for scaling information of LULC types and ESs across urban hierarchical levels. The analysis linked spatial scaling with complex socio-ecological urban systems and helped us understand the possible processes and mechanisms behind the scaling functions. The knowledge of the spatial scaling of LULC and ESs can be translated into sustainable landscape planning and policy practices, thus providing references for similar highly developed urban regions globally.
This study used the hierarchical scaling strategy to investigate how different LULC types and ESs change with the spatial extent of the three largest urban agglomerations in China. We aimed to achieve the following research objectives: (1) analyze the scaling relations of LULC and ESs across different urban hierarchical levels; (2) understand the underlying mechanisms behind the scaling relations of LULC and ESs; and (3) explore how the scaling relations of ESs respond to LULC and its policy implications.

Study area
The three largest urban agglomerations, BTH, YRD, and PRD ( Fig. 1), cover only 5% of Chinese territory (DUSNB 2018) but generate 40% of the total national urban impervious surface (Ma et al. 2019). They are the primary carriers of China's socioeconomic development and new urbanization processes (Fang and Yu 2017). In 2018, BTH, YRD, and PRD accounted for 8.1%, 11.0%, and 5.1% of the total national urban population and 9.5%, 16.7%, and 9.2% of the total gross domestic product (GDP), respectively (DUSNB 2018). In addition to their importance to China's economic and social development, these three urban agglomerations constitute the typical and core regions of prominent environmental issues in China (Zhang et al. 2017;Gao et al. 2019).
Administrative urban hierarchy is the most common and actual ecosystem management boundary for regional planning and policy implementation (Forgione et al. 2016;Nikodinoska et al. 2018;Xing et al. 2021). BTH, YRD, and PRD consist of mature and complete nested urban hierarchical levels that extend from the city proper, metropolitan region to the broader urban agglomeration (Ma et al. 2019). The lowest city proper level (the city proper of Beijing, Shanghai, or Guangzhou) belongs to a metropolitan region (Beijing, Shanghai, or Guangzhou metropolitan region), and is also embedded in the highest urban agglomeration level (BTH, YRD, or PRD urban agglomeration) (Fig. 1).

Data acquisition
For the three urban agglomerations, LULC data from 2018 with a spatial resolution of 30 9 30 m were downloaded from the Resource and Environment Science and Data Center (http://www.resdc.cn/). Digital elevation model data with a spatial resolution of 30 9 30 m were derived from the Geospatial Data Cloud (http://www.gscloud.cn/). Climate data, including the annual average precipitation and temperature in all meteorological stations, were collected from the China Meteorological Data Service Center (http://data.cma.cn/). Potential evapotranspiration data were collected from the Global Aridity and PET Database (https://cgiarcsi.community/data/globalaridity-and-pet-database/). Soil properties data, including bulk density, soil organic carbon, clay content, sand content, silt content, and soil depth, were obtained from the World Soil Information database (https://soilgrids.org/).

LULC maps in urban agglomerations
LULC maps in 2018 for the BTH, YRD, and PRD urban agglomerations were generated by researchers at the Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, through manual visual interpretation. The original data were Landsat 8 OLI_TIRS remote-sensing images that were downloaded from the United States Geological Survey (https://www.usgs.gov/). Sampling surveys in fields and high-resolution Google Earth images were used to verify the classification accuracy. The assessment results indicated reasonable classifications, with an overall classification accuracy of 95% and a kappa coefficient of 0.92. The raster data comprised six LULC types: (1) cropland, including irrigated cropland, dry cropland, and orchards; (2) forests, including natural forests, artificial forests, coppice, shrub forests, and sparse forests; (3) grassland, including all types of grasslands with coverage of more than 5%, pastures, and hay; (4) water bodies, including rivers, lakes, reservoirs, ponds, swamps, tidal lands, and wetlands; (5) developed land, including urban construction land, rural residential land, industrial land, traffic road, airport, and mining area; and (6) barren land, including sandy land, desert, saline land, bare land, bare rock, and other unutilized lands.

Quantifying multiple ES indicators in urban agglomerations
This study selected eight ES indicators: carbon storage, food production, water yield, air pollution removal, nitrogen retention, soil retention, habitat quality, and recreational opportunity. The selection was based on three criteria: (1) the indicators belong to the basic provisioning, regulating, and cultural services (Millennium Ecosystem Assessment 2005); (2) they have great significance to the sustainable development of three urban agglomerations (Lan et al. 2017;Zhang et al. 2017;Liu et al. 2018;Sun et al. 2018;Luo et al. 2020); and (3) they represent the information of water resources, food, climate, soil, habitat, and culture in urban agglomerations, and thus can comprehensively reflect the characteristics of complex urban ecosystems.
We used the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) model to quantify the carbon storage, water yield, nitrogen retention, soil retention, and habitat quality for the three urban agglomerations. Food production was estimated by using the Global Agro-Ecological Zones (GAEZ) model. The removal of PM 2.5 (atmospheric aerosol particles with a diameter less than 2.5 lm) was estimated on the basis of the adsorption capacity of vegetation to PM 2.5 pollutants. Recreational The concentric circle illustrations on the right represent the spatial extension of three urban agglomerations at the administrative hierarchy levels. The CP of Beijing/Shanghai/ Guangzhou refer to the city proper of Beijing/Shanghai/ Guangzhou; Beijing/Shanghai/Guangzhou refers to the Beijing/Shanghai/Guangzhou metropolitan region; BTH/YRD/ PRD refers to the Beijing-Tianjin-Hebei/Yangtze River Delta/ Pearl River Delta urban agglomeration opportunity was quantified by the Recreation opportunity Spectrum (ROS) model. Table 1 lists the ES indicators and assessment descriptions. In addition, the calculation process and detailed parameters are shown in Supplementary Information (SI).
Quantifying the scaling relations of LULC and ESs with respect to spatial extent ES supply depends on ecosystem structures and processes. Thus, LULC change has become a major driving force that impacts regional ecosystem change (Hasan et al. 2020). Both LULC compositions and configuration are expected to significantly impact ESs (Lei et al. 2021). Therefore, exploring how LULC change responds to an increased spatial extent is an important basis for understanding and explaining the spatial scaling effect of ESs. The form of scalograms (Wu 2004;Frazier 2016;Ma et al. 2019) was adopted in this study to reflect how different LULC types and ESs respond to changing spatial extents in the three urban agglomerations. A series of concentric circles with gradually expanding radii distance (Unit: km) can be used to reflect the changes of LULC and the characteristics of landscape patterns in different buffer zones Kuang et al. 2014;Ma et al. 2019) across urban hierarchical levels very well ( Fig. 1, Table 2). Thus, scalograms were constructed by plotting the proportion of LULC and the average ESs with respect to increasing spatial extents (Wu 2004).
Using SPSS 25.0, we adopted the curve estimation regression models (Jomnonkwao et al. 2020) to investigate the scaling relations. Compared with other statistical methods, the regression model is more suitable to express the curve fitting relationship between the spatial extent and quantitative indicator (Wu et al. 2004;Xu et al. 2020). It is very useful when predicting unknown scaling relations according to independent variables (Ma et al. 2019). In this study, the independent variable (X) was the radius of the concentric circle (i.e., spatial extent). The dependent variable (Y) was the proportion of LULC or the average ES value within the corresponding concentric circle area. We performed a check of the stationarity of data using the EViews 10.0 software (Liu et al. 2016) and checked the independence of residuals using SPSS 27.0 (Li et al. 2012). If the independent and dependent variables pass the stationarity test (P \ 0.05) and the residuals are independent, the regression method can be adopted to reflect their relationships. The P value and coefficient of determination (R 2 ) were used to indicate the goodness of fit of the curve (Massada and Radeloff, 2010). The regression model with P \ 0.001 and the largest R 2 was selected as the optimal one (Alhamad et al. 2011). Food production The food production potential is calculated by considering the limiting factors, such as water resource, soil properties, topographic conditions, cultivated land distribution, and management measures (Liu et al. 2015) Water yield The InVEST Water Yield model estimates the water of different landscape areas. It is based on the principle of water balance on a grid map (Sharp et al. 2020) PM 2.5 removal The PM 2.5 removal is calculated based on the adsorption of PM 2.5 per unit area for different LULC types

Nitrogen retention
The InVEST Nutrient Delivery Ratio model assesses nitrogen retention according to the nutrient pollutants removed by vegetation and soil in surface runoff (Sharp et al. 2020) Soil retention The InVEST Sediment Delivery model computes sediment retention based on the amount of soil loss and its delivery ratio (Hamel et al. 2015) Habitat quality The InVEST Habitat Quality model combines the information of LULC types and threats to biodiversity to produce habitat quality maps (Sharp et al. 2020) Recreational opportunity The ROS model is associated with the recreational index and accessibility index for reaching recreational locations (Lavorel et al. 2020)

Scaling relations of LULC for three urban agglomerations
For the three urban agglomerations, the proportion of developed land decreased as the spatial extent of the radius increased. In contrast, the proportions of forests and cropland increased with the expansion of the radius. The proportion of water bodies showed irregular fluctuations with the expansion of scale. In addition, the proportion of grassland in BTH increased when the radius expanded, while the proportions of grassland and barren land were exceedingly small in YRD and PRD (Figs. 2, 3). Overall, the proportions of cropland, forest, and developed land exhibited appropriate and predictable scaling laws at certain levels of urban hierarchy. The typical scaling curves for different urban hierarchical levels were classified into three rules: linear, sigmoidal, and exponential relationships. Among them, the exponential relationships were divided into exponential growth and exponential decay (Fig. A1). However, the scaling relations of water bodies and barren land were unpredictable at nearly all three urban hierarchical levels. The proportion of grassland showed predictable scaling relationships in BTH but was unpredictable in YRD and PRD (Table 3).
For the BTH, linear scaling relations were found at the city proper level for most of the land use types. Next, the general relationships changed to exponential and sigmoidal functions at the broader metropolitan and urban agglomeration levels. For YRD, most of the LULC types did not reveal predictable scaling relations and functions except for cropland and developed land. Both showed linear relationships at the city proper level and showed sigmoidal and exponential The center of concentric circles is the administrative center of a city proper. For each urban hierarchical level, the concentric circle with the maximum radius can cover the boundary of the corresponding level. BTH/YRD/PRD represents the Beijing-Tianjin-Hebei/ Yangtze River Delta/Pearl River Delta urban agglomeration, CP of Beijing/Shanghai/Guangzhou refers to the city proper of Beijing/ Shanghai/Guangzhou  Table 2 growth or exponential decay at higher urban hierarchical levels, respectively. For PRD, the cropland, forest, and developed land scaled in a predictable way. The proportion of cropland increased as exponential curves, while forest increased as sigmoidal curves. The proportion of developed land decreased in a linear relationship at the city proper level, while the relationship switched to exponential decay at the higher levels (Table 3).

Scaling relations of multiple ESs in BTH
For the spatial patterns of multiple ESs in BTH, the values for most of ESs increase with increasing radii extent, except for water yield and soil retention (Fig. 4). According to the scalograms of multiple ESs, with respect to increasing spatial extent in BTH, carbon storage, food production, water yield, PM 2.5 removal, nitrogen retention, and habitat quality However, soil retention did not show any simple and predictable scaling relations (Fig. 5f). In general, the two form functions, which were power law and exponential functions, fit the scaling relations of ESs better in the BTH agglomeration. The ESs showed pervasive power-law scaling relations at the city proper level while exhibiting exponential scaling relations at the higher metropolitan and urban agglomeration levels (Fig. 5). Specifically, the values of carbon storage, food production, PM 2.5 removal, and habitat quality increased consistently and followed a power-law function with a scaling exponent greater than 1 within the city proper of Beijing. As the spatial extent increased in the metropolitan region and urban agglomeration, the indicator values showed exponential functions with negative exponents (Fig. 5a, b, d, and g). For nitrogen retention, the scaling exponent was smaller than 1 in the power-law function within the city proper of Beijing. As the spatial extent of the radii increased beyond the city proper, the nitrogen retention values also showed exponential decay with a negative exponent and continued to increase until it eventually reaches an asymptote (Fig. 5e). In addition to these indicators, the water yield showed a decreasing trend with the expansion of spatial extent of the radii. The index values followed a power-law decline within the city proper of Beijing. They continued to decrease with exponential scaling relations at higher urban hierarchical levels (Fig. 5c). For recreational opportunities, the index values first decreased and then increased with the increasing spatial extent of the radii within the city proper. As the spatial extent increased, the index values followed exponential functions and continued to grow until they stabilized (Fig. 5h). Barren land --- The archetypes of scaling relations were shown in Fig. A1. For all scaling relation equations in Table 3, the relationships between the proportions of land use/land cover (LULC) types and the increasing spatial extents were significant (**P \ 0.001, R 2 [ 0.990). The dash in Table 3 means the scaling relation of the LULC type was unpredictable at the corresponding urban hierarchical level

Scaling relations of multiple ESs in YRD
For the spatial patterns of multiple ESs in YRD, the values for most of ESs become higher with increasing spatial extent, except for water yield, nitrogen retention, and soil retention (Fig. A2). According to the scalograms of multiple ESs of the increasing spatial extent in YRD, only food production and water yield displayed universal and predictable scaling relations at different urban hierarchical levels ( Fig. 6b and c). For carbon storage, PM 2.5 removal, and habitat quality, the index values followed predictable scaling relations within the Shanghai metropolitan region. As the spatial extent of the radii increased beyond the metropolitan region, the scaling relations became unpredictable (Fig. 6a, d and g). In addition, for nitrogen retention, soil retention, and recreational opportunities, there were no predictable scaling relations that could fit the curves (Fig. 6e, f and h). Specifically, food production values increased with a power-law scaling relationship at the city proper level, while shifting to exponential growth scaling relationships at the higher urban hierarchical levels (Fig. 6b). Water yield values decreased with a powerlaw scaling relation at the city proper level, followed by exponential decay scaling as the spatial extent of the radii covered the metropolitan region and urban agglomeration (Fig. 6c). For carbon storage, PM 2.5 removal, and habitat quality, the index values exhibited power-law scaling relations within the city proper of Shanghai. As the spatial extent of the radii further increased to the Shanghai metropolitan region, carbon storage still followed a power-law function, while PM 2.5 removal and habitat quality followed an exponential growth relationship. The values of these three ESs showed upward staircase-like curves when the spatial extent increased to urban agglomeration ( Fig. 6a, d, and g). In addition, the values of nitrogen retention did not change significantly at various urban hierarchical levels (Fig. 6e). The soil retention values become higher when spatial extent increased beyond the metropolitan region (Fig. 6f). The recreational opportunity index values showed an upward staircaselike curve with increasing spatial extent (Fig. 6h).

Scaling relations of multiple ESs in PRD
For the spatial patterns of multiple ESs in PRD, the values for most of ESs increase with increasing spatial extent of the radii, except for water yield and soil retention (Fig. A3). According to the scalograms of multiple ESs, with respect to increasing spatial extent of the radii in PRD, carbon storage, food production, water yield, PM 2.5 removal, nitrogen retention, and habitat quality all exhibited predictable scaling relationships at different urban hierarchical levels ( Fig. 7a-e and g). As with the BTH agglomeration, the scaling relationships of soil retention in PRD were unpredictable (Fig. 7f). Recreational opportunity showed a predictable scaling relation only at the urban agglomeration level (Fig. 7h). In general, the scaling relationship also presented two types of functions, which were power law and exponential relationships. Moreover, exponential relationships predominantly existed at higher urban hierarchical levels (Fig. 7).
Carbon storage and PM 2.5 removal values followed power-law functions within the Guangzhou metropolitan region. Thereafter, they showed exponential functions as the spatial extent of the radii further increased to the urban agglomeration ( Fig. 7a and d). The values of food production, nitrogen retention, and habitat quality followed power-law functions within the city proper and then increased exponentially as the spatial extent of the radii increased beyond the city proper (Fig. 7b, e, and g). For water yield, the index values showed a monotonously decreasing trend as the spatial extent increased. The index values followed the power-law functions within the Guangzhou metropolitan region and then showed the exponential function as the spatial extent of the radii further increased (Fig. 7c). For recreational opportunity, the index values began to follow an exponential function only when the spatial extent increased to the urban agglomeration level (Fig. 7h).

Discussion
How did the scaling relations of LULC change across different urban hierarchical levels and agglomerations? Different LULC types showed diversified distribution trends when the spatial extent of the radii increased. The developed lands exhibited decreasing trends from city proper outwards while croplands and forests showed an increasing trend. This condition was mainly because croplands and forests were substantially converted to developed land during the rapid process of urbanization in China, especially for large urban agglomerations (Liu et al. 2019;Kuang 2020b;Zhou et al. 2021). In contrast, water bodies and grasslands showed irregular and random characteristics with the increase of radii spatial extent. This phenomenon was predominantly related to natural geographical conditions (e.g., topography, climate, etc.) and regional planning (e.g., park layouts, the demarcation of nature reserves) (Yang et al. 2019b;Sun et al. 2020). In addition, barren lands had no c Fig. 6 Scalograms of multiple ecosystem service indicators with respect to increasing concentric circle radii at various urban hierarchical levels: the city proper of Shanghai, Shanghai metropolitan region, and the Yangtze River Delta (YRD). (a) Carbon storage; (b) Food production; (c) Water yield; (d) PM 2.5 removal; (e) Nitrogen retention; (f) Soil retention; (g) Habitat quality index; (h) Recreational opportunity index. The annotations are labeled on the Y-axis. For all regression curves and equations, the significance level of P-values \ 0.001 obvious predictable scaling relations, as their area within the three highly urbanized agglomerations was exceedingly small. Among all LULC types, developed land and cropland were the most predictable at three urban hierarchical levels (Table 3). Developed land and cropland are typically the most severely affected by human activities, especially in urban centers (D'Amour et al. 2016). Socioeconomic factors such as population agglomeration and gross domestic product growth predominantly affect the spatial pattern of developed land and cropland (Arcaute et al. 2015;Li et al. 2017;Gu 2019). The scaling relations of socioeconomic factors can be well approximated by common scaling laws. For example, since the population density and the Gross Domestic Product followed predictable functions in the three urban agglomerations (Fig. A4), the two LULC types were also predictable across radii extent. In a prior study, six world metropolitan regions with high urbanization levels, including Phoenix and Baltimore in the US; Santiago, Chile; London, England; Paris, France; and Berlin, Germany, exhibited general urban impervious surfaces scaling relations that were quite similar to the three Chinese metropolitan regions tested here: Beijing, Shanghai, and Guangzhou (Ma et al. 2019). Thus, the predictable scaling relationships identified in this study may provide information for the extrapolation of developed lands at certain urban hierarchical levels in similarly high developed megacities.
Different LULC types presented various scaling relationships at the different urban hierarchical levels, of which scaling was most predictable at the city proper level. For example, water bodies and barren land showed predictable scaling functions at the city proper level but were unpredictable at the broader metropolitan and urban agglomeration levels (Table 3). This situation was predominantly because the influencing factors of LULC distributions became more diverse with the gradual expansion of spatial amplitude (i.e., the addition of more natural, geological, and policy factors, such as topography, climate, and regional planning) (Verburg et al. 2003;Du et al. 2014). It is presumed that these factors made the scaling relationships more complex and increased the unpredictability. At the city proper level, linear relationships (i.e., power-law with scaling exponent = 1) existed for most of the LULC types. The results were relatively consistent with those of previous studies showing that developed lands followed power-law scaling relationships with increasing extent in highly urbanized areas (Ma et al. 2019). In urban centers, human production and living activities severely disturb the land spatial composition. Studies have proven that high population densities influenced the spatial patterns of urban developed lands within the city proper of the three regions (Ma et al. 2016a(Ma et al. , 2019. As the spatial extent of the radii increased to cover the broader metropolitan and urban hierarchical levels, population densities became smaller and more stable (Fig. A4). Thus, human disturbances gradually weakened beyond the suburbs Lan et al. 2021). Other complicated natural or policy factors are likely the dominant drivers at broader scales, which leads to the unpredictability of LULC metrics ( Fig. 3; Table 3).
Different urban agglomerations also presented various scaling relationships for LULC, among which the relations in BTH were the most predictable and those in YRD were the least predictable. This condition was mainly related to the following properties of each urban agglomeration: (1) the constraints of natural and topography conditions; (2) the urban expansion patterns and development modes; (3) the socioeconomic planning and ecological policies. BTH was a simple combined urban agglomeration, with Beijing as the core metropolitan area and Tianjin as the largest sub-core city (Fig. 2) (Fang and Yu 2017). The urban landscape pattern in the Beijing metropolitan region exhibited typical concentric ring structure characteristics along with the spatial distribution of the ring-roads (Wu et al. 2020a, b). Therefore, the concentric ring structure of LULC patterns in BTH were consistent with the concentric circles of expanding radii distance in the analysis, and thus made the scaling relations more predictable. However, the spatial patterns of different LULC presented complex and unpredictable scaling relations in the YRD. The c Fig. 7 Scalograms of multiple ecosystem service indicators with respect to increasing concentric circle radii at various urban hierarchical levels: the city proper of Guangzhou, Guangzhou metropolitan region, and the Pearl River Delta (PRD). (a) Carbon storage; (b) Food production; (c) Water yield; (d) PM 2.5 removal; (e) Nitrogen retention; (f) Soil retention; (g) Habitat quality index; (h) Recreational opportunity index. The annotations are labeled on the Y-axis. For all regression curves and equations, the significance level P-values \ 0.001 hypothesized reasons for this unpredictability are as follows: (1) the city proper of Shanghai can only expand westward due to topography constraints (Fig. 1); (2) the scattered sprawl and polycentric development in the rapid development stage of urbanization created a spatial polycentric mega-city region with several urban cores Tao et al. 2020), which led to the fragmentation of forest, croplands, and other landscape patches ; (3) the ''Development Planning for Yangtze River Delta Region'' in 2016 (www.gov.cn/) approved 26 prefecture-level cities as core cities in the YRD. The disparities in infrastructure construction, industrial planning, and eco-environmental protection projects among these core cities (Zhu et al. 2016) aggravated the unpredictability of LULC for the whole urban agglomeration area.
How did the scaling relations of multiple ESs change and respond to LULC across different urban hierarchical levels and agglomerations?
Most ES values increased with the increasing spatial extent of the radii, except for water yield, nitrogen retention, and soil retention in the three urban agglomerations. The average values of carbon storage, PM 2.5 removal, habitat quality, and recreational opportunity were lower in urban centers and became higher in broader metropolitan and urban agglomeration regions, due to the increase of ecological lands, such as forests, wetlands, and grasslands (Fig. 2). Studies have shown that ecological land per unit can provide more regulation and cultural service values, compared to artificial land (Luo et al. 2018;Baumeister et al. 2020). In contrast, the values of water yield were highest in urban centers. This situation was mainly because the reference evapotranspiration in urban areas was less than in suburban vegetated areas (Yang et al. 2019a;Benra et al. 2021). It is worth noting that soil retention changed irregularly with an increase in radii spatial extent. The phenomenon was likely largely determined by comprehensive factors, such as topography (Sun et al. 2014), rainfall intensity (Rodríguez-Caballero et al. 2013), and vegetation distribution (Korkanç and Dorum 2019).
For different types of ESs, provisioning services were most predictable. Food production and water yield showed predictable scaling relations at all hierarchical levels and urban agglomerations. For example, food production revealed power law and exponential growth relationships at all three urban agglomerations. Landscape structure (e.g., composition and configuration) is often considered to be related to ecological processes and ESs (Leitão and Ahern 2002). To be specific, the scaling relations of food production in this study might be associated with the predictable distribution of croplands. Both food production and croplands increased rapidly on smaller scales (power-law functions) and became stable at larger urban agglomeration scales (exponential functions). In addition, some types of regulating ESs only exhibited predictable functions at certain urban levels, such as for carbon storage and PM 2.5 removal within YRD (Fig. 6). When the spatial extent of the radii extended to the broader urban hierarchical levels, the scaling relationships of these two indicators transformed into complex staircase-like patterns. A previous study has shown that the scalograms of urban impervious surfaces exhibited scale breaks (change points) that corresponded roughly to the urban administrative levels (Wu and Li 2006). Similar to urban impervious surfaces, the breakpoint changes also existed in croplands, forests, and developed lands (Fig. 3f), which might cause the staircase-like patterns of carbon storage and PM 2.5 removal within YRD. In addition to the composition of LULC, the configuration of LULC also appears to be related to ESs across scales (Duarte et al. 2018). For example, the scaling relations of the landscape shape index (Fig. 8) were similar to those of habitat quality (Figs. 5g, 6g, and 7g) in three urban agglomerations. As the spatial extent of the radii increased, the landscape shape index for cropland ( Fig. 8a-c), forest ( Fig. 8d-f), and developed land (Fig. 8g-i) increased. More complex landscape shape index values were related to ecological processes such as material cycling, energy exchange, and species dispersal, and thus improved habitat quality and biodiversity (Dramstad et al. 1996;Karimi et al. 2021).
Two types of predictable scaling functions for ESs were observed at the different urban hierarchical levels, among which the power-laws predominantly expressed ES patterns at smaller levels while exponential relationships were more suitable for larger levels. Although the power-law function was found to be ubiquitous when representing the scaling relationship between ecological attributes and measurable scales (Newman 2005;Fisher et al. 2008;Spence 2009), it did not capture the scaling relations across all LULC and ESs at different urban hierarchical levels. Zhao and Liu (2014) found that a similar power-law described the relationship between critical scale resolution of the carbon cycle and spatial extent. They proposed that future studies should further investigate the compatibility of power-law with the scaling relations of ecological indicators. This study fills the research gap in the scaling of ESs and demonstrates that power-law scaling can be fitted well only at certain ranges or urban hierarchical levels. At the city proper level, landscape structures are often shaped by anthropogenic factors related to developed lands, such as demographic, economic, and traffic (Li et al. 2013;Xie et al. 2017;Ma et al. 2019). Consequently, most of the ESs revealed similar predictable power-law functions that were consistent with developed lands. As the spatial extent of the radii increased to capture the broader metropolitan or urban agglomeration regions, cropland, forest, grasslands, and water bodies accounted for larger proportions and were changed into dominant LULC types (Fig. 3). Natural factors, such as topography and climate conditions, which are macroscopic and stable, became the main driving forces in shaping landscape patterns and processes Smith et al. 2019). Thus, the ES patterns became stable and exhibited significant exponential relationships on a broader scale.
Different urban agglomerations also presented various predictable scaling relationships for ESs, Fig. 8 Scalograms of landscape shape index for different land use types with respect to increasing concentric circle radii at three urban agglomerations among which scaling relationships in the BTH were the most quantitatively predictable while scaling relationships in the YRD were unpredictable. The study showed that the more complex the spatial configuration of LULC, the greater the variability of ES patterns Yee et al. 2021) and the poorer the predictability of scaling relations. Therefore, the predictability of spatial scaling of ESs may be closely related to the complexity of LULC patterns. Compared with YRD and PRD, the shape of the BTH is more regular, and the LULC is distributed more evenly in all directions (Fig. 1). The urbanization of BTH was dominated by the ''urban edge-expansion mode (the newly developed patches expanded from the fringes of existing urban centers)'' (Yu and Zhou 2017), which may have resulted in the scaling relations for the ESs being more predictable. However, the scaling relations of ESs in YRD were difficult to predict for two main reasons: (1) the topographical constraints and polycentric urban sprawl in the rapid development stage of urbanization created a more complex urban shape and increased landscape fragmentation (Li and Zhou 2019;Tao et al. 2020), and thus intensified the spatial heterogeneity of ESs. (2) Large proportions of ESs were provided by forest and water bodies in the YRD; however, the scaling relationships of these two LULC types were unpredictable (Table 3), which may have also affected the predictability of spatial scaling of ESs.
How to integrate the scaling relations of ESs into land use planning and policy practice?
The scaling relations of ESs could reflect the actual spatial characteristics of the ecological environment at each urban hierarchical level of the urban agglomerations. We used a piecewise function to simulate the scaling curve of ESs, which can reflect the critical thresholds (change points) of various ES indicators (Momblanch et al. 2016) that were the turning points of the change rates of scaling relations for ESs (Fig. 9a). The range of change points represents the actual geographical boundary of the ES changes, which was not consistent with the urban administrative boundary. Fig. 9b and c illustrate the distance (variance) between the change point of ES scaling relations and the administrative boundary for multiple ESs in the three urban agglomerations. In the city proper, the change point values for most of the ESs were larger than the boundary value of the city proper of Beijing, while the change point values for all ESs were less than the boundary value of the city proper of Guangzhou. This suggests that the eco-environmental effects of urbanization of Beijing city proper go beyond the scale of administrative planning, and hence, the actual management extent of ESs should be broader than the administrative boundary. In contrast, the eco-environmental effects of urbanization for Guangzhou city proper did not extend to the administrative boundary, and the actual impact extent of ESs was smaller than the administrative boundary. At the metropolitan level, the actual impact extent of most of the ESs was broader than the administrative boundary in the Shanghai metropolitan region, while the actual ES boundaries were smaller in the Beijing and Guangzhou metropolitan regions. Therefore, the Shanghai metropolitan region should pay more attention to the impacts of urbanization on ESs. To effectively protect the ecological environment and natural resources, policy makers should incorporate the actual extent of ES management into sustainable land use planning, rather than merely relying on administrative boundaries (Hein et al. 2006;Brunet et al. 2018). The proposed ES boundary changes can provide a theoretical basis for the boundary determination of the actual ES management and landscape planning in urban agglomerations.
This study confirmed that the scaling of ESs was closely related to land use patterns in the three urban agglomerations. Thus, the supply of multiple ESs can be effectively improved by managing the spatial patterns of LULC (Mitchell et al. 2015). At the city proper level, the eco-environmental effects of urbanization in the city proper of Beijing had obviously gone beyond the administrative scale (Fig. 9b). This type of human-dominated urban center (Robinson et al. 2009) was the most vulnerable and the worstperforming area for multiple ESs (Fig. 5). It thus requires prioritization for land use optimization. At the city proper level, LULC patterns were more predictable (Table 3), and thus, landscape planning can promote ESs more predictably. However, urbanization expansion and LULC spatial patterns emerged from a myriad of policy decisions and planning processes (Batty 2008). To conserve ESs, more compact urban expansion forms should be encouraged, such as improving land use efficiency to maintain rational densification (Shoemaker et al. 2019;Wang et al. 2019a). For the city proper of Beijing, policymakers should implement land use regulations or strategies, such as land-use system reform and a new type of urbanization to reconcile the conflict between developed land and cropland . At the metropolitan level, the eco-environmental effects of urbanization in the Shanghai metropolitan region exceeded the administrative boundary (Fig. 9c). Regulating services, including PM 2.5 removal and soil retention performed poorly in the metropolitan region (Fig. 6). Conserving natural land cover, such as by setting an ''ecological red line'' to protect those green and blue spaces that are affected by urbanization most seriously (Lin and Li 2019) can effectively improve the regulating services (Wang et al. 2020). In addition, since the scaling relations of ESs were similar to landscape configuration metrics, decision-makers can adjust landscape structures (Rieb and Bennett 2020;Lei et al. 2021) to influence the ecological process, support landscape multifunctionality (Bolliger et al. 2011), and improve multiple regulating and cultural ESs (Cong et al. 2016;Lamy et al. 2016).

Limitations and future perspectives
Although the InVEST models contain uncertainties and do not fully reflect the actual indicator values (Pickard et al. 2017;Daneshi et al. 2021), previous Fig. 9 a Example of a scalogram for ecosystem service (ES) change points versus the administrative boundary points. The solid black dots represent the turning points for the ES change rates; the solid red dots represent the administrative boundary points. The ''d cp '' represents the distance between the first change point of ES scaling relations and the administrative boundary of the city proper; the ''d mr '' represents the distance between the second change point of ES scaling relations and the administrative boundary of the metropolitan region. The distance was calculated by the change point value minus the administrative boundary value. b The ''d cp '' for multiple ESs in Beijing-Tianjin-Hebei (BTH), Yangtze River Delta (YRD), and Pearl River Delta (PRD) urban agglomerations. CS = carbon storage; FP = food production; WY = water yield; PR = PM 2.5 removal; NR = nitrogen retention; SR = soil retention; HQ = habitat quality; RS = recreational opportunity. c The ''d mr '' for multiple ESs in BTH, YRD, and PRD urban agglomerations studies have confirmed that they performed well in certain regions (Redhead et al. 2016;Sun et al. 2018), which included urbanized areas similar to the ones in the current study. Thus, the results of this study in the three urban agglomerations are expected to be relatively reliable. The models also have been applied at multiple scales and have proven to be effective (Nelson et al. 2009;Grafius et al. 2016;Bagstad et al. 2018). Therefore, their use was justified in the current study to explore the scaling relations of ESs.
Limitations also existed in the design of the scalograms and statistical techniques for identifying the scaling relationships of LULC and ESs. In this study, the scalograms were created using the radii of circles rather than the true extent of land area analyzed. In order to test whether using the actual land cover extents might change the results, this study made a preliminary attempt with six typical indicators in three urban agglomerations. Results showed that when using the true extent of land area as an indicator for x-axis (Fig. A5), the trends of regression curves and the types of scaling functions for corresponding ESs were the same as using the radii of circles (Figs. 5, 6, and 7). Therefore, although uncertainties existed in the scalograms plotting, the results of this study are still reliable and can provide rules and references for the scaling relations of ESs. It is worth mentioning that the regression techniques we used were unable to avoid the autocorrelation issues that result because each successive radii contained the data information from the radii before it, which means each concentric circle is dependent on the prior circle to some extent. The corresponding planning and policy implications then had certain limitations as they were put forward based on the existing methods. Future research should increase the sampled database of urban agglomerations, test alternative statistical models that may be able to reduce the dependence, and further validate the scaling relations that were proposed in this study. More targeted policy implications also need to be put forward based on the findings.
The scaling relations of ESs were influenced by diverse factors. In addition to landscape composition metrics (i.e., LULC patterns), exploring the scaling relations of landscape configuration metrics also can help us better understand the underlying biophysical mechanisms of landscape patterns affecting the ES scaling relations. Future research still requires some simulation-based work to better understand how the shape of the landscape might affect the scaling functions and coefficient estimates. Natural and socio-economic indicators also affect the ES patterns. Previous research has mainly focused on exploring how these attributes change with different city-scales (where population is usually used as the proxy indicator) (Bettencourt and Lobo 2016;Zhao et al. 2018). These studies often adopt power-law functions to express the relationships between city attributes and scales (Wu 2004;Bettencourt et al. 2010;Bettencourt 2013). In the future, we should further compare different fitting curves derived from the scaling relations of the city scale indicator (e.g., population) or spatial extent.
In this study, we selected three highly developed national-level urban agglomerations in China to reflect the scaling relations of multiple LULC types and ESs. The three urban agglomerations exhibited many similar and consistent trends in the scaling relations. The significant findings and rules in this study can help provide paradigms for other similar megacities in China or other countries. Previous research also proved that the scaling relations of developed lands in other six world metropolitan regions were closely resembled those of Beijing, Shanghai, and Guangzhou (Ma et al. 2019). Future studies still must validate whether the scaling relations of different LULC types and ESs will hold in urban agglomerations with other development levels.

Conclusions
Exploring how LULC and ES patterns change spatially across multiple scales is helpful for promoting urban landscape sustainability. Based on LULC maps and ES assessment results, this study used the scalogram forms to explore how different LULC types and ES indicators respond to changing spatial extents in the three largest urban agglomerations of China. The results revealed that different LULC mainly exhibited three types of predictable scaling relations, linear, exponential, and sigmoidal relationships, among which developed land and cropland were the most predictable. LULC types at the city proper level were more predictable than the metropolitan and urban agglomeration levels. Most of the ES indicators increased when the spatial extent increased, but their scaling relations varied; provisioning services were the most predictable, while soil retention was unpredictable. The ESs predominantly presented two types of predictable scaling functions, among which ES patterns were mainly expressed as a power-law, particularly at the lower levels, while exponential relationships were applicable to the higher levels. Among the three urban agglomerations, BTH performed the best, while YRD performed the worst in terms of the predictability of the ES spatial scaling. Since the scaling relations of ESs might be deeply affected by the LULC patterns, implementing land use composition and configuring optimization strategies are conducive to ecological conservation planning. The scaling relations of ESs can provide a scientific basis for the boundary determination of urban land use management.