How to Detect Scale Effect of Ecosystem Services Supply ? A Comprehensive Insight from Xilinhot in Inner Mongolia , China

Spatial scale plays a crucial role in the assessment and management of ecosystem services (ES), yet explicit information for identifying and understanding the scale effect on ES supply remains limited. In an attempt to detect scale effect on ES supply from a comprehensive perspective, this study developed a framework for integrating scale effect in three aspects, including individual ES patterns, pairwise ES interactions, and ecosystem service bundles (ESB). The framework was tested in Xilinhot, a prairie landscape city of Inner Mongolia, at four different levels of spatial scale. The results indicated that, most ES showed a decreasing clustering at coarser scales in terms of spatial pattern. At the same time, coarser scales resulted in fewer trade-offs and stronger synergies between pairwise ES. The identification of ESB varied greatly with scale, and this change reflected in the composition of ES variables and spatial distribution of bundles. We attributed the scale effect of the above three aspects to differences in social-ecological factors and their driving mechanisms at different scales. This comprehensive framework could support local managers to coordinate the management of multiple ES at different scales.


Introduction
Ecosystem services (ES) are benefits that humans obtain directly or indirectly from ecosystems [1].ES can promote human well-being, and therefore, have attracted much attention of the general public and governments worldwide [2].Over the past 20 years, research has focused on quantifying, assessing, mapping, and managing ES [3,4].However, many questions remains to be answered about the relationships among multi ES across heterogeneous landscapes [5,6].
In the current study, spatial extent is defined as the overall size of the study area, while spatial grain is defined as the individual unit of observation or analysis [7].Differences in ecological patterns and processes, observed at different scales, are termed scale dependence or scale effect [8].Some researchers firstly proposed that the scale dependence of ES should be determined by four aspects, including supplier, beneficiary, measurement, and management [9][10][11][12][13].However, the scale issue is complex in the context ES assessments.Ecosystems have different hierarchical structures and the supply of ES depends on different social-ecological process.In addition, the unit scale of supply determines the output of ES.For example, crop resistance to pests, weeds, and disease can be provided at the genetic level, but the ES that vegetation provides in terms of regulating water flow is provided at the habitat and community levels [14].ES at various scales interact in different ways, which produce different mutual benefits and mutual constraints.At the same time, the understanding and demand of ES by the stakeholders at different spatial scales have their own emphasis, while management decisions exert their own influences on the spatial distribution of ES [15,16].In general, local residents are more concerned about primary services, such as food supply and aesthetic appreciation, which provide direct benefits [17,18].Meanwhile, policy makers at regional levels [19][20][21], national [22,23] and global levels [24] tend to promote ES with long-term benefits to social welfare, such as climate regulation.As the spatial pattern of ES often manifest differently at multi-level spatial scales, it is important to identify these differences in order to promote the sustainable management of ecosystems [25,26].
A principal challenge is that ES do not vary independently of each other and may have highly non-linear relationships [27].Usually, when people prefer to consume or maximize the benefits of one or more ES, they will affect the supply of other ES intentionally or unintentionally.Such relationship mechanisms among ES are generally defined as tradeoffs or synergies [25,28,29].To the best of our knowledge, most ES studies have focused on the spatial distribution of ES supply or pairwise ES interactions on a single space scale.Specific regular grid unit [30,31] or administrative division unit [18,20,21,[32][33][34] were widely applied as the research unit in these studies.However, ES assessments conducted at only a single scale may miss or distort the associations between ES [35].For example, efforts to strengthen food production services at a local scale may affect water quality at broader scales, potentially threatening the quality of aquatic habitats [36].Such trade-offs can only be observed when multiple scales are considered.Some researchers have attempted to focus on scale dependence by examining pairwise ES interactions.However, such attempts have been limited to the upscaling on the spatial scale of the same property (such as kilometer grid) [37][38][39], and several researchers even set up only two scale levels for simple comparison.Obviously, it is not yet possible to explain how trade-offs or synergies are related across different levels of scale.Besides this, multiple ES, usually more than three, can form bundles in the space due to their similar responses to social and ecological factors.The concept of Ecosystem Service Bundles (ESB) has been proposed with an operable framework to help in the search for general rules determining associations among multiple ES [28].Thus, the scale dependence of ES supply is multi-dimensional, which can be embodied or detected in the spatial pattern of individual ES, pairwise ES interactions, and ecosystem service bundles.However, most researchers have overlooked the scale dependence of ESB due to method applicability or data reliability [35,38].This research deficiency leads us to incomprehensive understanding of the scale effect on ES supply and hinders our ability to apply this knowledge to management decisions.
In this study, we thus chose to focus on the scale effect of ES supply and conducted a case study in Xilinhot, a pastoral city in Inner Mongolia steppe area of China.The landscape of this research area, dominated by the natural grassland ecosystem, is relatively homogeneous because local people carried out weak transformation to the original landscape.We did not put the research area in the urbanized area because some services supply in urban ecological systems are land-independent that do not compete with each other in the space [40], such as the industry sector (manufacture) and the tertiary sector (education).Usually, such ES are considered as not directly driven by the land use or land cover (LULC), which may bring troubles to spatially explicit quantization and mapping.This paper presents a framework to integrate and detect scale effect of ecosystem services spatial patterns and associations among them.Our analysis proceeded in the following three steps.Firstly, we mapped the actual or potential values of eight typical ecosystem services across LULC types.Next, we set up four level spatial scales, covering different grain sizes and properties, including

Study Area
Xilinhot is a prairie landscape city in central Inner Mongolia (43 • 02 -44 • 52 N and 115 • 18 -117 • 06 E) and covers a total area of 14,785 km 2 (Figure 1A,B).The terrain is hilly in the South with gently undulating plains in the North and the average altitude is 988.5 m.The study area is located in the arid and semi-arid climate zone, which is characterized by cold temperatures, frequent drought, and strong wind [41].The annual average temperature is stable at 0-3 • C, and July is the hottest month (21 • C).The average annual precipitation (PPT) is 295 mm, and much of the PPT occurs in July and August due to the large seasonal variability of this region.As the core area of XilinGol Grassland National Nature Reserve in China, the zonal soil is chestnut soil and the zonal vegetation is the temperate steppe dominated by xerophilous and bunch grass.The main constructive species and dominant species of steppe vegetation are Stipa grandis, S. kiriensis, Leymus chinensis, etc. Taking land use patterns of Xilinhot in 2010 as an example (Figure 1D), grasslands accounted for 88.63% of the total study area, followed by unused land (including sand land, bare land and saline land) (7.34%), farmland (2.43%), construction land (0.67%), open water (0.5%), and forests (including shrubs) (0.43%).

Study Area
Xilinhot is a prairie landscape city in central Inner Mongolia (43°02′-44°52′ N and 115°18′-117°06′ E) and covers a total area of 14,785 km 2 (Figure 1A,B).The terrain is hilly in the South with gently undulating plains in the North and the average altitude is 988.5 m.The study area is located in the arid and semi-arid climate zone, which is characterized by cold temperatures, frequent drought, and strong wind [41].The annual average temperature is stable at 0-3 °C, and July is the hottest month (21 °C).The average annual precipitation (PPT) is 295 mm, and much of the PPT occurs in July and August due to the large seasonal variability of this region.As the core area of XilinGol Grassland National Nature Reserve in China, the zonal soil is chestnut soil and the zonal vegetation is the temperate steppe dominated by xerophilous and bunch grass.The main constructive species and dominant species of steppe vegetation are Stipa grandis, S. kiriensis, Leymus chinensis, etc. Taking land use patterns of Xilinhot in 2010 as an example (Figure 1D), grasslands accounted for 88.63% of the total study area, followed by unused land (including sand land, bare land and saline land) (7.34%), farmland (2.43%), construction land (0.67%), open water (0.5%), and forests (including shrubs) (0.43%).

Quantification of ES and Data Sources
Specific features of the study region shaped our selection of the ES and ES assessment indicators used in this study.There is a fact that the study area is not only an important livestock production base in China, but also the priority implementation area of some ecological projects, including the Grain for Green Program (GFGP).During China's 12th Five-Year Plan, a social and economic development initiative lasting from 2011-2015, ecological restoration policies implemented in Inner Mongolia prairie area focused on returning farmland to grasslands, implementing rotational grazing and zoning, and prohibiting over-grazing.The study area lies within the Northern part of the Beijing-Tianjin region, an area prone to soil erosion by water and wind, where the government launched a comprehensive watershed management project (Beijing-Tianjin Sandstorm Source Control Project) to prevent grassland degradation and desertification [42].Therefore, our selection of ES was based on the following two considerations: (1) The indicators should be importance of the region's social economy (e.g., animal husbandry is a pillar industry in the region); and (2) the indicators should represent the special status as a barrier of ecologic safety in Northern China.
In the end, we selected eight ES indicators, covering a range of provisioning services (grass production, livestock density, and water yield), regulating services (soil conservation, sand fixation, and carbon storage), maintenance services (habitat quality), and cultural services (landscape aesthetics) using the Millennium Ecosystem Assessment (MEA) classification systems (Table 1).In this paper, we mapped the physical quantity for each ES at 2010 to make sure we can analyze the ES supply under the same time node.In order to map each ES in the spatially explicit solution, we use an evaluation model or production function to quantify each ES based on high spatial resolution data.Table 1 provides quantitative methods and references for all services.The detailed description of the quantification of services is given in Supplementary Materials.
The main data used in this study were summarized in Table 2. Based on the ArcGIS (Version 10.3), using tools including reprojection, rasterization, and resampling, all data were unified by a projection coordinate system and a spatial resolution of (30 m in this study), thus we established the standardized spatial database of research area.

Setting Space Scales of Multiple Levels
We had the choice of multiple spatial scale here due to three considerations, including the extent limit of the study area, the spatial resolution of the model input data, and the relationship between the ES and the relevant scales.We finally chose four scales to represent different grain sizes and ecological hierarchies, and their spatial distribution was shown in Figure 2. The individual point (IP) scale (1 km grid) is the most commonly used scale in current research, representing the smallest observation and analysis scale in many ES mapping studies [34,35,39].We chose a 1 km grid as the point scale because the original data of livestock density is based on a 1 km grid and the land use pattern of grassland is not as fragmented as in urban areas.The landscape patch (LP) scale (5 km grid) corresponds to the scale of plant microhabitats and this size is equivalent to the unit of pasture that each herdsman family can manage.The plant community (PC) scale is a completely natural, and is controlled by the biophysics properties of the dominant population.The sub-watershed (SW) scale represents a type of geomorphic scale and is currently the minimum scale used in ecological restoration projects for wind and water erosion in the study area [41][42][43].We did not add a spatial level of administrative division scale, such as a village or township scale, because administrative boundaries may fluctuate throughout recent years.

Setting Space Scales of Multiple Levels
We had the choice of multiple spatial scale here due to three considerations, including the extent limit of the study area, the spatial resolution of the model input data, and the relationship between the ES and the relevant scales.We finally chose four scales to represent different grain sizes and ecological hierarchies, and their spatial distribution was shown in Figure 2. The individual point (IP) scale (1 km grid) is the most commonly used scale in current research, representing the smallest observation and analysis scale in many ES mapping studies [34,35,39].We chose a 1 km grid as the point scale because the original data of livestock density is based on a 1 km grid and the land use pattern of grassland is not as fragmented as in urban areas.The landscape patch (LP) scale (5 km grid) corresponds to the scale of plant microhabitats and this size is equivalent to the unit of pasture that each herdsman family can manage.The plant community (PC) scale is a completely natural, and is controlled by the biophysics properties of the dominant population.The sub-watershed (SW) scale represents a type of geomorphic scale and is currently the minimum scale used in ecological restoration projects for wind and water erosion in the study area [41][42][43].We did not add a spatial level of administrative division scale, such as a village or township scale, because administrative boundaries may fluctuate throughout recent years.
Vector map layers at the IP scale (1 km grid) and LP scale (5 km grid) were both created using the Fishnet tool in ArcGIS.Vector range of the PC scale was represented by vegetation community regionalization data directly.The SW scale used the vector range of sub watershed, which was obtained by using the Soil and Water Assessment Tool (SWAT) based on DEM data.The conversion process of ES data from fine scale to coarse scale was based on the Zonal Statistics tool in ArcGIS, so we can extract the statistical values of all ES raster layers within the vector zones of four spatial scales.In addition, the average value method was used to extract raster values [44].

Data Analysis
After the mapping of the eight ES, we set the range of each ES raster data to a fixed scale (for example, 0-1 in this study) by the method of min-max normalization for the convenience of data management and analysis [20,21,45].The spatial clustering of each ES was determined using both global and local Moran's I with queen contiguity (i.e., common borders and corners) to visualize and compare spatial distribution patterns [39].Correlation analysis was performed on each pair of services in R statistical software [46].Due to the nonlinear and non-normal characteristics of geospatial data distribution, Spearman's rank correlation analyses and corresponding significance tests were implemented using the "Hmisc" R package.Positive and negative correlations, respectively, represented the relationship of synergy and trade-off between pairwise ES.The absolute value of correlation coefficient changed from small to large, which represented the internal Vector map layers at the IP scale (1 km grid) and LP scale (5 km grid) were both created using the Fishnet tool in ArcGIS.Vector range of the PC scale was represented by vegetation community regionalization data directly.The SW scale used the vector range of sub watershed, which was obtained by using the Soil and Water Assessment Tool (SWAT) based on DEM data.The conversion process of ES data from fine scale to coarse scale was based on the Zonal Statistics tool in ArcGIS, so we can extract the statistical values of all ES raster layers within the vector zones of four spatial scales.In addition, the average value method was used to extract raster values [44].

Data Analysis
After the mapping of the eight ES, we set the range of each ES raster data to a fixed scale (for example, 0-1 in this study) by the method of min-max normalization for the convenience of data management and analysis [20,21,45].The spatial clustering of each ES was determined using both global and local Moran's I with queen contiguity (i.e., common borders and corners) to visualize and compare spatial distribution patterns [39].Correlation analysis was performed on each pair of services in R statistical software [46].Due to the nonlinear and non-normal characteristics of geospatial data distribution, Spearman's rank correlation analyses and corresponding significance tests were implemented using the "Hmisc" R package.Positive and negative correlations, respectively, represented the relationship of synergy and trade-off between pairwise ES.The absolute value of correlation coefficient changed from small to large, which represented the internal relationship between paired ES from weak to strong.We used the "psych" R package to conduct visualization treatments on the scatter distribution and correlation coefficient matrices for the data.
In order to identify ESB accurately, we have adopted two steps, including Principal Component Analysis (PCA) and Self-Organizing Map (SOM).As a precursor to cluster analysis, PCA is perfect for our dataset with multiple ES variables because an approximate linear relationship was commonly found between different ES variables in this study.Based on the principle of dimensionality reduction and ordination, this method can capture the variability in the dataset by transforming spatially correlated ES into uncorrelated several bundles of ES [47].Subsequently, the main multivariate relationships among ES variables were quantified to assess whether multiple ES could co-exist in a certain ESB [48].The load map of principal components provided guidance for more stable clustering solutions.We applied the SOM method to identify the ecosystem service bundles on account of its robustness to abnormal values and strong non-linearity mapping capacity [49].Consequently, we used the "kohonen" R package to learn and train dataset to cluster locations of space cells objectively according to their similarity in multi-ES supply.In order to avoid the local optimal solution, we set the number of training to ten times the number of space cells, the initial learning rate and the final learning rate are 0.05 and 0.01, respectively.The algorithm of SOM was parametrized to build three to eight clusters at each scale, respectively, and we then used the Silhouette Width Index (SWI) to determine the optimal number of clusters [34].The "starplot" R function was used to visualize the provision of ecosystem services by different clusters.The optimal clustering results will be mapped to vector data of different scales to visualize the spatial distribution pattern of the ecosystem service bundles.

Spatial Distribution of Ecosystem Services Supply
While there are certain similarities in terms of the spatial characteristics of different ES (such as CS and HQ), each ES showed unique spatial characteristics in maps of ES supply distribution (Figure 3).For example, the spatial distribution of the GP supply increased from a low level in the West to a higher level in the East, while the spatial distribution of LD and WY gradually decreased from a high level in the Southeast to a lower level in the Northwest.In particular, the decrease in WY was the most obvious.
Sustainability 2018, 10, x FOR PEER REVIEW 7 of 20 relationship between paired ES from weak to strong.We used the "psych" R package to conduct visualization treatments on the scatter distribution and correlation coefficient matrices for the data.
In order to identify ESB accurately, we have adopted two steps, including Principal Component Analysis (PCA) and Self-Organizing Map (SOM).As a precursor to cluster analysis, PCA is perfect for our dataset with multiple ES variables because an approximate linear relationship was commonly found between different ES variables in this study.Based on the principle of dimensionality reduction and ordination, this method can capture the variability in the dataset by transforming spatially correlated ES into uncorrelated several bundles of ES [47].Subsequently, the main multivariate relationships among ES variables were quantified to assess whether multiple ES could co-exist in a certain ESB [48].The load map of principal components provided guidance for more stable clustering solutions.We applied the SOM method to identify the ecosystem service bundles on account of its robustness to abnormal values and strong non-linearity mapping capacity [49].Consequently, we used the "kohonen" R package to learn and train dataset to cluster locations of space cells objectively according to their similarity in multi-ES supply.In order to avoid the local optimal solution, we set the number of training to ten times the number of space cells, the initial learning rate and the final learning rate are 0.05 and 0.01, respectively.The algorithm of SOM was parametrized to build three to eight clusters at each scale, respectively, and we then used the Silhouette Width Index (SWI) to determine the optimal number of clusters [34].The "starplot" R function was used to visualize the provision of ecosystem services by different clusters.The optimal clustering results will be mapped to vector data of different scales to visualize the spatial distribution pattern of the ecosystem service bundles.

Spatial Distribution of Ecosystem Services Supply
While there are certain similarities in terms of the spatial characteristics of different ES (such as CS and HQ), each ES showed unique spatial characteristics in maps of ES supply distribution (Figure 3).For example, the spatial distribution of the GP supply increased from a low level in the West to a higher level in the East, while the spatial distribution of LD and WY gradually decreased from a high level in the Southeast to a lower level in the Northwest.In particular, the decrease in WY was the most obvious.

Anselin Local Moran's I
Anselin Local Moran's I was used to further explore statistically significant hot or cold spots and spatial outliers for all eight ES at the four different spatial scales (Figure 5).We found that, at the two finer scales (IP and LP scales), local spatial clustering patterns were dominated by High-High Cluster and Low-Low Cluster, while High-Low Outlier and Low-High Outlier associations were rare.However, there was a tendency that the High-Low Outlier and Low-High Outlier patterns increased sharply at coarser scales.This is because the smaller and isolated hot spots (High-High Cluster) or cold spots (Low-Low Cluster) are prone to be cut or filled at coarser scales [39].Nevertheless, at the SW scale, the hot spots of LD and WY have a large overlap area, which further validated our interpretation of livestock clustering.

Anselin Local Moran's I
Anselin Local Moran's I was used to further explore statistically significant hot or cold spots and spatial outliers for all eight ES at the four different spatial scales (Figure 5).We found that, at the two finer scales (IP and LP scales), local spatial clustering patterns were dominated by High-High Cluster and Low-Low Cluster, while High-Low Outlier and Low-High Outlier associations were rare.However, there was a tendency that the High-Low Outlier and Low-High Outlier patterns increased sharply at coarser scales.This is because the smaller and isolated hot spots (High-High Cluster) or cold spots (Low-Low Cluster) are prone to be cut or filled at coarser scales [39].Nevertheless, at the SW scale, the hot spots of LD and WY have a large overlap area, which further validated our interpretation of livestock clustering.

Significance of Trade-Offs and Synergies
We analyzed a total of 28 potential pairwise ES and their interactions among the eight ES (Figure 6).Apparently, as the scale became coarser, accompanied by the decreased number of pairwise ES that could still maintain significant correlation (p < 0.05).At the IP scale, 27 pairs were significantly correlated with each other, and this number only dropped to 26 at the LP scale.However, as the scale became the PC scale, the number rapidly dropped to only 18, and further decreased to 12 at the SW scale.Only one pair, CS-HQ, kept very significantly correlative with each other at all four scales (p < 0.001).This was because the results of significance tests were easily impacted by the sample size, which decreased exponentially from fine to coarse scales as the numbers of spatial units at the four scales were 16,172, 712, 86, and 22 respectively.At the same time, we found that the number of ES pairs with significant negative correlation was decreasing in the process of scaling up (p < 0.05).There was only a significant trade-off between the SC and LA services at the SW scale.

Significance of Trade-Offs and Synergies
We analyzed a total of 28 potential pairwise ES and their interactions among the eight ES (Figure 6).Apparently, as the scale became coarser, accompanied by the decreased number of pairwise ES that could still maintain significant correlation (p < 0.05).At the IP scale, 27 pairs were significantly correlated with each other, and this number only dropped to 26 at the LP scale.However, as the scale became the PC scale, the number rapidly dropped to only 18, and further decreased to 12 at the SW scale.Only one pair, CS-HQ, kept very significantly correlative with each other at all four scales (p < 0.001).This was because the results of significance tests were easily impacted by the sample size, which decreased exponentially from fine to coarse scales as the numbers of spatial units at the four scales were 16,172, 712, 86, and 22 respectively.At the same time, we found that the number of ES pairs with significant negative correlation was decreasing in the process of scaling up (p < 0.05).There was only a significant trade-off between the SC and LA services at the SW scale.

Strength of Trade-Offs and Synergies
We summarized nine styles of change trends in this study by which the relationship between pairwise ES has changed in the process of upscaling (Figure 7A-I).The vast majority of paired ES, which shows a positive correlation, the strength of correlation increased in the intervals with significance.However, for pairwise ES interactions showing negative correlations, the strength change of correlation has relative uncertainty in the intervals with significance including unidirectional decreases or increases, or an initial decrease followed by an increase.It is worth noting that WY-SC did not pass the significance test at the IP, LP, or SW scales, and it showed no correlation (|r| ≤ 0.1) at both IP and LP scales.Nevertheless, this pair showed a significantly moderate trade-off at the PC scale (−0.5 < r ≤ −0.3).Such results also verified that the relationship between pairwise ES may be easily missed at one single scale.However, although HQ-LA showed similar relationship changes, we generally thought that there was no correlation between these two services, because correlation coefficient at the IP scale was much smaller than 0.1.

Strength of Trade-Offs and Synergies
We summarized nine styles of change trends in this study by which the relationship between pairwise ES has changed in the process of upscaling (Figure 7A-I).The vast majority of paired ES, which shows a positive correlation, the strength of correlation increased in the intervals with significance.However, for pairwise ES interactions showing negative correlations, the strength change of correlation has relative uncertainty in the intervals with significance including unidirectional decreases or increases, or an initial decrease followed by an increase.It is worth noting that WY-SC did not pass the significance test at the IP, LP, or SW scales, and it showed no correlation (|r| ≤ 0.1) at both IP and LP scales.Nevertheless, this pair showed a significantly moderate trade-off at the PC scale (−0.5 < r ≤ −0.3).Such results also verified that the relationship between pairwise ES may be easily missed at one single scale.However, although HQ-LA showed similar relationship changes, we generally thought that there was no correlation between these two services, because correlation coefficient at the IP scale was much smaller than 0.1.

Figure 7.
Strength of correlations between pairs of ecosystem services varies with four different scales of space.The symbol at the bottom of the interval represents a significant difference; the green check marks represent significant correlation coefficients under this scale (p < 0.05) and red cross marks represent correlation coefficients with no significance.

Multivariate Relationships
As shown in Figure 8, most ES could be fully explained by the first two main component axes Dim1 and Dim2, while the interpretation axes of some specific ES (e.g., LA) could be found in higherorder principal components.However, in general, higher-order principal components were negligible because they only accounted for a small part of the variance and mainly contained random changes.There was a clear trend that, as spatial scale became coarser, the percentage of the total variance that could be explained by the first two principal components increased.At the IP scale, the first two principal components could explain 48.7% of the total variance.However, at the SW scale, the first two principal components could explain 74.3% of the total variance.However, the ES with the highest representativeness underwent major changes as the spatial scale changes.At the IP and LP scales, WY has the highest representativeness.Meanwhile the representativeness of SF and SC was weaker at these two scales because they were highly affected by vegetation and terrain conditions, which showed weaker spatial heterogeneity in the study area at the finer spatial scales (i.e., 1 km 2 or 5 km 2 cells).However, at the relatively coarse PC scale, SF and CS showed the highest representativeness.

Multivariate Relationships
As shown in Figure 8, most ES could be fully explained by the first two main component axes Dim1 and Dim2, while the interpretation axes of some specific ES (e.g., LA) could be found in higher-order principal components.However, in general, higher-order principal components were negligible because they only accounted for a small part of the variance and mainly contained random changes.There was a clear trend that, as spatial scale became coarser, the percentage of the total variance that could be explained by the first two principal components increased.At the IP scale, the first two principal components could explain 48.7% of the total variance.However, at the SW scale, the first two principal components could explain 74.3% of the total variance.However, the ES with the highest representativeness underwent major changes as the spatial scale changes.At the IP and LP scales, WY has the highest representativeness.Meanwhile the representativeness of SF and SC was weaker at these two scales because they were highly affected by vegetation and terrain conditions, which showed weaker spatial heterogeneity in the study area at the finer spatial scales (i.e., 1 km 2 or 5 km 2 cells).However, at the relatively coarse PC scale, SF and CS showed the highest representativeness.The more closely related ES variables are, the higher the degree of positive correlation is, the more likely they are to be clustered in the latter SOM.Additionally, the longer the ES variable arrow is, the farther from the origin, the more the variable is explained by the two principal components.As the length of some variable arrows is very close, we use the color gradient to distinguish the different COS2 values.

Spatially Explicit Bundles
Based on map visualization, our SOM results revealed that the process of ESB formation by multiple ES was clearly dependent on spatial scale.We found that the clustering results at the PC scale always maintained the highest value of SWI among all scales, which provided the best separation effect between the clusters (Figure 9).In terms of the optimal clustering solutions at each scale, except the number of four clusters was associated with the highest SWI value of 0.264 at the IP scale.In order to compare the combination and distribution of ESB at different scales more intuitively, we uniformly set the number of clusters to three.We found that, with the changing scale, the component and spatial distribution location of each class of ESB were markedly different (Figure 10).Taking the Class A bundle as an example, SF supply was relatively strong in class A at the IP and SW scales, while it was quite weak at the LP and PC scales.On the SW scale, Class A bundle could not be identified in the South of the study area, but Class A bundle was identified in the same location at the first three scales.Moreover, we can see that with the scale change, the area proportion of different classes of ESB has also changed.The more closely related ES variables are, the higher the degree of positive correlation is, the more likely they are to be clustered in the latter SOM.Additionally, the longer the ES variable arrow is, the farther from the origin, the more the variable is explained by the two principal components.As the length of some variable arrows is very close, we use the color gradient to distinguish the different COS2 values.

Spatially Explicit Bundles
Based on map visualization, our SOM results revealed that the process of ESB formation by multiple ES was clearly dependent on spatial scale.We found that the clustering results at the PC scale always maintained the highest value of SWI among all scales, which provided the best separation effect between the clusters (Figure 9).In terms of the optimal clustering solutions at each scale, except the number of four clusters was associated with the highest SWI value of 0.264 at the IP scale.In order to compare the combination and distribution of ESB at different scales more intuitively, we uniformly set the number of clusters to three.We found that, with the changing scale, the component and spatial distribution location of each class of ESB were markedly different (Figure 10).Taking the Class A bundle as an example, SF supply was relatively strong in class A at the IP and SW scales, while it was quite weak at the LP and PC scales.On the SW scale, Class A bundle could not be identified in the South of the study area, but Class A bundle was identified in the same location at the first three scales.Moreover, we can see that with the scale change, the area proportion of different classes of ESB has also changed.

The Scale Effect of Ecosystem Services Supply
In this study, the random distribution of several ES, including HQ and LA, could be seen at coarse spatial scales.Similar distributions were observed for the education service in the Yangtze River Delta region of China [20], and the tourism service in the Quebec region of Canada [35].We found that randomly distributed ES can be provided directly or influenced indirectly by semi-natural and artificial systems, which are often under intense human modification in large scale research units (such as the expansion of construction land and grazing intensity of pasture in this study).While previous findings shows that the spatial pattern of ES supply can be optimized through human configuration, we must be alert to that discrete distribution of ES hotspots may lead to spatial mismatches between supply and demand in local ecosystem management [52].A better

The Scale Effect of Ecosystem Services Supply
In this study, the random distribution of several ES, including HQ and LA, could be seen at coarse spatial scales.Similar distributions were observed for the education service in the Yangtze River Delta region of China [20], and the tourism service in the Quebec region of Canada [35].We found that randomly distributed ES can be provided directly or influenced indirectly by semi-natural and artificial systems, which are often under intense human modification in large scale research units (such as the expansion of construction land and grazing intensity of pasture in this study).While previous findings shows that the spatial pattern of ES supply can be optimized through human configuration, we must be alert to that discrete distribution of ES hotspots may lead to spatial mismatches between supply and demand in local ecosystem management [52].A better

The Scale Effect of Ecosystem Services Supply
In this study, the random distribution of several ES, including HQ and LA, could be seen at coarse spatial scales.Similar distributions were observed for the education service in the Yangtze River Delta region of China [20], and the tourism service in the Quebec region of Canada [35].We found that randomly distributed ES can be provided directly or influenced indirectly by semi-natural and artificial systems, which are often under intense human modification in large scale research units (such as the expansion of construction land and grazing intensity of pasture in this study).While previous findings shows that the spatial pattern of ES supply can be optimized through human configuration, we must be alert to that discrete distribution of ES hotspots may lead to spatial mismatches between supply and demand in local ecosystem management [52].A better understanding of the spatial clustering characteristics of ES could provide managers useful information to decide where to improve the supply of which ES at a particular scale.
When performing data integration of ES from a fine scale to a coarse scale, it may hinder local trade-offs between ES, especially when different ES exhibit spatial competition.Xu found that the ES of food supply and tourism generated trade-offs at fine spatial scales (≤7 km), but such trade-offs became synergies at a coarser spatial scale (>7 km) because the supply sources of tourism have increased at a coarser scale (water area is also considered a kind of tourism resource) [39].In our research, the same phenomenon was found between GP and SF, while this pair changed from trade-off (at the IP and LP scales) to synergy (at the SW scale).This is because large-scale units in the study can make the impact of meteorological factors such as wind speed and precipitation on these two services more significant, and the ability of high coverage vegetation weakening wind erosion can be brought into full play, while small-scale units may only form a local trade-off between these two services through topographic conditions.This explanation is based on the niche complementary effect in ecology in which a multi-species mixture, compared with a monoculture, can make full use one or more resources in the system to achieve higher productivity and benefits [20,38].Therefore, the units often include a greater number and type of ES at coarse scales.It is worth noting that there are most pairs of moderate and high trade-offs with significance at the PC scale (r ≤ −0.3), which might be caused by the spatial competition among different plant populations in grasslands.The relationships between ES, highly affected by the variability of vegetation characteristics, may be more easily observable in correlation analysis.
While we clustered three class of ESB, their degree of cluster separation identified at the four spatial scales was quite different.The SOM results obtained at the PC scale provided the highest SWI value (0.719), which was much larger than the corresponding values obtained at the other three scales.The PC scale also provided the highest degree of separation among clusters, indicating that the PC scale may provide a more reliable result for ESB identification in the study area.Since water is the most important factor in local development in arid regions, the PCA results obtained at the PC scale showed that the variance explained by WY was the highest among all eight ES variables.Therefore, in SOM, WY-dominated clusters (Class A Bundles) were most likely to be separated.ESB revealed that, in the typical steppe areas with arid and semi-arid climate, population distribution and grassland evolution, compared to LULC, played a more significant role in driving the spatial combination of multiple ES.Since the biophysical characteristics of different types of vegetation are significantly different, such as water consumption, evaporation, soil consolidation, and palatability to livestock [53][54][55].The coupling of multiple factors leads to the selection of scale is different from that of other regions with higher urbanization level.If only using the regular grid scale, we cannot accurately reveal associations among ES supply.Thus, a multi-scale analysis is necessary for us to obtain a more robust result in pairwise ES interactions and ES bundles.

Uncertainties in Scaling-Up Process
Uncertainties in ES modelling and underlying measurements have been widely discussed in many studies [56].However, scaling-up process or approaches can also bring uncertainty [57].Firstly, the problem of how to cluster ES from a fine scale to a coarse scale deserves wide attention.In this paper, we provided a solution based on GIS.The premise of this approach is to use polygon vector and zonal statistics to extract ES data at each scale, which is also a common sampling method in multiscale analysis.While some studies used random points to extract ES data at multiple scales, it remains questionable as to whether the results of random points at a coarser scale can reflect the overall condition of the 'units' due to the nugget effect.The nugget effect tends to originate from errors due to data with different spatial resolutions when we are sampling or resampling at different distances [58].However, heterogeneity and inherent variability of different ES dominated the different nugget effects, which remains to be understood and solved.
Secondly, the choice of spatial scales depends on the process of ES selection, while the resolution of ES quantification depends on the quantized model and input data, which limits the choice of the minimum scale.However, if multiple ES closely related to human activities are assessed in a large regional area (such as tourism service), it is often necessary to select the scales corresponding to multiple administrative regions, such as towns and cities [59].For example, the IP scale used in this study was a one km grid based on the following two reasons.First, it was constrained by the resolution of LD data.Second, a smaller scale would exponentially increase the amount of data (the 30 m scale would generate more than 200,000 grid data), which would create a heavy burden of computing.Therefore, the selection of one km as the individual point scale here was more of a compromise.However, the assessment resolution of some other services was less than one km.For example, the resolution of GP is 250 m and the resolution of CS is 30 m.Therefore, the setting of minimum scale to one km could easily lead to the loss of internal fluctuation details in some units with small resolution elements (LUCC and NDVI).This is the most easily overlooked source of uncertainty.However, with the increased resolutions in the data, such as land use, meteorology, and soil data, as well as further calibration of model parameters, the uncertainty will gradually decrease in future research.

A Framework for Integrating Patterns of Ecosystem Services and Associations among Them at Multiple Spatial Scales
The integration of the ES concept in decision-making has been increasingly promoted in the scientific literature, as well as in policy guidelines at different strategic levels [60].Most of the previous studies have focused on how to integrate ES into social values [61], or how to include ES and human welfare into the framework of management practice [18].However, few researchers have noticed the scale effect during evaluating the ES supply process.Based on our results, we have established a framework to integrate spatial patterns of ecosystem services and associations among them at multiple scales (Figure 11).
In this framework, we need to solve a core issue about how to detect the scale effect.In this study, we have coordinated the analysis routes in three areas.In terms of individual ES pattern, the results of cold and hot spot identification are relatively intuitive when the popular methods of Moran's I and Gi* statistics are used [62].We used a relatively simple correlation analysis to quantify pairwise ES interaction and its strength.Of course, the recently proposed statistical models, such as Root Mean Square Error (RMSE) [63,64] and Production Possibilities Frontier (PPF) [65], can also be used to describe such relationship more accurately, which depends on the actual situation of ES distribution in the study area.Many researchers have used cluster analyses to identify ESB, however, the algorithm of cluster analysis and the number of clusters needs to be subjectively determined.Therefore, the exploratory nature of cluster analysis makes it unsuitable for understanding the causality in the associations among multiple ES.It is necessary to perform PCA before clustering, so as to visualize major multivariate relationships.
Finally, social and ecological variables, which can affect ES supply at different scales, should be identified.Currently, the understanding regarding the differences among social and ecological variables at different scales, as well as the driving mechanisms of ES trade-offs, is scarce and should be urgently solved to characterize and explain the cross-scale correlation between these factors [5].Furthermore, we can introduce the policy objectives and stakeholder needs at multiple scales, so as to co-ordinate ecological scales and administrative scales in the multi-level management of ES [66].We will carry out relevant studies in the next step.
scientific literature, as well as in policy guidelines at different strategic levels [60].Most of the previous studies have focused on how to integrate ES into social values [61], or how to include ES and human welfare into the framework of management practice [18].However, few researchers have noticed the scale effect during evaluating the ES supply process.Based on our results, we have established a framework to integrate spatial patterns of ecosystem services and associations among them at multiple scales (Figure 11).

Conclusions
It is necessary to understand the effects of spatial scale on the measurement, analysis, and management of ES.We set up our study area in a city in Inner Mongolian dominated by natural grassland landscape.Based on the quantification and mapping of eight representative ecosystem services, we set up four levels of spatial scale to detect scale effect on ES supply.The results showed that scale effect of ES supply can be revealed in three aspects, including individual ES patterns, pairwise ES interactions, and ecosystem service bundles (ESB).This study proved that single-scale analysis can only capture the limited information of the ES supply pattern at a certain scale.Because the degree of heterogeneity of the spatial unit cannot be clarified, it may lead to distortion or omission of the relationship among certain ES.Our research has provided a framework to integrate patterns of ecosystem services and associations among them at multiple spatial scales.This framework is expected to be used in other regions or other types of ecosystems around the world, so that the trend of changes in this scale effect can be further tested and the results can provide an important reference for policy-makers to coordinate different scales.Our future work will focus on clarifying the driving mechanisms of tradeoffs and synergies at different scales.
(1) individual point (IP) scale; (2) landscape patch (LP) scale; (3) plant community (PC) scale; and (4) sub-watershed (SW) scale.Finally, we analyzed scale effects of ES supply from three aspects (1) individual ES patterns; (2) pairwise ES interactions; and (3) ecosystem service bundles.Our framework can help in obtaining a systematic understanding of on scale effect of ecosystem service supply, which may support the land use planning at local scale.Sustainability 2018, 10, 3654 3 of 21

Figure 1 .
Figure 1.(A) The location of the Inner Mongolia Autonomous Region in China; (B) the location of Xilinhot in the Inner Mongolia Autonomous Region; (C) typical grassland landscape in the study area; and (D) 2010 land use pattern for the study area."High coverage grassland" indicates that grass cover 50% or more of the land, "medium coverage grassland" indicates 20% to <50% coverage, and "low coverage grassland" indicates 5% to <20%).

