Assessment Method and Scale of Observation Influence Ecosystem Service Bundles

The understanding of relationships between ecosystem services and the appropriate spatial scales for their analysis and characterization represent opportunities for sustainable land management. Bundles have appeared as an integrated method to assess and visualize consistent associations among multiple ecosystem services. Most of the bundle assessments focused on a static framework at a specific spatial scale. Here, we addressed the effects of applying two cluster analyses (static and dynamic) for assessing bundles of ecosystem services across four different scales of observation (two administrative boundaries and two sizes of grids) over 13 years (from 2000 to 2013). We used the ecosystem services matrix to model and map the potential supply of seven ecosystem services in a case study system in the central high-Andean Puna of Peru. We developed a sensitivity analysis to test the robustness of the matrix. The differences between the configuration, spatial patterns, and historical trajectories of bundles were measured and compared. We focused on two hypotheses: first, bundles of ecosystem services are mainly affected by the method applied for assessing them; second, these bundles are influenced by the scale of observation over time. For the first hypothesis, the results suggested that the selection of a method for assessing bundles have inferences on the interactions with land-use change. The diverse implications to management on ecosystem services support that static and dynamic assessments can be complementary to obtain better contributions for decision-making. For the second hypothesis, our study showed that municipality and grid-scales kept similar sensitivity in capturing the aspects of ecosystem service bundles. Then, in favorable research conditions, we recommend the combination of a municipal and a fine-grid scale to assure robustness and successfully land-use planning processes.


Introduction
The ecosystem services (hereafter ES) concept-the benefits obtained from nature for human well-being [1]-has become an integrated framework in sustainability science [2]. The ES framework facilitates ecosystem conservation opportunities [3] and affords innovative and valuable data to help decision-making [4]. In that sense, ES research is a significant and rising field of research [3], gathering studies around the world that are largely focused on the assessment and the management of the state of ES [3,5]. Relationships between ES are an issue that has received increasing interest in the literature [6][7][8][9][10]. These review studies addressed the importance of the analysis of relationships between ES.
ES relationships vary over time [11][12][13][14] and depending on the scale of observation. For example, a situation of mutual enhancement among a pair of ES at the county level could become in an In our work, we studied the effects of applying two cluster analyses for assessing bundles of ES over time (years 2000, 2009, and 2013). We computed one (static) analysis with ES values and one (dynamic) with the amounts of changes in ES values at two times (∆ES, for short). We performed the assessments across four scales of observation: two administrative levels (provincial and municipal) and two grid resolutions (3 × 3 km and 0.25 × 0.25 km). To determine the differences between the results of each method, we measured and compared the configuration, spatial patterns, and historical trajectories of ES bundles. Additionally, a sensitivity analysis that simulated a scenario with changes in the scores of ES potential supply over time tested the inconsistencies of the ES matrix on the findings.

Study Site
The selected area is a section (12 provinces, 27,612 km 2 ) of a larger study site (24 provinces, 64,025 km 2 ) of previous studies [12,44] in the central high-Andean moist Puna (administrative departments of Junin, Huancavelica, and Ayacucho) (Figure 1). During the study period, these chosen provinces are characterized by high land-use change intensity [44], mainly due to farming expansion, agricultural de-intensification and deforestation [12]. Additionally, this study site has a high population density in the moist Puna, about 44 inhab./km 2 (2013), with the strongest urban development in the metropolitan areas of the two major cities, Huancayo and Ayacucho [54].
Land 2020, 9, x FOR PEER REVIEW 3 of 19 assessments across four scales of observation: two administrative levels (provincial and municipal) and two grid resolutions (3 × 3 km and 0.25 × 0.25 km). To determine the differences between the results of each method, we measured and compared the configuration, spatial patterns, and historical trajectories of ES bundles. Additionally, a sensitivity analysis that simulated a scenario with changes in the scores of ES potential supply over time tested the inconsistencies of the ES matrix on the findings.