Figure 1 .
Figure 1.(A) The location of the Inner Mongolia Autonomous Region in China; (B) the location of Xilinhot in the Inner Mongolia Autonomous Region; (C) typical grassland landscape in the study area; and (D) 2010 land use pattern for the study area."High coverage grassland" indicates that grass cover 50% or more of the land, "medium coverage grassland" indicates 20% to <50% coverage, and "low coverage grassland" indicates 5% to <20%).

Figure 2 .
Figure 2. Vector range and unit division of four different spatial scales analyzed in this study.As the 1 km grid for the individual point (IP) scale is relatively small, we show only the Northeast corner of the study area.

Figure 2 .
Figure 2. Vector range and unit division of four different spatial scales analyzed in this study.As the 1 km grid for the individual point (IP) scale is relatively small, we show only the Northeast corner of the study area.

Figure 3 .
Figure 3. Spatial distribution of eight ES across Xilinhot.Each color represents a different ES; the variation of color from deep to shallow indicates the change in ES supply from high to low.

Figure 3 . 21 3. 2 . 20 3. 2 .
Figure 3. Spatial distribution of eight ES across Xilinhot.Each color represents a different ES; the variation of color from deep to shallow indicates the change in ES supply from high to low.

Figure 4 .
Figure 4. Global Moran's I results for ES at four spatial scales.All eight ES passed the significance test at the IP, LP, and PC scale (p < 0.01), while four ES including SF, CS, HQ, and LA failed to pass the significance test at the SW scale (p < 0.01).