Study Site
The selected area is a section (12 provinces, 27,612 km 2 ) of a larger study site (24 provinces, 64,025 km 2 ) of previous studies [12,44] in the central high-Andean moist Puna (administrative departments of Junin, Huancavelica, and Ayacucho) ( Figure 1). During the study period, these chosen provinces are characterized by high land-use change intensity [44], mainly due to farming expansion, agricultural de-intensification and deforestation [12]. Additionally, this study site has a high population density in the moist Puna, about 44 inhab./km 2 (2013), with the strongest urban development in the metropolitan areas of the two major cities, Huancayo and Ayacucho [54]. This landscape is dominated by an expansion of livestock breeding in the upper lands and an increase in farming in the fertile lowlands. This is typical of many mountain agroecosystems across the world. Most of this territory is embedded within the Mantaro river basin, including ecosystem services associated with agricultural practices (crops and livestock provision, regulation of soil erosion and maintenance of soil quality), hydrological cycle (water purification and water flow regulation) and climate regulation. However, the main land use/land cover in the study area consists of natural grasslands (59%), shrublands (16%), and agricultural lands (15%) by 2013 (Table S1, Supplementary Materials). The diverse combination of them formed three groups of landscapes. A group of provinces (Acobamba, Huamanga, Huanta, and Vilcas Huaman) show a uniform distribution of the main land use/land cover (hereafter LULC) units (31%, 34%, and 31%, respectively). The second pattern, described by two provinces (Jauja and Chupaca), displays a territory mainly dominated by two LULC units (72% of natural grassland and 19% of farming areas). The third group of provinces (Angaraes, Cangallo, Concepcion, Huancavelica, Huancayo, and This landscape is dominated by an expansion of livestock breeding in the upper lands and an increase in farming in the fertile lowlands. This is typical of many mountain agroecosystems across the world. Most of this territory is embedded within the Mantaro river basin, including ecosystem services associated with agricultural practices (crops and livestock provision, regulation of soil erosion and maintenance of soil quality), hydrological cycle (water purification and water flow regulation) and climate regulation. However, the main land use/land cover in the study area consists of natural grasslands (59%), shrublands (16%), and agricultural lands (15%) by 2013 (Table S1, Supplementary Materials). The diverse combination of them formed three groups of landscapes. A group of provinces (Acobamba, Huamanga, Huanta, and Vilcas Huaman) show a uniform distribution of the main land use/land cover (hereafter LULC) units (31%, 34%, and 31%, respectively). The second pattern, described by two provinces (Jauja and Chupaca), displays a territory mainly dominated by two LULC units (72% of natural grassland and 19% of farming areas). The third group of provinces (Angaraes, Cangallo, Concepcion, Huancavelica, Huancayo, and Huaytara) discloses a landscape characterized by low farming development (14%), the highest high-Andean wetland coverage (3%) and high pasture extent (63%).

Land Use/Land Cover Data Sources
The study area is covered by 11 LULC types ( Figure 1). These categories, included in the Peruvian standardized nomenclature of the Corine Land Cover (CLC), were derived from three-time data sources: the map of high-Andean ecosystems in 2000 [45], the official flora cover map from 2009 [55], and the official flora cover map from 2013 [56]. According to the official sources, the maps were submitted to a verification and field survey procedure for improving the accuracy of the land use/land cover classification. However, the different geographical scales made necessary a generalization of the land use/land cover classes. Table S2 (see Supplementary Materials) shows the harmonization of the three-time step features to obtain a common legend of eleven LULC units. Moreover, the description of satellite images, mapping scale, minimum mapping area, and type of data of the three source maps are specified in Table S3 (see Supplementary Materials).

Ecosystem Services Potential Supply
The study is based on the capacity matrix that was done specifically for assessing the ES in the study area, developed by Madrigal-Martínez and Miralles i García [12] (subsequently referred to as the High-Andean Study). The assessment obtained the potential supply of 7 site-specific ES identified by the Common International Classification of Ecosystem Services [57]-two regulating services related to mediation of flows (regulation of soil erosion and water flow regulation); one ES related to filtration, sequestration, storage or accumulation by ecosystems (water purification); two services linked to the maintenance of physical, chemical, biological conditions (soil quality and global climate regulation) and, finally, two provisioning services related to nutrition (crops and reared animals).
To develop the ES matrix, the High-Andean Study consulted 43 experts to rank the ES potential supply associated with a specific LULC on a relative spatial scale, ranging from 0 (no relevant ES potential supply) up to 5 (very high ES potential supply). The experts were carefully selected to increase confidence according to their specific skills on ES and the moist Puna ecosystems. Additionally, the survey was thoroughly described individually, and they scored only the LULC/ES pairs that were sure in their judgments. Each response was collected and deprived of outliers using the interquartile range method. Then, a final score was computed using the mean. Furthermore, the potential supply of the LULC in provisioning services was achieved from official model results included in land planning instruments of the administrative departments under study. Table S4 (see Supplementary Materials) provides the maximum capacity of the eleven LULC categories to supply the seven ES. The ecosystem services were set as constant values assuming that land units are in good condition during the study period.

Scaling Method
The ES and ∆ES maps were derived from the matrix model and upscaled to four spatial scales: two administrative divisions (provincial and municipal) and two grids (coarse and fine). The four spatial scales were selected for their particular importance in spatial planning and ecosystem services mapping. The provincial level (~10 3 km 2 ) has a central role in the Peruvian planning system binding national and departmental directives with local interventions (Organic Law of Municipalities No. 27,972). The municipal level (~10 2 km 2 ) is where land-use management in urban areas and the countryside are made. Coarse-grid resolution (9 km 2 ) was chosen because it explores patterns of ecosystem services and approximates a locality. A fine-grid (0.25 km 2 ) was included because it is where individual land-use management and land-cover changes occur. This spatial scale was decided as the finest because the study maps are based on a geographical scale of 1:100 000, and following the Corine Land Cover approach and the official flora cover map from 2009 [55], this spatial resolution corresponds to the minimum mapping area. Both grid resolutions are important for planning green infrastructure to support human well-being.
The administrative areas were calculated using boundaries from the Peruvian National Institute of Informatics and Statistics. The 12 provincial units range from 750 to 6075 km 2 (with an average of 2301 km 2 ), whereas the 175 municipality units vary from 5 to 2176 km 2 (with an average of 158 km 2 ). On the other hand, the coarse-grid (3 × 3 km) and the fine-grid (0.5 × 0.5 km) resolutions were both generated using the Fishnet tool and the Geoprocessing tool in ArcGIS 10.3 [58]. The coarse-grid comprises 3019 cell units, while the fine-grid has 110,343 spatial units. The cells with at least 95% of their area within the boundaries of the study area were included.
After this, each of the four maps of spatial units was separately intersected with every LULC map of each year (2000, 2009, and 2013), obtaining 12 maps. Next, the ES matrix was applied on these 12 maps deriving 84 maps of ES potential supply. These potential supply maps were aggregated to their corresponding spatial resolution by using Equation (1): where ES ns is the potential supply of a given spatial unit s for a given ecosystem service n, ES i is the score assigned to a given LULC unit i, and A i is the area of that given LULC unit i within the given spatial unit n. S is the total area of the given spatial unit. Figure S5 (see Supplementary Materials) provides a graphical sample of the scaling method. Lastly, to obtain the upscaled ∆ES values over the two periods, from 2000 (t 1 ) to 2009 (t 2 ) and 2009 (t 2 ) to 2013 (t 3 ), Equation (2) was used: where ∆ES ns is the potential supply of a given spatial unit s for a given ecosystem service n of the final year t k+1 minus the potential supply of that given spatial unit s for the given ecosystem service n of the initial year t k .
The best number of clusters was determined using the "NbClust" R package [59] configured with the combination of "euclidean" distance measure, "kmeans" method, "alllong" index, and a significance value of 0.1 for Beale's index. This package was run (n = 4) with ES and ∆ES values at the provincial and municipal levels. The majority of indices proposed three clusters as the best number in all datasets. Bundle types were identified applying a k-means cluster analysis run with 10,000 iterations in R [60]. The k-means cluster analysis grouped the values in three specific combinations of ES based on their characteristics. For later comparisons, the bundles were named: bundle type 1, bundle type 2, and bundle type 3. Each bundle type was drawn using Excel 2015. The different aspects of bundles were analyzed with standard metrics (Table 1). Then, the results were compared to identify differences (effects) that can establish trends. To estimate the configuration metrics, Excel 2015 was used. The spatial patterns and historical trajectories were computed using ArcGIS 10.3 [58]. Table 1. Metrics (and their description) used for the achievement of the aspects of bundles.

Aspect
Metric Description

Configuration
True Diversity (Order 2) The diversity of a set of ES provided in a given bundle type is calculated as the effective number of ecosystem services based on Hill numbers [61,62]. For the "dynamic bundles", we used the absolute value of each amount of change in ES specified by a given bundle. This metric was included because it affords a stable, clearly understood, and sensitive overall similarity measure supporting cross-study assessments [11,62].
The sum of the absolute value of each ES (or ∆ES) specified by a given bundle type. The sum represents an overall level of the provisioning of services (or of the change in services). High absolute values thus indicate zones with a comparatively high supply of (or change in) multiple services, while low values indicate the opposite. This metric was included in the bundle analysis because policies are intended to protect the overall level of ES provision rather than, or in addition, to the provision of individual services.

Spatial patterns
Percentage of land The proportional abundance of a given bundle type in a given year or a given period across the study area. It is a landscape metric that acts as a proxy for change, thus allowing for the interpretation of spatial patterns over time and space. This metric measured the results of both cluster analyses.

Historical trajectories
Percentage of land change The proportion of land changing from one bundle in a year or period t to another in a year or period t+1 on the same spatial scale. This metric measured the results of both cluster analyses.

Sensitivity Analysis
A sensitivity analysis of the ES matrix was applied to test the robustness of the methodological approach. The analysis consisted of the development of a sensitivity scenario based on a four steps method adapted from the five common stages of a scenario development [63]. In the first step, the aim of the sensitivity analysis was defined-to test how changes in the scores of ES potential supply of the High-Andean Study matrix affects the results over time. In the second step, two key drivers and their trends that affected (positively or negatively) the potential supply of services were identified from interviews with five experts: climate change and technological improvement of agriculture and forestry.
In stage three, the scenario assumptions were deducted using the trends of the key drivers. These trends were simulated as a rate of positive/negative change (+/− 0.1 per year) on the ES values of the LULC units. Climate change had negative consequences on regulating services supplied by the following ecosystems: natural grasslands, shrublands, forests, glaciers, and high-Andean wetlands. On the contrary, well-managed farming enhanced regulating (erosion, water flow, and soil quality) and provisioning services of agricultural areas and reduced the pollution of rivers and lakes, recovering their functions of purifying water and flow control. Likewise, the technological improvement of forest plantations increased, regulating services (soil quality, control of soil erosion, water flow, and global climate regulation). The scores of ES for continuous urban fabric and sparsely vegetated areas stayed unaffected.
In stage four, with the simulated scores of ES, two new model matrices for 2009 and 2013 were generated (see Tables S6 and S7 in the accompanying Supplementary Materials), whereas, for 2000, that created by the High-Andean Study was used. From these matrices, the ES maps at the four spatial scales were derived running the scaling method defined in Section 2.4. Finally, the assessments of relationships between ES were performed following Section 2.5.

Static Cluster Analysis
The results of the two metrics used to evaluate the effects of the four spatial scales on the configuration of bundles showed similarities and disparities (Figure 2A). Regarding similarities, all the bundles provided an effective number of ES that ranged from 6.51 to 6.87. Concerning dissimilarities, most of the bundle types indicated disproportions among the abundance of ES. However, it showed a trend towards being higher for large spatial scales. Additionally, there was a trend of increasing of ES abundance from bundle type 1 to type 3 at each spatial resolution, but it had more similarities when the spatial scale increased. In that way, the provincial level was defined by the slight variation of ES values of the three bundle types. However, at the municipal level, type 3 was a multifunctional bundle, type 2 was a multifunctional agricultural bundle, and type 1 corresponded to an agriculture bundle. The coarse-grid scale mainly differed from the municipal in the bundle type 1 (agriculture and sparsely vegetated areas). However, at the fine-grid, the ES bundling showed a multifunctional bundle (type 3), an agriculture bundle (type 2), and an urban and sparsely vegetated area bundle (type 1).
The sensitivity analysis showed similarities between the effective number of ES provided by all the bundles, whereas the highest differences were detected among the abundance of bundles ( Figure 2B). The diversity and the abundance of ES provided in bundles type 3 and type 2 was similar at the four scales of observation, whereas in type 1, differed. Thus, type 3 was a bundle with the highest values of regulating services, and type 2 was a bundle with the highest values in crop and livestock services. However, type 1 at the provincial level kept similarities with type 2, whereas at the municipal and grid scales had the lowest values of ES defined by urban and sparsely vegetated areas.
The spatial distribution of bundles obtained from ES values showed higher similarities among the three smaller spatial scales ( Figure 3A). Thus, bundle type 3 dominated the territory (percentage of land >63%) over the three years. Nevertheless, the agricultural bundle had higher correspondences between grid-scales. At the provincial level, the three types of bundles were more evenly distributed ( Figure 3A). The sensitivity analysis showed that the similarities between the spatial distribution of bundles followed a trend towards being higher for small spatial scales ( Figure 3B). Then, at the municipal level and the two grid-scales, bundles kept fair spatial consistency across time, especially for types 2 and 3. On the contrary, at the provincial level, the territory was defined by a bundle type each year.
The analysis of historical trajectories showed that the bundle provided by any given land changed through time at each spatial scale but followed a decreasing trend from large to small (Table S8, Supplementary Materials). During the total study period, at the provincial level, 68% followed any trajectory of change, whereas this change was 30% at the municipal level. In the same way, the coarse-grid and fine-grid showed inferior variations of 24% and 14%, respectively. Furthermore, there was a second trend towards a higher number of transitions for fine spatial scales. These two trends were confirmed by the sensitivity analysis (Table S9,

Dynamic Cluster Analysis
The analysis of the configuration of bundles at the four spatial scales presented similar measures of the effective number of ES changes (that ranged from 5.30 to 6.07), but differences in most of the N values ( Figure 4A). Only bundle type 2 did not manifest these dissimilarities, since describing a territory without land-use change at the four spatial resolutions, remaining with similar and lowest N (almost 0). On the contrary, the N values specified by bundles type 1 and type 3 decreased when the spatial scale increased. In this regard, bundle type 1 revealed an increasing pattern from larger to smaller spatial scales, that detected the reduction in regulating services, and the increase in provisioning ES. However, bundle type 3 specified a trend of increase in provisioning services and a decrease in regulating.

Dynamic Cluster Analysis
The analysis of the configuration of bundles at the four spatial scales presented similar measures of the effective number of ES changes (that ranged from 5.30 to 6.07), but differences in most of the N values ( Figure 4A). Only bundle type 2 did not manifest these dissimilarities, since describing a territory without land-use change at the four spatial resolutions, remaining with similar and lowest N (almost 0). On the contrary, the N values specified by bundles type 1 and type 3 decreased when the spatial scale increased. In this regard, bundle type 1 revealed an increasing pattern from larger to smaller spatial scales, that detected the reduction in regulating services, and the increase in provisioning ES. However, bundle type 3 specified a trend of increase in provisioning services and a decrease in regulating.
For the sensitivity analysis, Figure 4B shows the similarities and the differences between the configuration of bundles across the four spatial scales. Similarities of the 2 D metric are found for types 1 and 3, whereas type 2 showed higher differences across the four spatial scales. On the other hand, the N metric showed that for each bundle type, grid-scales had higher similarities between them and the municipality level. Furthermore, bundles type 1 and type 3 showed a consistent configuration of positive values of provisioning services and negative of regulating, whereas type 2 differed at the provincial level in the regulating services. Thus, bundles showed higher similarities among the three smaller spatial scales. For the sensitivity analysis, Figure 4B shows the similarities and the differences between the configuration of bundles across the four spatial scales. Similarities of the 2 D metric are found for types 1 and 3, whereas type 2 showed higher differences across the four spatial scales. On the other hand, the N metric showed that for each bundle type, grid-scales had higher similarities between them and the municipality level. Furthermore, bundles type 1 and type 3 showed a consistent configuration of positive values of provisioning services and negative of regulating, whereas type 2 differed at the provincial level in the regulating services. Thus, bundles showed higher similarities among the three smaller spatial scales.
The spatial distribution of bundles across the two smaller spatial scales displayed a consistent pattern that began to be less evident at the provincial level ( Figure 5A). In that sense, at the municipal level and on the two grid-scales, the territory seemed dominated by bundle type 2 (percentage of land >84%), whereas this percentage high declined at the provincial level. Likewise, the sensitivity analysis indicated fair robustness between municipal and grid-scales ( Figure 5B). However, there were minor areas with changes in ES supply only detected at grid resolutions. The spatial distribution of bundles across the two smaller spatial scales displayed a consistent pattern that began to be less evident at the provincial level ( Figure 5A). In that sense, at the municipal level and on the two grid-scales, the territory seemed dominated by bundle type 2 (percentage of land >84%), whereas this percentage high declined at the provincial level. Likewise, the sensitivity analysis indicated fair robustness between municipal and grid-scales ( Figure 5B). However, there were minor areas with changes in ES supply only detected at grid resolutions.
Historical trajectories of bundles achieved with ΔES values showed that the land that changed from one to another differed among spatial scales but was higher (52%) at the provincial level than at smaller levels (municipal: 24%; coarse-grid: 16%; fine-grid: 13%) (Table S10, Supplementary Materials). These transitions uncovered four main trajectories at all the spatial scales, and two more only found at the grid scales. Likewise, the sensitivity analysis showed that the proportion of land changing from one bundle to another was higher at the provincial level, and the number of trajectories was higher as the spatial scale decreased (Table S11, Supplementary Materials). Historical trajectories of bundles achieved with ∆ES values showed that the land that changed from one to another differed among spatial scales but was higher (52%) at the provincial level than at smaller levels (municipal: 24%; coarse-grid: 16%; fine-grid: 13%) (Table S10, Supplementary Materials). These transitions uncovered four main trajectories at all the spatial scales, and two more only found at the grid scales. Likewise, the sensitivity analysis showed that the proportion of land changing from one bundle to another was higher at the provincial level, and the number of trajectories was higher as the spatial scale decreased (Table S11, Supplementary Materials).

Discussion
In our study, the ES matrix contributes to the assessment of relationships between ES, applying two different methods (static and dynamic) across four scales of observation over time. At the spatial scale level, it revealed several findings consistent with those found by comparable biophysical assessment [16]. We analyzed the differences between each assessment method by comparing the results of standard metrics at each spatial scale over time. Subsequently, we discuss the main findings of the study validated by the sensitivity analysis (Table 2) and organized as scale and assessment method effects that might have implications on ES management.

Discussion
In our study, the ES matrix contributes to the assessment of relationships between ES, applying two different methods (static and dynamic) across four scales of observation over time. At the spatial scale level, it revealed several findings consistent with those found by comparable biophysical assessment [16]. We analyzed the differences between each assessment method by comparing the results of standard metrics at each spatial scale over time. Subsequently, we discuss the main findings of the study validated by the sensitivity analysis (Table 2) and organized as scale and assessment method effects that might have implications on ES management. Table 2. Scale and assessment method effects on bundles of ecosystem services.

Assessment Method Effect Spatial Scale Effect
• Configuration: disagreement in the direction of the relationships between multiple ES. • Spatial patterns: static cluster analysis captured only a snapshot of ES bundles at different years, whereas cluster analysis with ∆ES values displayed dynamics of ES bundles.
• Configuration: static cluster analysis displayed a trend towards more similarities among bundle types for large spatial scales, whereas dynamic cluster analysis showed a similar trend of positive and negative change in the ES supply at the three smaller spatial scales. • Spatial patterns: static cluster analysis suggested higher similarities between bundles at the municipal level and the two grid-scales, whereas dynamic cluster analysis showed some consistency across spatial scales.

•
Historical trajectories: both cluster analyses detected: (1) a trend towards a high percentage of land change for large spatial scales, and (2) a trend towards a high number of trajectories for fine spatial scales.

Effects of Different Cluster Assessments on Bundles of ES
Depending on the cluster assessment, we found relationships between multiple ES that shifted in different ways. This finding agrees with previous work that also confirmed that the chosen method influences the result [6,20,42]. In that sense, in our study, "static bundles" suggested a positive spatial co-occurrence among the seven ES. On the contrary, "dynamic bundles" proposed a negative relationship between provisioning and regulating services. The synergy detected with the static assessment shows an opportunity to enhance multiple ES simultaneously. However, it missed the trade-off between regulating and provisioning services, and it could represent an unexpected loss of success for ES management. In fact, it implicates missing opportunities for win-win solutions that involve investments in conservation, restoration, and sustainable ecosystem use [64].
The spatial distribution of bundles captured by each cluster assessment showed differences. Thus, ES values displayed a landscape characterized by bundles with a specific diversity and abundance of ecosystem services supply at each time-step. On the other hand, ∆ES values addressed the dynamics of ES bundles over the two time-periods. This last interpretation may facilitate the understanding of the instabilities that produce the temporal dynamics on ecosystems since trends expose whether there has been a change and the specified event that caused it [65]. This finding concerning "dynamic bundles" is consistent with previous research for the knowledge of land-changes dynamics [12].

Effects of Different Scales of Observation on Bundles of ES
The static assessment of bundles suggested that the configuration followed a trend towards more similarities at large spatial scales (Figure 2). This effect may explain that large spatial units follow a multifunctional landscape allowing relationships between ES to concur in synergy. It is understandable because the impacts of management actions at a fine-scale may be insignificant at a larger spatial scale if the land-use type affected is scarce, which is related to the capacity to capture local heterogeneity. Thus, the relationships between ES are conditioned by the geographical size of any single land-use change in the spatial unit. Consequently, at the grid scales, bundle types were more specialized according to one LULC unit (this was evident at the fine-grid scale). However, the provincial level provided a comparative abundance of ES because they were characterized by a similar combination of land-units. This similarity indicates comparable levels of land-use diversity that produces akin multifunctionality at large spatial units. Although multifunctionality is location related [66], this effect is observed in previous work of ES bundles across different administrative levels [16,24]. For instance, this generalization of the configuration can be inconvenient when we need to identify areas of highest/lowest supply of ES (hotspots/coldspots) for spatial prioritization or designing green infrastructure. For, as has been observed in our study site, the bundles of small size only persist across grid-scales. It implies a loss of bundle diversity when we upscale, which agrees with Zen et al. [67]. Then, large scales (dramatically at the provincial level) may fail to observe determinant factors and their influence on the sustainability of the ecosystems and their services. It reinforces the assumption that the increase in the spatial scale of observation brings a homogenization of the landscape [68], and only the mainland changes are significant [12].
At the three smaller spatial scales, bundles showed a similar configuration of positive and negative change in ES supply ( Figure 4B), reflecting higher accuracy with the rate of change established by the different drivers (climate change, and technological improvement of agriculture and forestry). Needless to say, these bundles offer a basic view of the dynamic of ES that may help in planning win-win solutions. However, this basic picture depends on the size of the spatial unit, since it determines the intensity of drivers of change. In our study, as large as the spatial scale was, the land-use change impacts were more buffered. Although the provincial bundles detailed many similarities with the smaller scales of observation, the contrasts involve caution when using this spatial scale for the management of ecosystem services.
Static cluster analysis suggested high similarities between the spatial distribution of bundles at the municipal level and the two grid-scales. Consequently, it manifested fair robustness across the three smaller spatial scales, which differed with Raudsepp-Hearne and Peterson [16]. It may be related to the Andean study area, which is a landscape with ecosystem services more evenly distributed, and some amount of each ES facilitating multifunctionality can be found at the municipality level. Thus, the variation of bundling across a territory depends on the spatial heterogeneity of services since spatial homogeneity uncovers the same type of bundle across spatial scales. This diversity of findings recommends that researchers and decision-makers should be aware of the size and the heterogeneity of the spatial units to improve the aims of ES analyses [69]. Even though many times, there are limitations related to data scarcity or availability, which impede the research from being conducted optimally. We agree with previous research that considering at least two spatial scales should assure robustness [70,71], but we suggest a fine-grid scale and the municipality level. A fine-scale is important to show specific spots at local level that give a better panorama for well-informed planning decisions, whereas, at the municipality level is where political decisions are made and socioeconomic data are available. However, it is worth emphasizing that our study shows sufficient consistency between the municipal scale and the grid-scales.
The spatial distribution of bundles resulting from ∆ES values revealed some consistency across spatial scales. However, bundling generalization was more evident as the scale of observation increased. This effect produces homogeneity at broad resolutions that can lead to shape a territory with similar land-use change intensity and overlooking fine-grained information needed for spatial conservation planning [72]. In our study site, at the provincial level, that generalization obscures changes in ecosystem services at lower levels that may be of importance for planning and management solutions. However, Madrigal-Martínez and Miralles i García [12] showed that, in research conditions of data scarcity, it is possible to address knowledge about land-change dynamics affecting ES that may help for policy and planning purposes at the provincial level.
For historical trajectories of bundles, both cluster analyses indicated that the area providing any given bundle changes higher at broad spatial scales over time. It implies that objects (land-units) within a large spatial unit are strongly associated, and a substantial change in one of them affects the total, whereas minor and static zones are overlooked. In our study area, this was more evident at the provincial level, in which the variation in ES supply of a given province was due to changes only in a few land-units. It is a consequence of upscaling that has direct impacts on the intensity of land-use change affecting ES. Low intense land-use change is not significant at broad scales [12]. In that sense, only at the grid-scales minor land-use changes that configured small size bundles were detected. This effect was detected in both cluster analyses and showed a trend towards a high number of trajectories for fine spatial scales. For example, we observe that bundles characterized by an increase in regulating services at grid-scales disappear at large (municipal and provincial). It reveals that changes at larger spatial scales have a buffer effect, whereas, at the fine-scales, bundles are more sensitive to temporal changes shaped by the direct local-scale drivers. This finding supports the assumption that knowledge of local contexts of ES is policy-relevant since their changes in values and demand are finer observable over time [73,74]. Therefore, the assessment of the spatial extension under the influence of drivers could help with the understanding of the stability of ES provision, endorsing robustness for the development of sustainable management and conservation strategies.

Methodological Limitations
In this study, the analyses presented should be understood as using the best existing data of an acceptable quality to admit a robust demonstration. Even so, the method (ES matrix) brings potential limitations to the study, and technical and thematic uncertainties [75]. In that sense, we highlight that the capacity matrix simplifies landscape functionality producing uncertainties in the quantification of ES (e.g., regulating services). It is due to this that some ES are not only dependent on the presence of certain land use/land cover types but also their spatial configuration. Moreover, management actions on each land-use may affect ES flow differently (specially in provisioning services), and this effect could be measured vaguely for the matrix. Another limitation lies in that the reduced and diverse data sources of land use/land cover classes made a generalization of the landscape necessary, which could influence the bundles that emerge at larger spatial scales. In fact, a more precise number of land use/land cover classes could result in the reconfiguration of bundles [69]. Additionally, in ES matrix models, the multifunctionality is strongly dependent on the number of services provided by the different land use/land cover types [76]. On the other hand, when data at a fine-scale were summarized at the administrative levels (aggregation effect), they could cause a loss of information [77]. Finally, the data source (the map of high-Andean ecosystems) has a vague delimitation for two land units (agricultural areas and forest plantation), comprehending them in only one land-use category (Areas modified by human action). However, we considered this limitation of minor importance because this aspect was clarified using the land-use types from the two official flora cover maps.

Conclusions
We developed a study that addressed the effects of different cluster methods for assessing bundles of ES across different scales of observation over time, using an example in the high-Andean moist Puna. We aimed to detect the differences in applying two cluster analyses-for ES values and ∆ES values-and the effects of different scales of observation-two administrative levels and two grid resolutions-on ES bundles over time. To address these objectives, we investigated two hypotheses: (1) bundles of ES differ on the method applied for assessing them; (2) these bundles are affected by the scale of observation. Our analysis uncovered consistent differences suggesting that the selection of a method for assessing bundles of ES might define the results, and the scale of observation influenced them.
"Static" bundles suggested synergies between provisioning and regulating services, whereas "dynamic" indicated negative relationships. Then, the assumption of a general pattern of trade-offs between these groups of services needs to be analyzed in detail [27,42]. The diverse interpretations found in our study suggest that both assessment methods have implications for management of ES, and both can be complementary to obtain better contributions for decision-making. However, if research objectives are focused on the understanding of the instabilities that produce the temporal dynamics on ecosystems, we recommend the assessment of "dynamic" bundles since these are more sensitive to changes of the different drivers across spatial scales. Moreover, any spatial scale can be eligible, but large administrative levels need caution.
The differences addressed over time showed confident generalization to advise the pros and cons of which spatial scale to use. The municipality level showed sufficient consistency with grid-scales, which may be enough to guide policy, as other studies highlighted [16,28]. However, for spatial conservation, the fine-grid scale could be needed to visualize small patch sizes. Then, as a rule, resulting from the study, ES bundles at grid scales characterized by a high level of dispersion and small patch size disappear or are imperceptible at administrative levels. Indeed, at heterogeneous landscapes, bundling becomes complex, whereas bundles are very similar across different spatial scales on homogeneous landscapes. In that sense, bundles at administrative levels tend to describe landscape multifunctionality, whereas fine-grained resolutions define more specialized bundles.
Finally, we have shown that the ES matrix and standard metrics give guidance to show the implications of choosing a method and a scale of observation in bundle assessment. To the best of our knowledge, this is the first study in which such a comprehensive step by step framework comparing "dynamic" and "static" bundles of ES has been developed. Bearing in mind the potential of bundles to support decision-making, the results might help the choice of bundling methods during the design of research projects. Our findings fill the knowledge gap on relationships between multiple ES utilizing cluster techniques robustly. Future studies should focus on a much more exhaustive list of ES. Additionally, more research is required to assess bundles at different spatial extensions and on landscapes with diverse levels of spatial heterogeneity.