Figure 4 .
Figure 4. Global Moran's I results for ES at four spatial scales.All eight ES passed the significance test at the IP, LP, and PC scale (p < 0.01), while four ES including SF, CS, HQ, and LA failed to pass the significance test at the SW scale (p < 0.01).

Figure 5 .
Figure 5. Local patterns of spatial associations for ES at four different spatial scales.All eight services passed the significance test at four scales (p < 0.05).

Figure 5 .
Figure 5. Local patterns of spatial associations for ES at four different spatial scales.All eight services passed the significance test at four scales (p < 0.05).

Figure 6 .
Figure 6.Spearman's correlation coefficients for pairwise ES interactions at four different spatial scales.Positive correlations are in the blue spectrum and negative correlations are in the red spectrum.The significance level is shown as stars on the upper right side of each correlation coefficient (No * means p ≥ 0.05; * means p < 0.05; ** means p < 0.01; *** means p < 0.001).

Figure 6 .
Figure 6.Spearman's correlation coefficients for pairwise ES interactions at four different spatial scales.Positive correlations are in the blue spectrum and negative correlations are in the red spectrum.The significance level is shown as stars on the upper right side of each correlation coefficient (No * means p ≥ 0.05; * means p < 0.05; ** means p < 0.01; *** means p < 0.001).

Figure 7 .
Figure 7. Strength of correlations between pairs of ecosystem services varies with four different scales of space.The symbol at the bottom of the interval represents a significant difference; the green check marks represent significant correlation coefficients under this scale (p < 0.05) and red cross marks represent correlation coefficients with no significance.

Figure 8 .
Figure 8. Correlation circles of principle component analysis for four different scales of space.It can be interpreted as follow:The more closely related ES variables are, the higher the degree of positive correlation is, the more likely they are to be clustered in the latter SOM.Additionally, the longer the ES variable arrow is, the farther from the origin, the more the variable is explained by the two principal components.As the length of some variable arrows is very close, we use the color gradient to distinguish the different COS2 values.

Figure 8 .
Figure 8. Correlation circles of principle component analysis for four different scales of space.It can be interpreted as follow:The more closely related ES variables are, the higher the degree of positive correlation is, the more likely they are to be clustered in the latter SOM.Additionally, the longer the ES variable arrow is, the farther from the origin, the more the variable is explained by the two principal components.As the length of some variable arrows is very close, we use the color gradient to distinguish the different COS2 values.

Figure 9 .
Figure 9. Changes in the Silhouette Width Index of different SOM clusters from three to eight at four different spatial scales.

Figure 10 .
Figure 10.Ecosystem service bundles identified at four different spatial scales.The rose diagram on the left represents the composition of ES in each ESB, and the rectangle on the right side of different colors represents three classes of ESB under each scale.The bundle with the highest WY supply capacity was labeled as A, the bundle with the highest CS supply capacity was labeled as B, and the bundle with the lowest multi-functional capacity at different scales was labeled as C.

Figure 9 . 20 Figure 9 .
Figure 9. Changes in the Silhouette Width Index of different SOM clusters from three to eight at four different spatial scales.

Figure 10 .
Figure 10.Ecosystem service bundles identified at four different spatial scales.The rose diagram on the left represents the composition of ES in each ESB, and the rectangle on the right side of different colors represents three classes of ESB under each scale.The bundle with the highest WY supply capacity was labeled as A, the bundle with the highest CS supply capacity was labeled as B, and the bundle with the lowest multi-functional capacity at different scales was labeled as C.

Figure 10 .
Figure 10.Ecosystem service bundles identified at four different spatial scales.The rose diagram on the left represents the composition of ES in each ESB, and the rectangle on the right side of different colors represents three classes of ESB under each scale.The bundle with the highest WY supply capacity was labeled as A, the bundle with the highest CS supply capacity was labeled as B, and the bundle with the lowest multi-functional capacity at different scales was labeled as C.

Figure 11 .
Figure 11.A framework for integrating patterns of ecosystem services and associations among them at multiple spatial scales.

Figure 11 .
Figure 11.A framework for integrating patterns of ecosystem services and associations among them at multiple spatial scales.

Table 1 .
Ecosystem services analyzed in this study.

Table 2 .
Data sources used in this study.
Cold and Arid Regions Science Data Center at Lanzhou (http://westdc.westgis.ac.cn/data/)Climate data Temperature and precipitation data from 13 meteorological stations within the study area, as interpolated by ANUSPLIN at 30 m spatial resolution China Meteorological Sharing Service System (http://data.cma.cn/data/)